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

// example_precision_showdown.cpp
// "double vs Rational vs Float" — Hilbert matrix precision showdown
//
// The Hilbert matrix H(i,j) = 1/(i+j+1) is one of the most ill-conditioned
// matrices in numerical analysis. This demo solves Hx = b for increasing n
// and compares the accuracy of double, Float, and Rational.
//
// Build:
//   cd build && cmake --build . --config Release --target example-precision-showdown

#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/core/mp/Float.hpp>
#include <math/core/mp/Rational.hpp>
#include <math/linalg/decomposition.hpp>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <chrono>

using namespace calx;
using namespace calx::algorithms;

// Build n x n Hilbert matrix
template<typename T>
Matrix<T> hilbert(int n) {
    Matrix<T> H(n, n);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            H(i, j) = T(1) / T(i + j + 1);
    return H;
}

// Build Hilbert matrix for Rational (no T(int) needed)
Matrix<Rational> hilbert_rational(int n) {
    Matrix<Rational> H(n, n);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < n; j++)
            H(i, j) = Rational(1, i + j + 1);
    return H;
}

// Right-hand side: b = H * [1, 1, ..., 1]^T  so the exact solution is x = [1, 1, ..., 1]
template<typename T>
Vector<T> make_rhs(const Matrix<T>& H, int n) {
    Vector<T> b(n);
    for (int i = 0; i < n; i++) {
        T s = T(0);
        for (int j = 0; j < n; j++) s = s + H(i, j);
        b[i] = s;
    }
    return b;
}

Vector<Rational> make_rhs_rational(const Matrix<Rational>& H, int n) {
    Vector<Rational> b(n);
    for (int i = 0; i < n; i++) {
        Rational s;
        for (int j = 0; j < n; j++) s = s + H(i, j);
        b[i] = s;
    }
    return b;
}

int main() {
    std::cout << "Hilbert Matrix Precision Showdown" << std::endl;
    std::cout << "==================================" << std::endl;
    std::cout << "Solve Hx = b where the exact solution is x = (1, 1, ..., 1)" << std::endl;
    std::cout << "Error = max |x_i - 1|" << std::endl;
    std::cout << std::endl;

    std::cout << std::setw(4) << "n"
              << std::setw(18) << "double"
              << std::setw(18) << "Float(200d)"
              << std::setw(18) << "Rational"
              << std::setw(16) << "cond(H) approx" << std::endl;
    std::cout << std::string(74, '-') << std::endl;

    for (int n = 1; n <= 20; n++) {
        // --- double ---
        double err_double = 0;
        {
            auto H = hilbert<double>(n);
            auto b = make_rhs(H, n);
            try {
                auto x = lu_solve(H, b);
                for (int i = 0; i < n; i++)
                    err_double = std::max(err_double, std::abs(x[i] - 1.0));
            } catch (...) {
                err_double = std::numeric_limits<double>::infinity();
            }
        }

        // --- Float (200 digits ~ 673 bits) ---
        double err_float = 0;
        {
            Float::setDefaultPrecision(200);
            auto H = hilbert<Float>(n);
            auto b = make_rhs(H, n);
            try {
                auto x = lu_solve(H, b);
                for (int i = 0; i < n; i++) {
                    Float diff = calx::abs(x[i] - Float(1));
                    err_float = std::max(err_float, diff.toDouble());
                }
            } catch (...) {
                err_float = std::numeric_limits<double>::infinity();
            }
        }

        // --- Rational (exact) ---
        double err_rational = 0;
        {
            auto H = hilbert_rational(n);
            auto b = make_rhs_rational(H, n);
            auto x = gaussian_elimination(H, b, PivotStrategy::Partial);
            for (int i = 0; i < n; i++) {
                if (x[i] != Rational(1)) {
                    err_rational = 1.0; // should never happen
                }
            }
        }

        // Approximate condition number (product of diagonal of U / min diagonal)
        double log_cond = 0;
        {
            auto H = hilbert<double>(n);
            try {
                auto [LU, piv] = lu_decomposition(H);
                double max_d = 0, min_d = 1e300;
                for (int i = 0; i < n; i++) {
                    double d = std::abs(LU(i, i));
                    if (d > max_d) max_d = d;
                    if (d > 0 && d < min_d) min_d = d;
                }
                log_cond = std::log10(max_d / min_d) * n * 0.5; // rough estimate
            } catch (...) {}
        }

        // Format output
        auto fmt_err = [](double e) -> std::string {
            char buf[32];
            if (e == 0.0)
                snprintf(buf, sizeof(buf), "0 (exact)");
            else if (e >= 1.0)
                snprintf(buf, sizeof(buf), "FAILED");
            else if (std::isinf(e))
                snprintf(buf, sizeof(buf), "SINGULAR");
            else
                snprintf(buf, sizeof(buf), "%.2e", e);
            return buf;
        };

        std::cout << std::setw(4) << n
                  << std::setw(18) << fmt_err(err_double)
                  << std::setw(18) << fmt_err(err_float)
                  << std::setw(18) << fmt_err(err_rational);
        if (log_cond > 0)
            std::cout << std::setw(12) << "~10^" << (int)log_cond;
        std::cout << std::endl;
    }

    std::cout << std::endl;
    std::cout << "Observations:" << std::endl;
    std::cout << "  - double loses all accuracy around n=13 (condition number ~10^17)" << std::endl;
    std::cout << "  - Float(200 digits) maintains accuracy up to n=20+" << std::endl;
    std::cout << "  - Rational always gives the exact answer (error = 0)" << std::endl;

    // --- Bonus: Show the exact inverse of H(5) ---
    std::cout << std::endl;
    std::cout << "Bonus: Exact inverse of H(5) (all entries are integers!):" << std::endl;
    {
        auto H = hilbert_rational(5);
        auto Hinv = lu_inverse(H);
        for (int i = 0; i < 5; i++) {
            for (int j = 0; j < 5; j++)
                std::cout << std::setw(9) << Hinv(i, j);
            std::cout << std::endl;
        }
    }

    return 0;
}