First ten-digit prime in the consecutive digits of e

prime_in_e.cpp

/*
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;
}