sangi

コード例

Int — 多倍長整数

#include <math/core/mp/Int.hpp>
#include <iostream>
using namespace sangi;

int main() {
    Int a("123456789012345678901234567890");
    Int b("987654321098765432109876543210");

    std::cout << "a + b = " << a + b << std::endl;
    std::cout << "a * b = " << a * b << std::endl;
    std::cout << "gcd   = " << gcd(a, b) << std::endl;

    // 素数判定 (Miller-Rabin)
    Int p("170141183460469231731687303715884105727");  // 2^127 - 1
    std::cout << p << " is prime: " << std::boolalpha
              << isProbablePrime(p) << std::endl;

    // NaN/Infinity 安全伝播
    Int zero;
    Int nan = Int(1) / zero;
    std::cout << "1/0 = " << nan << std::endl;  // NaN
}

Float — 多倍長浮動小数点

#include <math/core/mp/Float.hpp>
#include <iostream>
using namespace sangi;

int main() {
    int prec = 50;  // 50 桁精度
    Float::setDefaultPrecision(prec);

    // 数学定数 (スレッドローカルキャッシュ、2 回目以降は即座に返る)
    Float pi = Float::pi(prec);    // π を 50 桁で計算 (Chudnovsky)

    // 超越関数
    Float y = exp(pi, prec);       // e^π
    Float z = log(y, prec);        // log(e^π) = π

    std::cout << "exp(pi) = " << y.toDecimalString(prec) << std::endl;
    std::cout << "log(exp(pi)) = " << z.toDecimalString(prec) << std::endl;

    // 三角関数
    Float s = sin(pi / Float(6, prec), prec);  // sin(π/6) = 0.5
    std::cout << "sin(pi/6) = " << s.toDecimalString(20) << std::endl;
}

Rational — 有理数

#include <math/core/mp/Rational.hpp>
#include <iostream>
using namespace sangi;

int main() {
    Rational a(1, 3);   // 1/3
    Rational b(1, 6);   // 1/6

    Rational sum = a + b;
    std::cout << "1/3 + 1/6 = " << sum << std::endl;  // 1/2 (自動約分)

    Rational product = a * b;
    std::cout << "1/3 * 1/6 = " << product << std::endl;  // 1/18

    // 多倍長整数との組み合わせ
    Rational big(Int("1000000000000000000"), Int("3"));
    std::cout << big << std::endl;
}

Complex — 複素数

#include <math/core/Complex.hpp>
#include <iostream>
using namespace sangi;
using namespace sangi::literals;

int main() {
    Complex<double> z(3.0, 4.0);
    std::cout << "z     = " << z << std::endl;       // (3, 4)
    std::cout << "|z|   = " << abs(z) << std::endl;   // 5
    std::cout << "arg   = " << arg(z) << std::endl;   // 0.9273
    std::cout << "conj  = " << conj(z) << std::endl;  // (3, -4)

    // Euler's formula: e^(i*pi) = -1
    auto euler = exp(Complex<double>(0, std::numbers::pi));
    std::cout << "e^(i*pi) = " << euler << std::endl;  // (-1, 0)

    // User-defined literal
    auto w = 2.0 + 3.0_i;
    std::cout << "w*w = " << w * w << std::endl;      // (-5, 12)
}

Vector — 数値ベクトル

#include <math/core/vector.hpp>
#include <iostream>
using namespace sangi;

int main() {
    Vector<double> a{1.0, 2.0, 3.0};
    Vector<double> b{4.0, 5.0, 6.0};

    // 式テンプレート: 中間オブジェクトなしで評価
    Vector<double> c = 2.0 * a + b;
    std::cout << "dot(a, b) = " << dot(a, b) << std::endl;  // 32
    std::cout << "norm(a)   = " << norm(a) << std::endl;     // 3.74166

    // 3D 外積 (StaticVector)
    StaticVector<double, 3> i{1, 0, 0}, j{0, 1, 0};
    auto k = cross(i, j);  // [0, 0, 1]
}

Matrix — 行列

#include <math/core/matrix.hpp>
#include <math/core/mp/Rational.hpp>
#include <math/linalg/decomposition.hpp>
#include <iostream>
using namespace sangi;
using namespace sangi::algorithms;

int main() {
    // 行列演算
    Matrix<double> A({{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}});
    Vector<double> b({8, -11, -3});
    auto x = lu_solve(A, b);        // x = {2, 3, -1}

    // operator^: 転置, 逆行列, 冪乗
    auto At   = A ^ 'T';     // 転置
    auto Ainv = A ^ (-1);    // 逆行列
    auto A3   = A ^ 3;       // A*A*A

    // Rational で厳密な逆行列
    Matrix<Rational> H(3, 3);
    for (int i = 0; i < 3; i++)
        for (int j = 0; j < 3; j++)
            H(i, j) = Rational(1, i + j + 1);  // Hilbert 行列
    auto Hinv = lu_inverse(H);               // 全要素が整数 (厳密)
    auto I = H * Hinv;                       // 厳密に単位行列
}

LinAlg — 線形代数

#include <math/core/matrix.hpp>
#include <math/linalg/decomposition.hpp>
#include <math/linalg/eigenvalues.hpp>
#include <iostream>
using namespace sangi;
using namespace sangi::algorithms;

int main() {
    Matrix<double> A({{4, 1, 1}, {1, 3, 0}, {1, 0, 2}});

    // SVD 分解
    auto [U, sigma, Vt] = svd_decomposition(A);

    // 対称固有値
    auto [evals, evecs] = eigen_symmetric(A);

    // ガウス消去法 (ピボット戦略選択)
    Vector<double> b({6, 4, 3});
    auto x = gaussian_elimination(A, b, PivotStrategy::Full);

    // Bareiss: 整数行列の行列式 (除算なし)
    Matrix<int> M({{1, 2, 3}, {4, 5, 6}, {7, 8, 10}});
    auto det = bareiss_determinant(M);  // -3
}

Polynomial — 多項式

#include <math/core/polynomial.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // x^2 - 1 = (x-1)(x+1)
    Polynomial<double> p({-1, 0, 1});  // -1 + 0x + x^2
    Polynomial<double> q({1, 1});      //  1 + x

    auto [quot, rem] = p.divmod(q);
    std::cout << "quot = " << quot << std::endl;  // -1 + x
    std::cout << "rem  = " << rem << std::endl;   // 0

    // 微分・積分
    auto dp = p.derivative();    // 2x
    auto ip = p.integral();      // -x + x^3/3

    // 求根 (全複素根)
    auto roots = findAllRoots(p);  // {-1, 1}
}

FFT — 高速フーリエ変換

#include <math/fft/fft.hpp>
#include <math/fft/fft_utils.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // 複素 FFT (in-place)
    std::vector<Complex<double>> data = {1, 2, 3, 4};
    FFT<Complex<double>>::fft(data);   // 順変換
    FFT<Complex<double>>::ifft(data);  // 逆変換 (元に戻る)

    // 実数 FFT (N → N/2+1 の複素スペクトル)
    std::vector<double> signal = {1, 0, -1, 0, 1, 0, -1, 0};
    auto spectrum = RealFFT<double>::fft(signal);
    auto amp = amplitude_spectrum<double>(spectrum);

    // 多項式乗算 (畳み込み)
    std::vector<double> a = {1, 2, 3}, b = {4, 5};
    auto c = RealFFT<double>::convolve(a, b);  // {4, 13, 22, 15}
}

Roots — 求根アルゴリズム

#include <math/roots/roots.hpp>
#include <math/core/polynomial.hpp>
#include <iostream>
using namespace sangi;

int main() {
    // Brent 法: f(x) = x^3 - 2 の根 (= ∛2)
    auto f = [](double x) { return x * x * x - 2.0; };
    double root = brentRoot(f, 1.0, 2.0);
    std::cout << "cbrt(2) = " << root << std::endl;  // 1.25992...

    // Newton 法: f(x) = sin(x) - 0.5 の根
    auto g  = [](double x) { return std::sin(x) - 0.5; };
    auto dg = [](double x) { return std::cos(x); };
    double r2 = newtonRoot(g, dg, 0.5);
    std::cout << "arcsin(0.5) = " << r2 << std::endl;  // π/6

    // 多項式の全複素根 (Aberth-Ehrlich)
    Polynomial<double> p({6, -5, 1});  // x^2 - 5x + 6 = (x-2)(x-3)
    auto roots = findAllRoots(p);
    for (auto& z : roots)
        std::cout << z << std::endl;  // 2, 3
}

Optimization — 最適化

#include <math/optimization/optimization_nd.hpp>
#include <iostream>
using namespace sangi;

int main() {
    using V = Vector<double>;
    using M = Matrix<double>;

    // Rosenbrock 関数の BFGS 最小化
    auto f  = [](const V& x) { return (1-x[0])*(1-x[0]) + 100*(x[1]-x[0]*x[0])*(x[1]-x[0]*x[0]); };
    auto df = [](const V& x) {
        return V{-2*(1-x[0]) - 400*x[0]*(x[1]-x[0]*x[0]),
                 200*(x[1]-x[0]*x[0])};
    };
    auto result = bfgs_minimize<double, V, M>(f, df, V{-1, -1});
    std::cout << "min at " << result.point << std::endl;   // [1, 1]
    std::cout << "f(min) = " << result.value << std::endl;  // ~0
}