sangi

Examples

Int — Arbitrary-Precision Integer

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

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

    // NaN/Infinity safe propagation
    Int zero;
    Int nan = Int(1) / zero;
    std::cout << "1/0 = " << nan << std::endl;  // NaN
}

Float — Arbitrary-Precision Floating-Point

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

int main() {
    int prec = 50;  // 50-digit precision
    Float::setDefaultPrecision(prec);

    // Mathematical constants (thread-local cached; instant on subsequent calls)
    Float pi = Float::pi(prec);    // π to 50 digits (Chudnovsky)

    // Transcendental functions
    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;

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

Rational — Arbitrary-Precision Rational Number

#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 (auto-reduced)

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

    // Combine with arbitrary-precision integers
    Rational big(Int("1000000000000000000"), Int("3"));
    std::cout << big << std::endl;
}

Complex — Complex Numbers

#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 — Numeric 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};

    // Expression templates: no temporaries until assignment
    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 cross product (StaticVector)
    StaticVector<double, 3> i{1, 0, 0}, j{0, 1, 0};
    auto k = cross(i, j);  // [0, 0, 1]
}

Matrix — Dense 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 operations
    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^: transpose, inverse, power
    auto At   = A ^ 'T';     // transpose
    auto Ainv = A ^ (-1);    // inverse
    auto A3   = A ^ 3;       // A*A*A

    // Exact inverse with 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 matrix
    auto Hinv = lu_inverse(H);               // all integer entries (exact)
    auto I = H * Hinv;                       // exactly identity
}

LinAlg — Linear Algebra

#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 decomposition
    auto [U, sigma, Vt] = svd_decomposition(A);

    // Symmetric eigenvalues
    auto [evals, evecs] = eigen_symmetric(A);

    // Gaussian elimination (pivot strategy selection)
    Vector<double> b({6, 4, 3});
    auto x = gaussian_elimination(A, b, PivotStrategy::Full);

    // Bareiss: integer determinant (division-free)
    Matrix<int> M({{1, 2, 3}, {4, 5, 6}, {7, 8, 10}});
    auto det = bareiss_determinant(M);  // -3
}

Polynomial — 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

    // Differentiation and integration
    auto dp = p.derivative();    // 2x
    auto ip = p.integral();      // -x + x^3/3

    // Root finding (all complex roots)
    auto roots = findAllRoots(p);  // {-1, 1}
}

FFT — Fast Fourier Transform

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

int main() {
    // Complex FFT (in-place)
    std::vector<Complex<double>> data = {1, 2, 3, 4};
    FFT<Complex<double>>::fft(data);   // forward
    FFT<Complex<double>>::ifft(data);  // inverse (recovers original)

    // Real FFT (N -> N/2+1 complex spectrum)
    std::vector<double> signal = {1, 0, -1, 0, 1, 0, -1, 0};
    auto spectrum = RealFFT<double>::fft(signal);
    auto amp = amplitude_spectrum<double>(spectrum);

    // Polynomial multiplication (convolution)
    std::vector<double> a = {1, 2, 3}, b = {4, 5};
    auto c = RealFFT<double>::convolve(a, b);  // {4, 13, 22, 15}
}

Roots — Root-Finding Algorithms

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

int main() {
    // Brent's method: root of 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's method: root of 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

    // All complex roots of a polynomial (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 — Optimization

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

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

    // BFGS minimization of the Rosenbrock function
    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
}