A production-quality, highly optimized mathematical library implementing algorithms for Number Theory, Algebra, and Probability & Statistics in modern C.
- Number Theory: GCD, LCM, Euler's Totient, Carmichael's Lambda, Miller-Rabin primality testing, modular arithmetic, prime factorization
- Algebra: Polynomial operations with Horner's method, matrix arithmetic, determinants, matrix inverse, linear equation solving
- Probability & Statistics: xoshiro256** RNG, 6 statistical distributions, descriptive statistics, special functions
Core algorithms for computational number theory (math_number_theory.c):
| Function | Description | Complexity |
|---|---|---|
| math_gcd() | Greatest Common Divisor (Euclidean) | O(log min(a,b)) |
| math_lcm() | Least Common Multiple | O(log min(a,b)) |
| math_extended_gcd() | Bézout coefficients for modular inverse | O(log min(a,b)) |
| math_euler_totient() | Euler's φ(n) - count of coprimes | O(√n) |
| math_carmichael_lambda() | Carmichael's λ(n) function | O(√n) |
| math_isqrt() | Integer square root (Newton's method) | O(log n) |
| math_mod_exp() | Modular exponentiation (square-and-multiply) | O(log exp) |
| math_mod_inverse() | Modular multiplicative inverse | O(log n) |
| math_is_prime_miller_rabin() | Probabilistic primality test | O(k log³ n) |
| math_factorize() | Prime factorization | O(√n) |
| math_divisor_count() | Number of divisors τ(n) | O(√n) |
| math_divisor_sum() | Sum of divisors σ(n) | O(√n) |
| math_mobius() | Möbius function μ(n) | O(√n) |
Key Optimizations:
- Bit operations for powers of 2
- __uint128_t for overflow-free multiplication when available
- Multiplicative properties for divisor functions
- Deterministic witnesses for Miller-Rabin primality testing
Polynomial and matrix operations (math_algebra.c):
Polynomial Operations:
- Evaluation using Horner's method - O(n) with minimal multiplications
- Addition, subtraction, multiplication
- Symbolic differentiation
Matrix Operations:
- Basic operations: add, subtract, multiply, transpose, scalar multiply
- Determinant via LU decomposition - O(n³)
- Matrix inverse via Gauss-Jordan elimination - O(n³)
- Linear system solver: Ax = b via Gaussian elimination with partial pivoting - O(n³)
- Trace computation
Key Optimizations:
- Horner's method reduces polynomial multiplications by 50%
- Cache-friendly matrix access patterns
- Partial pivoting for numerical stability
- LU decomposition for efficient determinant computation
Comprehensive statistical toolkit (math_probability.c):
Random Number Generation:
- xoshiro256** algorithm - period 2^256 - 1
- Uniform distribution with 53-bit precision
- Range-limited generation for doubles and integers
- Splitmix64 initialization for good state distribution
Statistical Distributions:
- Normal/Gaussian - Box-Muller transform
- Exponential - Inverse transform method
- Gamma - Marsaglia & Tsang method
- Bernoulli - Single trial
- Binomial - n trials with probability p
- Poisson - Knuth algorithm (small λ) + ratio of uniforms (large λ)
Descriptive Statistics:
- Mean, variance, standard deviation
- Median, percentiles (with interpolation)
- Skewness (third standardized moment)
- Kurtosis (excess kurtosis)
- Covariance, correlation (Pearson)
Probability Functions:
- Binomial PMF using log-gamma
- Poisson PMF using log-gamma
- Normal PDF and CDF
Special Functions:
- Gamma function - Lanczos approximation
- Log-gamma for numerical stability
- Beta function
- Error function (erf) - Abramowitz & Stegun approximation
- Complementary error function (erfc)
Combinatorics:
- Factorial (n! for n <= 20 with overflow protection)
- Binomial coefficients C(n,k) with symmetry optimization
- Permutations P(n,k)
#include "math_lib.h"
#include <stdio.h>
int main(void) {
// GCD and LCM
uint64_t gcd = math_gcd(48, 18);
printf("gcd(48, 18) = %lu\n", gcd); // Output: 6
// Euler's Totient Function
uint64_t phi = math_euler_totient(100);
printf("φ(100) = %lu\n", phi); // Output: 40
// Modular Exponentiation
uint64_t result = math_mod_exp(2, 100, 13);
printf("2^100 mod 13 = %lu\n", result); // Output: 4
// Primality Testing
bool is_prime = math_is_prime_miller_rabin(1009, 10);
printf("Is 1009 prime? %s\n", is_prime ? "Yes" : "No"); // Output: Yes
// Prime Factorization
factorization_t factors;
math_factorize(360, &factors);
printf("360 = ");
for (uint32_t i = 0; i < factors.count; i++) {
printf("%lu^%u ", factors.factors[i].prime, factors.factors[i].exponent);
}
printf("\n"); // Output: 360 = 2^3 3^2 5^1
return 0;
}#include "math_lib.h"
#include <stdio.h>
int main(void) {
// Polynomial: p(x) = 2 + 3x + x²
double coeffs[] = {2.0, 3.0, 1.0};
double value = math_poly_eval_horner(coeffs, 2, 2.0);
printf("p(2) = %.2f\n", value); // Output: 12.00
// Create a 2x2 matrix
matrix_t* A = math_matrix_create(2, 2);
math_matrix_set(A, 0, 0, 1.0);
math_matrix_set(A, 0, 1, 2.0);
math_matrix_set(A, 1, 0, 3.0);
math_matrix_set(A, 1, 1, 4.0);
// Calculate determinant
double det = math_matrix_determinant(A);
printf("det(A) = %.2f\n", det); // Output: -2.00
// Calculate inverse
matrix_t* A_inv = math_matrix_inverse(A);
if (A_inv) {
printf("A is invertible\n");
math_matrix_destroy(A_inv);
}
// Solve linear system Ax = b
double b[] = {5.0, 11.0};
double x[2];
linear_solve_status_t status = math_solve_linear(A, b, x);
if (status == LINEAR_SOLVE_OK) {
printf("Solution: x = %.2f, y = %.2f\n", x[0], x[1]);
}
math_matrix_destroy(A);
return 0;
}#include "math_lib.h"
#include <stdio.h>
int main(void) {
// Initialize RNG
rng_state_t rng;
math_rng_init(&rng, 42);
// Generate random numbers
printf("Uniform [0,1): %.4f\n", math_rng_uniform(&rng));
printf("Integer [1,6]: %ld\n", math_rng_int_range(&rng, 1, 6));
// Sample from distributions
double normal = math_dist_normal(&rng, 100.0, 15.0);
printf("Normal(100, 15): %.2f\n", normal);
uint32_t poisson = math_dist_poisson(&rng, 5.0);
printf("Poisson(5): %u\n", poisson);
// Statistical analysis
double data[] = {2.0, 4.0, 6.0, 8.0, 10.0};
size_t n = 5;
printf("Mean: %.2f\n", math_stats_mean(data, n));
printf("Std Dev: %.2f\n", math_stats_stddev(data, n));
// Special functions
double gamma = math_gamma_function(5.0);
printf("Γ(5) = %.2f\n", gamma); // Output: 24.00
// Combinatorics
uint64_t binom = math_binomial_coeff(10, 3);
printf("C(10,3) = %lu\n", binom); // Output: 120
return 0;
}// Basic number theory
uint64_t math_gcd(uint64_t a, uint64_t b);
uint64_t math_lcm(uint64_t a, uint64_t b);
int64_t math_extended_gcd(int64_t a, int64_t b, int64_t* x, int64_t* y);
// Arithmetic functions
uint64_t math_euler_totient(uint64_t n);
uint64_t math_carmichael_lambda(uint64_t n);
uint64_t math_isqrt(uint64_t n);
// Modular arithmetic
int64_t math_mod(int64_t x, int64_t m);
uint64_t math_mod_add(uint64_t a, uint64_t b, uint64_t m);
uint64_t math_mod_sub(uint64_t a, uint64_t b, uint64_t m);
uint64_t math_mod_mul(uint64_t a, uint64_t b, uint64_t m);
uint64_t math_mod_exp(uint64_t base, uint64_t exp, uint64_t m);
uint64_t math_mod_inverse(uint64_t a, uint64_t m);
// Primality and factorization
bool math_is_prime_simple(uint64_t n);
bool math_is_prime_miller_rabin(uint64_t n, uint32_t iterations);
void math_factorize(uint64_t n, factorization_t* result);
// Divisor functions
uint64_t math_divisor_count(uint64_t n);
uint64_t math_divisor_sum(uint64_t n);
int32_t math_mobius(uint64_t n);// Polynomial operations
polynomial_t* math_poly_create(uint32_t degree);
void math_poly_destroy(polynomial_t* poly);
double math_poly_eval(const polynomial_t* poly, double x);
double math_poly_eval_horner(const double* coeffs, uint32_t degree, double x);
polynomial_t* math_poly_add(const polynomial_t* a, const polynomial_t* b);
polynomial_t* math_poly_sub(const polynomial_t* a, const polynomial_t* b);
polynomial_t* math_poly_mul(const polynomial_t* a, const polynomial_t* b);
polynomial_t* math_poly_derivative(const polynomial_t* poly);
// Matrix operations
matrix_t* math_matrix_create(uint32_t rows, uint32_t cols);
void math_matrix_destroy(matrix_t* mat);
double math_matrix_get(const matrix_t* mat, uint32_t row, uint32_t col);
void math_matrix_set(matrix_t* mat, uint32_t row, uint32_t col, double value);
matrix_t* math_matrix_add(const matrix_t* a, const matrix_t* b);
matrix_t* math_matrix_sub(const matrix_t* a, const matrix_t* b);
matrix_t* math_matrix_mul(const matrix_t* a, const matrix_t* b);
matrix_t* math_matrix_scalar_mul(const matrix_t* mat, double scalar);
matrix_t* math_matrix_transpose(const matrix_t* mat);
double math_matrix_determinant(const matrix_t* mat);
double math_matrix_trace(const matrix_t* mat);
matrix_t* math_matrix_inverse(const matrix_t* mat);
linear_solve_status_t math_solve_linear(const matrix_t* A, const double* b, double* x);// Random number generation
void math_rng_init(rng_state_t* rng, uint64_t seed);
uint64_t math_rng_next(rng_state_t* rng);
double math_rng_uniform(rng_state_t* rng);
double math_rng_uniform_range(rng_state_t* rng, double min, double max);
int64_t math_rng_int_range(rng_state_t* rng, int64_t min, int64_t max);
// Statistical distributions
double math_dist_normal(rng_state_t* rng, double mean, double stddev);
double math_dist_exponential(rng_state_t* rng, double lambda);
double math_dist_gamma(rng_state_t* rng, double shape, double scale);
bool math_dist_bernoulli(rng_state_t* rng, double p);
uint32_t math_dist_binomial(rng_state_t* rng, uint32_t n, double p);
uint32_t math_dist_poisson(rng_state_t* rng, double lambda);
// Descriptive statistics
double math_stats_mean(const double* data, size_t n);
double math_stats_variance(const double* data, size_t n);
double math_stats_stddev(const double* data, size_t n);
double math_stats_median(double* data, size_t n);
double math_stats_percentile(double* data, size_t n, double p);
double math_stats_skewness(const double* data, size_t n);
double math_stats_kurtosis(const double* data, size_t n);
double math_stats_covariance(const double* x, const double* y, size_t n);
double math_stats_correlation(const double* x, const double* y, size_t n);
// Probability functions
double math_pmf_binomial(uint32_t k, uint32_t n, double p);
double math_pmf_poisson(uint32_t k, double lambda);
double math_pdf_normal(double x, double mean, double stddev);
double math_cdf_normal(double x, double mean, double stddev);
// Special functions
double math_gamma_function(double x);
double math_log_gamma(double x);
double math_beta_function(double a, double b);
double math_erf(double x);
double math_erfc(double x);
// Combinatorics
uint64_t math_factorial(uint32_t n);
uint64_t math_binomial_coeff(uint32_t n, uint32_t k);
uint64_t math_permutation(uint32_t n, uint32_t k);| Operation | Time Complexity | Space Complexity |
|---|---|---|
| GCD (Euclidean) | O(log min(a,b)) | O(1) |
| Euler's φ(n) | O(√n) | O(1) |
| Modular exponentiation | O(log exp) | O(1) |
| Miller-Rabin primality | O(k log³ n) | O(1) |
| Prime factorization | O(√n) | O(log n) |
| Polynomial evaluation | O(degree) | O(1) |
| Matrix multiplication | O(n³) | O(n²) |
| Matrix determinant | O(n³) | O(n²) |
| Linear system solve | O(n³) | O(n²) |
| RNG (xoshiro256**) | O(1) | O(1) |
| Normal distribution | O(1) amortized | O(1) |
- Bit Operations: Powers of 2 handled with shifts and masks
- Multiplicative Properties: Divisor functions computed from prime factorization
- Horner's Method: Polynomial evaluation with minimal multiplications
- Cache-Friendly: Matrix operations optimized for memory access patterns
- Early Termination: Algorithms exit as soon as result is determined
- Overflow Protection: Safe arithmetic for modular operations
- Knuth, "The Art of Computer Programming, Vol. 2: Seminumerical Algorithms"
- Bach & Shallit, "Algorithmic Number Theory"
- Golub & Van Loan, "Matrix Computations"
- Vigna, "It Is High Time We Let Go of the Mersenne Twister" (xoshiro256**)
- Marsaglia & Tsang, "A Fast, Easily Implemented Method for Sampling from Decreasing or Symmetric Unimodal Density Functions"
- Box & Muller, "A Note on the Generation of Random Normal Deviates"
- Lanczos, "A Precision Approximation of the Gamma Function"
- Abramowitz & Stegun, "Handbook of Mathematical Functions"
H.Overman
All rights reserved. This software is provided for educational and research purposes.
Email: opsec.ee@pm.me
For bug reports, feature requests, or questions about the library.