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

#include <math/core/matrix.hpp>
#include <math/linalg/solvers.hpp>
#include <cmath>
#include <iomanip>
#include <iostream>

using namespace sangi;

void print_mat(const Matrix<double>& m) {
    std::cout << '{';
    for (std::size_t i = 0; i < m.rows(); ++i) {
        std::cout << (i ? ", {" : "{");
        for (std::size_t j = 0; j < m.cols(); ++j) {
            if (j) std::cout << ", ";
            std::cout << m(i, j);
        }
        std::cout << '}';
    }
    std::cout << '}';
}

template <typename V>
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);

    // ---- arithmetic ----------------------------------------------
    {
        std::cout << "[arithmetic]\n";
        Matrix<double> A = {{1, 2}, {3, 4}};
        Matrix<double> B = {{5, 6}, {7, 8}};
        Matrix<double> C = A * B;
        std::cout << "  A*B = "; print_mat(C); std::cout << "\n";

        Vector<double> v({1.0, 2.0});
        Vector<double> w = A * v;
        std::cout << "  A*v = "; print_vec(w); std::cout << "\n";

        Matrix<double> At = A ^ 'T';
        std::cout << "  A^T = "; print_mat(At); std::cout << "\n";

        Matrix<double> Ainv = A ^ (-1);
        std::cout << "  A^-1 = "; print_mat(Ainv); std::cout << "\n";

        Matrix<double> A3 = A ^ 3;
        std::cout << "  A^3 = "; print_mat(A3); std::cout << "\n";
    }

    // ---- determinant / inverse / trace / norm --------------------
    {
        std::cout << "[determinant/inverse/trace/norm]\n";
        Matrix<double> A = {{1, 2}, {3, 4}};
        std::cout << "  det(A) = " << A.determinant() << "\n";
        std::cout << "  trace(A) = " << A.trace() << "\n";
        std::cout << "  norm_F(A) = " << A.norm() << "\n";
        Matrix<double> Ainv = A.inverse();
        Matrix<double> I = A * Ainv;
        std::cout << "  A*A^-1 = "; print_mat(I); std::cout << "\n";
    }

    // ---- identity factory -----------------------------------------
    {
        std::cout << "[identity factory]\n";
        Matrix<double> I3 = Matrix<double>::identity(3);
        std::cout << "  identity(3) = "; print_mat(I3); std::cout << "\n";
    }

    // ---- solve via Matrix.solve() (newly added member) -----------
    {
        std::cout << "[Matrix.solve]\n";
        // A x = b, where A = [[4, 1], [1, 3]] and b = [9, 8]
        Matrix<double> A = {{4, 1}, {1, 3}};
        Vector<double> b({9.0, 8.0});
        auto x_opt = A.solve(b);
        if (x_opt) {
            std::cout << "  x = "; print_vec(*x_opt); std::cout << "\n";
            Vector<double> r = A * (*x_opt);
            std::cout << "  A*x = "; print_vec(r); std::cout << " (should be b)\n";
        } else {
            std::cout << "  singular system\n";
        }
    }

    return 0;
}