Rational Demo — The Power of Arbitrary-Precision Rationals

// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// example_rational.cpp
// Arbitrary-precision rational demo: 4 scenarios showcasing Rational
//
// Build:
//   cl /std:c++latest /EHsc /O2 /MD /utf-8 /I<calx>/include examples\example_rational.cpp
//      /link calx_rational.lib calx_float.lib calx_int.lib /MACHINE:X64

#include <math/core/mp/Rational.hpp>
#include <math/core/mp/Rational/RationalConstants.hpp>
#include <iostream>
#include <iomanip>
#include <string>
#include <cmath>
#include <vector>

using namespace calx;

// ============================================================================
// Demo 1: Floating-Point Trap vs Exact Rational
//   — Common computations where double fails
// ============================================================================
static void demo_float_trap() {
    std::cout << "============================================================\n";
    std::cout << " Demo 1: Floating-Point Trap vs Exact Rational\n";
    std::cout << "============================================================\n\n";

    // 0.1 + 0.2 == 0.3 ?
    {
        double da = 0.1, db = 0.2, dc = 0.3;
        Rational ra(1, 10), rb(2, 10), rc(3, 10);

        std::cout << "  [double]   0.1 + 0.2 == 0.3 ? "
                  << ((da + db == dc) ? "true" : "FALSE") << std::endl;
        std::cout << std::setprecision(20);
        std::cout << "             0.1 + 0.2 = " << (da + db) << std::endl;
        std::cout << "             0.3       = " << dc << std::endl;

        std::cout << "  [Rational] 1/10 + 2/10 == 3/10 ? "
                  << (((ra + rb) <=> rc) == 0 ? "true" : "FALSE") << std::endl;
        std::cout << "             1/10 + 2/10 = " << (ra + rb) << "\n\n";
    }

    // H_100 via forward sum vs reverse sum vs exact rational
    {
        constexpr int N = 100;

        // double: forward sum (1 -> N)
        double fwd = 0.0;
        for (int k = 1; k <= N; ++k) fwd += 1.0 / k;

        // double: reverse sum (N -> 1)
        double bwd = 0.0;
        for (int k = N; k >= 1; --k) bwd += 1.0 / k;

        // Rational: exact value
        Rational exact = harmonicNumber(N);

        std::cout << "  Harmonic number H_" << N << " = 1 + 1/2 + 1/3 + ... + 1/" << N << ":\n";
        std::cout << std::setprecision(17);
        std::cout << "    double (forward sum): " << fwd << std::endl;
        std::cout << "    double (reverse sum): " << bwd << std::endl;
        std::cout << "    Rational (exact)    : " << exact.toDecimal(17) << std::endl;
        std::cout << "    Digits in exact numerator:   " << exact.numerator().toString().size() << " digits\n";
        std::cout << "    Digits in exact denominator: " << exact.denominator().toString().size() << " digits\n";
    }
    std::cout << std::endl;
}

// ============================================================================
// Demo 2: Continued Fractions — Best Approximations of Pi
//   — Why is 355/113 so remarkably accurate?
// ============================================================================
static void demo_continued_fraction() {
    std::cout << "============================================================\n";
    std::cout << " Demo 2: Continued Fractions — Best Approximations of Pi\n";
    std::cout << "============================================================\n\n";

    // Continued fraction expansion of pi: [3; 7, 15, 1, 292, 1, 1, 1, 2, ...]
    std::vector<Int> pi_cf = {
        Int(3), Int(7), Int(15), Int(1), Int(292),
        Int(1), Int(1), Int(1), Int(2), Int(1),
        Int(3), Int(1), Int(14)
    };

    // Reference value of pi (sufficient digits for comparison)
    const std::string pi_str = "3.14159265358979323846264338327950288";

    std::cout << "  pi = [3; 7, 15, 1, 292, 1, 1, 1, 2, 1, 3, 1, 14, ...]\n\n";
    std::cout << std::setw(6) << "Terms" << "  "
              << std::setw(20) << "Approximation" << "  "
              << std::setw(24) << "Decimal value" << "  "
              << "Correct digits" << std::endl;
    std::cout << std::string(75, '-') << std::endl;

    for (size_t depth = 1; depth <= pi_cf.size(); ++depth) {
        // Reconstruct rational from first depth terms of the continued fraction
        std::vector<Int> partial(pi_cf.begin(), pi_cf.begin() + depth);
        Rational approx = Rational::fromContinuedFraction(partial);

        // Decimal expansion
        std::string dec = approx.toDecimal(20);

        // Count correct digits (after decimal point)
        int correct = 0;
        size_t start = 2;  // after "3."
        for (size_t i = start; i < std::min(dec.size(), pi_str.size()); ++i) {
            if (dec[i] == pi_str[i]) ++correct;
            else break;
        }

        // Fraction display
        std::string frac = approx.numerator().toString() + "/"
                         + approx.denominator().toString();

        std::cout << std::setw(6) << depth << "  "
                  << std::setw(20) << frac << "  "
                  << std::setw(24) << dec << "  "
                  << correct << " digits" << std::endl;
    }

    std::cout << "\n  Note: Adding term [292] dramatically improves precision from 355/113.\n";
    std::cout << "  The larger the quotient, the better the preceding approximation.\n";
    std::cout << "  355/113 achieves 6 decimal digits with just a 3-digit denominator —\n";
    std::cout << "  discovered by Zu Chongzhi in 5th-century China.\n";
    std::cout << std::endl;
}

// ============================================================================
// Demo 3: Harmonic Numbers & Bernoulli Numbers
//   — Number-theoretic gems computed with arbitrary-precision rationals
// ============================================================================
static void demo_harmonic_bernoulli() {
    std::cout << "============================================================\n";
    std::cout << " Demo 3: Harmonic Numbers & Bernoulli Numbers\n";
    std::cout << "============================================================\n\n";

    // Harmonic numbers
    std::cout << "  --- Harmonic numbers H_n (exact values) ---\n";
    for (int n : {1, 5, 10, 20, 50}) {
        Rational Hn = harmonicNumber(n);
        std::cout << "  H_" << std::setw(2) << n << " = "
                  << std::setw(8) << Hn.toDecimal(12) << "...";
        if (n <= 10) {
            std::cout << "  = " << Hn;
        }
        std::cout << std::endl;
    }

    // Bernoulli numbers
    std::cout << "\n  --- Bernoulli numbers B_n (exact values) ---\n";
    std::cout << "  (B_1 = -1/2; B_n = 0 for odd n > 1)\n\n";
    for (int n : {0, 1, 2, 4, 6, 8, 10, 12, 20, 30}) {
        Rational Bn = bernoulli(n);
        std::cout << "  B_" << std::setw(2) << n << " = " << Bn << std::endl;
    }

    std::cout << "\n  The numerator and denominator of B_30 are large integers:\n";
    Rational B30 = bernoulli(30);
    std::cout << "    Numerator:   " << B30.numerator() << std::endl;
    std::cout << "    Denominator: " << B30.denominator() << std::endl;

    std::cout << std::endl;
}

// ============================================================================
// Demo 4: Stern-Brocot Tree & Farey Sequence
//   — Exploring the beautiful structure of rational numbers
// ============================================================================
static void demo_stern_brocot_farey() {
    std::cout << "============================================================\n";
    std::cout << " Demo 4: Stern-Brocot Tree & Farey Sequence\n";
    std::cout << "============================================================\n\n";

    // Stern-Brocot tree: enumerate via nextMinimal
    std::cout << "  --- Stern-Brocot order (first 20 by ascending height) ---\n  ";
    {
        Rational x(1, 1);
        std::cout << x;
        for (int i = 1; i < 20; ++i) {
            x = nextMinimal(x);
            std::cout << ", " << x;
        }
        std::cout << ", ...\n\n";
    }

    // Calkin-Wilf sequence
    std::cout << "  --- Calkin-Wilf sequence (first 15) ---\n  ";
    {
        Rational x(1, 1);
        std::cout << x;
        for (int i = 1; i < 15; ++i) {
            x = nextCalkinWilf(x);
            std::cout << ", " << x;
        }
        std::cout << ", ...\n\n";
    }

    // Farey sequence F_7
    std::cout << "  --- Farey sequence F_7 (irreducible fractions with denominator <= 7) ---\n  ";
    {
        std::vector<Rational> farey;
        farey.push_back(Rational(0, 1));
        farey.push_back(Rational(1, 7));

        // Enumerate 0 < p/q <= 1, q <= 7 using nextMinimal
        Rational x(0, 1);
        farey.clear();
        farey.push_back(x);
        for (;;) {
            auto [left, right] = fareyNeighbors(x, Int(7));
            if (right > Rational(1)) break;
            farey.push_back(right);
            x = right;
        }

        for (size_t i = 0; i < farey.size(); ++i) {
            if (i > 0) std::cout << ", ";
            std::cout << farey[i];
        }
        std::cout << "\n";
        std::cout << "  Total: " << farey.size() << " irreducible fractions\n\n";
    }

    // simplestBetween
    std::cout << "  --- simplestBetween: simplest rational between two given rationals ---\n";
    {
        struct { Rational a, b; } cases[] = {
            {Rational(3, 7),  Rational(2, 5)},
            {Rational(1, 3),  Rational(1, 2)},
            {Rational(7, 10), Rational(5, 7)},
            {Rational(99, 100), Rational(100, 101)},
        };
        for (auto& [a, b] : cases) {
            Rational lo = (a < b) ? a : b;
            Rational hi = (a < b) ? b : a;
            Rational s = simplestBetween(lo, hi);
            std::cout << "    simplestBetween(" << lo << ", " << hi << ") = " << s
                      << std::endl;
        }
    }

    std::cout << std::endl;
}

// ============================================================================
// main
// ============================================================================
int main() {
    std::cout << "###########################################################\n";
    std::cout << "#  calx Rational Demo — Arbitrary-Precision Rationals      #\n";
    std::cout << "###########################################################\n\n";

    demo_float_trap();
    demo_continued_fraction();
    demo_harmonic_bernoulli();
    demo_stern_brocot_farey();

    std::cout << "###########################################################\n";
    std::cout << "#  All demos completed                                     #\n";
    std::cout << "###########################################################\n";

    return 0;
}