Help us shape the future of Go! Take the 2026 Go Developer Survey

special

package module
v0.1.0 Latest Latest
Warning

This package is not in the latest version of its module.

 Go to latest Published: Oct 11, 2026 License: BSD-3-Clause Imports: 3 Imported by: 0

README

special

Special mathematical functions in pure Go (no cgo, no assembly): the gamma family, error functions, incomplete gamma and beta functions and their inverses, Bessel functions, zeta functions, combinatorics, logistic helpers and complete elliptic integrals — with an element-wise layer over go-ndarray.

It is the numerical foundation of the go-numerics statistics library: every distribution's cdf and quantile function is computed with these functions, so their accuracy is the product. Each function is measured in ULPs against mpmath at high precision, over its whole domain — tails, zeros, poles, huge and tiny arguments, subnormals, NaN and ±Inf — and SciPy is measured on exactly the same points.

go get github.com/go-numerics/special

package main

import (
	"fmt"

	"github.com/go-numerics/special"
)

func main() {
	fmt.Println(special.Gamma(4.5))              // 11.631728396567448
	fmt.Println(special.Erfc(10))                // 2.088487583762545e-45
	fmt.Println(special.Ndtr(-30))               // 4.906713927148187e-198
	fmt.Println(special.Ndtri(1e-300))           // -37.0470962993612
	fmt.Println(special.GammaP(2.5, 1.5))        // regularised lower incomplete gamma
	fmt.Println(special.BetaIncInv(2, 3, 0.975)) // a beta quantile
	fmt.Println(special.Jv(2.5, 10))             // Bessel function of real order
}

Arrays (github.com/go-numerics/special/nd): any real dtype in, Float32 kept as Float32, everything else Float64; NumPy broadcasting; large arrays are split across cores.

x, _ := ndarray.Linspace(-5, 5, 1_000_001)
p, _ := nd.Ndtr(x)

Accuracy

Measured on about 250 000 points in 419 regions over 62 functions, this package's maximum error is at or below SciPy's in 418 regions (the exception, K1 at small x, is a tie at 1.05 ulp). Typical maxima:

this package SciPy 1.18.1
Gamma, LogGamma, Beta (all regions) 0.5–1 ulp 3–3e15 ulp
Digamma near its negative zeros 31 ulp 2.5e17 ulp
Erfc to the subnormal range 0.5 ulp 2e13 ulp
Ndtr lower tail, x in [-38.5, -10) 0.5 ulp 1.8e12 ulp
Ndtri, ErfInv, ErfcInv 0.5–0.6 ulp 1.4–4.4 ulp (wrong at subnormals)
GammaP / GammaQ, a in [1e4, 1e15] 0.6 ulp 8.9e15 ulp
BetaInc, a and b in [100, 1e4) 9 ulp 3487 ulp
J0 / Y0 next to their zeros 0.5 ulp 3e10 / 2e12 ulp
Jv at x in [1e4, 1e300] 1.7 ulp 5e18 ulp

The one area not solved: inverses of the incomplete beta function with a or b far below 1e-3 (both packages can be wrong there).

Full method, sources of every algorithm, negative control and limitations: docs/accuracy.md; every region: docs/accuracy-table.md.

Speed

Accuracy has a price: single-threaded, a Go loop is about 3x slower than SciPy's C for the gamma and error-function families and 10–15x for the Bessel functions of real order and the incomplete beta function. The nd functions spread large arrays over the cores and beat SciPy for 21 of 31 functions measured (typically 1.5–3x on a 16-core host; e.g. Gamma 6.4 vs 13.4 ns/element, Ndtr 4.7 vs 9.0, GammaP 37 vs 137).

Method and full table: docs/perf.md.

Functions

family functions
gamma Gamma, LogGamma, LogGammaSign, GammaSign, Digamma, Polygamma, Beta, LogBeta
incomplete gamma / beta GammaP, GammaQ, GammaPInv, GammaQInv, BetaInc, BetaIncC, BetaIncInv, BetaIncCInv
error function, normal Erf, Erfc, Erfcx, ErfInv, ErfcInv, Ndtr, Ndtri, LogNdtr
combinatorics Comb, Perm, Binom, Factorial, Factorial2, LogFactorial
logistic, log helpers Expit, LogExpit, Logit, Log1pExp, Xlogy, Xlog1py, LogSumExp, Softmax, LogSoftmax
zeta Zeta, HurwitzZeta
Bessel J0, J1, Jn, Jv, Y0, Y1, Yn, Yv, I0, I1, In, Iv, K0, K1, Kn, Kv, I0e, I1e, Ive, K0e, K1e, Kve, JnZeros
elliptic EllipK, EllipKm1, EllipE

Conventions: functions return NaN outside their domain and ±Inf at poles, as the IEEE 754 operations do; nothing panics on a float argument. Edge values follow SciPy's where SciPy's are mathematically sound; the exceptions are documented at each function.

Portability

Pure Go; CI runs the full accuracy suite on amd64, arm64, 386, arm (32-bit), riscv64, loong64, ppc64le and s390x, compiles the tests for every go tool dist list target, and holds a 100% statement-coverage gate.

One finding along the way: on amd64, Go's math.Log (assembly) returns about -709.09 for every subnormal argument (it should go down to -744.44); this package never calls it on a value that can be subnormal.

License

BSD 3-Clause. Ported algorithms keep their attribution in place: Cephes (Stephen L. Moshier, distributed under BSD by SciPy) and Boost.Math (Boost Software License 1.0).

Documentation

Overview

Package special implements special mathematical functions in pure Go (no cgo): the gamma family, error functions, incomplete gamma and beta functions and their inverses, Bessel functions, zeta functions and combinatorial and logistic helpers.

Accuracy is the product. Every function is measured in ULPs against mpmath evaluated at high precision over its whole domain, tails, poles and special values included; the fixtures and the measured table are in the repository (testdata/fixtures, docs/accuracy.md). Each algorithm cites its source in a comment at its definition.

Conventions: a function returns NaN outside its domain and ±Inf at a pole, as the IEEE 754 operations do; no function panics on a float argument.

Example 
package main

import (
	"fmt"

	"github.com/go-numerics/special"
)

func main() {
	fmt.Println(special.Gamma(4.5))
	fmt.Println(special.Erfc(10))              // 2.088487583762545e-45, no underflow to 0
	fmt.Println(special.Ndtr(-30))             // the far tail keeps full relative accuracy
	fmt.Println(special.Ndtri(1e-300))         // and so does its inverse
	fmt.Println(special.LogGamma(1 + 0x1p-30)) // log Γ near its zero at 1
}

Output:
11.631728396567448
2.088487583762545e-45
4.906713927148187e-198
-37.0470962993612
-5.375739784311044e-10

Index

Examples

Constants

This section is empty.

Variables

This section is empty.

Functions

func Beta 

func Beta(a, b float64) float64

Beta returns the beta function B(a, b) = Γ(a)Γ(b)/Γ(a+b).

Special values: NaN if either argument is NaN; +Inf when a or b is 0 or a negative integer and a+b is not a non-positive integer; 0 when a+b is a non-positive integer and neither a nor b is; overflow to ±Inf and underflow to 0 as the true value dictates.

Algorithm: for positive arguments below 10 the three gamma values come from Gamma's double-double core and the rounding of a+b is corrected to first order (Γ(s+e) = Γ(s)(1 + ψ(s)e)); otherwise B = exp(log B) with log B from LogBeta's cancellation-free double-double evaluation.

func BetaInc 

func BetaInc(a, b, x float64) float64

BetaInc returns the regularised incomplete beta function I_x(a, b) = (1/B(a, b)) ∫₀ˣ t^(a-1) (1-t)^(b-1) dt (DLMF 8.17.2).

Domain: a >= 0, b >= 0, 0 <= x <= 1, not both a and b zero or infinite. Special values: I_0 = 0 and I_1 = 1; for a = 0 (or b = +Inf), I_x = 1 for x > 0; for b = 0 (or a = +Inf), I_x = 0 for x < 1; NaN for an argument outside the domain or NaN.

Algorithm: DiDonato & Morris' TOMS 708 (BRATIO) with double-double prefixes; see the comment at the top of incbeta.go.

func BetaIncC 

func BetaIncC(a, b, x float64) float64

BetaIncC returns the complement of the regularised incomplete beta function, 1 - I_x(a, b) = I_(1-x)(b, a), computed directly where it is the smaller of the two so that it keeps its relative accuracy for x near 1.

Domain and special values: as BetaInc, with the complementary values.

Algorithm: as BetaInc.

func BetaIncCInv 

func BetaIncCInv(a, b, q float64) float64

BetaIncCInv returns the inverse of the complement of the regularised incomplete beta function: the x in [0, 1] with 1 - I_x(a, b) = q.

Domain and special values: as BetaIncInv, with BetaIncCInv(a, b, 0) = 1 and BetaIncCInv(a, b, 1) = 0.

Algorithm: as BetaIncInv, with q as the exact target of the upper tail (so a small q keeps full relative accuracy in 1 - x).

func BetaIncInv 

func BetaIncInv(a, b, p float64) float64

BetaIncInv returns the inverse of the regularised incomplete beta function: the x in [0, 1] with I_x(a, b) = p.

Domain: a > 0, b > 0, 0 <= p <= 1. Special values: BetaIncInv(a, b, 0) = 0 and BetaIncInv(a, b, 1) = 1 for any a, b >= 0 that are not NaN, with a = 0 or b = 0 giving NaN (SciPy's conventions); NaN for p outside [0, 1], a or b negative, infinite (for 0 < p < 1) or NaN. A root closer to 0 than the smallest subnormal returns 0, one closer to 1 than half an ulp of 1 returns 1.

Algorithm: safeguarded Newton iteration in log x or log(1-x) on the smaller tail; see the comment at the top of incbeta_inv.go.

func Binom 

func Binom(n, k float64) float64

Binom returns the generalised binomial coefficient Γ(n+1) / (Γ(k+1) Γ(n-k+1)) for real n and k, as SciPy's binom.

For integer k with 0 <= k <= 64 it is the product Π (n-k+i)/i in double-double (any real n, including negative); when n+1 > 0 and n-k+1 > 0 it is 1/((n+1) B(n-k+1, k+1)) with the beta function's cancellation-free evaluation; otherwise the logarithms of the three gamma functions (each in double-double, at double-double arguments) are combined before a single exponential, with their signs. It is 0 where 1/Γ(k+1) or 1/Γ(n-k+1) vanishes (a non-positive integer argument), NaN for NaN arguments or a pole of Γ(n+1) not cancelled by the denominator.

func Comb 

func Comb(n, k int64) float64

Comb returns the binomial coefficient C(n, k) = n! / (k! (n-k)!), the number of k-element subsets of an n-set, as a float64: 0 when k < 0, k > n or n < 0 (as SciPy's comb with exact=False).

Algorithm: for min(k, n-k) <= 64 the product Π (n-k+i)/i in double-double, so every result below 2^53 is exact and larger ones are correctly rounded; beyond, 1/((n+1) B(k+1, n-k+1)) with the beta function's cancellation-free double-double logarithm (LogBeta), which avoids the catastrophic cancellation of log n! - log k! - log (n-k)! for large n.

func Digamma 

func Digamma(x float64) float64

Digamma returns the digamma function ψ(x) = d/dx log Γ(x).

Special values: ψ(+Inf) = +Inf; ψ(+0) = -Inf, ψ(-0) = +Inf; NaN at the poles x = -1, -2, ..., at -Inf and for NaN.

Algorithm: on [1, 2] ψ(x) = (x - x0) D(x-1) with x0 the positive zero in double-double and D fitted to mpmath (the form of Boost.Math's digamma_imp_1_2); (0, 10) is reduced to [1, 2] by ψ(x+1) = ψ(x) + 1/x with the reciprocals in double-double; x >= 10 uses the asymptotic series (DLMF 5.11.2). For -40 < x < 0, ψ(x) = ψ(t) - Σ 1/(x+k) with t = x+n in [0, 1) and exact x+k, summed in double-double; below, the reflection ψ(x) = ψ(-x) - 1/x - π cot(πx) (DLMF 5.5.4) with exact reduction of πx.

func EllipE 

func EllipE(m float64) float64

EllipE returns the complete elliptic integral of the second kind E(m) = ∫₀^{π/2} sqrt(1 - m sin²θ) dθ, parameter m = k², for m <= 1.

Special values: E(1) = 1, E(-Inf) = +Inf, NaN for m > 1 or NaN.

func EllipK 

func EllipK(m float64) float64

EllipK returns the complete elliptic integral of the first kind K(m) = ∫₀^{π/2} dθ / sqrt(1 - m sin²θ), parameter m = k² (SciPy's ellipk convention), for m <= 1.

Special values: K(1) = +Inf, K(-Inf) = 0, NaN for m > 1 or NaN. Near m = 1 use EllipKm1 with p = 1-m to keep the digits m itself cannot hold.

func EllipKm1 

func EllipKm1(p float64) float64

EllipKm1 returns K(1 - p) for p >= 0 (SciPy's ellipkm1): the form to use near m = 1, where p carries digits that 1-p would lose.

func Erf 

func Erf(x float64) float64

Erf returns the error function erf(x) = 2/sqrt(pi) ∫₀ˣ exp(-t²) dt.

Special values: erf(±0) = ±0, erf(±Inf) = ±1, erf(NaN) = NaN.

Algorithm: for |x| <= 1/2, x E(x²) with E fitted to mpmath; beyond, 1 - erfc(|x|) in double-double with erfc from erfcMant.

func ErfInv 

func ErfInv(y float64) float64

ErfInv returns the inverse error function: ErfInv(Erf(x)) = x.

Special values: ErfInv(±1) = ±Inf, ErfInv(±0) = ±0; NaN for |y| > 1 or NaN.

Algorithm: Boost.Math's rational initial approximations (erf_inv_imp) followed by one Halley step whose residual erf(x) - y (or erfc(x) - (1-y) for |y| > 1/2, with 1-y exact) is computed in double-double, which takes the result to about half an ulp.

func Erfc 

func Erfc(x float64) float64

Erfc returns the complementary error function erfc(x) = 1 - erf(x).

Special values: erfc(+Inf) = 0, erfc(-Inf) = 2, erfc(NaN) = NaN; erfc underflows to 0 for x > 27.2.

Algorithm: for x >= 1/2, exp(-x²) erfcx(x) with x² in double-double and erfcx fitted to mpmath, multiplied before the single rounding; for |x| < 1/2, 1 - erf(x) in double-double; for x <= -1/2, 2 - erfc(-x).

func ErfcInv 

func ErfcInv(q float64) float64

ErfcInv returns the inverse complementary error function: ErfcInv(Erfc(x)) = x, for 0 <= q <= 2.

Special values: ErfcInv(0) = +Inf, ErfcInv(2) = -Inf, ErfcInv(1) = 0; NaN outside [0, 2] or for NaN.

Algorithm: as ErfInv, refining against erfc(x) - q directly, so that small q (down to the smallest subnormal) keep full relative accuracy.

func Erfcx 

func Erfcx(x float64) float64

Erfcx returns the scaled complementary error function erfcx(x) = exp(x²) erfc(x).

Special values: erfcx(+Inf) = 0, erfcx(-Inf) = +Inf, erfcx(NaN) = NaN; erfcx overflows to +Inf for x < -26.628.

Algorithm: for x >= 1/2 the fitted erfcx directly (no exponential at all); for |x| < 1/2, exp(x²)(1 - erf(x)) in double-double; for x <= -1/2, 2 exp(x²) - erfcx(-x).

func Expit 

func Expit(x float64) float64

Expit returns the logistic sigmoid 1/(1 + exp(-x)).

The exponential of the smaller magnitude side is formed as a double-double and the quotient taken in double-double, so the result is within about half an ulp everywhere, including the lower tail where it is exp(x) to full relative accuracy. Expit(-Inf) = 0, Expit(+Inf) = 1, NaN for NaN.

func Factorial 

func Factorial(n float64) float64

Factorial returns n! = Γ(n+1). Integers 0..170 come from a table of correctly rounded values; other non-negative n use the gamma function (with n+1 never rounded: Γ(1+n) from the core for n < 1, n Γ(n) above).

Special values: NaN for NaN; 0 for n < 0 (the convention of SciPy's factorial); +Inf for n > 170.6243769563027.

func Factorial2 

func Factorial2(n int) float64

Factorial2 returns the double factorial n!! = n (n-2) (n-4) ..., ending at 1 or 2. 0!! = (-1)!! = 1 by the usual convention (SciPy's factorial2 returns 0 for -1); n < -1 gives 0, as in SciPy. The product is carried in double-double and rounded once, so the result is correctly rounded (up to the 2^-100 double-double error) and exact whenever n!! < 2^53; it overflows to +Inf for n > 300.

func Gamma 

func Gamma(x float64) float64

Gamma returns the gamma function Γ(x) = ∫₀^∞ t^(x-1) e^(-t) dt, continued to negative x by Γ(x) = Γ(x+1)/x.

Special values: Γ(+Inf) = +Inf; Γ(±0) = ±Inf; Γ(x) = NaN for x a negative integer, x = -Inf or NaN; Γ overflows to +Inf for x > 171.6243769563027 and underflows to ±0 for x < -184.

Algorithm: on t in [-0.5, 1], Γ(1+t) = 1 + t(t-1)F(t) with F a piecewise Chebyshev interpolant fitted to mpmath (testdata/gen/fit.py); other x in (0, 10) are reduced to it by the recurrence Γ(x+1) = xΓ(x), whose factors x-k are exact and whose product is carried in double-double; x >= 10 uses Stirling's series (DLMF 5.11.1) with (x-1/2)log x - x formed in double-double and exponentiated by dd.Exp; x < -34 uses the reflection formula Γ(x) = -π / (x sin(πx) Γ(-x)) with sin(πx) reduced exactly. Integers are a table of correctly rounded factorials. Measured error against mpmath: docs/accuracy.md.

func GammaP 

func GammaP(a, x float64) float64

GammaP returns the regularised lower incomplete gamma function P(a, x) = γ(a, x)/Γ(a) = (1/Γ(a)) ∫₀ˣ t^(a-1) e^(-t) dt (DLMF 8.2.4).

Domain: a >= 0, x >= 0. Special values: P(a, 0) = 0 for a > 0; P(a, +Inf) = 1; P(0, x) = 1 for x > 0; P(+Inf, x) = 0 for finite x; NaN for a < 0, x < 0, a = x = 0, a = x = +Inf, or a NaN argument.

Algorithm: see the comment at the top of incgamma.go (DiDonato & Morris / Boost.Math method selection, Temme's expansion, double-double prefix).

func GammaPInv 

func GammaPInv(a, p float64) float64

GammaPInv returns the inverse of the regularised lower incomplete gamma function: the x >= 0 with P(a, x) = p.

Domain: a >= 0, 0 <= p <= 1. Special values: GammaPInv(a, 0) = 0, GammaPInv(a, 1) = +Inf; for a = 0 or a = +Inf only those two are defined (NaN otherwise); NaN for a < 0, p outside [0, 1] or a NaN argument. A root below the smallest subnormal returns 0.

Algorithm: DiDonato & Morris' initial approximation (via Boost.Math's find_inverse_gamma) refined by a safeguarded Newton iteration in log x; see the comment at the top of incgamma_inv.go.

func GammaQ 

func GammaQ(a, x float64) float64

GammaQ returns the regularised upper incomplete gamma function Q(a, x) = Γ(a, x)/Γ(a) = 1 - P(a, x) (DLMF 8.2.4), computed directly (not as 1 - P) where it is the smaller of the two, so it keeps its relative accuracy in the upper tail down to the subnormal range.

Domain and special values: as GammaP, with Q = 1 - P at the edges (Q(a, 0) = 1, Q(a, +Inf) = 0, Q(0, x) = 0, Q(+Inf, x) = 1).

Algorithm: as GammaP.

func GammaQInv 

func GammaQInv(a, q float64) float64

GammaQInv returns the inverse of the regularised upper incomplete gamma function: the x >= 0 with Q(a, x) = q.

Domain: a >= 0, 0 <= q <= 1. Special values: GammaQInv(a, 0) = +Inf, GammaQInv(a, 1) = 0; for a = 0 or a = +Inf only those two are defined (NaN otherwise); NaN for a < 0, q outside [0, 1] or a NaN argument.

Algorithm: as GammaPInv, matching Q(x) = q directly when q <= 1/2 (so small q keep full relative accuracy in x).

func GammaSign 

func GammaSign(x float64) float64

GammaSign returns the sign of Γ(x): +1 or -1, ±1 at ±0 (following the sign of the infinity Γ returns there), and NaN at the poles x = -1, -2, ..., at -Inf and for NaN.

func HurwitzZeta 

func HurwitzZeta(s, q float64) float64

HurwitzZeta returns the Hurwitz zeta function ζ(s, q) = Σ_{k>=0} (k+q)^-s for s > 1.

Special values: +Inf at s = 1 and at q = 0, -1, -2, ...; NaN for s < 1, for NaN arguments and for q < 0 with non-integer s (where (k+q)^-s is not real); ζ(s, +Inf) = 0.

Algorithm: Euler–Maclaurin summation in double-double (hurwitzPos). For q < 0 and integer s, the negative terms (q+k)^-s, k < ceil(-q), are summed explicitly with their signs and the rest by hurwitzPos.

func I0 

func I0(x float64) float64

I0 returns the modified Bessel function of the first kind of order zero, I0(x) = sum (x^2/4)^k / (k!)^2, for real x (even). I0(0) = 1, I0(±Inf) = +Inf, I0(NaN) = NaN; overflows to +Inf for |x| > ~713.99.

Algorithm: Boost.Math bessel_i0.hpp double-precision minimax approximations (polynomial in x^2/4 below 7.75, in 1/x above).

func I0e 

func I0e(x float64) float64

I0e returns I0(x) e^-|x|, defined for all real x; I0e(±Inf) = 0.

func I1 

func I1(x float64) float64

I1 returns the modified Bessel function of the first kind of order one, I1(x) = (x/2) sum (x^2/4)^k / (k! (k+1)!), for real x (odd). I1(0) = 0, I1(±Inf) = ±Inf, I1(NaN) = NaN; overflows for |x| > ~713.98.

Algorithm: Boost.Math bessel_i1.hpp double-precision minimax approximations.

func I1e 

func I1e(x float64) float64

I1e returns I1(x) e^-|x| (odd in x); I1e(±Inf) = ±0.

func In 

func In(n int, x float64) float64

In returns the modified Bessel function of the first kind of integer order n, Iv(float64(n), x).

func Iv 

func Iv(v, x float64) float64

Iv returns the modified Bessel function of the first kind of real order v, I_v(x) = sum_k (x/2)^(2k+v) / (k! Gamma(k+v+1)) (DLMF 10.25.2). For non-integer v it is real only for x >= 0 (NaN for x < 0); for integer v, I_v(-x) = (-1)^v I_v(x) and I_{-n} = I_n. Iv(v, 0) is 1 for v = 0, 0 for v > 0 or integer v, and ±Inf (the sign of 1/Gamma(v+1)) for non-integer v < 0; Iv(v, +Inf) = +Inf; Iv(+Inf, x) = 0, Iv(-Inf, x) = NaN; overflow gives +Inf.

Algorithm: Boost.Math bessel_ik.hpp (Temme's series and Steed's CF2 for K, CF1 and the Wronskian for I, power series for small x), carried exponentially scaled, with the large-x expansion DLMF 10.40.1 and the reflection I_{-v} = I_v + (2/pi) sin(v pi) K_v (DLMF 10.27.2). I0 and I1 serve the orders 0 and 1.

func Ive 

func Ive(v, x float64) float64

Ive returns Iv(v, x) e^-|x|, with the conventions of Iv; Ive(v, ±Inf) is ±0 (the limit), never NaN, and the result does not overflow.

Algorithm: as Iv, computed scaled throughout.

func J0 

func J0(x float64) float64

J0 returns the Bessel function of the first kind of order zero, 

J0(x) = (1/pi) * integral_0^pi cos(x sin t) dt,

for real x (even in x). J0(0) = 1, J0(±Inf) = 0 (the limit; SciPy returns NaN), J0(NaN) = NaN.

Algorithm: the power series for |x| <= 1.5; a Taylor series about the nearest zero of J0 or J1 (tabulated in double-double, coefficients from Bessel's equation; Harrison 2009) up to 100; Hankel's expansion in modulus-phase form with the phase in double-double and exact (Payne-Hanek) reduction beyond. The error stays within a few ulps of the value, also next to the zeros.

func J1 

func J1(x float64) float64

J1 returns the Bessel function of the first kind of order one, J1(x) = (1/pi) * integral_0^pi cos(t - x sin t) dt, for real x (odd in x). J1(0) = 0, J1(±Inf) = 0, J1(NaN) = NaN.

Algorithm: as J0 (power series for |x| <= 2, Taylor series about the zeros of J1 and of its derivative up to 100, modulus-phase Hankel expansion beyond).

func Jn 

func Jn(n int, x float64) float64

Jn returns the Bessel function of the first kind of integer order n, for any real x: J_{-n}(x) = (-1)^n J_n(x) and J_n(-x) = (-1)^n J_n(x). Jn(n, 0) = 0 for n != 0, Jn(n, ±Inf) = 0, Jn(n, NaN) = NaN.

Algorithm: J0 and J1 for n = 0, 1; otherwise as Jv: the power series for small x, the double-double modulus-phase Hankel expansion for large x, and Steed's method in double-double between (at the order n itself for x >= n, after backward recurrence from n for x < n). The forward recurrence from J0 and J1 used by Cephes and Boost is not used: it loses the relative accuracy next to the zeros of J_n.

func JnZeros 

func JnZeros(n int, nt int) []float64

JnZeros returns the first nt positive zeros of the Bessel function J_n, in increasing order. Since J_{-n} = (-1)^n J_n, a negative n gives the zeros of J_|n|. nt <= 0 gives an empty result. (x = 0, a zero of J_n for n != 0, is not counted, as in SciPy's jn_zeros.)

Algorithm: an initial estimate from McMahon's expansion (DLMF 10.21.19) for n = 0, and from Olver's uniform expansion (A&S 9.5.22, with the Airy zeros of A&S 10.4.94) for n >= 1, as in Boost.Math bessel_jy_zero.hpp; then Halley's iteration on Jn, with dJ_n/dx = J_{n-1} - (n/x) J_n and the second derivative from Bessel's equation, run to convergence.

func Jv 

func Jv(v, x float64) float64

Jv returns the Bessel function of the first kind of real order v, 

J_v(x) = sum_k (-1)^k (x/2)^(2k+v) / (k! Gamma(k+v+1))   (DLMF 10.2.2).

For integer v it is Jn(int(v), x). For non-integer v the function is real only for x >= 0: Jv(v, x < 0) = NaN. Jv(v, 0) is 0 for v > 0 and ±Inf (the sign of 1/Gamma(v+1)) for non-integer v < 0; Jv(v, +Inf) = 0; Jv(+Inf, x) = 0 and Jv(-Inf, x) = NaN; NaN in, NaN out.

Algorithm: after Boost.Math bessel_jy.hpp (N. M. Temme, J. Comput. Phys. 21, 1976, for Y at x <= 1/4; Steed's method with the continued fractions CF1 and CF2 of A. R. Barnett et al., Comput. Phys. Commun. 8, 1974, above; the power series for small x), carried in double-double; the large-x expansion in modulus-phase form (DLMF 10.18.17-18) with the phase in double-double and Payne-Hanek reduction. Negative orders use Steed's method directly when x >= |v|, else the reflection formula J_{-nu} = cos(nu pi) J_nu - sin(nu pi) Y_nu with exact sinPi/cosPi. Orders beyond 1e7 are supported only where the power series or the large-x expansion applies; elsewhere the result is NaN. The cost of Steed's method grows like x, so above x = 1e6 the large-x expansion is used even next to a zero.

func K0 

func K0(x float64) float64

K0 returns the modified Bessel function of the second kind of order zero, K0(x) = integral_0^Inf cos(x sinh t) dt, for x > 0. K0(0) = +Inf, K0(+Inf) = 0, K0(x < 0) = NaN, K0(NaN) = NaN; underflows to 0 beyond x ~ 745.

Algorithm: Boost.Math bessel_k0.hpp double-precision minimax approximations (log-term form on (0, 1], rational in 1/x above).

func K0e 

func K0e(x float64) float64

K0e returns K0(x) e^x for x > 0 (NaN for x < 0, +Inf at 0, 0 at +Inf).

func K1 

func K1(x float64) float64

K1 returns the modified Bessel function of the second kind of order one for x > 0. K1(0) = +Inf, K1(+Inf) = 0, K1(x < 0) = NaN, K1(NaN) = NaN; overflows to +Inf for x below about 5.6e-309 and underflows beyond ~745.

Algorithm: Boost.Math bessel_k1.hpp double-precision minimax approximations.

func K1e 

func K1e(x float64) float64

K1e returns K1(x) e^x for x > 0 (NaN for x < 0, +Inf at 0, 0 at +Inf).

func Kn 

func Kn(n int, x float64) float64

Kn returns the modified Bessel function of the second kind of integer order n, Kv(float64(n), x).

func Kv 

func Kv(v, x float64) float64

Kv returns the modified Bessel function of the second kind of real order v, K_v(x) = (pi/2) (I_{-v}(x) - I_v(x)) / sin(v pi) (DLMF 10.27.4), even in v, defined for x > 0. Kv(v, 0) = +Inf, Kv(v, +Inf) = 0, Kv(v, x < 0) = NaN, Kv(±Inf, x) = +Inf; NaN in, NaN out. Overflow gives +Inf and underflow 0.

Algorithm: K0/K1 for v = 0, 1; otherwise Temme's series (x <= 2) or Steed's CF2 (in double-double) for the fractional order and forward recurrence in double-double (Boost.Math bessel_ik.hpp), or the large-x expansion DLMF 10.40.2 where it converges.

func Kve 

func Kve(v, x float64) float64

Kve returns Kv(v, x) e^x, with the conventions of Kv (Kve(v, +Inf) = 0).

Algorithm: as Kv, computed scaled throughout.

func Log1pExp 

func Log1pExp(x float64) float64

Log1pExp returns log(1 + exp(x)), the softplus function, accurate in both tails: x + log1p(exp(-x)) for x > 0 (no overflow) and log1p(exp(x)) with exp(x) in double-double for x <= 0 (full relative accuracy as it tends to exp(x)). Log1pExp(-Inf) = 0, Log1pExp(+Inf) = +Inf, NaN for NaN.

func LogBeta 

func LogBeta(a, b float64) float64

LogBeta returns log|B(a, b)|, the logarithm of the absolute value of the beta function B(a, b) = Γ(a)Γ(b)/Γ(a+b).

Special values: NaN if either argument is NaN; +Inf when a or b is 0 or a negative integer (and a+b is not); -Inf when B(a, b) = 0 (a+b a non-positive integer, a and b not).

Algorithm: for positive arguments with max(a, b) >= 10, log Γ(b) - log Γ(a+b) is evaluated as a single cancellation-free expression from Stirling's series (logGammaDiff) and added to log Γ(a) in double-double; otherwise the three log-gammas are combined in double-double with the rounding of a+b corrected to first order. Negative arguments go through LogGamma with the same correction.

func LogExpit 

func LogExpit(x float64) float64

LogExpit returns log(Expit(x)) = -log(1 + exp(-x)) without underflow for large negative x (where it is x - log1p(exp(x)) ~ x) and with full relative accuracy for large positive x (where it is -exp(-x)).

func LogFactorial 

func LogFactorial(n float64) float64

LogFactorial returns log(n!) = log Γ(n+1) for n >= 0 without forming n+1 (which would round for non-integer n). NaN for n < 0 or NaN; +Inf for n = +Inf.

func LogGamma 

func LogGamma(x float64) float64

LogGamma returns log|Γ(x)|, the natural logarithm of the absolute value of the gamma function.

Special values: LogGamma(±Inf) = +Inf; LogGamma(x) = +Inf at the poles x = 0, -1, -2, ...; NaN for NaN; +Inf for x > 2.556348e305 (overflow).

Algorithm: near the zeros at x = 1 and x = 2 the result is log1p of Γ(1+t)-1 computed from the same fitted core as Gamma (gamma1pm1), so the relative accuracy holds as the result goes to zero; on (2.5, 10) it is the logarithm of Gamma's double-double value; from 10 Stirling's series (DLMF 5.11.1) is summed in double-double. For x < 0 it is log|Γ(1+t)| - log|t Π(x+k)| with exact factors for x > -34, and the reflection formula log π - log|x sin πx| - log Γ(-x) below.

func LogGammaSign 

func LogGammaSign(x float64) (float64, int)

LogGammaSign returns log|Γ(x)| and the sign of Γ(x) (+1 or -1), like math.Lgamma. At a pole the sign is +1, except for -0 where it is -1.

func LogNdtr 

func LogNdtr(x float64) float64

LogNdtr returns log Φ(x), the logarithm of the standard normal cdf, accurate in both tails: for x -> -Inf it is -x²/2 - log(-x) - ... with no underflow, for large x it is -Φ(-x) with full relative accuracy.

Special values: LogNdtr(-Inf) = -Inf, LogNdtr(+Inf) = 0, LogNdtr(NaN) = NaN.

Algorithm: for x <= -1/sqrt(2), log erfcx(z) - z² - log 2 with z = -x/sqrt(2) in double-double; for x >= 1/sqrt(2), log1p(-Φ(-x)) with Φ(-x) in double-double; in between, the logarithm of Φ in double-double.

func LogSoftmax 

func LogSoftmax(dst, x []float64) []float64

LogSoftmax writes x_i - LogSumExp(x) into dst (allocated when nil or too short) and returns it.

func LogSumExp 

func LogSumExp(x []float64) float64

LogSumExp returns log(Σ exp(x_i)) without overflow or underflow: the maximum m is factored out and log1p of Σ_{i≠argmax} exp(x_i - m) is taken in double-double (each x_i - m exact, each exponential to about 2^-54), so the result is within about half an ulp even when one term dominates. It returns -Inf for an empty slice or when every x_i is -Inf, +Inf if any x_i is +Inf, NaN if any is NaN.

func Logit 

func Logit(p float64) float64

Logit returns the log-odds log(p/(1-p)), the inverse of Expit, for 0 <= p <= 1: -Inf at 0, +Inf at 1, NaN outside or for NaN. 1-p and the quotient are carried in double-double and the logarithm keeps its relative accuracy as p/(1-p) approaches 1, so the result is accurate around p = 1/2 where it crosses zero.

func Ndtr 

func Ndtr(x float64) float64

Ndtr returns the standard normal cumulative distribution function Φ(x) = (1/sqrt(2π)) ∫_{-∞}^x exp(-t²/2) dt = erfc(-x/sqrt 2)/2.

Special values: Φ(-Inf) = 0, Φ(+Inf) = 1, Φ(NaN) = NaN.

Algorithm: the argument -x/sqrt(2) is carried in double-double into erfc (erfcMant), so the lower tail keeps full relative accuracy down to the subnormal range; the rounding of x/sqrt(2) that limits implementations computing erfc(-x*M_SQRT1_2) (an error growing like x² ulp) does not occur.

func Ndtri 

func Ndtri(p float64) float64

Ndtri returns the inverse of the standard normal cdf: Ndtri(Ndtr(x)) = x, for 0 <= p <= 1.

Special values: Ndtri(0) = -Inf, Ndtri(1) = +Inf, Ndtri(1/2) = 0; NaN outside [0, 1] or for NaN.

Algorithm: the guess -sqrt(2) erfcinv(2p) from Boost.Math's approximations, refined by one Halley step on Φ(x) - p with the relative residual from Ndtr's double-double evaluation (so p in the far lower tail, down to subnormals, keeps full relative accuracy). For p > 1/2, Ndtri(p) = -Ndtri(1-p) with 1-p exact.

func Perm 

func Perm(n, k int64) float64

Perm returns the number of k-permutations of n, n!/(n-k)! = n (n-1) ... (n-k+1), as a float64: 0 when k < 0, k > n or n < 0 (as SciPy's perm with exact=False).

Algorithm: the product in double-double for k <= 64 (exact below 2^53, correctly rounded above); beyond, exp(log Γ(n+1) - log Γ(n-k+1)) with the difference formed without cancellation (logGammaDiff) when n-k+1 is at least 10 and k <= n-k+1.

func Polygamma 

func Polygamma(n int, x float64) float64

Polygamma returns the polygamma function ψ⁽ⁿ⁾(x), the n-th derivative of the digamma function; Polygamma(0, x) = Digamma(x).

Special values: NaN for n < 0, NaN or x = -Inf; ψ⁽ⁿ⁾(+Inf) = 0 (n >= 1); at the poles x = 0, -1, -2, ... the value is (-1)^(n+1) Inf.

Algorithm: ψ⁽ⁿ⁾(x) = (-1)^(n+1) n! ζ(n+1, x) (DLMF 25.11.12), with ζ from the double-double Euler–Maclaurin summation and n! in double-double, the two combined as scaled values so that the product overflows only when the result does.

func Softmax 

func Softmax(dst, x []float64) []float64

Softmax writes exp(x_i) / Σ_j exp(x_j) into dst (allocated when nil or too short) and returns it; the maximum is subtracted first and the sum is accumulated in double-double. An empty x gives an empty result; if x holds NaN or +Inf the results are NaN (as SciPy's softmax).

func Xlog1py 

func Xlog1py(x, y float64) float64

Xlog1py returns x*log1p(y), defined as 0 when x = 0 and y is not NaN, as SciPy's xlog1py, with log1p(y) in double-double: one rounding.

func Xlogy 

func Xlogy(x, y float64) float64

Xlogy returns x*log(y), defined as 0 when x = 0 and y is not NaN (so that 0*log(0) terms of entropies vanish), as SciPy's xlogy. The product is formed with log(y) in double-double: one rounding.

func Y0 

func Y0(x float64) float64

Y0 returns the Bessel function of the second kind of order zero, Y0(x) = (2/pi)(log(x/2) + gamma) J0(x) + ..., defined for x > 0. Y0(0) = -Inf, Y0(+Inf) = 0, Y0(x < 0) = NaN, Y0(NaN) = NaN.

Algorithm: the series DLMF 10.8.2 for x <= 0.6; a Taylor series about the zeros of Y0 and Y1 up to 100 (the point chosen by |x - z|/z, since the logarithmic singularity at 0 bounds the radius of convergence); the modulus-phase Hankel expansion beyond.

func Y1 

func Y1(x float64) float64

Y1 returns the Bessel function of the second kind of order one, defined for x > 0. Y1(0) = -Inf, Y1(+Inf) = 0, Y1(x < 0) = NaN, Y1(NaN) = NaN; Y1(x) overflows to -Inf for x below about 3.5e-309.

Algorithm: as Y0 (series DLMF 10.8.1 for x <= 1.5).

func Yn 

func Yn(n int, x float64) float64

Yn returns the Bessel function of the second kind of integer order n for x >= 0, with Y_{-n}(x) = (-1)^n Y_n(x). Yn(n, 0) = -Inf (+Inf for odd negative n), Yn(n, x < 0) = NaN, Yn(n, +Inf) = 0, Yn(n, NaN) = NaN; large negative values overflow to -Inf.

Algorithm: Y0 and Y1 for n = 0, 1; otherwise as Yv (Temme's series and forward recurrence for small x, Steed's method in double-double, the modulus-phase Hankel expansion for large x).

func Yv 

func Yv(v, x float64) float64

Yv returns the Bessel function of the second kind of real order v, Y_v(x) = (J_v(x) cos(v pi) - J_{-v}(x)) / sin(v pi) (DLMF 10.2.3), and its limit for integer v (Yn). It is real only for x >= 0: Yv(v, x < 0) = NaN. Yv(v, 0) = -Inf for v > 0 or v = 0; for v < 0 it is the limit -Inf*cos(v pi) (0 at half-integers, where Y_{-n-1/2} = ±J_{n+1/2}); Yv(v, +Inf) = 0; Yv(+Inf, x > 0) = -Inf, Yv(-Inf, x) = NaN; NaN in, NaN out. Large negative values overflow to -Inf.

Algorithm: as Jv; for small x and non-integer v, the two power series of the definition where the second is negligible.

func Zeta 

func Zeta(s float64) float64

Zeta returns the Riemann zeta function ζ(s) for real s.

Special values: ζ(1) = +Inf, ζ(+Inf) = 1, ζ(0) = -1/2, ζ(-2k) = 0 for integers k >= 1; NaN for NaN and -Inf. ζ overflows to ±Inf for s < -260 (roughly).

Algorithm: on [-1, 6], ζ(s) = 1/(s-1) + R(s) with R entire and fitted to mpmath (zz_zetaPieces_coeffs.go), so the pole is exact; on (6, 64), 1 + g(s) 2^-s with g fitted to mpmath; for s < -1 the reflection formula ζ(s) = 2^s π^(s-1) sin(πs/2) Γ(1-s) ζ(1-s) (DLMF 25.4.1) with the powers and the gamma function combined in double-double before one exponential and sin(πs/2) reduced exactly, so the trivial zeros are exact and nearby values keep their relative accuracy.

Types

This section is empty.

Source Files

View all Source files

Directories

Path Synopsis
internal
dd
Package dd is a minimal double-double ("float128-lite") arithmetic used where a special function must carry more than 53 bits through an intermediate step: the argument of an exponential, a long product, a cancellation.
 Package dd is a minimal double-double ("float128-lite") arithmetic used where a special function must carry more than 53 bits through an intermediate step: the argument of an exponential, a long product, a cancellation.
Package nd applies the special functions of package special element-wise to go-ndarray arrays (github.com/go-ndarray/ndarray).
 Package nd applies the special functions of package special element-wise to go-ndarray arrays (github.com/go-ndarray/ndarray).

Jump to

Keyboard shortcuts

? : This menu
/ : Search site
f or F : Jump to
y or Y : Canonical URL
go.dev uses cookies from Google to deliver and enhance the quality of its services and to analyze traffic. Learn more.