sml-specfun

Real-valued special functions in pure Standard ML — the analytic substrate that statistics and probability lean on, in one vetted, dependency-free suite: the gamma family (gamma, lgamma, beta, digamma/psi), the error functions (erf, erfc, and the inverse erfInv), and the regularized incomplete gamma (gammaIncP/gammaIncQ) and beta (betaInc) functions. No FFI, no external dependencies, and deterministic — every kernel uses only basis-library arithmetic and behaves identically under both MLton and Poly/ML.

=== sml-specfun demo =========================================

Gamma family
  gamma(1/2)  [= sqrt pi]    = 1.772453850906
  gamma(5)    [= 4! = 24]    = 24.000000000000
  lgamma(100) [= ln 99!]     = 359.134205369575
  beta(2,3)   [= 1/12]       = 0.083333333333
  digamma(1)  [= -euler]     = ~0.577215664901

Error functions
  erf(1)                     = 0.842700792950
  erfc(2)                    = 0.004677734981
  erfc(5)  [far tail]        = 0.000000000001537
  erfInv(0.8427007929)       = 1.000000000000

Regularized incomplete functions
  P(3,5)  lower gamma        = 0.875347980517
  Q(3,5)  upper gamma        = 0.124652019483
  I_0.5(2,3) inc. beta       = 0.687500000000

=============================================================

Generated by examples/demo.sml (make example). Every printed real is formatted with Real.fmt (StringCvt.FIX ...), which always emits a decimal point, so this output is byte-identical run to run and across MLton and Poly/ML (we never use Real.toString, whose formatting differs between the two). The leading ~ on digamma(1) is SML's negative sign.

Status

• 65 assertions, green on MLton and Poly/ML (N passed, 0 failed), with byte-identical output.
• Basis-library only — lib/.../specfun.{sig,sml} has no dependencies.
• Deterministic across compilers: every series / continued-fraction / Newton kernel iterates to a fixed value tolerance (never a wall-clock or step budget that could diverge between compilers), using only +,-,*,/ and the Basis Math primitives.
• This generalizes the Abramowitz & Stegun erf approximation currently inlined in sml-stats into a reusable suite, unblocking proper distribution CDFs (Student-t, chi-squared, F, gamma, beta) that can share one vetted implementation.

Install

With smlpkg:

smlpkg add github.com/sjqtentacles/sml-specfun
smlpkg sync

Include the library MLB from your own build (basis-only, no vendored deps):

local
  $(SML_LIB)/basis/basis.mlb
  lib/github.com/sjqtentacles/sml-specfun/sources.mlb (via smlpkg)
in
  ...
end

This brings structure Specfun into scope.

Quick start

val gp  = Specfun.gamma 0.5            (* sqrt pi = 1.7724538509... *)
val lg  = Specfun.lgamma 100.0         (* ln 99!  = 359.1342053695... *)
val b   = Specfun.beta (2.0, 3.0)      (* 1/12 *)
val e   = Specfun.erf 1.0              (* 0.8427007929... *)
val tail = Specfun.erfc 5.0            (* 1.537e-12, accurate in the tail *)
val xi  = Specfun.erfInv 0.95          (* inverse error function *)

(* regularized incomplete functions -- the building blocks for CDFs *)
val p   = Specfun.gammaIncP (3.0, 5.0) (* lower P(a,x) *)
val q   = Specfun.gammaIncQ (3.0, 5.0) (* upper Q(a,x), P + Q = 1 *)
val ix  = Specfun.betaInc (2.0, 3.0, 0.5)  (* regularized I_x(a,b) *)

API (signature SPECFUN)

val eps : real                              (* internal convergence tolerance *)

(* gamma family *)
val gamma   : real -> real                  (* Lanczos, reflection for x < 1/2 *)
val lgamma  : real -> real                  (* ln |G(x)|, stable for large x   *)
val beta    : real * real -> real           (* B(a,b) = G(a)G(b)/G(a+b)        *)
val lbeta   : real * real -> real
val digamma : real -> real                  (* psi = d/dx ln G(x)              *)
val psi     : real -> real                  (* alias for digamma              *)

(* error functions *)
val erf     : real -> real
val erfc    : real -> real                  (* 1 - erf, accurate in the tail   *)
val erfInv  : real -> real                  (* inverse of erf on (-1, 1)       *)

(* regularized incomplete functions *)
val gammaIncP : real * real -> real         (* lower P(a,x), a > 0, x >= 0     *)
val gammaIncQ : real * real -> real         (* upper Q(a,x), P + Q = 1         *)
val betaInc   : real * real * real -> real  (* I_x(a,b), a,b > 0, 0 <= x <= 1  *)

Methods & conventions

• gamma / lgamma — the Lanczos approximation (g = 7, 9-term series). gamma uses the reflection formula G(x)G(1−x) = π / sin(πx) for x < 1/2; lgamma returns ln |G(x)| and stays finite where gamma would overflow (e.g. lgamma 100 = ln 99!).
• beta / lbeta — lbeta(a,b) = lgamma a + lgamma b − lgamma(a+b), and beta = exp ∘ lbeta.
• digamma / psi — the asymptotic series, pushing the argument up with the recurrence ψ(x) = ψ(x+1) − 1/x until it is large, with reflection for the left half-plane. Poles at non-positive integers raise Domain.
• erf / erfc — expressed through the regularized incomplete gamma: erf x = sign(x)·P(½, x²) and erfc x = Q(½, x²) for x ≥ 0, so the far tail of erfc stays accurate instead of the catastrophic 1 − erf x.
• erfInv — Winitzki's closed-form initial guess refined by Newton's method on erf, iterated to tolerance.
• gammaIncP / gammaIncQ — a power series for x < a+1 and a Lentz-evaluated continued fraction for x ≥ a+1, the textbook split. Both raise Domain for a ≤ 0 or x < 0.
• betaInc — the prefactor xᵃ(1−x)ᵇ / (a·B(a,b)) times a Lentz continued fraction, taking the I_x(a,b) = 1 − I_{1−x}(b,a) branch for fast convergence. Raises Domain outside a,b > 0, 0 ≤ x ≤ 1.

Build & test

make test        # MLton: build + run the suite
make test-poly   # Poly/ML: use-and-run the suite
make all-tests   # both
make example     # print the deterministic reference table above
make clean

Both compilers run the same strict-TDD suite (test/), every comparison made through an explicit epsilon (never string-comparing reals), against known reference values:

• gamma — gamma(½) = √π, gamma(5) = 24, half-integer values, the reflection branch (gamma(−½) = −2√π), and the recurrence G(x+1) = x·G(x).
• lgamma — lgamma(½) = ln√π, lgamma(10) = ln 9!, lgamma(100) = ln 99!.
• beta / digamma — B(2,3) = 1/12, B(½,½) = π; ψ(1) = −γ, ψ(½) = −γ − 2 ln 2, ψ(10), and the recurrence ψ(x+1) = ψ(x) + 1/x.
• erf / erfc — erf(1) = 0.8427007929…, erfc(2) = 0.0046777349…, the far tail erfc(5) ≈ 1.537e-12, oddness of erf, and erf x + erfc x = 1.
• erfInv — erfInv(0.5), and the round-trips erfInv(erf x) = x, erf(erfInv y) = y.
• incomplete gamma/beta — closed forms for integer a (P(1,x) = 1 − e⁻ˣ, P(3,5)), P + Q = 1, the cross-check erf(x) = P(½, x²), and incomplete-beta closed forms (I_x(1,1) = x, I_{0.5}(2,3) = 0.6875, I_x(½,½) = (2/π)·arcsin√x) plus the symmetry I_x(a,b) = 1 − I_{1−x}(b,a).

Poly/ML note

CI builds Poly/ML 5.9.1 from source rather than using the Ubuntu package (Poly/ML 5.7.1), whose X86 code generator crashes (asGenReg raised while compiling) on heavy real-arithmetic code. See .github/workflows/ci.yml.

License

MIT — see LICENSE.