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

// example_int_demo.cpp
// Multi-precision integer demo: showcasing Int's power in 5 scenarios

#include <math/core/mp/Int.hpp>
#include <math/core/mp/Int/IntGCD.hpp>
#include <math/core/mp/Int/IntSqrt.hpp>
#include <math/core/mp/Int/IntModular.hpp>
#include <math/core/mp/Int/IntPrime.hpp>
#include <math/core/mp/Int/IntCombinatorics.hpp>
#include <math/core/mp/Int/IntFactorization.hpp>
#include <iostream>
#include <iomanip>
#include <string>
#include <chrono>
#include <vector>

using namespace calx;

// ============================================================================
// Demo 1: Fibonacci Numbers -- The World of Giant Integers
//   -- F(1000) has 209 digits, F(10000) has 2092 digits
// ============================================================================
static void demo_fibonacci() {
    std::cout << "============================================================\n";
    std::cout << " Demo 1: Fibonacci Numbers -- The World of Giant Integers\n";
    std::cout << "============================================================\n\n";

    // Small values
    std::cout << "  F(10)  = ";
    {
        Int a(0), b(1);
        for (int i = 2; i <= 10; ++i) {
            Int t = a + b;
            a = std::move(b);
            b = std::move(t);
        }
        std::cout << b << std::endl;
    }

    // F(100)
    std::cout << "  F(100) = ";
    {
        Int a(0), b(1);
        for (int i = 2; i <= 100; ++i) {
            Int t = a + b;
            a = std::move(b);
            b = std::move(t);
        }
        std::cout << b << std::endl;
    }

    // F(1000)
    std::cout << "\n  F(1000) = ";
    {
        Int a(0), b(1);
        for (int i = 2; i <= 1000; ++i) {
            Int t = a + b;
            a = std::move(b);
            b = std::move(t);
        }
        std::string s = b.toString();
        // First 50 digits + ... + last 20 digits
        if (s.size() > 80) {
            std::cout << s.substr(0, 50) << "..." << s.substr(s.size() - 20)
                      << "\n  (" << s.size() << " digits)\n";
        } else {
            std::cout << s << "\n";
        }
    }

    // F(10000) -- timed
    {
        auto start = std::chrono::high_resolution_clock::now();
        Int a(0), b(1);
        for (int i = 2; i <= 10000; ++i) {
            Int t = a + b;
            a = std::move(b);
            b = std::move(t);
        }
        auto end = std::chrono::high_resolution_clock::now();
        auto us = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        std::string s = b.toString();
        std::cout << "\n  F(10000) = " << s.substr(0, 50) << "..."
                  << s.substr(s.size() - 20)
                  << "\n  (" << s.size() << " digits, " << us << " us)\n";
    }

    // Property of Fibonacci numbers: gcd(F(m), F(n)) = F(gcd(m, n))
    std::cout << "\n  Property: gcd(F(m), F(n)) = F(gcd(m, n))\n";
    {
        // Compute F(12) and F(18)
        auto fib = [](int n) -> Int {
            Int a(0), b(1);
            for (int i = 2; i <= n; ++i) {
                Int t = a + b;
                a = std::move(b);
                b = std::move(t);
            }
            return n == 0 ? a : b;
        };

        Int f12 = fib(12);
        Int f18 = fib(18);
        Int g = gcd(f12, f18);
        Int f6 = fib(6);

        std::cout << "    F(12) = " << f12 << ", F(18) = " << f18 << "\n";
        std::cout << "    gcd(F(12), F(18)) = " << g << "\n";
        std::cout << "    F(gcd(12,18)) = F(6) = " << f6 << "\n";
        std::cout << "    Match: " << ((g <=> f6) == 0 ? "YES" : "NO") << "\n";
    }
    std::cout << std::endl;
}

// ============================================================================
// Demo 2: RSA Encryption -- Mini Demo
//   -- Encryption and decryption using large primes, modular exponentiation, and modular inverse
// ============================================================================
static void demo_rsa() {
    std::cout << "============================================================\n";
    std::cout << " Demo 2: RSA Encryption -- Mini Demo\n";
    std::cout << "============================================================\n\n";

    // Using Mersenne primes (verified in Demo 3)
    // M_89 = 2^89 - 1  (27-digit prime)
    // M_107 = 2^107 - 1 (33-digit prime)
    Int p = pow(Int(2), 89u) - Int(1);
    Int q = pow(Int(2), 107u) - Int(1);

    std::cout << "  p = M_89  = " << p << " (" << p.toString().size() << " digits)\n";
    std::cout << "  q = M_107 = " << q << " (" << q.toString().size() << " digits)\n";

    // Primality check
    std::cout << "  p is prime: " << (isProbablePrime(p) ? "YES" : "NO") << "\n";
    std::cout << "  q is prime: " << (isProbablePrime(q) ? "YES" : "NO") << "\n\n";

    // n = p * q
    Int n = p * q;
    std::cout << "  n = p * q\n";
    std::cout << "    = " << n << "\n";
    std::cout << "    (" << n.toString().size() << " digits)\n\n";

    // phi(n) = (p-1)(q-1)
    Int phi = (p - Int(1)) * (q - Int(1));

    // Public exponent e = 65537
    Int e(65537);
    std::cout << "  Public key: (n, e=" << e << ")\n";

    // Private exponent d = e^-1 mod phi(n)
    Int d = IntModular::inverseMod(e, phi);
    std::cout << "  Private key: d = " << d.toString().substr(0, 40) << "...\n";
    std::cout << "    (" << d.toString().size() << " digits)\n\n";

    // Encryption: encrypt message m
    Int m(123456789);
    std::cout << "  Plaintext:  m = " << m << "\n";

    // c = m^e mod n
    Int c = powMod(m, e, n);
    std::cout << "  Ciphertext: c = " << c << "\n";

    // Decryption: m' = c^d mod n
    Int m_dec = powMod(c, d, n);
    std::cout << "  Decrypted:  m'= " << m_dec << "\n";
    std::cout << "  Match:      " << ((m <=> m_dec) == 0 ? "YES" : "NO") << "\n";

    // A slightly larger message
    Int m2("314159265358979323846264338327950288");
    std::cout << "\n  Plaintext:  m = " << m2 << "\n";
    Int c2 = powMod(m2, e, n);
    std::cout << "  Ciphertext: c = " << c2 << "\n";
    Int m2_dec = powMod(c2, d, n);
    std::cout << "  Decrypted:  m'= " << m2_dec << "\n";
    std::cout << "  Match:      " << ((m2 <=> m2_dec) == 0 ? "YES" : "NO") << "\n";

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

// ============================================================================
// Demo 3: Mersenne Prime Verification
//   -- Verify special values of p where 2^p - 1 is prime
// ============================================================================
static void demo_mersenne() {
    std::cout << "============================================================\n";
    std::cout << " Demo 3: Mersenne Primes -- 2^p - 1\n";
    std::cout << "============================================================\n\n";

    // Known Mersenne prime exponents
    int exponents[] = {2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127};
    // Non-primes (for comparison)
    int non_mersenne[] = {11, 23, 29, 37, 41, 43};

    std::cout << "  Mersenne primes M_p = 2^p - 1:\n\n";
    std::cout << std::setw(6) << "p" << "  "
              << std::setw(8) << "prime?" << "  "
              << std::setw(10) << "digits" << "  "
              << "M_p" << std::endl;
    std::cout << std::string(70, '-') << std::endl;

    for (int p : exponents) {
        Int mp = pow(Int(2), static_cast<unsigned>(p)) - Int(1);
        bool is_prime = isProbablePrime(mp);
        std::string s = mp.toString();
        std::string display = (s.size() <= 40) ? s : s.substr(0, 30) + "..." + s.substr(s.size() - 5);

        std::cout << std::setw(6) << p << "  "
                  << std::setw(8) << (is_prime ? "YES" : "no") << "  "
                  << std::setw(10) << s.size() << "  "
                  << display << std::endl;
    }

    std::cout << "\n  Examples of p where 2^p - 1 is composite (with factorization):\n";
    for (int p : non_mersenne) {
        Int mp = pow(Int(2), static_cast<unsigned>(p)) - Int(1);
        std::cout << "    M_" << p << " = " << mp << " = ";
        auto factors = IntFactorization::factorize(mp);
        for (size_t i = 0; i < factors.size(); ++i) {
            if (i > 0) std::cout << " x ";
            std::cout << factors[i].first;
            if (factors[i].second > 1) std::cout << "^" << factors[i].second;
        }
        std::cout << "\n";
    }

    std::cout << "\n  M_127 is a 39-digit prime -- verified by Lucas in 1876 by hand.\n";
    std::cout << "  The largest prime ever verified by a human before computers.\n";
    std::cout << std::endl;
}

// ============================================================================
// Demo 4: Large Factorials & Binomial Coefficients
//   -- Number of digits in 1000!, exact value of C(100,50)
// ============================================================================
static void demo_factorial_binomial() {
    std::cout << "============================================================\n";
    std::cout << " Demo 4: Factorial & Binomial Coefficients\n";
    std::cout << "============================================================\n\n";

    // Factorial
    std::cout << "  --- Factorial n! ---\n";
    for (int n : {10, 20, 50, 100, 1000}) {
        auto start = std::chrono::high_resolution_clock::now();
        Int fact = IntCombinatorics::factorial(Int(n));
        auto end = std::chrono::high_resolution_clock::now();
        auto us = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();

        std::string s = fact.toString();
        if (s.size() <= 50) {
            std::cout << "  " << std::setw(4) << n << "! = " << s << "\n";
        } else {
            std::cout << "  " << std::setw(4) << n << "! = " << s.substr(0, 30)
                      << "..." << s.substr(s.size() - 10) << " ("
                      << s.size() << " digits, " << us << " us)\n";
        }
    }

    // Number of trailing zeros in 1000!
    std::cout << "\n  Number of trailing zeros in 1000!: ";
    {
        Int fact1000 = IntCombinatorics::factorial(Int(1000));
        auto [q, count] = removeFactor(fact1000, Int(10));
        std::cout << count << "\n";
        std::cout << "  (theoretical: floor(1000/5)+floor(1000/25)+floor(1000/125)+floor(1000/625) = "
                  << (200 + 40 + 8 + 1) << ")\n";
    }

    // Binomial coefficients
    std::cout << "\n  --- Binomial Coefficient C(n, k) ---\n";
    {
        struct { int n, k; } cases[] = {{10, 5}, {20, 10}, {50, 25}, {100, 50}};
        for (auto& [n, k] : cases) {
            Int c = IntCombinatorics::binomial(Int(n), Int(k));
            std::string s = c.toString();
            std::cout << "  C(" << std::setw(3) << n << ", " << std::setw(2) << k << ") = ";
            if (s.size() <= 40) {
                std::cout << s;
            } else {
                std::cout << s.substr(0, 30) << "... (" << s.size() << " digits)";
            }
            std::cout << "\n";
        }
    }

    // Pascal's triangle
    std::cout << "\n  Pascal's Triangle (n=0..10):\n";
    for (int n = 0; n <= 10; ++n) {
        std::cout << std::string(static_cast<size_t>(30 - n * 3), ' ');
        for (int k = 0; k <= n; ++k) {
            Int c = IntCombinatorics::binomial(Int(n), Int(k));
            std::cout << std::setw(5) << c << " ";
        }
        std::cout << "\n";
    }

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

// ============================================================================
// Demo 5: Extended GCD & Modular Exponentiation
//   -- Solutions to linear Diophantine equations and modular powers of large numbers
// ============================================================================
static void demo_number_theory() {
    std::cout << "============================================================\n";
    std::cout << " Demo 5: Extended GCD & Modular Exponentiation\n";
    std::cout << "============================================================\n\n";

    // Extended Euclidean algorithm: find x, y such that ax + by = gcd(a, b)
    std::cout << "  --- Extended GCD: ax + by = gcd(a, b) ---\n\n";
    {
        struct { int a, b; } cases[] = {{35, 15}, {240, 46}, {1071, 462}};
        for (auto& [a, b] : cases) {
            Int x, y;
            Int g = IntGCD::extendedGcd(Int(a), Int(b), x, y);
            std::cout << "    " << a << " * (" << x << ") + "
                      << b << " * (" << y << ") = " << g << "\n";
            // Verification
            Int check = Int(a) * x + Int(b) * y;
            std::cout << "    Verify: " << check << " = gcd(" << a << ", " << b << ") "
                      << ((check <=> g) == 0 ? "OK" : "NG") << "\n\n";
        }
    }

    // Large modular exponentiation
    std::cout << "  --- Modular Exponentiation: a^n mod m ---\n\n";
    {
        // 2^1000 mod (10^9 + 7)
        Int mod(1000000007);
        Int result = powMod(Int(2), Int(1000), mod);
        std::cout << "    2^1000 mod (10^9+7) = " << result << "\n";

        // Fermat's little theorem: a^(p-1) = 1 (mod p)
        Int p(1000000007);
        Int a(314159265);
        Int fermat = powMod(a, p - Int(1), p);
        std::cout << "\n    Fermat's little theorem: a^(p-1) mod p = 1\n";
        std::cout << "    " << a << "^(10^9+6) mod (10^9+7) = " << fermat << "\n";
        std::cout << "    (= 1: " << ((fermat <=> Int(1)) == 0 ? "YES" : "NO") << ")\n";
    }

    // Integer square root of large numbers
    std::cout << "\n  --- Integer Square Root ---\n\n";
    {
        // sqrt(10^100) = 10^50
        Int n = pow(Int(10), 100u);
        Int s = IntSqrt::sqrt(n);
        std::cout << "    isqrt(10^100) = " << s << "\n";
        std::cout << "    (= 10^50: " << ((s <=> pow(Int(10), 50u)) == 0 ? "YES" : "NO") << ")\n";

        // sqrt(2^256)
        Int n2 = pow(Int(2), 256u);
        Int s2 = IntSqrt::sqrt(n2);
        std::cout << "\n    isqrt(2^256) = " << s2 << "\n";
        std::cout << "    (= 2^128: " << ((s2 <=> pow(Int(2), 128u)) == 0 ? "YES" : "NO") << ")\n";

        // Perfect square test
        Int perfect = Int(12345678) * Int(12345678);
        Int root;
        bool is_sq = IntSqrt::isSquare(perfect, &root);
        std::cout << "\n    12345678^2 = " << perfect << "\n";
        std::cout << "    isSquare: " << (is_sq ? "YES" : "NO")
                  << ", sqrt = " << root << "\n";
    }

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

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

    demo_fibonacci();
    demo_rsa();
    demo_mersenne();
    demo_factorial_binomial();
    demo_number_theory();

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

    return 0;
}