diff --git a/lib/std/math/distributions.c3 b/lib/std/math/distributions.c3
new file mode 100644
index 000000000..62963f02d
--- /dev/null
+++ b/lib/std/math/distributions.c3
@@ -0,0 +1,589 @@
+// Copyright (c) 2026 Koni Marti. All rights reserved.
+// Use of this source code is governed by the MIT license.
+<*
+ This module provides a comprehensive set of continuous and discrete
+ probability distributions with support for PDF, CDF, inverse CDF,
+ random sampling, mean, and variance calculations.
+*>
+module std::math::distributions;
+import std::math::random;
+
+// Distribution interface defining common statistical operations
+interface Distribution
+{
+	<* Calculate the mean (expected value) of the distribution *>
+	fn double mean();
+
+	<* Calculate the variance of the distribution *>
+	fn double variance();
+}
+
+interface ContinuousDistribution : Distribution
+{
+	 <* Probability density function (PDF) *>
+	fn double pdf(double x);
+
+	<* Cumulative distribution function (CDF) *>
+	fn double cdf(double x);
+
+	<* Inverse cumulative distribution function (quantile function) *>
+	fn double quantile(double p);
+
+	<* Generate a random sample from the distribution *>
+	fn double random(Random rand);
+}
+
+interface DiscreteDistribution : Distribution
+{
+	<* Probability mass function (PMF) *>
+	fn double pmf(int k);
+
+	<* Cumulative distribution function (CDF) *>
+	fn double cdf(int k);
+
+	<* Inverse cumulative distribution function *>
+	fn int quantile(double p);
+
+	<* Generate a random sample from the distribution *>
+	fn int random(Random rand);
+}
+
+<*
+ Uniform distribution over [a, b]
+ *>
+struct UniformDist (ContinuousDistribution)
+{
+	double a;
+	double b;
+}
+
+<*
+ @require b > a : "Upper bound must be greater than lower bound."
+*>
+fn UniformDist uniform(double a, double b)
+{
+	return (UniformDist){ a, b };
+}
+
+fn double UniformDist.mean(&self) @dynamic
+{
+	return (self.a + self.b) / 2.0;
+}
+
+fn double UniformDist.variance(&self) @dynamic
+{
+	double range = self.b - self.a;
+	return range * range / 12.0;
+}
+
+fn double UniformDist.pdf(&self, double x) @dynamic
+{
+	if (x < self.a || x > self.b) return 0;
+	return 1.0 / (self.b - self.a);
+}
+
+fn double UniformDist.cdf(&self, double x) @dynamic
+{
+	if (x < self.a) return 0.0;
+	if (x > self.b) return 1.0;
+	return (x - self.a) / (self.b - self.a);
+}
+
+<*
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+*>
+fn double UniformDist.quantile(&self, double p) @dynamic
+{
+	return self.a + p * (self.b - self.a);
+}
+
+fn double UniformDist.random(&self, Random rand) @dynamic
+{
+	return self.a + random::next_double(rand) * (self.b - self.a);
+}
+
+<*
+ Normal (Gaussian) distribution
+ *>
+struct NormalDist (ContinuousDistribution)
+{
+	double mu;
+	double sigma;
+}
+
+<*
+ @require sigma > 0.0 : "Standard deviation must be positive"
+*>
+fn NormalDist normal(double mu = 0.0, double sigma = 1.0)
+{
+	return (NormalDist){ mu, sigma };
+}
+
+fn double NormalDist.mean(&self) @dynamic
+{
+	return self.mu;
+}
+
+fn double NormalDist.variance(&self) @dynamic
+{
+	return self.sigma * self.sigma;
+}
+
+fn double NormalDist.pdf(&self, double x) @dynamic
+{
+	double z = (x - self.mu) / self.sigma;
+	return math::exp(-0.5 * z * z) / (self.sigma * math::sqrt(2.0 * math::PI));
+}
+
+fn double NormalDist.cdf(&self, double x) @dynamic
+{
+	double z = (x - self.mu) / self.sigma;
+	return math::clamp(0.5 * (1.0 + math::erf(z / math::SQRT2)), 0.0, 1.0);
+}
+
+<*
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+*>
+fn double NormalDist.quantile(&self, double p) @dynamic
+{
+	double z = inverse_erf(2.0 * p - 1.0) * math::SQRT2;
+	return self.mu + self.sigma * z;
+}
+
+fn double NormalDist.random(&self, Random rand) @dynamic
+{
+	// Box-Muller transform.
+	double u1 = random::next_double(rand);
+	double u2 = random::next_double(rand);
+	double z = math::sqrt(-2.0 * math::ln(u1)) * math::cos(2.0 * math::PI * u2);
+	return self.mu + self.sigma * z;
+}
+
+<*
+ Exponential distribution
+ *>
+struct ExponentialDist (ContinuousDistribution)
+{
+	double lambda;
+}
+
+<*
+ @require lambda > 0.0 : "Rate parameter must be positive."
+*>
+fn ExponentialDist exponential(double lambda = 1.0)
+{
+	return (ExponentialDist){ lambda };
+}
+
+<*
+ @require self.lambda > 0.0 : "Rate parameter must be positive."
+*>
+fn double ExponentialDist.mean(&self) @dynamic
+{
+	return 1.0 / self.lambda;
+}
+
+<*
+ @require self.lambda > 0.0 : "Rate parameter must be positive."
+*>
+fn double ExponentialDist.variance(&self) @dynamic
+{
+	return 1.0 / (self.lambda * self.lambda);
+}
+
+fn double ExponentialDist.pdf(&self, double x) @dynamic
+{
+	if (x < 0.0) return 0.0;
+	return self.lambda * math::exp(-self.lambda * x);
+}
+
+fn double ExponentialDist.cdf(&self, double x) @dynamic
+{
+	if (x < 0.0) return 0.0;
+	return math::clamp(1.0 - math::exp(-self.lambda * x), 0.0, 1.0);
+}
+
+<*
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+ @require self.lambda > 0.0 : "Rate parameter must be positive."
+*>
+fn double ExponentialDist.quantile(&self, double p) @dynamic
+{
+	return -math::ln(1.0 - p) / self.lambda;
+}
+
+<*
+ @require self.lambda > 0.0 : "Rate parameter must be positive."
+*>
+fn double ExponentialDist.random(&self, Random rand) @dynamic
+{
+	return -math::ln(1.0 - random::next_double(rand)) / self.lambda;
+}
+
+<*
+ Student's t-distribution
+ *>
+struct TDist (ContinuousDistribution)
+{
+	double df;
+}
+
+<*
+ @require df > 0.0 : "Degrees of freedom must be positive."
+*>
+fn TDist t_distribution(double df)
+{
+	return (TDist){ df };
+}
+
+fn double TDist.mean(&self) @dynamic
+{
+	if (self.df <= 1.0) return double.nan;
+	return 0.0;
+}
+
+fn double TDist.variance(&self) @dynamic
+{
+	if (self.df <= 1.0) return double.nan;
+	if (self.df <= 2.0) return double.inf;
+	return self.df / (self.df - 2.0);
+}
+
+fn double TDist.pdf(&self, double x) @dynamic
+{
+	double v = self.df;
+	double coef = math::tgamma((v + 1.0) / 2.0) /
+		(math::sqrt(v * math::PI) * math::tgamma(v / 2.0));
+	return coef * math::pow(1.0 + x * x / v, -(v + 1.0) / 2.0);
+}
+
+fn double TDist.cdf(&self, double x) @dynamic
+{
+	double v = self.df;
+	if (x == 0.0) return 0.5;
+
+	double t = v / (v + x * x);
+	double a = v / 2.0;
+	double b = 0.5;
+
+	// Using regularized incomplete beta function.
+	double beta_cdf = incomplete_beta(t, a, b);
+
+	double p = x >= 0.0 ? 1.0 - 0.5 * beta_cdf : 0.5 * beta_cdf;
+	return math::clamp(p, 0.0, 1.0);
+}
+
+<*
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+*>
+fn double TDist.quantile(&self, double p) @dynamic
+{
+	if (p == 0.5) return 0.0;
+	double x = (p < 0.5) ? -1.0 : 1.0;
+	return newton_raphson(self, x, p) ?? double.nan;
+}
+
+fn double TDist.random(&self, Random rand) @dynamic
+{
+	// Generate using relationship with normal and chi-squared
+	NormalDist std_normal = normal(0.0, 1.0);
+	double z = std_normal.random(rand);
+	double v = chi_squared_sample(self.df, rand);
+	return z / math::sqrt(v / self.df);
+}
+
+<*
+ F-distribution
+*>
+struct FDist (ContinuousDistribution)
+{
+	double d1;
+	double d2;
+}
+
+<*
+ @require d1 > 0.0 && d2 > 0.0 : "Degrees of freedom must be positive."
+*>
+fn FDist f_distribution(double d1, double d2)
+{
+	return (FDist){ d1, d2 };
+}
+
+fn double FDist.mean(&self) @dynamic
+{
+	if (self.d2 <= 2.0) return double.nan;
+	return self.d2 / (self.d2 - 2.0);
+}
+
+fn double FDist.variance(&self) @dynamic
+{
+	if (self.d2 <= 4.0) return double.nan;
+	double d1 = self.d1;
+	double d2 = self.d2;
+	return 2.0 * d2 * d2 * (d1 + d2 - 2.0) /
+		(d1 * (d2 - 2.0) * (d2 - 2.0) * (d2 - 4.0));
+}
+
+fn double FDist.pdf(&self, double x) @dynamic
+{
+	if (x < 0.0) return 0.0;
+
+	double d1 = self.d1;
+	double d2 = self.d2;
+
+	double num = math::pow(d1 * x, d1) * math::pow(d2, d2);
+	double denom = math::pow(d1 * x + d2, d1 + d2);
+	double beta_term = x * beta_function(d1 / 2.0, d2 / 2.0);
+
+	return math::sqrt(num / denom) / beta_term;
+}
+
+
+fn double FDist.cdf(&self, double x) @dynamic
+{
+	if (x <= 0.0) return 0.0;
+
+	double d1 = self.d1;
+	double d2 = self.d2;
+	double t = d1 * x / (d1 * x + d2);
+
+	double p = incomplete_beta(t, d1 / 2.0, d2 / 2.0);
+	return math::clamp(p, 0.0, 1.0);
+}
+
+<*
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+*>
+fn double FDist.quantile(&self, double p) @dynamic
+{
+	return find_quantile(self, 0.0, 1000.0, p);
+}
+
+fn double FDist.random(&self, Random rand) @dynamic
+{
+	// Generate using ratio of chi-squared variables.
+	double u1 = chi_squared_sample(self.d1, rand);
+	double u2 = chi_squared_sample(self.d2, rand);
+	return (u1 / self.d1) / (u2 / self.d2);
+}
+
+<*
+ Chi-squared distribution
+ *>
+struct ChiSquaredDist (ContinuousDistribution)
+{
+	double k;
+}
+
+<*
+ @require k > 0.0 : "Degrees of freedom must be positive"
+*>
+fn ChiSquaredDist chi_squared(double k)
+{
+	return (ChiSquaredDist){ k };
+}
+
+fn double ChiSquaredDist.mean(&self) @dynamic
+{
+	return self.k;
+}
+
+fn double ChiSquaredDist.variance(&self) @dynamic
+{
+	return 2.0 * self.k;
+}
+
+fn double ChiSquaredDist.pdf(&self, double x) @dynamic
+{
+	if (x < 0.0) return 0.0;
+	if (x == 0.0 && self.k < 2.0) return double.inf;
+	if (x == 0.0) return 0.0;
+
+	double k = self.k;
+	return math::pow(x, k / 2.0 - 1.0) * math::exp(-x / 2.0) /
+		   (math::pow(2.0, k / 2.0) * math::tgamma(k / 2.0));
+}
+
+fn double ChiSquaredDist.cdf(&self, double x) @dynamic
+{
+	if (x <= 0.0) return 0.0;
+	double p = lower_incomplete_gamma(self.k / 2.0, x / 2.0);
+	return math::clamp(p, 0.0, 1.0);
+}
+
+<*
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+*>
+fn double ChiSquaredDist.quantile(&self, double p) @dynamic
+{
+	double low = 0.0;
+	double high = self.k + 10.0 * math::sqrt(2.0 * self.k);
+	return find_quantile(self, low, high, p);
+}
+
+fn double ChiSquaredDist.random(&self, Random rand) @dynamic
+{
+	return chi_squared_sample(self.k, rand);
+}
+
+<*
+ Binomial distribution
+ *>
+struct BinomialDist (DiscreteDistribution)
+{
+	int n;
+	double p;
+}
+
+<*
+ @require n >= 0 : "Number of trials must be non-negative."
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+*>
+fn BinomialDist binomial(int n, double p)
+{
+	return (BinomialDist){ n, p };
+}
+
+fn double BinomialDist.mean(&self) @dynamic
+{
+	return (double)self.n * self.p;
+}
+
+fn double BinomialDist.variance(&self) @dynamic
+{
+	return (double)self.n * self.p * (1.0 - self.p);
+}
+
+fn double BinomialDist.pmf(&self, int k) @dynamic
+{
+	if (k < 0 || k > self.n) return 0.0;
+	return binomial_coefficient(self.n, k) *
+		math::pow(self.p, (double)k) *
+		math::pow(1.0 - self.p, (double)(self.n - k));
+}
+
+fn double BinomialDist.cdf(&self, int k) @dynamic
+{
+	if (k < 0) return 0.0;
+	if (k >= self.n) return 1.0;
+
+	double sum = 0.0;
+	for (int i = 0; i <= k; i++)
+	{
+		sum += self.pmf(i);
+	}
+	return sum;
+}
+
+<*
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+*>
+fn int BinomialDist.quantile(&self, double p) @dynamic
+{
+	double cumulative = 0.0;
+	for (int k = 0; k <= self.n; k++)
+	{
+		cumulative += self.pmf(k);
+		if (cumulative >= p) return k;
+	}
+	return self.n;
+}
+
+fn int BinomialDist.random(&self, Random rand) @dynamic
+{
+	// Generate using Bernoulli trials.
+	int successes = 0;
+	for (int i = 0; i < self.n; i++)
+	{
+		if (random::next_double(rand) < self.p)
+		{
+			successes++;
+		}
+	}
+	return successes;
+}
+
+<*
+ Poisson distribution
+ *>
+struct PoissonDist (DiscreteDistribution)
+{
+	double lambda;
+}
+
+<*
+ @require lambda > 0.0 : "Rate parameter must be positive."
+*>
+fn PoissonDist poisson(double lambda)
+{
+	return (PoissonDist){ lambda };
+}
+
+fn double PoissonDist.mean(&self) @dynamic
+{
+	return self.lambda;
+}
+
+fn double PoissonDist.variance(&self) @dynamic
+{
+	return self.lambda;
+}
+
+fn double PoissonDist.pmf(&self, int k) @dynamic
+{
+	if (k < 0) return 0.0;
+	return math::exp(-self.lambda + (double)k * math::ln(self.lambda) - ln_factorial(k));
+}
+
+fn double PoissonDist.cdf(&self, int k) @dynamic
+{
+	if (k < 0) return 0.0;
+
+	double sum = 0.0;
+	for (int i = 0; i <= k; i++)
+	{
+		sum += self.pmf(i);
+	}
+	return sum;
+}
+
+<*
+ @require p >= 0.0 && p <= 1.0 : "Probability must be between 0 and 1."
+*>
+fn int PoissonDist.quantile(&self, double p) @dynamic
+{
+	double cumulative = 0.0;
+	int k = 0;
+	while (cumulative < p)
+	{
+		cumulative += self.pmf(k);
+		if (cumulative >= p) return k;
+		k++;
+		if (k > 1_000_000) break; // Safety limit
+	}
+	return k;
+}
+
+fn int PoissonDist.random(&self, Random rand) @dynamic
+{
+	// Knuth's algorithm for small lambda.
+	if (self.lambda < 30.0)
+	{
+		double l = math::exp(-self.lambda);
+		int k = 0;
+		double p = 1.0;
+		do
+		{
+			k++;
+			p *= random::next_double(rand);
+		} while (p > l);
+		return (k - 1);
+	}
+	else
+	{
+		// Use normal approximation for large lambda
+		NormalDist approx = normal(self.lambda, math::sqrt(self.lambda));
+		return (int)math::max(0.0, math::round(approx.random(rand)));
+	}
+}
+
diff --git a/lib/std/math/distributions_private.c3 b/lib/std/math/distributions_private.c3
new file mode 100644
index 000000000..904ab7d08
--- /dev/null
+++ b/lib/std/math/distributions_private.c3
@@ -0,0 +1,407 @@
+// Copyright (c) 2026 Koni Marti. All rights reserved.
+// Use of this source code is governed by the MIT license.
+module std::math::distributions @private;
+import std::math::random;
+
+struct ConvergenceControl
+{
+	usz max_iter;
+	double epsilon;
+}
+
+const DEFAULT_CONV = (ConvergenceControl){ 600, 1e-12 };
+const RELAXED_CONV = (ConvergenceControl){ 600, 1e-6 };
+
+faultdef NOT_CONVERGED;
+
+<*
+ Compute binomial coefficient C(n, k)
+*>
+fn double binomial_coefficient(int n, int k)
+{
+	if (k < 0 || k > n) return 0.0;
+	if (k == 0 || k == n) return 1.0;
+
+	// Use symmetry.
+	if (k > n - k) k = n - k;
+
+	double result = 1.0;
+	for (int i = 0; i < k; i++)
+	{
+		result *= (double)(n - i);
+		result /= (double)(i + 1);
+	}
+	return result;
+}
+
+<*
+ Natural logarithm of factorial.
+*>
+fn double ln_factorial(int n)
+{
+	if (n < 0) return double.nan;
+	if (n <= 1) return 0.0;
+	return math::lgamma((double)(n + 1));
+}
+
+<*
+ Beta function B(a, b)
+ @require a > 0
+ @require b > 0
+*>
+fn double beta_function(double a, double b)
+{
+	return math::exp(math::lgamma(a) + math::lgamma(b) - math::lgamma(a + b));
+}
+
+<*
+ Regularized incomplete beta function I_x(a, b)
+ Based on https://github.com/codeplea/incbeta/blob/master/incbeta.c
+*>
+fn double incomplete_beta(double x, double a, double b, ConvergenceControl conv = DEFAULT_CONV)
+{
+	if (x < 0.0 || x > 1.0) return double.nan;
+	if (x == 0.0) return 0.0;
+	if (x == 1.0) return 1.0;
+
+	if (x > (a + 1.0)/(a + b + 2.0)) return 1.0 - incomplete_beta(1.0 - x, b,a);
+
+	const double TINY = 1e-30;
+
+	// Find the first part before the continued fraction.
+	double lbeta_ab = math::lgamma(a) + math::lgamma(b) - math::lgamma(a + b);
+	double front = math::exp(math::ln(x) * a + math::ln(1.0 - x) * b - lbeta_ab) / a;
+
+	// Use Lentz's algorithm to evaluate the continued fraction.
+	double f = 1.0;
+	double c = 1.0;
+	double d = 0.0;
+
+	usz m;
+	for (usz i = 0; i <= conv.max_iter; ++i)
+	{
+		m = i/2;
+
+		double numerator;
+
+		if (i == 0)
+		{
+			numerator = 1.0;
+		}
+		else if (i % 2 == 0)
+		{
+			numerator = (m * (b - m) * x) / ((a + 2.0 * m - 1.0) * (a + 2.0 * m));
+		}
+		else
+		{
+			numerator = -((a + m) * (a + b + m) * x) / ((a + 2.0 * m) * (a + 2.0 * m + 1));
+		}
+
+		d = 1.0 + numerator * d;
+		if (math::abs(d) < TINY) d = TINY;
+		d = 1.0 / d;
+
+		c = 1.0 + numerator / c;
+		if (math::abs(c) < TINY) c = TINY;
+
+		double cd = c*d;
+		f *= cd;
+
+		if (math::abs(1.0 - cd) < conv.epsilon) return front * (f - 1.0);
+	}
+
+	return double.nan; // Not convered.
+}
+
+<*
+ Calculates the p-th quantile for a continuous distribution using bisection.
+ @return? NOT_CONVERGED
+*>
+fn double? bisection_search(ContinuousDistribution dist, double low, double high, double p,
+	ConvergenceControl conv = DEFAULT_CONV)
+{
+	// Expand upper bound if needed.
+	while (dist.cdf(high) < p) high *= 2.0;
+
+	// Bisection search.
+	for (usz i = 0; i < conv.max_iter; i++)
+	{
+		double mid = (low + high) * 0.5;
+		if (high - low < conv.epsilon) return mid;
+
+		if (dist.cdf(mid) < p)
+		{
+			low = mid;
+		}
+		else
+		{
+			high = mid;
+		}
+	}
+	return NOT_CONVERGED~;
+}
+
+<*
+ Calculates the p-th quantile for a continuous distribution using Newton-Raphson.
+ @return? NOT_CONVERGED
+*>
+fn double? newton_raphson(ContinuousDistribution dist, double x, double p,
+	ConvergenceControl conv = DEFAULT_CONV)
+{
+	double delta, pdf;
+	for (usz i = 0; i < conv.max_iter; i++)
+	{
+		pdf = dist.pdf(x);
+		if (pdf < 1e-300) break;
+
+		delta = (dist.cdf(x) - p) / pdf;
+		x -= delta;
+
+		if (math::abs(delta) < conv.epsilon) return x;
+	}
+	return NOT_CONVERGED~;
+}
+
+<*
+ Calculates the p-th quantile for a continuous distribution.
+ @require p >= 0.0 && p <= 1.0
+ @require low < high
+*>
+fn double find_quantile(ContinuousDistribution dist, double low, double high, double p)
+{
+	double mid = bisection_search(dist, low, high, p, RELAXED_CONV) ?? (low + high) * 0.5;
+	return newton_raphson(dist, mid, p, DEFAULT_CONV) ?? mid;
+}
+
+<*
+ Generate a chi-squared random sample.
+ @param k : "Degrees of freedom"
+ @require k > 0.0
+*>
+fn double chi_squared_sample(double k, Random rand)
+{
+	// Sum of k squared standard normals.
+	NormalDist std_normal = normal(0.0, 1.0);
+	double sum = 0.0;
+
+	int k_int = (int)k;
+	for (int i = 0; i < k_int; i++)
+	{
+		double z = std_normal.random(rand);
+		sum += z * z;
+	}
+
+	// Handle fractional degrees of freedom.
+	double frac = k - (double)k_int;
+	if (frac > 0.0)
+	{
+		double z = std_normal.random(rand);
+		sum += frac * z * z;
+	}
+
+	return sum;
+}
+
+<*
+ Inverse of the error function (math::erf).
+ Based on Golang's math.Erfinv.
+*>
+fn double inverse_erf(double x)
+{
+	if (x < -1 || x > 1)
+	{
+		return double.nan;
+	}
+	else if (x == 1.0)
+	{
+		return double.inf;
+	}
+	else if (x == -1.0)
+	{
+		return -double.inf;
+	}
+
+	const double LN2 = 6.931471805599453094172321214581e-1;
+
+	const double A0 = 1.1975323115670912564578e0;
+	const double A1 = 4.7072688112383978012285e1;
+	const double A2 = 6.9706266534389598238465e2;
+	const double A3 = 4.8548868893843886794648e3;
+	const double A4 = 1.6235862515167575384252e4;
+	const double A5 = 2.3782041382114385731252e4;
+	const double A6 = 1.1819493347062294404278e4;
+	const double A7 = 8.8709406962545514830200e2;
+
+	const double B0 = 1.0000000000000000000e0;
+	const double B1 = 4.2313330701600911252e1;
+	const double B2 = 6.8718700749205790830e2;
+	const double B3 = 5.3941960214247511077e3;
+	const double B4 = 2.1213794301586595867e4;
+	const double B5 = 3.9307895800092710610e4;
+	const double B6 = 2.8729085735721942674e4;
+	const double B7 = 5.2264952788528545610e3;
+
+	const double C0 = 1.42343711074968357734e0;
+	const double C1 = 4.63033784615654529590e0;
+	const double C2 = 5.76949722146069140550e0;
+	const double C3 = 3.64784832476320460504e0;
+	const double C4 = 1.27045825245236838258e0;
+	const double C5 = 2.41780725177450611770e-1;
+	const double C6 = 2.27238449892691845833e-2;
+	const double C7 = 7.74545014278341407640e-4;
+
+	const double D0 = 1.4142135623730950488016887e0;
+	const double D1 = 2.9036514445419946173133295e0;
+	const double D2 = 2.3707661626024532365971225e0;
+	const double D3 = 9.7547832001787427186894837e-1;
+	const double D4 = 2.0945065210512749128288442e-1;
+	const double D5 = 2.1494160384252876777097297e-2;
+	const double D6 = 7.7441459065157709165577218e-4;
+	const double D7 = 1.4859850019840355905497876e-9;
+
+	const double E0 = 6.65790464350110377720e0;
+	const double E1 = 5.46378491116411436990e0;
+	const double E2 = 1.78482653991729133580e0;
+	const double E3 = 2.96560571828504891230e-1;
+	const double E4 = 2.65321895265761230930e-2;
+	const double E5 = 1.24266094738807843860e-3;
+	const double E6 = 2.71155556874348757815e-5;
+	const double E7 = 2.01033439929228813265e-7;
+
+	const double F0 = 1.414213562373095048801689e0;
+	const double F1 = 8.482908416595164588112026e-1;
+	const double F2 = 1.936480946950659106176712e-1;
+	const double F3 = 2.103693768272068968719679e-2;
+	const double F4 = 1.112800997078859844711555e-3;
+	const double F5 = 2.611088405080593625138020e-5;
+	const double F6 = 2.010321207683943062279931e-7;
+	const double F7 = 2.891024605872965461538222e-15;
+
+	double sign = 1.0;
+	if (x < 0)
+	{
+		x = -x;
+		sign = -1.0;
+	}
+
+	double ans;
+	if (x <= 0.85)
+	{
+		double r = 0.180625 - 0.25 * x *x;
+		double z1 = ((((((A7 * r + A6) * r + A5) * r + A4) * r + A3) * r + A2) * r + A1) * r + A0;
+		double z2 = ((((((B7 * r + B6) * r + B5) * r + B4) * r + B3) * r + B2) * r + B1) * r + B0;
+		ans = (x * z1) / z2;
+	}
+	else
+	{
+		double z1, z2;
+		double r = math::sqrt(LN2 - math::ln(1.0 - x));
+		if (r <= 5.0)
+		{
+			r -= 1.6;
+			z1 = ((((((C7 * r + C6) * r + C5) * r + C4) * r + C3) * r + C2) * r + C1) * r + C0;
+			z2 = ((((((D7 * r + D6) * r + D5) * r + D4) * r + D3) * r + D2) * r + D1) * r + D0;
+		}
+		else
+		{
+			r -= 5.0;
+			z1 = ((((((E7 * r + E6) * r + E5) * r + E4) * r + E3) * r + E2) * r + E1) * r + E0;
+			z2 = ((((((F7 * r + F6) * r + F5) * r + F4) * r + F3) * r + F2) * r + F1) * r + F0;
+		}
+		ans = z1 / z2;
+	}
+	return sign * ans;
+}
+
+<*
+ Regularized Lower incomplete gamma function.
+ Returns nan when not converged.
+ @param s : "Shape parameter"
+ @param x : "Upper limit of integration"
+ @require s > 0.0  : "s must be positive."
+ @require x >= 0.0 : "x must be non-negative."
+*>
+fn double lower_incomplete_gamma(double s, double x)
+{
+	if (x == 0.0) return 0.0;
+	if (x == double.inf) return 1.0;
+
+	// Use series expansion for x < s+1
+	if (x < s + 1.0)
+	{
+	    return incomplete_gamma_series_expansion(s, x) ?? double.nan;
+	}
+	else
+	{
+		return 1.0 - incomplete_gamma_continued_fraction(s, x) ?? double.nan;
+	}
+}
+
+<*
+ Lower incomplete gamma series expansion.
+ @return? NOT_CONVERGED
+*>
+fn double? incomplete_gamma_series_expansion(double s, double x, ConvergenceControl conv = DEFAULT_CONV)
+{
+	double lnpre = s * math::ln(x) - x - math::lgamma(s + 1.0);
+	if (lnpre < -708.0) return 0.0;   // result underflows to zero
+
+	double term = 1.0;
+	double sum  = 1.0;
+	double ap   = s;
+
+	for (int n = 1; n <= conv.max_iter; n++)
+	{
+		ap   += 1.0;
+		term *= x / ap;
+		sum  += term;
+		if (math::abs(term) < math::abs(sum) * conv.epsilon)
+		{
+		    return math::exp(lnpre) * sum;
+		}
+	}
+
+ 	// Non-convergence.
+	return NOT_CONVERGED~;
+}
+
+<*
+ Modified Lentz continued fraction.
+ @return? NOT_CONVERGED
+*>
+fn double? incomplete_gamma_continued_fraction(double s, double x, ConvergenceControl conv = DEFAULT_CONV)
+{
+	double lnpre = s * math::ln(x) - x - math::lgamma(s);
+
+	// Lentz initialisation: fpmin guards against division by zero.
+	double fpmin = 1e-300;
+
+	double b = x + 1.0 - s;
+	double c = 1.0 / fpmin;
+	double d = 1.0 / b;
+	double h = d;
+
+	for (int i = 1; i <= conv.max_iter; i++)
+	{
+		double an = (double)i * (s - (double)i);
+		b += 2.0;
+
+		d = an * d + b;
+		if (math::abs(d) < fpmin) d = fpmin;
+
+		c = b + an / c;
+		if (math::abs(c) < fpmin) c = fpmin;
+
+		d = 1.0 / d;
+		double delta = d * c;
+		h *= delta;
+
+		if (math::abs(delta - 1.0) < conv.epsilon)
+		{
+			return math::exp(lnpre) * h;
+		}
+	}
+
+ 	// Non-convergence.
+	return NOT_CONVERGED~;
+}
diff --git a/lib/std/math/math.c3 b/lib/std/math/math.c3
index 36d4768b6..8bbd7ffb3 100644
--- a/lib/std/math/math.c3
+++ b/lib/std/math/math.c3
@@ -124,6 +124,42 @@ macro sign(x)
 	$endif
 }
 
+<*
+ @require values::@is_int(x) || values::@is_float(x) : "Expected an integer or floating point value"
+*>
+macro erf(x)
+{
+	$if $typeof(x) == float:
+		return _erff(x);
+	$else
+		return _erf(x);
+	$endif
+}
+
+<*
+ @require values::@is_int(x) || values::@is_float(x) : "Expected an integer or floating point value"
+*>
+macro tgamma(x)
+{
+	$if $typeof(x) == float:
+		return _tgammaf(x);
+	$else
+		return _tgamma(x);
+	$endif
+}
+
+<*
+ @require values::@is_int(x) || values::@is_float(x) : "Expected an integer or floating point value"
+*>
+macro lgamma(x)
+{
+	$if $typeof(x) == float:
+		return _lgammaf(x);
+	$else
+		return _lgamma(x);
+	$endif
+}
+
 <*
  @require values::@is_int(x) || values::@is_float(x) : "Expected an integer or floating point value"
  @require values::@is_int(y) || values::@is_float(y) : "Expected an integer or floating point value"
@@ -1088,6 +1124,14 @@ extern fn float _acoshf(float x) @MathLibc("acoshf");
 extern fn float _asinhf(float x) @MathLibc("asinhf");
 extern fn float _atanhf(float x) @MathLibc("atanhf");
 
+extern fn double _erf(double x) @MathLibc("erf");
+extern fn double _lgamma(double x) @MathLibc("lgamma");
+extern fn double _tgamma(double x) @MathLibc("tgamma");
+
+extern fn float _erff(float x) @MathLibc("erf");
+extern fn float _lgammaf(float x) @MathLibc("lgammaf");
+extern fn float _tgammaf(float x) @MathLibc("tgammaf");
+
 
 fn double _frexp(double x, int* e)
 {
diff --git a/lib/std/math/math_nolibc/erf.c3 b/lib/std/math/math_nolibc/erf.c3
new file mode 100644
index 000000000..649924f34
--- /dev/null
+++ b/lib/std/math/math_nolibc/erf.c3
@@ -0,0 +1,77 @@
+module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH));
+
+<*
+ Error function for float.
+*>
+fn float _erff(float x)
+{
+	// Handle special cases.
+	if (x == 0.0) return 0.0;
+	if (x == float.inf) return 1.0;
+	if (x == -float.inf) return -1.0;
+
+	// Use symmetry: erf(-x) = -erf(x)
+	float sign = 1.0;
+	if (x < 0.0)
+	{
+	    x = -x;
+	    sign = -1.0;
+	}
+
+	const float P1 =  0.4067420160f;
+	const float P2 =  0.0072279182f;
+	const float A1 =  0.3168798904f;
+	const float A2 = -0.1383293141f;
+	const float A3 =  1.0868083034f;
+	const float A4 = -1.1169415512f;
+	const float A5 =  1.2064490307f;
+	const float A6 = -0.3931277152f;
+	const float A7 =  0.0382613542f;
+
+	float t = 1.0f / (1.0f + P1 * x + P2 * x * x);
+	float t2 = t * t;
+	float sum = A1 * t + A2 * t2 + A3 * t2 * t + A4 * t2 * t2;
+	sum += A5 * t2 * t2 * t + A6 * t2 * t2 * t2 + A7 * t2 * t2 * t2 * t;
+
+	// For intermediate x, use continued fraction.
+	return sign * (1.0f - sum * math::exp(-x * x));
+}
+
+<*
+ Error function for double.
+*>
+fn double _erf(double x)
+{
+	// Handle special cases.
+	if (x == 0.0) return 0.0;
+	if (x == double.inf) return 1.0;
+	if (x == -double.inf) return -1.0;
+
+	// Use symmetry: erf(-x) = -erf(x)
+	double sign = 1.0;
+	if (x < 0.0)
+	{
+	    x = -x;
+	    sign = -1.0;
+	}
+
+	const double P1 =  0.406742016006509;
+	const double P2 =  0.0072279182302319;
+	const double A1 =  0.316879890481381;
+	const double A2 = -0.138329314150635;
+	const double A3 =  1.08680830347054;
+	const double A4 = -1.11694155120396;
+	const double A5 =  1.20644903073232;
+	const double A6 = -0.393127715207728;
+	const double A7 =  0.0382613542530727;
+
+	double t = 1.0 / (1.0 + P1 * x + P2 * x * x);
+	double t2 = t * t;
+	double sum = A1 * t + A2 * t2 + A3 * t2 * t + A4 * t2 * t2;
+	sum += A5 * t2 * t2 * t + A6 * t2 * t2 * t2 + A7 * t2 * t2 * t2 * t;
+
+	// For intermediate x, use continued fraction.
+	return sign * (1.0 - sum * math::exp(-x * x));
+}
+
+
diff --git a/lib/std/math/math_nolibc/gamma.c3 b/lib/std/math/math_nolibc/gamma.c3
new file mode 100644
index 000000000..b8f62708b
--- /dev/null
+++ b/lib/std/math/math_nolibc/gamma.c3
@@ -0,0 +1,70 @@
+module std::math::nolibc @if(env::NO_LIBC || $feature(C3_MATH));
+
+// Constants
+const double SQRT_PI = 1.7724538509055160272981674833411451;
+const double LN_SQRT_2PI = 0.9189385332046727417803297364056176;
+
+<*
+ Natural logarithm of the gamma function using the Lanczos approximation.
+ @require x > 0.0 : "x must be positive."
+*>
+fn double lgamma(double x)
+{
+	// Lanczos approximation coefficients (g = 7, n = 9)
+	const double[*] LANCZOS_COEF = {
+	    0.99999999999980993,
+	    676.5203681218851,
+	    -1259.1392167224028,
+	    771.32342877765313,
+	    -176.61502916214059,
+	    12.507343278686905,
+	    -0.13857109526572012,
+	    9.9843695780195716e-6,
+	    1.5056327351493116e-7
+	};
+	const double LANCZOS_G = 7.0;
+
+	// For small x, use the reflection formula.
+	if (x < 0.5)
+	{
+	    return math::ln(math::PI) - math::ln(math::sin(math::PI * x)) - lgamma(1.0 - x);
+	}
+
+	// Shift to use approximation around x >= 1.5
+	x -= 1.0;
+
+	double base = x + LANCZOS_G + 0.5;
+	double sum = LANCZOS_COEF[0];
+
+	for (int i = 1; i < 9; i++)
+	{
+	    sum += LANCZOS_COEF[i] / (x + (double)i);
+	}
+
+	return LN_SQRT_2PI + math::ln(sum) + (x + 0.5) * math::ln(base) - base;
+}
+
+<*
+ Gamma function.
+ Valid for x > 0 and some negative non-integer values.
+*>
+fn double tgamma(double x)
+{
+	// Handle special cases.
+	if (x == 0.0) return double.inf;
+
+	// Check for negative integers (poles).
+	if (x < 0.0 && x == math::floor(x))
+	{
+	    return double.nan;
+	}
+
+	// For positive values, use exp(lgamma(x)).
+	if (x > 0.0)
+	{
+	    return math::exp(lgamma(x));
+	}
+
+	// For negative non-integer values, use the reflection formula.
+	return math::PI / (math::sin(math::PI * x) * tgamma(1.0 - x));
+}
diff --git a/releasenotes.md b/releasenotes.md
index 6781a2be8..ab8021956 100644
--- a/releasenotes.md
+++ b/releasenotes.md
@@ -28,6 +28,7 @@
 - Add Xorshiro128++.
 - Add single-byte code page support (DOS/OEM, Windows/ANSI, and ISO/IEC 8859).
 - Add `array::even`, `array::odd`, and `array::unlace` macros. #2892
+- Add discrete and continuous distributions in `std::math::distributions`.
 
 ### Fixes
 - Add error message if directory with output file name already exists
diff --git a/test/unit/stdlib/math/distributions.c3 b/test/unit/stdlib/math/distributions.c3
new file mode 100644
index 000000000..b981cc414
--- /dev/null
+++ b/test/unit/stdlib/math/distributions.c3
@@ -0,0 +1,530 @@
+// Copyright (c) 2026 Koni Marti. All rights reserved.
+// Use of this source code is governed by the MIT license.
+module std::math::distributions_test;
+import std::io, std::math, std::math::distributions @public;
+
+const double TOLERANCE = 1e-10;
+const double RELAXED_TOLERANCE = 1e-6;
+const double VERY_RELAXED_TOLERANCE = 1e-3;
+
+macro approx(double left, double right, double tol, String $msg = "")
+{
+	assert(math::is_approx(left, right, tol), "%s != %s (eps: %g): %s", left, right, tol, $msg);
+}
+
+fn void test_uniform_mean() @test
+{
+	UniformDist dist = distributions::uniform(0.0, 10.0);
+	double mean = dist.mean();
+	approx(mean, 5.0, TOLERANCE, "Uniform mean should be (a+b)/2");
+}
+
+fn void test_uniform_variance() @test
+{
+	UniformDist dist = distributions::uniform(0.0, 10.0);
+	double variance = dist.variance();
+	double expected = 100.0 / 12.0;
+	approx(variance, expected, TOLERANCE, "Uniform variance should be (b-a)²/12");
+}
+
+fn void test_uniform_pdf() @test
+{
+	UniformDist dist = distributions::uniform(0.0, 10.0);
+
+	// PDF should be constant in range
+	approx(dist.pdf(5.0), 0.1, TOLERANCE, "PDF at midpoint");
+	approx(dist.pdf(0.0), 0.1, TOLERANCE, "PDF at lower bound");
+	approx(dist.pdf(10.0), 0.1, TOLERANCE, "PDF at upper bound");
+
+	// PDF should be 0 outside range
+	approx(dist.pdf(-1.0), 0.0, TOLERANCE, "PDF below range");
+	approx(dist.pdf(11.0), 0.0, TOLERANCE, "PDF above range");
+}
+
+fn void test_uniform_cdf() @test
+{
+	UniformDist dist = distributions::uniform(0.0, 10.0);
+
+	approx(dist.cdf(-1.0), 0.0, TOLERANCE, "CDF below range");
+	approx(dist.cdf(0.0), 0.0, TOLERANCE, "CDF at lower bound");
+	approx(dist.cdf(5.0), 0.5, TOLERANCE, "CDF at midpoint");
+	approx(dist.cdf(7.5), 0.75, TOLERANCE, "CDF at 75%");
+	approx(dist.cdf(10.0), 1.0, TOLERANCE, "CDF at upper bound");
+	approx(dist.cdf(11.0), 1.0, TOLERANCE, "CDF above range");
+}
+
+fn void test_uniform_quantile() @test
+{
+	UniformDist dist = distributions::uniform(0.0, 10.0);
+
+	approx(dist.quantile(0.0), 0.0, TOLERANCE);
+	approx(dist.quantile(0.25), 2.5, TOLERANCE);
+	approx(dist.quantile(0.5), 5.0, TOLERANCE);
+	approx(dist.quantile(0.75), 7.5, TOLERANCE);
+	approx(dist.quantile(1.0), 10.0, TOLERANCE);
+}
+
+fn void test_uniform_round_trip() @test
+{
+	UniformDist dist = distributions::uniform(0.0, 10.0);
+
+	for (double x = 1.0; x <= 9.0; x += 2.0)
+	{
+		double p = dist.cdf(x);
+		double x_recovered = dist.quantile(p);
+		approx(x_recovered, x, TOLERANCE, "CDF/inverse_CDF round-trip");
+	}
+}
+
+fn void test_uniform_random_samples() @test
+{
+	UniformDist dist = distributions::uniform(0.0, 10.0);
+
+	DefaultRandom r;
+	random::seed_entropy(&r);
+
+	double sum;
+	int n_samples = 10_000;
+
+	for (int i = 0; i < n_samples; i++)
+	{
+		double sample = dist.random(&r);
+		sum += sample;
+		assert(sample >= 0.0 && sample <= 10.0, "Random sample in range");
+	}
+	approx(sum/n_samples, 5.000, 0.1);
+}
+
+fn void test_normal_mean_variance() @test
+{
+	NormalDist dist = distributions::normal(10.0, 2.0);
+	approx(dist.mean(), 10.0, TOLERANCE);
+	approx(dist.variance(), 4.0, TOLERANCE);
+}
+
+fn void test_normal_pdf_at_mean() @test
+{
+	NormalDist std_normal = distributions::normal(0.0, 1.0);
+	double pdf_0 = std_normal.pdf(0.0);
+	double expected = 1.0 / math::sqrt(2.0 * math::PI);
+	approx(pdf_0, expected, TOLERANCE, "PDF at mean");
+}
+
+fn void test_normal_pdf_symmetry() @test
+{
+	NormalDist std_normal = distributions::normal(0.0, 1.0);
+	approx(std_normal.pdf(-1.0), std_normal.pdf(1.0), TOLERANCE, "PDF symmetry");
+	approx(std_normal.pdf(-2.0), std_normal.pdf(2.0), TOLERANCE, "PDF symmetry at 2");
+}
+
+fn void test_normal_cdf_known_values() @test
+{
+	NormalDist std_normal = distributions::normal(0.0, 1.0);
+
+	approx(std_normal.cdf(0.0),   0.500, TOLERANCE, "CDF at mean");
+	approx(std_normal.cdf(-1.96), 0.025, VERY_RELAXED_TOLERANCE, "CDF at -1.96");
+	approx(std_normal.cdf(1.96),  0.975, VERY_RELAXED_TOLERANCE, "CDF at 1.96");
+	approx(std_normal.cdf(-2.58), 0.005, VERY_RELAXED_TOLERANCE, "CDF at -2.58");
+	approx(std_normal.cdf(2.58),  0.995, VERY_RELAXED_TOLERANCE, "CDF at 2.58");
+}
+
+fn void test_normal_quantile() @test
+{
+	NormalDist std_normal = distributions::normal(0.0, 1.0);
+
+	approx(std_normal.quantile(0.5), 0.0, TOLERANCE, "Median");
+
+	double z_975 = std_normal.quantile(0.975);
+	approx(z_975, 1.96, 0.02, "97.5th percentile");
+
+	double z_025 = std_normal.quantile(0.025);
+	approx(z_025, -1.96, 0.02, "2.5th percentile");
+}
+
+fn void test_normal_round_trip() @test
+{
+	NormalDist std_normal = distributions::normal(0.0, 1.0);
+
+	for (double x = -3.0; x <= 3.0; x += 1.0)
+	{
+		double p = std_normal.cdf(x);
+		double x_recovered = std_normal.quantile(p);
+		approx(x_recovered, x, TOLERANCE, "Normal round-trip");
+	}
+}
+
+fn void test_normal_random_samples() @test
+{
+	NormalDist std_normal = distributions::normal(0.0, 1.0);
+
+	double sum = 0.0;
+	double sum_sq = 0.0;
+	int n_samples = 10_000;
+
+	DefaultRandom r;
+	random::seed_entropy(&r);
+
+	for (int i = 0; i < n_samples; i++)
+	{
+		double sample = std_normal.random(&r);
+		sum += sample;
+		sum_sq += sample * sample;
+	}
+
+	double sample_mean = sum / (double)n_samples;
+	double sample_var = sum_sq / (double)n_samples - sample_mean * sample_mean;
+
+	approx(sample_mean, 0.0, 0.1, "Sample mean ~0");
+	approx(sample_var, 1.0, 0.1, "Sample variance ~1");
+}
+
+fn void test_normal_custom_parameters() @test
+{
+	NormalDist custom = distributions::normal(100.0, 15.0);
+
+	approx(custom.mean(), 100.0, TOLERANCE);
+	approx(custom.variance(), 225.0, TOLERANCE);
+	approx(custom.cdf(100.0), 0.5, TOLERANCE, "CDF at mean");
+}
+
+fn void test_exponential_mean_variance() @test
+{
+	ExponentialDist dist = distributions::exponential(2.0);
+	approx(dist.mean(), 0.5, TOLERANCE);
+	approx(dist.variance(), 0.25, TOLERANCE);
+}
+
+fn void test_exponential_pdf() @test
+{
+	ExponentialDist dist = distributions::exponential(2.0);
+
+	approx(dist.pdf(0.0), 2.0, TOLERANCE, "PDF at 0");
+	approx(dist.pdf(-1.0), 0.0, TOLERANCE, "PDF for negative x");
+
+	double pdf_1 = dist.pdf(1.0);
+	double expected = 2.0 * math::exp(-2.0);
+	approx(pdf_1, expected, TOLERANCE, "PDF at 1");
+}
+
+fn void test_exponential_cdf() @test
+{
+	ExponentialDist dist = distributions::exponential(2.0);
+
+	approx(dist.cdf(0.0), 0.0, TOLERANCE, "CDF at 0");
+	approx(dist.cdf(-1.0), 0.0, TOLERANCE, "CDF for negative x");
+
+	double cdf_mean = dist.cdf(dist.mean());
+	approx(cdf_mean, 1.0 - math::exp(-1.0), TOLERANCE, "CDF at mean");
+}
+
+fn void test_exponential_quantile() @test
+{
+	ExponentialDist dist = distributions::exponential(2.0);
+
+	approx(dist.quantile(0.0), 0.0, TOLERANCE);
+
+	double median = dist.quantile(0.5);
+	double expected_median = math::ln(2.0) / 2.0;
+	approx(median, expected_median, TOLERANCE, "Median");
+}
+
+fn void test_exponential_round_trip() @test
+{
+	ExponentialDist dist = distributions::exponential(2.0);
+
+	for (double p = 0.1; p <= 0.9; p += 0.1)
+	{
+		double x = dist.quantile(p);
+		double p_recovered = dist.cdf(x);
+		approx(p_recovered, p, TOLERANCE, "Exponential round-trip");
+	}
+}
+
+fn void test_exponential_random_samples() @test
+{
+	ExponentialDist dist = distributions::exponential(2.0);
+
+	DefaultRandom r;
+	random::seed_entropy(&r);
+
+	for (int i = 0; i < 100; i++)
+	{
+		double sample = dist.random(&r);
+		assert(sample >= 0.0, "Random sample non-negative");
+	}
+}
+
+fn void test_t_mean_variance() @test
+{
+	TDist dist = distributions::t_distribution(10.0);
+	approx(dist.mean(), 0.0, TOLERANCE);
+	approx(dist.variance(), 10.0/8.0, TOLERANCE);
+}
+
+fn void test_t_pdf_symmetry() @test
+{
+	TDist dist = distributions::t_distribution(10.0);
+	approx(dist.pdf(-1.0), dist.pdf(1.0), TOLERANCE, "PDF symmetry");
+	approx(dist.pdf(-2.0), dist.pdf(2.0), TOLERANCE, "PDF symmetry at 2");
+}
+
+fn void test_t_cdf_symmetry() @test
+{
+	TDist dist = distributions::t_distribution(10.0);
+
+	approx(dist.cdf(0.0), 0.5, TOLERANCE, "CDF at 0");
+
+	double cdf_1 = dist.cdf(1.0);
+	double cdf_neg_1 = dist.cdf(-1.0);
+	approx(cdf_1 + cdf_neg_1, 1.0, TOLERANCE, "CDF symmetry");
+}
+
+fn void test_t_quantile() @test
+{
+	TDist dist = distributions::t_distribution(10.0);
+
+	approx(dist.quantile(0.5), 0.0, RELAXED_TOLERANCE, "Median is 0");
+
+	double upper = dist.quantile(0.975);
+	double lower = dist.quantile(0.025);
+	approx(upper, -lower, RELAXED_TOLERANCE, "Inverse CDF symmetry");
+}
+
+
+fn void test_t_random_samples() @test
+{
+	TDist dist = distributions::t_distribution(10.0);
+
+	double sum = 0.0;
+	int n_samples = 10_000;
+
+	DefaultRandom r;
+	random::seed_entropy(&r);
+
+	for (int i = 0; i < n_samples; i++)
+	{
+		sum += dist.random(&r);
+	}
+
+	double sample_mean = sum / (double)n_samples;
+	approx(sample_mean, 0.0, 0.1, "Sample mean ~0");
+}
+
+fn void test_f_mean() @test
+{
+	FDist dist =distributions::f_distribution(5.0, 10.0);
+	double expected_mean = 10.0 / 8.0;
+	approx(dist.mean(), expected_mean, TOLERANCE);
+}
+
+fn void test_f_pdf() @test
+{
+	// FIXME: add tests for pdf
+	FDist dist =distributions::f_distribution(5.0, 10.0);
+	approx(dist.pdf(0.5), 0.687607, RELAXED_TOLERANCE, "PDF at x=0.5 for F(5,10) is wrong");
+	approx(dist.pdf(10.0), 0.000478163, RELAXED_TOLERANCE, "PDF at x=10.0 for F(5,10) is wrong");
+}
+
+fn void test_f_cdf() @test
+{
+	FDist dist =distributions::f_distribution(5.0, 10.0);
+	approx(dist.cdf(0.0), 0.0, TOLERANCE, "CDF at 0");
+	approx(dist.cdf(-1.0), 0.0, TOLERANCE, "CDF for negative x");
+}
+
+fn void test_f_quantile() @test
+{
+	FDist dist =distributions::f_distribution(5.0, 10.0);
+	double x_median = dist.quantile(0.5);
+	assert(x_median > 0.0, "Median positive");
+
+	double x_25 = dist.quantile(0.25);
+	double x_75 = dist.quantile(0.75);
+	assert(x_75 > x_25, "Inverse CDF monotonic");
+}
+
+fn void test_f_random_samples() @test
+{
+	FDist dist =distributions::f_distribution(5.0, 10.0);
+
+ 	DefaultRandom r;
+	random::seed_entropy(&r);
+
+	for (int i = 0; i < 100; i++)
+	{
+		double sample = dist.random(&r);
+		assert(sample > 0.0, "Random sample positive");
+	}
+}
+
+fn void test_chi_squared_mean_variance() @test
+{
+	ChiSquaredDist dist = distributions::chi_squared(5.0);
+	approx(dist.mean(), 5.0, TOLERANCE);
+	approx(dist.variance(), 10.0, TOLERANCE);
+}
+
+fn void test_chi_squared_pdf() @test
+{
+	ChiSquaredDist dist = distributions::chi_squared(5.0);
+	approx(dist.pdf(-1.0), 0.0, TOLERANCE, "PDF for negative x");
+	approx(dist.pdf(1.145), 0.091910, RELAXED_TOLERANCE);
+	approx(dist.pdf(5.000), 0.122042, RELAXED_TOLERANCE);
+}
+
+fn void test_chi_squared_cdf() @test
+{
+	ChiSquaredDist dist = distributions::chi_squared(1.0);
+	approx(dist.cdf(1.0), 0.682689, RELAXED_TOLERANCE, "");
+
+	dist.k = 5;
+	approx(dist.cdf(5.0), 0.5841, VERY_RELAXED_TOLERANCE);
+}
+
+fn void test_chi_squared_quantile() @test
+{
+	ChiSquaredDist dist = distributions::chi_squared(5.0);
+	approx(dist.quantile(0.95), 11.0705, 0.1, "95th percentile");
+}
+
+fn void test_chi_squared_round_trip() @test
+{
+	ChiSquaredDist dist = distributions::chi_squared(5.0);
+	for (double p = 0.1; p <= 0.9; p += 0.2)
+	{
+		double x = dist.quantile(p);
+		assert(x > 0.0, "Inverse CDF positive");
+		double p_recovered = dist.cdf(x);
+		approx(p_recovered, p, RELAXED_TOLERANCE, "Chi-squared round-trip");
+	}
+}
+
+fn void test_chi_squared_random_samples() @test
+{
+	ChiSquaredDist dist = distributions::chi_squared(5.0);
+
+ 	DefaultRandom r;
+	random::seed_entropy(&r);
+
+	for (int i = 0; i < 100; i++)
+	{
+		double sample = dist.random(&r);
+		assert(sample >= 0.0, "Random sample non-negative");
+	}
+}
+
+fn void test_binomial_mean_variance() @test
+{
+	BinomialDist dist = distributions::binomial(10, 0.5);
+	approx(dist.mean(), 5.0, TOLERANCE);
+	approx(dist.variance(), 2.5, TOLERANCE);
+}
+
+fn void test_binomial_pmf_symmetry() @test
+{
+	BinomialDist dist = distributions::binomial(10, 0.5);
+	for (int k = 0; k <= 5; k++)
+	{
+		double pmf_k = dist.pmf(k);
+		double pmf_10_minus_k = dist.pmf(10 - k);
+		approx(pmf_k, pmf_10_minus_k, TOLERANCE, "PMF symmetry");
+	}
+}
+
+fn void test_binomial_pmf_sums_to_one() @test
+{
+	BinomialDist dist = distributions::binomial(10, 0.5);
+	double sum = 0.0;
+	for (int k = 0; k <= 10; k++)
+	{
+		sum += dist.pmf(k);
+	}
+	approx(sum, 1.0, TOLERANCE, "PMF sums to 1");
+}
+
+fn void test_binomial_cdf() @test
+{
+	BinomialDist dist = distributions::binomial(10, 0.5);
+	approx(dist.cdf(-1), 0.0, TOLERANCE, "CDF below 0");
+	approx(dist.cdf(10), 1.0, TOLERANCE, "CDF at n");
+	approx(dist.cdf(5), 0.623, VERY_RELAXED_TOLERANCE, "CDF at mean");
+}
+
+fn void test_binomial_random_samples() @test
+{
+	BinomialDist dist = distributions::binomial(10, 0.5);
+
+	DefaultRandom r;
+	random::seed_entropy(&r);
+
+	for (int i = 0; i < 100; i++)
+	{
+		double sample = dist.random(&r);
+		assert(sample >= 0.0 && sample <= 10.0, "Random in range");
+	}
+}
+
+fn void test_poisson_mean_variance() @test
+{
+	PoissonDist dist = distributions::poisson(3.0);
+	approx(dist.mean(), 3.0, TOLERANCE);
+	approx(dist.variance(), 3.0, TOLERANCE);
+}
+
+fn void test_poisson_pmf_sums_to_one() @test
+{
+	PoissonDist dist = distributions::poisson(3.0);
+	double sum = 0.0;
+	for (int k = 0; k <= 20; k++)
+	{
+		sum += dist.pmf(k);
+	}
+	assert(sum > 0.99, "PMF sums to ~1");
+}
+
+fn void test_poisson_pmf_at_zero() @test
+{
+	PoissonDist dist = distributions::poisson(3.0);
+	double pmf_0 = dist.pmf(0);
+	double expected = math::exp(-3.0);
+	approx(pmf_0, expected, TOLERANCE, "PMF at 0");
+}
+
+fn void test_poisson_cdf() @test
+{
+	PoissonDist dist = distributions::poisson(3.0);
+	approx(dist.cdf(-1), 0.0, TOLERANCE, "CDF for negative k");
+}
+
+fn void test_poisson_random_samples() @test
+{
+	PoissonDist dist = distributions::poisson(3.0);
+
+	DefaultRandom r;
+	random::seed_entropy(&r);
+
+	for (int i = 0; i < 100; i++)
+	{
+		double sample = dist.random(&r);
+		assert(sample >= 0.0, "Random non-negative");
+	}
+}
+
+fn void test_poisson_large_lambda() @test
+{
+	PoissonDist large = distributions::poisson(50.0);
+	approx(large.mean(), 50.0, TOLERANCE, "Large lambda mean");
+}
+
+fn void test_beta_function() @test
+{
+	assert(distributions::beta_function(2.5, 1.5) == distributions::beta_function(1.5, 2.5), "Beta function not symmetrical.");
+	assert(math::is_approx(distributions::beta_function(3.0, 2.0), 0.08333333333, 1e-6), "Beta(3,2) is wrong");
+	assert(math::is_approx(distributions::beta_function(4.0, 1.0), 1.0/4.0, 1e-6), "Beta(4,1) is wrong");
+}
+
+fn void test_lower_incomplete_gamma() @test
+{
+	approx(distributions::lower_incomplete_gamma(0.5,0.5), 0.682689, 1e-4, "");
+}
+