// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later
//
// example_doc_linalg_samples.cpp
// Verification harness for samples in api/LinAlg.html (ja+en) "Examples".

#include <math/linalg/decomposition.hpp>
#include <math/linalg/solvers.hpp>
#include <math/core/mp/Rational.hpp>
#include <iomanip>
#include <iostream>

using namespace sangi;
using namespace sangi::algorithms;

template <typename V>
static void print_vec(const V& v) {
    std::cout << "[";
    for (std::size_t i = 0; i < v.size(); ++i) {
        if (i) std::cout << ", ";
        std::cout << v[i];
    }
    std::cout << "]";
}

int main() {
    std::cout << std::setprecision(15);

    // ---- 1. Solve a linear system via LU decomposition -----------
    {
        std::cout << "[1. LU decomp + solve]\n";
        Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
        Vector<double> b({8.0, -11.0, -3.0});
        auto x = solve(A, b);
        std::cout << "  x = "; print_vec(x); std::cout << "  (expected [2, 3, -1])\n";

        auto det = lu_determinant(A);
        std::cout << "  det(A) = " << det << "\n";

        // For the docs: print the return value of lu_decomposition directly
        auto [LU, perm] = lu_decomposition(A);
        std::cout << "  LU =\n";
        for (std::size_t i = 0; i < LU.rows(); ++i) {
            std::cout << "    ";
            for (std::size_t j = 0; j < LU.cols(); ++j) std::cout << LU(i, j) << " ";
            std::cout << "\n";
        }
        std::cout << "  perm = "; print_vec(perm); std::cout << "\n";
    }

    // ---- 2. SVD low-rank approximation ---------------------------
    {
        std::cout << "[2. SVD low-rank]\n";
        Matrix<double> A = {{1, 2, 3}, {4, 5, 6}, {7, 8, 9}};
        auto [U, sigma, Vt] = svd_decomposition(A);
        std::cout << "  Singular values = ";
        for (std::size_t i = 0; i < sigma.size(); ++i) std::cout << sigma[i] << " ";
        std::cout << "\n  (expected ≈ 16.85, 1.07, 0; rank = 2)\n";

        Matrix<double> U1(U.rows(), 1), Vt1(1, Vt.cols());
        for (std::size_t i = 0; i < U.rows(); ++i) U1(i, 0) = U(i, 0);
        for (std::size_t j = 0; j < Vt.cols(); ++j) Vt1(0, j) = Vt(0, j);
        Vector<double> s1({sigma[0]});
        auto A1 = reconstruct_matrix_from_svd(U1, s1, Vt1);
        std::cout << "  rank-1 approx (0,0) = " << A1(0, 0)
                  << "   (rank-1 ≈ 1.74; full ≈ 1)\n";
    }

    // ---- 3. gaussian_elimination + Full pivot --------------------
    {
        std::cout << "[3. Gauss elim with PivotStrategy]\n";
        Matrix<double> A = {
            {1e-18, 1.0, 1.0},
            {1.0,   1.0, 2.0},
            {1.0,   2.0, 1.0}
        };
        Vector<double> b({2.0, 4.0, 4.0});
        auto x_partial = gaussian_elimination(A, b, PivotStrategy::Partial);
        auto x_full    = gaussian_elimination(A, b, PivotStrategy::Full);
        std::cout << "  Partial: "; print_vec(x_partial); std::cout << "\n";
        std::cout << "  Full:    "; print_vec(x_full); std::cout << "\n";
    }

    // ---- 4. Hilbert exact solve via Rational ---------------------
    {
        std::cout << "[4. Hilbert exact Rational]\n";
        const std::size_t n = 4;
        Matrix<Rational> H(n, n);
        for (std::size_t i = 0; i < n; ++i)
            for (std::size_t j = 0; j < n; ++j)
                H(i, j) = Rational(1, static_cast<int>(i + j + 1));

        Vector<Rational> bb(n);
        for (std::size_t i = 0; i < n; ++i) {
            Rational s(0);
            for (std::size_t j = 0; j < n; ++j) s = s + H(i, j);
            bb[i] = s;
        }
        auto x = gaussian_elimination(H, bb);
        std::cout << "  x = [";
        for (std::size_t i = 0; i < n; ++i) {
            if (i) std::cout << ", ";
            std::cout << x[i].toString();
        }
        std::cout << "]   (expected [1, 1, 1, 1] exact)\n";
    }

    return 0;
}