SIGMA — Statistical Inference, Generation, Moments & Analysis

Zero-dependency, header-only, SIMD-accelerated statistical computing in C11.

     ____ ___ ____ __  __    _
    / ___|_ _/ ___|  \/  |  / \
    \___ \| | |  _| |\/| | / _ \
     ___) | | |_| | |  | |/ ___ \
    |____/___\____|_|  |_/_/   \_\

SIGMA is a header-only C11 library for high-performance statistical computing, featuring probability distribution functions (PDF, CDF, quantile, random), a pluggable SIMD-accelerated PRNG, summary statistics with online algorithms, and zero external dependencies.

Features

• 15+ Probability Distributions - Normal, Exponential, Gamma, Chi-squared, Beta, Cauchy, Logistic, Laplace, Log-normal, Weibull, Uniform, Student's t, F, Binomial, Poisson
• Four functions per distribution - d (density), p (CDF), q (quantile), r (random)
• SIMD-Accelerated PRNG - xoshiro256** with AVX2 (4-wide) and NEON (2-wide) parallel streams
• Pluggable by design - Swap memory allocators and PRNGs without changing code
• HPC Best Practices - restrict, alignas, branch-free selection, prefetch hints
• Online Statistics - Welford's algorithm, mergeable moments (parallel reduction)
• Buffer-fill APIs - Fill entire arrays with random variates in one call
• Zero dependencies - Just C11 and <math.h>

Quick Start

#include <sigma/sigma.h>
#include <stdio.h>

int main(void) {
    // Initialize RNG with seed 42
    sigma_state_t rng;
    sigma_rng_init(&rng, 42);

    // Single values
    double x = sigma_rnorm_d(&rng, 0.0, 1.0);   // N(0,1)
    double p = sigma_pnorm_d(1.96, 0.0, 1.0);   // ~0.975
    double q = sigma_qnorm_d(0.975, 0.0, 1.0);  // ~1.96

    // Fill buffer with 1M normals
    double buf[1024 * 1024];
    sigma_rnorm_fill_d(&rng, 0.0, 1.0, buf, 1024 * 1024);

    return 0;
}

Building

SIGMA is header-only. Just add the include/ directory to your compiler's include path:

gcc -I./include -O3 -mavx2 -mfma myprogram.c -o myprogram -lm

CMake

mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release ..
cmake --build . --parallel

Build Options

Option	Default	Description
SIGMA_BUILD_TESTS	ON	Build test suite
SIGMA_BUILD_BENCHMARKS	ON	Build benchmarks
SIGMA_BUILD_EXAMPLES	ON	Build example programs
SIGMA_ENABLE_AVX2	ON	Enable AVX2 SIMD (x86_64)
SIGMA_ENABLE_NEON	OFF	Enable ARM NEON
SIGMA_ENABLE_OPENMP	OFF	Enable OpenMP parallelization

Distributions

Distribution	PDF	CDF	Quantile	Random	Fill
Normal (norm)	dnorm	pnorm	qnorm	rnorm	rnorm_fill
Exponential (exp)	dexp	pexp	qexp	rexp	rexp_fill
Gamma (gamma)	dgamma	pgamma	qgamma	rgamma	rgamma_fill
Chi-squared (chisq)	dchisq	pchisq	qchisq	rchisq	rchisq_fill
Beta (beta)	dbeta	pbeta	qbeta	rbeta	rbeta_fill
Cauchy (cauchy)	dcauchy	pcauchy	qcauchy	rcauchy	rcauchy_fill
Logistic (logis)	dlogis	plogis	qlogis	rlogis	rlogis_fill
Laplace (laplace)	dlaplace	plaplace	qlaplace	rlaplace	rlaplace_fill
Log-normal (lnorm)	dlnorm	plnorm	qlnorm	rlnorm	rlnorm_fill
Weibull (weibull)	dweibull	pweibull	qweibull	rweibull	rweibull_fill
Uniform (unif)	dunif	punif	qunif	runif	runif_fill
Student's t (t)	dt	pt	qt	rt	rt_fill
F (f)	df	pf	qf	rf	rf_fill
Binomial (binom)	dbinom	pbinom	qbinom	rbinom	-
Poisson (pois)	dpois	ppois	qpois	rpois	-

Naming follows the R convention: sigma_{d,p,q,r}<name>_d(...)

PRNG (xoshiro256**)

The default PRNG is xoshiro256** (Blackman & Vigna, 2018), a fast high-quality generator with:

• Period: 2^256 - 1
• Jump ahead: 2^128 steps via sigma_rng_jump()
• Seeding: splitmix64 for high-quality state initialization
• SIMD: 4-wide AVX2, 2-wide NEON (aarch64)

Seeding

sigma_state_t rng;
sigma_rng_init(&rng, 42);           // Seed with 42
sigma_rng_init(&rng, time(NULL));   // Random seed

Scalar vs Buffer

// Scalar - one at a time
uint64_t u = sigma_rng_next_u64(&rng);
double   d = sigma_rng_next_d(&rng);     // [0,1)
double   n = sigma_rng_next_normal(&rng); // N(0,1)

// Buffer fill - SIMD accelerated
double buf[1000000];
sigma_rng_fill_d(&rng, buf, 1000000);       // Uniform [0,1)
sigma_rng_fill_normal(&rng, buf, 1000000);  // Standard normal

Pluggable Allocator

Pass a custom allocator for zero-allocation hot paths:

// Stack arena - no malloc!
unsigned char arena_buf[16 * 1024 * 1024];
sigma_arena_t arena = SIGMA_ARENA_INIT(arena_buf);
sigma_allocator_t alloc = sigma_arena_wrap(&arena);

// Use with stats computation
sigma_full_stats_t st;
sigma_stats_compute_d(data, n, &st, &alloc);

// Reset arena instantly
sigma_arena_reset(&arena);

Summary Statistics

// Batch statistics
double mean = sigma_mean_d(data, n);
double var  = sigma_var_d(data, n);    // Sample variance
double std  = sigma_std_d(data, n);
double skew = sigma_skew_d(data, n);
double kurt = sigma_kurt_d(data, n);

// Online (single-pass, mergeable for parallel reduction)
sigma_online_stats_t s;
sigma_online_init(&s);
for (int i = 0; i < n; i++) sigma_online_push(&s, data[i]);
double v = sigma_online_var(&s);

// Merge two online stats (for parallel map-reduce)
sigma_online_merge(&stats_a, &stats_b);

// Full report
sigma_full_stats_t st;
sigma_stats_compute_d(data, n, &st, NULL);
// st.mean, st.median, st.std, st.skew, st.kurt, st.min, st.max

Examples

Example	Description
ex_basic	PDF/CDF/quantile for standard normal
ex_rng	Random variate generation, buffer fills
ex_stats	Summary statistics with SIMD
ex_monte_carlo	European option pricing via Monte Carlo
ex_bayesian	Beta-Bernoulli coin flip inference
ex_custom_allocator	Arena allocator for zero-allocation paths

Run examples:

cd build
./examples/ex_monte_carlo
./examples/ex_bayesian

Testing

cd build
ctest --output-on-failure

The test suite covers:

• Known-value tests for all distributions
• Quantile roundtrip consistency (q(p(x)) = x)
• Parameter validation (NaN for invalid inputs)
• Edge cases (boundary probabilities, extreme parameters)
• RNG reproducibility and statistical properties
• Online vs batch statistic consistency
• Large-mean cancellation regression tests
• Allocator correctness

Benchmarks

cd build
./benchmarks/sigma_bench

Architecture

include/sigma/
    sigma.h       - Umbrella header (include this)
    core.h        - Platform detection, SIMD macros, types
    alloc.h       - Pluggable allocator (arena, pool, default)
    rng.h         - xoshiro256** PRNG with SIMD fill
    simd_math.h   - Vectorized exp/log for AVX2/NEON
    special.h     - Gamma, beta, erf, incomplete functions
    dists.h       - All distribution families (d/p/q/r)
    stats.h       - Summary statistics with online algorithms

Performance Notes

• SIMD width: AVX2 processes 4 doubles in parallel; NEON (aarch64) processes 2
• Best throughput: Use buffer-fill APIs (*_fill_*) rather than scalar loops
• Memory alignment: aligned_alloc(32, ...) gives ~10-15% speedup on AVX2
• Arena allocators: Eliminate malloc overhead in tight loops
• Online stats: Single-pass, numerically stable, mergeable for parallel computation

Porting from R / statslib

R / statslib	SIGMA
dnorm(x, mu, sigma)	sigma_dnorm_d(x, mu, sigma)
pnorm(x, mu, sigma)	sigma_pnorm_d(x, mu, sigma)
qnorm(p, mu, sigma)	sigma_qnorm_d(p, mu, sigma)
rnorm(n, mu, sigma)	sigma_rnorm_fill_d(&rng, mu, sigma, buf, n)
mean(x)	sigma_mean_d(x, n)
var(x)	sigma_var_d(x, n)
sd(x)	sigma_std_d(x, n)

License

MIT License. See LICENSE for details.

Acknowledgments

• Blackman & Vigna for xoshiro256**
• Marsaglia for the polar method (normal generation)
• Marsaglia & Tsang for the gamma sampler
• Acklam for the rational approximation inverse normal CDF
• Keith O'Hara's statslib for the API design inspiration