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

// example_float_demo.cpp
// Multi-precision floating-point demo: showcasing Float's power in 5 scenarios

#include <math/core/mp/Float.hpp>
#include <math/core/mp/Float/FloatMath.hpp>
#include <iostream>
#include <iomanip>
#include <string>
#include <chrono>
#include <cmath>

using namespace calx;

// ============================================================================
// Demo 1: Mathematical Constants to 100 Digits
//   -- Compute pi, e, gamma, log2, sqrt(2), Catalan to high precision
// ============================================================================
static void demo_constants() {
    std::cout << "============================================================\n";
    std::cout << " Demo 1: Mathematical Constants -- 100 Digits\n";
    std::cout << "============================================================\n\n";

    int prec = 110;  // compute with extra margin
    int show = 100;

    struct {
        const char* name;
        const char* symbol;
        Float (*func)(int);
    } constants[] = {
        // Compute dependencies first since pi internally caches other constants
        {"sqrt 2",  "sqrt(2)",  Float::sqrt2},
        {"log 2",   "ln 2",     Float::log2},
        {"log 10",  "ln 10",    Float::log10},
        {"e",       "e",        Float::e},
        {"gamma",   "gamma",    Float::euler},
        {"catalan", "G",        Float::catalan},
        {"pi",      "pi",       Float::pi},
    };

    int idx = 0;
    for (auto& [name, symbol, func] : constants) {
        // Use a different precision for each constant to avoid cache hits
        int p = prec + idx++;
        auto start = std::chrono::high_resolution_clock::now();
        Float val = func(p);
        auto end = std::chrono::high_resolution_clock::now();
        auto us = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        std::string s = val.toDecimalString(show);
        // Line break every 60 digits
        std::cout << "  " << std::setw(8) << std::left << symbol << " = ";
        for (size_t i = 0; i < s.size(); i += 60) {
            if (i > 0) std::cout << std::string(13, ' ');
            size_t len = std::min<size_t>(60, s.size() - i);
            std::cout << s.substr(i, len) << "\n";
        }
        std::cout << std::string(13, ' ') << "(" << us << " us)\n\n";
    }
}

// ============================================================================
// Demo 2: Pi to 10,000 Digits (Chudnovsky)
//   -- Demonstrating computation speed
// ============================================================================
static void demo_pi_digits() {
    std::cout << "============================================================\n";
    std::cout << " Demo 2: Pi -- 1,000 and 10,000 Digits (Chudnovsky)\n";
    std::cout << "============================================================\n\n";

    // Compare computation speed at different digit counts
    int benchmarks[] = {1000, 10000, 100000, 1000000};
    for (int digits : benchmarks) {
        auto start = std::chrono::high_resolution_clock::now();
        Float pi = Float::pi(digits);
        auto end = std::chrono::high_resolution_clock::now();
        auto ms = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        if (digits == 1000) {
            // Show first 200 digits for 1000-digit computation
            std::string s = pi.toDecimalString(200);
            std::cout << "  pi (1,000 digits, first 200 shown):\n\n";
            for (size_t i = 0; i < s.size(); i += 60) {
                size_t len = std::min<size_t>(60, s.size() - i);
                std::cout << "    " << s.substr(i, len) << "\n";
            }
            std::cout << "    ...\n";
        }

        if (ms >= 1000) {
            std::cout << "  " << std::setw(10) << digits << " digits: "
                      << (ms / 1000) << " ms\n";
        } else {
            std::cout << "  " << std::setw(10) << digits << " digits: "
                      << ms << " us\n";
        }
    }

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

// ============================================================================
// Demo 3: Ramanujan's "Almost Integer"
//   -- e^(pi*sqrt(163)) is extremely close to an integer
// ============================================================================
static void demo_ramanujan() {
    std::cout << "============================================================\n";
    std::cout << " Demo 3: Ramanujan's \"Almost Integer\"\n";
    std::cout << "============================================================\n\n";

    int prec = 60;

    // e^(pi*sqrt(163)) ~ 262537412640768744 (extremely close to an integer)
    Float pi_val = Float::pi(prec);
    Float sqrt163 = sqrt(Float(163), prec);
    Float result = exp(pi_val * sqrt163, prec);

    std::string s = result.toDecimalString(40);

    std::cout << "  e^(pi * sqrt(163)) -- almost an integer:\n\n";
    std::cout << "    pi       = " << pi_val.toDecimalString(30) << "...\n";
    std::cout << "    sqrt(163)= " << sqrt163.toDecimalString(30) << "...\n";
    std::cout << "\n    e^(pi * sqrt(163)) =\n";
    std::cout << "      " << s << "\n";

    std::cout << "\n  The fractional part is 0.99999999999925...\n";
    std::cout << "  The difference from integer 262537412640768744 is about 7.5e-13.\n";
    std::cout << "  This phenomenon is explained by j-invariants and complex multiplication.\n";

    // Other "almost integers" from Heegner numbers
    std::cout << "\n  Other Heegner numbers:\n";
    for (int d : {19, 43, 67, 163}) {
        Float sq = sqrt(Float(d), prec);
        Float val = exp(pi_val * sq, prec);
        std::string vs = val.toDecimalString(20);
        std::cout << "    e^(pi*sqrt(" << std::setw(3) << d << ")) = " << vs << "\n";
    }

    std::cout << "\n  Larger d gives values closer to integers.\n";
    std::cout << std::endl;
}

// ============================================================================
// Demo 4: Breaking the Precision Barrier
//   -- Accurately compute what double cannot handle, using multi-precision
// ============================================================================
static void demo_precision_limit() {
    std::cout << "============================================================\n";
    std::cout << " Demo 4: Breaking the Precision Barrier\n";
    std::cout << "============================================================\n\n";

    // Catastrophic cancellation: (1 + x) - 1  (when x is extremely small)
    std::cout << "  --- Catastrophic Cancellation ---\n\n";
    {
        // x = 1e-15
        double xd = 1e-15;
        double rd = (1.0 + xd) - 1.0;

        int prec = 60;
        Float one(1);
        one.setPrecision(prec);
        Float xf = pow(Float(10), -15, prec);
        Float rf = (one + xf) - one;

        std::cout << "    x = 1e-15\n";
        std::cout << "    (1 + x) - 1:\n";
        std::cout << std::setprecision(20);
        std::cout << "      double: " << rd << "\n";
        std::cout << "      Float:  " << rf.toDecimalString(30) << "\n";
        std::cout << "      exact:  0.000000000000001\n\n";
    }

    // Basel problem: sum(1/k^2, k=1..N) -> pi^2/6
    std::cout << "  --- Basel Problem: sum(1/k^2) -> pi^2/6 ---\n\n";
    {
        int prec = 50;
        Float pi_val = Float::pi(prec);
        Float exact = sqr(pi_val, prec) / Float(6);

        // double: sum N=100000
        double sum_d = 0.0;
        for (int k = 1; k <= 100000; ++k)
            sum_d += 1.0 / (static_cast<double>(k) * k);

        // Float: sum N=1000
        Float sum_f = Float::zero();
        for (int k = 1; k <= 1000; ++k) {
            Float kf(k);
            sum_f = sum_f + Float(1) / (kf * kf);
        }

        std::cout << "    pi^2 / 6 (exact) = " << exact.toDecimalString(30) << "\n";
        std::cout << std::setprecision(15);
        std::cout << "    sum(1/k^2, k=1..100000) [double] = " << sum_d << "\n";
        std::cout << "    sum(1/k^2, k=1..1000)   [Float]  = " << sum_f.toDecimalString(30) << "\n";
        std::cout << "    (finite sum is less than exact value -- shows slow convergence)\n";
    }

    // High-precision verification of e^(ln x) = x
    std::cout << "\n  --- Identity Verification: exp(log(x)) = x ---\n\n";
    {
        int prec = 100;
        Float x("3.14159265358979323846264338327950288419716939937510");
        x.setPrecision(prec);

        Float lx = log(x, prec);
        Float elx = exp(lx, prec);

        std::cout << "    x              = " << x.toDecimalString(50) << "\n";
        std::cout << "    exp(log(x))    = " << elx.toDecimalString(50) << "\n";

        Float diff = abs(x - elx);
        std::cout << "    |x - exp(log(x))| = " << diff.toDecimalString(50) << "\n";
        std::cout << "    (computed at 50-digit precision, so internal error is below 10^-50)\n";
    }

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

// ============================================================================
// Demo 5: Independent Verification of Pi via Machin-type Formulas
//   -- pi/4 = 4*arctan(1/5) - arctan(1/239)
// ============================================================================
static void demo_machin() {
    std::cout << "============================================================\n";
    std::cout << " Demo 5: Machin's Formula -- Independent Verification of Pi\n";
    std::cout << "============================================================\n\n";

    int prec = 150;
    int show = 100;

    // Set precision before division (dividing at default precision loses digits)
    auto makeRecip = [&](int d) {
        Float num(1); num.setPrecision(prec);
        Float den(d); den.setPrecision(prec);
        return num / den;
    };

    // Machin's formula: pi/4 = 4*arctan(1/5) - arctan(1/239)
    std::cout << "  Machin's formula: pi/4 = 4*arctan(1/5) - arctan(1/239)\n\n";
    {
        Float one_5th = makeRecip(5);
        Float one_239th = makeRecip(239);

        auto start = std::chrono::high_resolution_clock::now();
        Float atan_1_5 = atan(one_5th, prec);
        Float atan_1_239 = atan(one_239th, prec);
        Float pi_machin = (atan_1_5 * Float(4) - atan_1_239) * Float(4);
        auto end = std::chrono::high_resolution_clock::now();
        auto us = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        auto start2 = std::chrono::high_resolution_clock::now();
        Float pi_chud = Float::pi(prec + 10);  // avoid cache hit
        auto end2 = std::chrono::high_resolution_clock::now();
        auto us2 = std::chrono::duration_cast<std::chrono::microseconds>(end2 - start2).count();

        std::cout << "    Machin:     " << pi_machin.toDecimalString(show) << "\n";
        std::cout << "    Chudnovsky: " << pi_chud.toDecimalString(show) << "\n";

        Float diff = abs(pi_machin - pi_chud);
        std::cout << "\n    Diff: " << diff.toDecimalString(show) << "\n";
        std::cout << "    (Machin: " << us << " us, Chudnovsky: " << us2 << " us)\n";
    }

    // Stormer's formula: pi/4 = 6*arctan(1/8) + 2*arctan(1/57) + arctan(1/239)
    std::cout << "\n  Stormer's formula: pi/4 = 6*atan(1/8) + 2*atan(1/57) + atan(1/239)\n\n";
    {
        auto start = std::chrono::high_resolution_clock::now();
        Float a8 = atan(makeRecip(8), prec);
        Float a57 = atan(makeRecip(57), prec);
        Float a239 = atan(makeRecip(239), prec);
        Float pi_stormer = (a8 * Float(6) + a57 * Float(2) + a239) * Float(4);
        auto end = std::chrono::high_resolution_clock::now();
        auto us = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        Float pi_chud = Float::pi(prec + 20);  // avoid cache hit

        std::cout << "    Stormer:    " << pi_stormer.toDecimalString(show) << "\n";
        std::cout << "    Chudnovsky: " << pi_chud.toDecimalString(show) << "\n";

        Float diff = abs(pi_stormer - pi_chud);
        std::cout << "\n    Diff: " << diff.toDecimalString(show) << "\n";
        std::cout << "    (Stormer: " << us << " us)\n";
    }

    std::cout << "\n  Obtaining identical results from different algorithms is\n";
    std::cout << "  strong evidence of correctness for calx's multi-precision operations (exp, log, atan).\n";
    std::cout << std::endl;
}

// ============================================================================
// main
// ============================================================================
int main() {
    std::cout << "###########################################################\n";
    std::cout << "#  calx Float Demo -- Power of Multi-Precision Floats      #\n";
    std::cout << "###########################################################\n\n";

    demo_constants();
    demo_pi_digits();
    demo_ramanujan();
    demo_precision_limit();
    demo_machin();

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

    return 0;
}