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

special

package module
v0.2.0 Latest Latest
Warning

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

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

README

go-numerics/special

special

Docs pkg.go.dev License

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.

🌐 Website · 📚 Documentation

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

Each row gives the Go function and, in brackets, the scipy.special function it is measured against — the pairing comes from the harness (bench/scipy_server.py and testdata/gen/fam_*.py), not from transliterating a name, which is why Log1pExp is scipy's softplus and LogGamma is gammaln. 63 of the 67 carry such a pairing; — marks the four that do not, and the notes under the table say what they are. readme_scipy_test.go reads the harness and fails if this column drifts from it.

family functions (scipy equivalent)
gamma Gamma (gamma), LogGamma (gammaln), LogGammaSign (—), GammaSign (gammasgn), Digamma (digamma), Polygamma (polygamma), Beta (beta), LogBeta (betaln)
incomplete gamma / beta GammaP (gammainc), GammaQ (gammaincc), GammaPInv (gammaincinv), GammaQInv (gammainccinv), BetaInc (betainc), BetaIncC (betaincc), BetaIncInv (betaincinv), BetaIncCInv (betainccinv)
error function, normal Erf (erf), Erfc (erfc), Erfcx (erfcx), ErfInv (erfinv), ErfcInv (erfcinv), Ndtr (ndtr), Ndtri (ndtri), LogNdtr (log_ndtr)
combinatorics Comb (comb), Perm (perm), Binom (binom), Factorial (factorial), Factorial2 (factorial2), LogFactorial (gammaln(n+1))
logistic, log helpers Expit (expit), LogExpit (log_expit), Logit (logit), Log1pExp (softplus), Xlogy (xlogy), Xlog1py (xlog1py), LogSumExp (logsumexp), Softmax (softmax), LogSoftmax (log_softmax)
zeta Zeta (zeta), HurwitzZeta (zeta)
Bessel J0 (j0), J1 (j1), Jn (jv), Jv (jv), Y0 (y0), Y1 (y1), Yn (yn), Yv (yv), I0 (i0), I1 (i1), In (—), Iv (iv), K0 (k0), K1 (k1), Kn (kn), Kv (kv), I0e (i0e), I1e (i1e), Ive (ive), K0e (k0e), K1e (k1e), Kve (kve), JnZeros (jn_zeros)
elliptic EllipK (ellipk), EllipKm1 (ellipkm1), EllipE (ellipe)

The four without a pairing, and why:

  • GammaSign and Kn — scipy does have gammasgn and kn; this package simply measures neither against them. The bracketed name is scipy's, not a measurement.
  • In — scipy has no in, and cannot: it is a Python keyword. Iv at integer order is the equivalent, and In(n, x) is Iv(float64(n), x). Kn/Kv and Jn/Jv stand in the same relation, which is why Jn is bracketed jv: scipy's jn is a deprecated alias of jv.
  • LogGammaSign — no scipy counterpart in any spelling; scipy splits the two results across gammaln and gammasgn.

The 20 …To slice kernels — GammaTo, ErfTo, XlogyTo and the rest — are the batch form of the function they are named after and compute the same thing, so they share its scipy equivalent rather than earning rows here. One exception is worth knowing: TrigammaTo is exported while Trigamma is not, so that kernel has no scalar form in the public API. It is psi'(x), which scipy spells polygamma(1, x) — the same function Polygamma already provides.

Zeta and HurwitzZeta share one scipy name, zeta(x, q=None), which takes the second argument optionally. LogFactorial is bracketed as a call rather than a name because the function is gammaln with its argument shifted, which is what the harness measures it against.

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 BetaTo added in v0.2.0 

func BetaTo(ζdst, ζa, ζb []float64, ζsa, ζsb int)

BetaTo sets dst[i] = Beta(a[i*sa], b[i*sb]) for i < len(dst) (a stride of 0 repeats a single value), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 DigammaTo added in v0.2.0 

func DigammaTo(ζdst, ζx []float64)

DigammaTo sets dst[i] = Digamma(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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: for |y| <= 1/2, y R(y²) with R fitted to mpmath; beyond, erfcinv(1 - |y|) (1 - |y| exact) as X(sqrt(-log(1 - |y|))) with X fitted to mpmath on pieces, the logarithm and the square root carried beyond float64; about half an ulp.

func ErfInvTo added in v0.2.0 

func ErfInvTo(ζdst, ζx []float64)

ErfInvTo sets dst[i] = ErfInv(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

func ErfTo added in v0.2.0 

func ErfTo(ζdst, ζx []float64)

ErfTo sets dst[i] = Erf(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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, with X evaluated at sqrt(-log q) directly, so that small q (down to the smallest subnormal) keep full relative accuracy.

func ErfcInvTo added in v0.2.0 

func ErfcInvTo(ζdst, ζx []float64)

ErfcInvTo sets dst[i] = ErfcInv(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

func ErfcTo added in v0.2.0 

func ErfcTo(ζdst, ζx []float64)

ErfcTo sets dst[i] = Erfc(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 ErfcxTo added in v0.2.0 

func ErfcxTo(ζdst, ζx []float64)

ErfcxTo sets dst[i] = Erfcx(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

func Expit 

func Expit(x float64) float64

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

e = exp(-|x|) comes from the float64 kernel as hi + lo (about 2^-68 relative); 1 + e is formed exactly (Fast2Sum) and the quotient 1/(1+e) or e/(1+e) is corrected by its FMA residual, so the result is within about half an ulp everywhere, including the lower tail where it is exp(x) to full relative accuracy. Below -660, where the low part of e would be subnormal, the double-double path is used. Expit(-Inf) = 0, Expit(+Inf) = 1, NaN for NaN.

func ExpitTo added in v0.2.0 

func ExpitTo(ζdst, ζx []float64)

ExpitTo sets dst[i] = Expit(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 Halley 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 GammaTo added in v0.2.0 

func GammaTo(ζdst, ζx []float64)

GammaTo sets dst[i] = Gamma(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 HurwitzZetaTo added in v0.2.0 

func HurwitzZetaTo(ζdst, ζa, ζb []float64, ζsa, ζsb int)

HurwitzZetaTo sets dst[i] = HurwitzZeta(a[i*sa], b[i*sb]) for i < len(dst) (a stride of 0 repeats a single value), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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: the power series below 1/8; up to 512, e^x times I0(x) e^-x from degree-11 interpolants on 8 pieces per octave of x (zz_besselI0eOct_coeffs.go, fitted to mpmath), with e^x from dd.ExpMant and the product rounded once; DLMF 10.40.1 beyond.

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: below 100, one polynomial per cell about the zeros of J0, of J1 and the midpoints between them, with the value and the derivative at the centre in double-double (Harrison 2009); 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 (cell polynomials 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's log-term form on (0, 1]; up to 512, e^-x times K0(x) e^x from degree-11 interpolants on 8 pieces per octave (fitted to mpmath), the product rounded once; DLMF 10.40.2 beyond.

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)) for x <= 0 (full relative accuracy as it tends to exp(x)). e = exp(-|x|) comes from the float64 kernel as hi + lo, w = 1 + e is formed exactly and log w = log(w.hi) + w.lo/w.hi from the log kernel (2^-68 absolute); x is added with TwoSum before the single final rounding. Log1pExp(-Inf) = 0, Log1pExp(+Inf) = +Inf, NaN for NaN.

func Log1pExpTo added in v0.2.0 

func Log1pExpTo(ζdst, ζx []float64)

Log1pExpTo sets dst[i] = Log1pExp(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 LogBetaTo added in v0.2.0 

func LogBetaTo(ζdst, ζa, ζb []float64, ζsa, ζsb int)

LogBetaTo sets dst[i] = LogBeta(a[i*sa], b[i*sb]) for i < len(dst) (a stride of 0 repeats a single value), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 LogGammaTo added in v0.2.0 

func LogGammaTo(ζdst, ζx []float64)

LogGammaTo sets dst[i] = LogGamma(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 LogNdtrTo added in v0.2.0 

func LogNdtrTo(ζdst, ζx []float64)

LogNdtrTo sets dst[i] = LogNdtr(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 is formed exactly (TwoSum), the quotient p/(1-p) is carried as hi + lo through its FMA residual, and log(hi) comes from the log kernel, whose error is relative to hi - 1 near 1: the result is accurate around p = 1/2 where it crosses zero.

func LogitTo added in v0.2.0 

func LogitTo(ζdst, ζx []float64)

LogitTo sets dst[i] = Logit(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 NdtrTo added in v0.2.0 

func NdtrTo(ζdst, ζx []float64)

NdtrTo sets dst[i] = Ndtr(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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: -sqrt(2) erfcinv(2p) with erfcinv from ErfcInv's fitted approximations as hi + lo and sqrt(2) as a double-double, multiplied before the single rounding (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 NdtriTo added in v0.2.0 

func NdtriTo(ζdst, ζx []float64)

NdtriTo sets dst[i] = Ndtri(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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 TrigammaTo added in v0.2.0 

func TrigammaTo(ζdst, ζx []float64)

TrigammaTo sets dst[i] = Polygamma(1, x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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. log(y) comes from the log kernel as hi + lo (2^-68 absolute) and the product is rounded once by an FMA. The kernel is written out here (the same operations as logHL, without its subnormal branch): at a few ns per call, the spill of a call under Go's ABI is a tenth of the cost.

func XlogyTo added in v0.2.0 

func XlogyTo(ζdst, ζa, ζb []float64, ζsa, ζsb int)

XlogyTo sets dst[i] = Xlogy(a[i*sa], b[i*sb]) for i < len(dst) (a stride of 0 repeats a single value), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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; cell polynomials about the zeros of Y0 and Y1 and the midpoints between them up to 100 (cells split further near the logarithmic singularity at 0); 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.

func ZetaTo added in v0.2.0 

func ZetaTo(ζdst, ζx []float64)

ZetaTo sets dst[i] = Zeta(x[i]) for i < len(dst) (x at least as long), bit-identical to the scalar function, in one loop with no call per element: the loop of the nd functions.

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.