Skip to content

Instantly share code, notes, and snippets.

@charmoniumQ
Last active December 4, 2018 01:35
Show Gist options
  • Save charmoniumQ/8698f1adfcffa53bd89efc4dbb10df88 to your computer and use it in GitHub Desktop.
Save charmoniumQ/8698f1adfcffa53bd89efc4dbb10df88 to your computer and use it in GitHub Desktop.
First ten-digit prime in the consecutive digits of e
/*
clang++ -g prime_in_e.cpp -std=c++11 -lgmp -lgmpxx -Wall -Werror -g -o test && ./test
*/
#include <iostream>
#include <vector>
#include <gmp.h>
#include <gmpxx.h>
#include <string>
#include <sstream>
#include <tuple>
#include <algorithm>
using namespace std;
typedef mpz_class integer;
typedef mpq_class rational;
typedef unsigned short int digit;
typedef unsigned short int small_int;
typedef unsigned long long tinteger; // true integer, fast-on-this-machine integer, can store 10-digits
vector<small_int> e_cfrac(small_int n);
rational cfrac2ratio(vector<small_int> digits);
tuple<integer, rational> int_frac_part(rational r);
vector<digit> ratio2digits(rational r, unsigned long n_digits);
string digits2str(vector<digit> digits);
tuple<integer, rational> int_frac_part(rational r);
tuple<integer, integer> num_den(rational r);
vector<digit> e_approx(small_int n);
small_int log10(rational);
tuple<small_int, tinteger> pow2_decomp(tinteger x);
vector<small_int> e_cfrac(small_int n) {
// generates n terms of https://oeis.org/A003417
// whichi s a continued fraction for e
n = n * 3;
vector<small_int> terms (n);
terms[0] = 2;
terms[1] = 1;
terms[2] = 2;
for (small_int i = 3; i < n; i += 3) {
terms[i] = 1;
terms[i+1] = 1;
terms[i+2] = terms[i-1] + 2;
}
return terms;
}
rational cfrac2ratio(vector<small_int> terms) {
// turns the truncated continued fraction into a rational
rational approx (1);
for (auto term = terms.crbegin(); term != terms.crend(); ++term) {
approx = rational(*term) + rational(1) / approx;
}
return approx;
}
tuple<integer, integer> num_den(rational r) {
integer num, den;
mpq_get_num(num.get_mpz_t(), r.get_mpq_t());
mpq_get_den(den.get_mpz_t(), r.get_mpq_t());
return tie(num, den);
}
vector<digit> ratio2digits(rational r, unsigned long n_digits) {
vector<digit> digits (n_digits);
integer int_part;
rational frac_part;
tie(int_part, frac_part) = int_frac_part(r);
integer num, den;
tie(num, den) = num_den(frac_part);
num *= integer(10);
//cout << r.get_str() << " = " << int_part.get_str() << " + " << frac_part.get_str() << endl;
for (tinteger i = 0; i < n_digits; ++i) {
integer dig = num / den;
// cout << num.get_str() << " / " << den.get_str() << " = " << dig.get_str() << endl;
digits[i] = static_cast<digit>(dig.get_si());
num = (num - dig * den) * 10;
}
return digits;
}
string digits2str(vector<digit> digits) {
ostringstream o;
for (const digit digit : digits) {
o << digit;
}
return o.str();
}
vector<digit> e_approx(small_int n) {
// returns the significant digits of an approximation of e based on the
// first n terms of the continued fraction
rational e1 = cfrac2ratio(e_cfrac(n));
rational e2 = cfrac2ratio(e_cfrac(n+1));
rational err = abs(e1 - e2);
digit sig_figs = -log10(err) - 2;
return ratio2digits(cfrac2ratio(e_cfrac(n)), sig_figs);
}
void test1() {
cout << "e = 2." + digits2str(e_approx(100)) << endl;
}
small_int log10(rational r) {
// computes log10 rounded away from zero
if (r == 1) {
return 0;
} else if (r > 1) {
for (small_int i = 0; ; ++i) {
r /= 10;
if (r <= 1) {
return i;
}
}
} else if (0 < r && r < 1) {
for (small_int i = 0; ; ++i) {
r *= 10;
if (r >= 1) {
return -i;
}
}
} else {
throw "r must be positive";
}
}
tinteger mod_exp(tinteger base, tinteger exponent, tinteger modulus) {
// computes base ^ exponent (mod modulus)
// https://en.wikipedia.org/wiki/Modular_exponentiation#Right-to-left_binary_method
if (modulus == 1) { return 0; }
tinteger result = 1;
base = base % modulus;
while (exponent > 0) {
if (exponent % 2 == 1) {
result = (result * base) % modulus;
}
exponent = exponent >> 1;
base = (base * base) % modulus;
}
return result;
}
tuple<small_int, tinteger> pow2_decomp(tinteger x) {
// returns a tuple of i and d where d is odd and 2^i * d = x
for (int i = 0; ; ++i) {
if ((x & 1) == 1) {
return tie(i, x);
}
x >>= 1;
}
}
// bool MR_is_prime(tinteger n) {
// cout << "testing " << n << endl;
// // https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
// unsigned s;
// tinteger d;
// tie(s, d) = pow2_decomp(n - 1);
// cout << "s, d, n = " << s << ", " << d << ", " << n << endl;
// cout << "2**s*d == n - 1" << endl;
// // https://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test#Deterministic_variants
// const small_int RABIN_TESTS = 4;
// tinteger as[RABIN_TESTS] = {2, 13, 23, 1662803};
// for (small_int i = 0; i < RABIN_TESTS; ++i) {
// tinteger a = as[i];
// cout << "a = " << a << endl;
// if (mod_exp(a, d, n) == 1) {
// cout << "a**d % n == 1" << endl;
// continue; // probable prime, do another iteration
// }
// for (small_int r = 0; r < s; ++r) {
// if (mod_exp(a, (1 << r)*d, n) == -1) {
// cout << "r = " << r << endl;
// cout << "a**(2**r*d) % n == -1 where r = " << r << endl;
// cout << "2**s*d == n - 1" << endl;
// goto continue_outer; // probable prime, do another iteration
// }
// }
// cout << a << " witnesses that " << n << " is composite" << endl;
// return false; // gotten this far without break-ing
// continue_outer: continue;
// }
// return true; // gotten this far without return-ing
// }
bool my_is_prime(tinteger n) {
return mpz_probab_prime_p(integer(long(n)).get_mpz_t(), 30) != 0;
}
void test2() {
// http://www.bigprimes.net/archive/prime/70000/
const int N_PRIMES = 100;
tinteger primes[N_PRIMES] = {
122947973, 122948509, 122948929, 122949353,
122947991, 122948513, 122948933, 122949371,
122948017, 122948519, 122948941, 122949373,
122948057, 122948533, 122948963, 122949377,
122948117, 122948537, 122948983, 122949391,
122948143, 122948557, 122948993, 122949419,
122948167, 122948561, 122948999, 122949439,
122948173, 122948659, 122949011, 122949461,
122948213, 122948681, 122949031, 122949469,
122948227, 122948689, 122949041, 122949521,
122948257, 122948713, 122949061, 122949523,
122948269, 122948723, 122949077, 122949553,
122948311, 122948729, 122949089, 122949581,
122948327, 122948747, 122949121, 122949601,
122948339, 122948773, 122949157, 122949653,
122948341, 122948779, 122949163, 122949667,
122948359, 122948789, 122949173, 122949679,
122948377, 122948803, 122949179, 122949707,
122948383, 122948849, 122949199, 122949713,
122948401, 122948857, 122949251, 122949721,
122948431, 122948887, 122949257, 122949727,
122948443, 122948899, 122949269, 122949733,
122948453, 122948909, 122949289, 122949793,
122948461, 122948921, 122949319, 122949811,
122948473, 122948927, 122949331, 122949823,
};
for (tinteger i = primes[0]; i <= primes[N_PRIMES - 1]; ++i) {
bool is_prime = find(primes, primes + N_PRIMES, i) != primes + N_PRIMES;
if (is_prime != my_is_prime(i)) {
cout << "Prime test failed for " << i << (is_prime ? " (prime)" : " (not prime)") << endl;
}
}
}
void test3() {
vector<digit> digits = e_approx(100);
for (small_int i = 0; i < digits.size() - 10; ++i) {
tinteger n = 0;
for (small_int j = 0; j < 10; ++j) {
n = n * 10 + digits[i + j];
}
if (my_is_prime(n)) {
cout << "beggining at the " << i << "th digit after the decimal" << endl;
cout << n << endl;
break;
}
}
}
int main() {
test3();
return 0;
}
tuple<integer, rational> int_frac_part(rational r) {
integer num, den;
tie(num, den) = num_den(r);
integer int_part = num / den;
integer remainder = num - den * int_part;
rational frac_part = rational(remainder) / rational(den);
return tie(int_part, frac_part);
}
rational abs(rational r) {
rational res;
mpq_abs(res.get_mpq_t(), r.get_mpq_t());
return res;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment