Matrix デモ — 行列演算の活用

#include <iostream>
#include <math/core/matrix.hpp>
#include <math/core/vector.hpp>
#include <math/linalg/decomposition.hpp>

using namespace calx;
using namespace calx::algorithms;

int main() {
    // 3×3 連立方程式  Ax = b
    Matrix<double> A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
    Vector<double> b = {8, -11, -3};

    // 方法1: LU 分解
    auto x1 = lu_solve(A, b);
    std::cout << "LU solve:       x = {";
    for (std::size_t i = 0; i < x1.size(); ++i)
        std::cout << (i ? ", " : "") << x1[i];
    std::cout << "}" << std::endl;

    // 方法2: ガウス消去法 (部分ピボット)
    auto x2 = gaussian_elimination(A, b, PivotStrategy::Partial);
    std::cout << "Gauss partial:  x = {";
    for (std::size_t i = 0; i < x2.size(); ++i)
        std::cout << (i ? ", " : "") << x2[i];
    std::cout << "}" << std::endl;

    // 方法3: ガウス消去法 (完全ピボット)
    auto x3 = gaussian_elimination(A, b, PivotStrategy::Full);
    std::cout << "Gauss full:     x = {";
    for (std::size_t i = 0; i < x3.size(); ++i)
        std::cout << (i ? ", " : "") << x3[i];
    std::cout << "}" << std::endl;

    // 検算: Ax = b
    Vector<double> r = A * x1 - b;
    std::cout << "Residual norm:  " << r.norm() << std::endl;
}

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

using namespace calx;
using namespace calx::algorithms;

int main() {
    // 5x4 行列 (ランク 2 構造 + 微小ノイズ)
    // A ≈ 3*u1*v1^T + 2*u2*v2^T
    Matrix<double> A = {
        { 3.01,  0.00,  6.00,  3.00},
        {10.00,  2.00, 11.98,  8.00},
        { 4.00,  2.00,  0.00,  2.00},
        {13.00,  2.02, 18.00, 11.00},
        {11.00,  4.00,  6.00,  6.99}
    };

    auto [U, sigma, Vt] = svd_decomposition(A);

    std::cout << "Singular values:" << std::endl;
    for (std::size_t i = 0; i < sigma.size(); i++)
        std::cout << "  sigma[" << i << "] = " << std::fixed
                  << std::setprecision(4) << sigma[i] << std::endl;

    // ランク 2 近似: 上位 2 個の特異値のみ残す
    std::size_t k = 2;
    Vector<double> s_trunc(sigma.size(), 0.0);
    for (std::size_t i = 0; i < k; i++) s_trunc[i] = sigma[i];
    auto A_approx = reconstruct_matrix_from_svd(U, s_trunc, Vt);

    // 近似誤差
    double err = 0;
    for (std::size_t i = 0; i < A.rows(); i++)
        for (std::size_t j = 0; j < A.cols(); j++)
            err += (A(i,j) - A_approx(i,j)) * (A(i,j) - A_approx(i,j));
    std::cout << "\nRank-2 approximation error (Frobenius): "
              << std::scientific << std::setprecision(2) << std::sqrt(err) << std::endl;
    std::cout << "Compression: " << A.rows() * A.cols() << " -> "
              << k * (A.rows() + A.cols() + 1) << " values" << std::endl;
}

// example_precision_showdown.cpp
// 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 <limits>

using namespace calx;
using namespace calx::algorithms;

// Hilbert matrix H(i,j) = 1/(i+j+1)
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;
}

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;
}

int main() {
    std::cout << std::setw(4) << "n"
              << std::setw(18) << "double"
              << std::setw(18) << "Float(200d)"
              << std::setw(18) << "Rational" << std::endl;

    for (int n = 1; n <= 20; n++) {
        // true solution: x1 = (1, 1, ..., 1)

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

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

        // --- Rational (exact) ---
        double err_r = 0;
        {
            auto H = hilbert_rational(n);
            Vector<Rational> x1(n, Rational(1));
            Vector<Rational> b = H * x1;    // b = H * x1
            auto x = gaussian_elimination(H, b, PivotStrategy::Partial);
            for (int i = 0; i < n; i++)
                if (x[i] != Rational(1)) err_r = 1.0;
        }

        // output
        auto fmt = [](double e) -> std::string {
            char buf[32];
            if (e == 0) return "0 (exact)";
            if (e >= 1.0) return "FAILED";
            if (std::isinf(e)) return "SINGULAR";
            snprintf(buf, sizeof(buf), "%.2e", e);
            return buf;
        };
        std::cout << std::setw(4) << n
                  << std::setw(18) << fmt(err_d)
                  << std::setw(18) << fmt(err_f)
                  << std::setw(18) << fmt(err_r) << std::endl;
    }
}