Rational デモ — 多倍長有理数の威力

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

// example_rational.cpp
// 多倍長有理数デモ: Rational の威力を4つのシナリオで実演
//
// 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: 浮動小数点の罠 vs 有理数の厳密性
//   — 身近な計算で double が失敗する例
// ============================================================================
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";
    }

    // Σ(1/k) を前方加算 vs 後方加算 vs 有理数
    {
        constexpr int N = 100;

        // double: 前方加算 (1→N)
        double fwd = 0.0;
        for (int k = 1; k <= N; ++k) fwd += 1.0 / k;

        // double: 後方加算 (N→1)
        double bwd = 0.0;
        for (int k = N; k >= 1; --k) bwd += 1.0 / k;

        // Rational: 厳密値
        Rational exact = harmonicNumber(N);

        std::cout << "  調和数 H_" << N << " = 1 + 1/2 + 1/3 + ... + 1/" << N << ":\n";
        std::cout << std::setprecision(17);
        std::cout << "    double (前方加算): " << fwd << std::endl;
        std::cout << "    double (後方加算): " << bwd << std::endl;
        std::cout << "    Rational (厳密) : " << exact.toDecimal(17) << std::endl;
        std::cout << "    厳密値の分子の桁数: " << exact.numerator().toString().size() << " 桁\n";
        std::cout << "    厳密値の分母の桁数: " << exact.denominator().toString().size() << " 桁\n";
    }
    std::cout << std::endl;
}

// ============================================================================
// Demo 2: 連分数展開で円周率の最良近似を探索
//   — なぜ 355/113 は驚くほど正確なのか
// ============================================================================
static void demo_continued_fraction() {
    std::cout << "============================================================\n";
    std::cout << " Demo 2: Continued Fractions — Best Approximations of Pi\n";
    std::cout << "============================================================\n\n";

    // π の連分数展開 [3; 7, 15, 1, 292, 1, 1, 1, 2, ...]
    // 最初の数項を手動で設定（Rational は有限精度なので）
    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)
    };

    // π の近似値（比較用、十分な桁数）
    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) << "項数" << "  "
              << std::setw(20) << "近似分数" << "  "
              << std::setw(24) << "小数値" << "  "
              << "正確な桁数" << std::endl;
    std::cout << std::string(75, '-') << std::endl;

    for (size_t depth = 1; depth <= pi_cf.size(); ++depth) {
        // 連分数を depth 項まで切って有理数に復元
        std::vector<Int> partial(pi_cf.begin(), pi_cf.begin() + depth);
        Rational approx = Rational::fromContinuedFraction(partial);

        // 小数展開
        std::string dec = approx.toDecimal(20);

        // 正確な桁数を数える（小数点以降）
        int correct = 0;
        size_t start = 2;  // "3." の後ろ
        for (size_t i = start; i < std::min(dec.size(), pi_str.size()); ++i) {
            if (dec[i] == pi_str[i]) ++correct;
            else break;
        }

        // 分数表示
        std::string frac = approx.numerator().toString() + "/"
                         + approx.denominator().toString();

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

    std::cout << "\n  注目: 項 [292] の追加で 355/113 → 精度が一気に向上。\n";
    std::cout << "  292 が大きいほど直前の近似が「良い」ことを意味します。\n";
    std::cout << "  355/113 は分母3桁で小数6桁正確 — 古代中国の祖沖之が発見。\n";
    std::cout << std::endl;
}

// ============================================================================
// Demo 3: 調和数 H_n の厳密値と Bernoulli 数
//   — 数論の宝石を多倍長有理数で計算
// ============================================================================
static void demo_harmonic_bernoulli() {
    std::cout << "============================================================\n";
    std::cout << " Demo 3: Harmonic Numbers & Bernoulli Numbers\n";
    std::cout << "============================================================\n\n";

    // 調和数
    std::cout << "  --- 調和数 H_n (厳密値) ---\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 数
    std::cout << "\n  --- Bernoulli 数 B_n (厳密値) ---\n";
    std::cout << "  (B_1 = -1/2, 奇数 n>1 では B_n = 0)\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  B_30 の分子・分母は巨大な整数:\n";
    Rational B30 = bernoulli(30);
    std::cout << "    分子: " << B30.numerator() << std::endl;
    std::cout << "    分母: " << B30.denominator() << std::endl;

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

// ============================================================================
// Demo 4: Stern-Brocot 木と Farey 数列
//   — 有理数の美しい構造を探索
// ============================================================================
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 木: nextMinimal で順に列挙
    std::cout << "  --- Stern-Brocot 順 (高さ昇順で最初の20個) ---\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 列
    std::cout << "  --- Calkin-Wilf 列 (最初の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 数列 F_7
    std::cout << "  --- Farey 数列 F_7 (分母 <= 7 の既約分数) ---\n  ";
    {
        // 0/1 から始めて mediant で構築
        std::vector<Rational> farey;
        farey.push_back(Rational(0, 1));
        farey.push_back(Rational(1, 7));

        // nextMinimal を使って 0 < p/q <= 1, q <= 7 を列挙
        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 << "  合計 " << farey.size() << " 個の既約分数\n\n";
    }

    // simplestBetween
    std::cout << "  --- simplestBetween: 2つの有理数の間の最も単純な有理数 ---\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 — 多倍長有理数の威力               #\n";
    std::cout << "###########################################################\n\n";

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

    std::cout << "###########################################################\n";
    std::cout << "#  全デモ完了                                             #\n";
    std::cout << "###########################################################\n";

    return 0;
}