Integrate BACKBEAT SDK and resolve KACHING license validation

Major integrations and fixes:
- Added BACKBEAT SDK integration for P2P operation timing
- Implemented beat-aware status tracking for distributed operations
- Added Docker secrets support for secure license management
- Resolved KACHING license validation via HTTPS/TLS
- Updated docker-compose configuration for clean stack deployment
- Disabled rollback policies to prevent deployment failures
- Added license credential storage (CHORUS-DEV-MULTI-001)

Technical improvements:
- BACKBEAT P2P operation tracking with phase management
- Enhanced configuration system with file-based secrets
- Improved error handling for license validation
- Clean separation of KACHING and CHORUS deployment stacks

🤖 Generated with [Claude Code](https://claude.ai/code)

Co-Authored-By: Claude <noreply@anthropic.com>

vendor/gonum.org/v1/gonum/mathext/README.md

@@ -0,0 +1,6 @@
+# mathext
+
+[![go.dev reference](https://pkg.go.dev/badge/gonum.org/v1/gonum/mathext)](https://pkg.go.dev/gonum.org/v1/gonum/mathext)
+[![GoDoc](https://godocs.io/gonum.org/v1/gonum/mathext?status.svg)](https://godocs.io/gonum.org/v1/gonum/mathext)
+
+Package mathext implements basic elementary functions not included in the Go standard library.

vendor/gonum.org/v1/gonum/mathext/airy.go

@@ -0,0 +1,41 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import "gonum.org/v1/gonum/mathext/internal/amos"
+
+// AiryAi returns the value of the Airy function at z. The Airy function here,
+// Ai(z), is one of the two linearly independent solutions to
+//
+//	y′′ - y*z = 0.
+//
+// See http://mathworld.wolfram.com/AiryFunctions.html for more detailed information.
+func AiryAi(z complex128) complex128 {
+	// id specifies the order of the derivative to compute,
+	// 0 for the function itself and 1 for the derivative.
+	// kode specifies the scaling option. See the function
+	// documentation for the exact behavior.
+	id := 0
+	kode := 1
+	air, aii, _, _ := amos.Zairy(real(z), imag(z), id, kode)
+	return complex(air, aii)
+}
+
+// AiryAiDeriv returns the value of the derivative of the Airy function at z. The
+// Airy function here, Ai(z), is one of the two linearly independent solutions to
+//
+//	y′′ - y*z = 0.
+//
+// See http://mathworld.wolfram.com/AiryFunctions.html for more detailed information.
+func AiryAiDeriv(z complex128) complex128 {
+	// id specifies the order of the derivative to compute,
+	// 0 for the function itself and 1 for the derivative.
+	// kode specifies the scaling option. See the function
+	// documentation for the exact behavior.
+	id := 1
+	kode := 1
+	air, aii, _, _ := amos.Zairy(real(z), imag(z), id, kode)
+	return complex(air, aii)
+}

vendor/gonum.org/v1/gonum/mathext/beta.go

@@ -0,0 +1,40 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import "gonum.org/v1/gonum/mathext/internal/gonum"
+
+// Beta returns the value of the complete beta function B(a, b). It is defined as
+//
+//	Γ(a)Γ(b) / Γ(a+b)
+//
+// Special cases are:
+//
+//	B(a,b) returns NaN if a or b is Inf
+//	B(a,b) returns NaN if a and b are 0
+//	B(a,b) returns NaN if a or b is NaN
+//	B(a,b) returns NaN if a or b is < 0
+//	B(a,b) returns +Inf if a xor b is 0.
+//
+// See http://mathworld.wolfram.com/BetaFunction.html for more detailed informations.
+func Beta(a, b float64) float64 {
+	return gonum.Beta(a, b)
+}
+
+// Lbeta returns the natural logarithm of the complete beta function B(a,b).
+// Lbeta is defined as:
+//
+//	Ln(Γ(a)Γ(b)/Γ(a+b))
+//
+// Special cases are:
+//
+//	Lbeta(a,b) returns NaN if a or b is Inf
+//	Lbeta(a,b) returns NaN if a and b are 0
+//	Lbeta(a,b) returns NaN if a or b is NaN
+//	Lbeta(a,b) returns NaN if a or b is < 0
+//	Lbeta(a,b) returns +Inf if a xor b is 0.
+func Lbeta(a, b float64) float64 {
+	return gonum.Lbeta(a, b)
+}

vendor/gonum.org/v1/gonum/mathext/betainc.go

@@ -0,0 +1,33 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import "gonum.org/v1/gonum/mathext/internal/cephes"
+
+// RegIncBeta returns the value of the regularized incomplete beta function
+// I(x;a,b). It is defined as
+//
+//	I(x;a,b) = B(x;a,b) / B(a,b)
+//	         = Γ(a+b) / (Γ(a)*Γ(b)) * int_0^x u^(a-1) * (1-u)^(b-1) du.
+//
+// The domain of definition is 0 <= x <= 1, and the parameters a and b must be positive.
+// For other values of x, a, and b RegIncBeta will panic.
+func RegIncBeta(a, b float64, x float64) float64 {
+	return cephes.Incbet(a, b, x)
+}
+
+// InvRegIncBeta computes the inverse of the regularized incomplete beta function.
+// It returns the x for which
+//
+//	y = I(x;a,b)
+//
+// The domain of definition is 0 <= y <= 1, and the parameters a and b must be
+// positive. For other values of x, a, and b InvRegIncBeta will panic.
+func InvRegIncBeta(a, b float64, y float64) float64 {
+	if y < 0 || 1 < y {
+		panic("mathext: parameter out of range")
+	}
+	return cephes.Incbi(a, b, y)
+}

vendor/gonum.org/v1/gonum/mathext/digamma.go

@@ -0,0 +1,45 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import (
+	"math"
+)
+
+// Digamma returns the logorithmic derivative of the gamma function at x.
+//
+//	ψ(x) = d/dx (Ln (Γ(x)).
+func Digamma(x float64) float64 {
+	// This is adapted from
+	// http://web.science.mq.edu.au/~mjohnson/code/digamma.c
+	var result float64
+	switch {
+	case math.IsNaN(x), math.IsInf(x, 1):
+		return x
+	case math.IsInf(x, -1):
+		return math.NaN()
+	case x == 0:
+		return math.Copysign(math.Inf(1), -x)
+	case x < 0:
+		if x == math.Floor(x) {
+			return math.NaN()
+		}
+		// Reflection formula, http://dlmf.nist.gov/5.5#E4
+		_, r := math.Modf(x)
+		result = -math.Pi / math.Tan(math.Pi*r)
+		x = 1 - x
+	}
+	for ; x < 7; x++ {
+		// Recurrence relation, http://dlmf.nist.gov/5.5#E2
+		result -= 1 / x
+	}
+	x -= 0.5
+	xx := 1 / x
+	xx2 := xx * xx
+	xx4 := xx2 * xx2
+	// Asymptotic expansion, http://dlmf.nist.gov/5.11#E2
+	result += math.Log(x) + (1.0/24.0)*xx2 - (7.0/960.0)*xx4 + (31.0/8064.0)*xx4*xx2 - (127.0/30720.0)*xx4*xx4
+	return result
+}

vendor/gonum.org/v1/gonum/mathext/doc.go

@@ -0,0 +1,7 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package mathext implements special math functions not implemented by the
+// Go standard library.
+package mathext // import "gonum.org/v1/gonum/mathext"

vendor/gonum.org/v1/gonum/mathext/ell_carlson.go

@@ -0,0 +1,168 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import (
+	"math"
+)
+
+// EllipticRF computes the symmetric elliptic integral R_F(x,y,z):
+//
+//	R_F(x,y,z) = (1/2)\int_{0}^{\infty}{1/s(t)} dt,
+//	s(t) = \sqrt{(t+x)(t+y)(t+z)}.
+//
+// The arguments x, y, z must satisfy the following conditions, otherwise the function returns math.NaN():
+//
+//	0 ≤ x,y,z ≤ upper,
+//	lower ≤ x+y,y+z,z+x,
+//
+// where:
+//
+//	lower = 5/(2^1022) = 1.112536929253601e-307,
+//	upper = (2^1022)/5 = 8.988465674311580e+306.
+//
+// The definition of the symmetric elliptic integral R_F can be found in NIST
+// Digital Library of Mathematical Functions (http://dlmf.nist.gov/19.16.E1).
+func EllipticRF(x, y, z float64) float64 {
+	// The original Fortran code was published as Algorithm 577 in ACM TOMS (http://doi.org/10.1145/355958.355970).
+	// This code is also available as a part of SLATEC Common Mathematical Library (http://netlib.org/slatec/index.html). Later, Carlson described
+	// an improved version in http://dx.doi.org/10.1007/BF02198293 (also available at https://arxiv.org/abs/math/9409227).
+	const (
+		lower = 5.0 / (1 << 256) / (1 << 256) / (1 << 256) / (1 << 254) // 5*2^-1022
+		upper = 1 / lower
+		tol   = 1.2674918778210762260320167734407048051023273568443e-02 // (3ε)^(1/8)
+	)
+	if x < 0 || y < 0 || z < 0 || math.IsNaN(x) || math.IsNaN(y) || math.IsNaN(z) {
+		return math.NaN()
+	}
+	if upper < x || upper < y || upper < z {
+		return math.NaN()
+	}
+	if x+y < lower || y+z < lower || z+x < lower {
+		return math.NaN()
+	}
+
+	A0 := (x + y + z) / 3
+	An := A0
+	Q := math.Max(math.Max(math.Abs(A0-x), math.Abs(A0-y)), math.Abs(A0-z)) / tol
+	xn, yn, zn := x, y, z
+	mul := 1.0
+
+	for Q >= mul*math.Abs(An) {
+		xnsqrt, ynsqrt, znsqrt := math.Sqrt(xn), math.Sqrt(yn), math.Sqrt(zn)
+		lambda := xnsqrt*ynsqrt + ynsqrt*znsqrt + znsqrt*xnsqrt
+		An = (An + lambda) * 0.25
+		xn = (xn + lambda) * 0.25
+		yn = (yn + lambda) * 0.25
+		zn = (zn + lambda) * 0.25
+		mul *= 4
+	}
+
+	X := (A0 - x) / (mul * An)
+	Y := (A0 - y) / (mul * An)
+	Z := -(X + Y)
+	E2 := X*Y - Z*Z
+	E3 := X * Y * Z
+
+	// http://dlmf.nist.gov/19.36.E1
+	return (1 - 1/10.0*E2 + 1/14.0*E3 + 1/24.0*E2*E2 - 3/44.0*E2*E3 - 5/208.0*E2*E2*E2 + 3/104.0*E3*E3 + 1/16.0*E2*E2*E3) / math.Sqrt(An)
+}
+
+// EllipticRD computes the symmetric elliptic integral R_D(x,y,z):
+//
+//	R_D(x,y,z) = (1/2)\int_{0}^{\infty}{1/(s(t)(t+z))} dt,
+//	s(t) = \sqrt{(t+x)(t+y)(t+z)}.
+//
+// The arguments x, y, z must satisfy the following conditions, otherwise the function returns math.NaN():
+//
+//	0 ≤ x,y ≤ upper,
+//	lower ≤ z ≤ upper,
+//	lower ≤ x+y,
+//
+// where:
+//
+//	lower = (5/(2^1022))^(1/3) = 4.809554074311679e-103,
+//	upper = ((2^1022)/5)^(1/3) = 2.079194837087086e+102.
+//
+// The definition of the symmetric elliptic integral R_D can be found in NIST
+// Digital Library of Mathematical Functions (http://dlmf.nist.gov/19.16.E5).
+func EllipticRD(x, y, z float64) float64 {
+	// The original Fortran code was published as Algorithm 577 in ACM TOMS (http://doi.org/10.1145/355958.355970).
+	// This code is also available as a part of SLATEC Common Mathematical Library (http://netlib.org/slatec/index.html). Later, Carlson described
+	// an improved version in http://dx.doi.org/10.1007/BF02198293 (also available at https://arxiv.org/abs/math/9409227).
+	const (
+		lower = 4.8095540743116787026618007863123676393525016818363e-103 // (5*2^-1022)^(1/3)
+		upper = 1 / lower
+		tol   = 9.0351169339315770474760122547068324993857488849382e-03 // (ε/5)^(1/8)
+	)
+	if x < 0 || y < 0 || math.IsNaN(x) || math.IsNaN(y) || math.IsNaN(z) {
+		return math.NaN()
+	}
+	if upper < x || upper < y || upper < z {
+		return math.NaN()
+	}
+	if x+y < lower || z < lower {
+		return math.NaN()
+	}
+
+	A0 := (x + y + 3*z) / 5
+	An := A0
+	Q := math.Max(math.Max(math.Abs(A0-x), math.Abs(A0-y)), math.Abs(A0-z)) / tol
+	xn, yn, zn := x, y, z
+	mul, s := 1.0, 0.0
+
+	for Q >= mul*math.Abs(An) {
+		xnsqrt, ynsqrt, znsqrt := math.Sqrt(xn), math.Sqrt(yn), math.Sqrt(zn)
+		lambda := xnsqrt*ynsqrt + ynsqrt*znsqrt + znsqrt*xnsqrt
+		s += 1 / (mul * znsqrt * (zn + lambda))
+		An = (An + lambda) * 0.25
+		xn = (xn + lambda) * 0.25
+		yn = (yn + lambda) * 0.25
+		zn = (zn + lambda) * 0.25
+		mul *= 4
+	}
+
+	X := (A0 - x) / (mul * An)
+	Y := (A0 - y) / (mul * An)
+	Z := -(X + Y) / 3
+	E2 := X*Y - 6*Z*Z
+	E3 := (3*X*Y - 8*Z*Z) * Z
+	E4 := 3 * (X*Y - Z*Z) * Z * Z
+	E5 := X * Y * Z * Z * Z
+
+	// http://dlmf.nist.gov/19.36.E2
+	return (1-3/14.0*E2+1/6.0*E3+9/88.0*E2*E2-3/22.0*E4-9/52.0*E2*E3+3/26.0*E5-1/16.0*E2*E2*E2+3/40.0*E3*E3+3/20.0*E2*E4+45/272.0*E2*E2*E3-9/68.0*(E3*E4+E2*E5))/(mul*An*math.Sqrt(An)) + 3*s
+}
+
+// EllipticF computes the Legendre's elliptic integral of the 1st kind F(phi,m), 0≤m<1:
+//
+//	F(\phi,m) = \int_{0}^{\phi} 1 / \sqrt{1-m\sin^2(\theta)} d\theta
+//
+// Legendre's elliptic integrals can be expressed as symmetric elliptic integrals, in this case:
+//
+//	F(\phi,m) = \sin\phi R_F(\cos^2\phi,1-m\sin^2\phi,1)
+//
+// The definition of F(phi,k) where k=sqrt(m) can be found in NIST Digital Library of Mathematical
+// Functions (http://dlmf.nist.gov/19.2.E4).
+func EllipticF(phi, m float64) float64 {
+	s, c := math.Sincos(phi)
+	return s * EllipticRF(c*c, 1-m*s*s, 1)
+}
+
+// EllipticE computes the Legendre's elliptic integral of the 2nd kind E(phi,m), 0≤m<1:
+//
+//	E(\phi,m) = \int_{0}^{\phi} \sqrt{1-m\sin^2(\theta)} d\theta
+//
+// Legendre's elliptic integrals can be expressed as symmetric elliptic integrals, in this case:
+//
+//	E(\phi,m) = \sin\phi R_F(\cos^2\phi,1-m\sin^2\phi,1)-(m/3)\sin^3\phi R_D(\cos^2\phi,1-m\sin^2\phi,1)
+//
+// The definition of E(phi,k) where k=sqrt(m) can be found in NIST Digital Library of Mathematical
+// Functions (http://dlmf.nist.gov/19.2.E5).
+func EllipticE(phi, m float64) float64 {
+	s, c := math.Sincos(phi)
+	x, y := c*c, 1-m*s*s
+	return s * (EllipticRF(x, y, 1) - (m/3)*s*s*EllipticRD(x, y, 1))
+}

vendor/gonum.org/v1/gonum/mathext/ell_complete.go

@@ -0,0 +1,355 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import (
+	"math"
+)
+
+// CompleteK computes the complete elliptic integral of the 1st kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1].
+//
+//	K(m) = \int_{0}^{π/2} 1/{\sqrt{1-m{\sin^2θ}}} dθ
+func CompleteK(m float64) float64 {
+	// Reference:
+	// Toshio Fukushima, Precise and fast computation of complete elliptic integrals
+	// by piecewise minimax rational function approximation,
+	// Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76.
+	// https://doi.org/10.1016/j.cam.2014.12.038
+	// Original Fortran code available at:
+	// https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation
+	if m < 0 || 1 < m || math.IsNaN(m) {
+		return math.NaN()
+	}
+
+	mc := 1 - m
+
+	if mc > 0.592990 {
+		t := 2.45694208987494165*mc - 1.45694208987494165
+		t2 := t * t
+		p := ((3703.75266375099019 + t2*(2744.82029097576810+t2*36.2381612593459565)) + t*(5462.47093231923466+t2*(543.839017382099411+t2*0.393188651542789784)))
+		q := ((2077.94377067058435 + t2*(1959.05960044399275+t2*43.5464368440078942)) + t*(3398.00069767755460+t2*(472.794455487539279+t2)))
+		return p / q
+	}
+	if mc > 0.350756 {
+		t := 4.12823963605439369*mc - 1.44800482178389491
+		t2 := t * t
+		p := ((4264.28203103974630 + t2*(3214.59187442783167+t2*43.2589626155454993)) + t*(6341.90978213264024+t2*(642.790566685354573+t2*0.475223892294445943)))
+		q := ((2125.06914237062279 + t2*(2006.03187933518870+t2*44.1848041560412224)) + t*(3479.95663350926514+t2*(482.900172581418890+t2)))
+		return p / q
+	}
+	if mc > 0.206924 {
+		t := 6.95255575949719117*mc - 1.43865064797819679
+		t2 := t * t
+		p := ((4870.25402224986382 + t2*(3738.29369283392307+t2*51.3609902253065926)) + t*(7307.18826377416591+t2*(754.928587580583704+t2*0.571948962277566451)))
+		q := ((2172.51745704102287 + t2*(2056.13612019430497+t2*44.9026847057686146)) + t*(3565.04737778032566+t2*(493.962405117599400+t2)))
+		return p / q
+	}
+	if mc > 0.121734 {
+		t := 11.7384669562155183*mc - 1.42897053644793990
+		t2 := t * t
+		p := ((5514.8512729127464 + t2*(4313.60788246750934+t2*60.598720224393536)) + t*(8350.4595896779631+t2*(880.27903031894216+t2*0.68504458747933773)))
+		q := ((2218.41682813309737 + t2*(2107.97379949034285+t2*45.6911096775045314)) + t*(3650.41829123846319+t2*(505.74295207655096+t2)))
+		return p / q
+	}
+	if mc > 0.071412 {
+		t := 19.8720241643813839*mc - 1.41910098962680339
+		t2 := t * t
+		p := ((6188.8743957372448 + t2*(4935.41351498551527+t2*70.981049144472361)) + t*(9459.3331440432847+t2*(1018.21910476032105+t2*0.81599895108245948)))
+		q := ((2260.73112539748448 + t2*(2159.68721749761492+t2*46.5298955058476510)) + t*(3732.66955095581621+t2*(517.86964191812384+t2)))
+		return p / q
+	}
+	if mc > 0.041770 {
+		t := 33.7359152553808785*mc - 1.40914918021725929
+		t2 := t * t
+		p := ((6879.5170681289562 + t2*(5594.8381504799829+t2*82.452856129147838)) + t*(10615.0836403687221+t2*(1167.26108955935542+t2*0.96592719058503951)))
+		q := ((2296.88303450660439 + t2*(2208.74949754945558+t2*47.3844470709989137)) + t*(3807.37745652028212+t2*(529.79651353072921+t2)))
+		return p / q
+	}
+	if mc > 0.024360 {
+		t := 57.4382538770821367*mc - 1.39919586444572085
+		t2 := t * t
+		p := ((7570.6827538712100 + t2*(6279.2661370014890+t2*94.886883830605940)) + t*(11792.9392624454532+t2*(1325.01058966228180+t2*1.13537029594409690)))
+		q := ((2324.04824540459984 + t2*(2252.22250562615338+t2*48.2089280211559345)) + t*(3869.56755306385732+t2*(540.85752251676412+t2)))
+		return p / q
+	}
+	if mc > 0.014165 {
+		t := 98.0872976949485042*mc - 1.38940657184894556
+		t2 := t * t
+		p := ((8247.2601660137746 + t2*(6974.7495213178613+t2*108.098282908839979)) + t*(12967.7060124572914+t2*(1488.54008220335966+t2*1.32411616748380686)))
+		q := ((2340.47337508405427 + t2*(2287.70677154700516+t2*48.9575432570382154)) + t*(3915.63324533769906+t2*(550.45072377717361+t2)))
+		return p / q
+	}
+	if mc > 0.008213 {
+		t := 168.010752688172043*mc - 1.37987231182795699
+		t2 := t * t
+		p := ((8894.2961573611293 + t2*(7666.5611739483371+t2*121.863474964652041)) + t*(14113.7038749808951+t2*(1654.60731579994159+t2*1.53112170837206117)))
+		q := ((2344.88618943372377 + t2*(2313.28396270968662+t2*49.5906602613891184)) + t*(3942.81065054556536+t2*(558.07615380622169+t2)))
+		return p / q
+	}
+	if mc > 0 {
+		t := 1.0 - 121.758188238159016*mc
+		p := -math.Log(mc*0.0625) * (34813.4518336350547 + t*(235.767716637974271+t*0.199792723884069485)) / (69483.5736412906324 + t*(614.265044703187382+t))
+		q := -mc * (9382.53386835986099 + t*(51.6478985993381223+t*0.00410754154682816898)) / (37327.7262507318317 + t*(408.017247271148538+t))
+		return p + q
+	}
+
+	return math.Inf(1)
+}
+
+// CompleteE computes the complete elliptic integral of the 2nd kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1].
+//
+//	E(m) = \int_{0}^{π/2} {\sqrt{1-m{\sin^2θ}}} dθ
+func CompleteE(m float64) float64 {
+	// Reference:
+	// Toshio Fukushima, Precise and fast computation of complete elliptic integrals
+	// by piecewise minimax rational function approximation,
+	// Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76.
+	// https://doi.org/10.1016/j.cam.2014.12.038
+	// Original Fortran code available at:
+	// https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation
+	if m < 0 || 1 < m || math.IsNaN(m) {
+		return math.NaN()
+	}
+
+	mc := 1 - m
+
+	if mc > 0.566638 {
+		t := 2.30753965506897236*mc - 1.30753965506897236
+		t2 := t * t
+		p := ((19702.2363352671642 + t2*(18177.1879313824040+t2*409.975559128654710)) + t*(31904.1559574281609+t2*(4362.94760768571862+t2*10.3244775335024885)))
+		q := ((14241.2135819448616 + t2*(10266.4884503526076+t2*117.162100771599098)) + t*(20909.9899599927367+t2*(1934.86289070792954+t2)))
+		return p / q
+	}
+	if mc > 0.315153 {
+		t := 3.97638030101198879*mc - 1.25316818100483130
+		t2 := t * t
+		p := ((16317.0721393008221 + t2*(15129.4009798463159+t2*326.113727011739428)) + t*(26627.8852140835023+t2*(3574.15857605556033+t2*7.93163724081373477)))
+		q := ((13047.1505096551210 + t2*(9964.25173735060361+t2*117.670514069579649)) + t*(19753.5762165922376+t2*(1918.72232033637537+t2)))
+		return p / q
+	}
+	if mc > 0.171355 {
+		t := 6.95419964116329852*mc - 1.19163687951153702
+		t2 := t * t
+		p := ((13577.3850240991520 + t2*(12871.9137872656293+t2*263.964361648520708)) + t*(22545.4744699553993+t2*(3000.74575264868572+t2*6.08522443139677663)))
+		q := ((11717.3306408059832 + t2*(9619.40382323874064+t2*118.690522739531267)) + t*(18431.1264424290258+t2*(1904.06010727307491+t2)))
+		return p / q
+	}
+	if mc > 0.090670 {
+		t := 12.3938774245522712*mc - 1.12375286608415443
+		t2 := t * t
+		p := ((11307.9485341543712 + t2*(11208.6068472959372+t2*219.253495956962613)) + t*(19328.6173704569489+t2*(2596.54874477084334+t2*4.66931143174036616)))
+		q := ((10307.6837501971393 + t2*(9241.7604666150102+t2*120.498555754227847)) + t*(16982.2450249024383+t2*(1893.41905403040679+t2)))
+		return p / q
+	}
+	if mc > 0.046453 {
+		t := 22.6157360291290680*mc - 1.05056878576113260
+		t2 := t * t
+		p := ((9383.1490856819874 + t2*(9977.2498973537718+t2*188.618148076418837)) + t*(16718.9730458676860+t2*(2323.49987246555537+t2*3.59313532204509922)))
+		q := ((8877.1964704758383 + t2*(8840.2771293410661+t2*123.422125687316355)) + t*(15450.0537230364062+t2*(1889.13672102820913+t2)))
+		return p / q
+	}
+	if mc > 0.022912 {
+		t := 42.4790790535661187*mc - 0.973280659275306911
+		t2 := t * t
+		p := ((7719.1171817802054 + t2*(9045.3996063894006+t2*169.386557799782496)) + t*(14521.7363804934985+t2*(2149.92068078627829+t2*2.78515570453129137)))
+		q := ((7479.7539074698012 + t2*(8420.3848818926324+t2*127.802109608726363)) + t*(13874.4978011497847+t2*(1892.69753150329759+t2)))
+		return p / q
+	}
+	if mc > 0.010809 {
+		t := 82.6241427745187144*mc - 0.893084359249772784
+		t2 := t * t
+		p := ((6261.6095608987273 + t2*(8304.3265605809870+t2*159.371262600702237)) + t*(12593.0874916293982+t2*(2048.68391263416822+t2*2.18867046462858104)))
+		q := ((6156.4532048239501 + t2*(7979.7435857665227+t2*133.911640385965187)) + t*(12283.8373999680518+t2*(1903.60556312663537+t2)))
+		return p / q
+	}
+	if mc > 0.004841 {
+		t := 167.560321715817694*mc - 0.811159517426273458
+		t2 := t * t
+		p := ((4978.06146583586728 + t2*(7664.6703673290453+t2*156.689647694892782)) + t*(10831.7178150656694+t2*(1995.66437151562090+t2*1.75859085945198570)))
+		q := ((4935.56743322938333 + t2*(7506.8028283118051+t2*141.854303920116856)) + t*(10694.5510113880077+t2*(1918.38517009740321+t2)))
+		return p / q
+	}
+	if mc > 0 {
+		t := 1.0 - 206.568890725056806*mc
+		p := -mc * math.Log(mc*0.0625) * (41566.6612602868736 + t*(154.034981522913482+t*0.0618072471798575991)) / (165964.442527585615 + t*(917.589668642251803+t))
+		q := (132232.803956682877 + t*(353.375480007017643-t*1.40105837312528026)) / (132393.665743088043 + t*(192.112635228732532-t))
+		return p + q
+	}
+
+	return 1
+}
+
+// CompleteB computes an associate complete elliptic integral of the 2nd kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1].
+//
+//	B(m) = \int_{0}^{π/2} {\cos^2θ} / {\sqrt{1-m{\sin^2θ}}} dθ
+func CompleteB(m float64) float64 {
+	// Reference:
+	// Toshio Fukushima, Precise and fast computation of complete elliptic integrals
+	// by piecewise minimax rational function approximation,
+	// Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76.
+	// https://doi.org/10.1016/j.cam.2014.12.038
+	// Original Fortran code available at:
+	// https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation
+	if m < 0 || 1 < m || math.IsNaN(m) {
+		return math.NaN()
+	}
+
+	mc := 1 - m
+
+	if mc > 0.555073 {
+		t := 2.24755971204264969*mc - 1.24755971204264969
+		t2 := t * t
+		p := ((2030.25011505956379 + t2*(1727.60635612511943+t2*25.0715510300422010)) + t*(3223.16236100954529+t2*(361.164121995173076+t2*0.280355207707726826)))
+		q := ((2420.64907902774675 + t2*(2327.48464880306840+t2*47.9870997057202318)) + t*(4034.28168313496638+t2*(549.234220839203960+t2)))
+		return p / q
+	}
+	if mc > 0.302367 {
+		t := 3.95716761770595079*mc - 1.19651690106289522
+		t2 := t * t
+		p := ((2209.26925068374373 + t2*(1981.37862223307242+t2*29.7612810087709299)) + t*(3606.58475322372526+t2*(422.693774742063054+t2*0.334623999861181980)))
+		q := ((2499.57898767250755 + t2*(2467.63998386656941+t2*50.0198090806651216)) + t*(4236.30953048456334+t2*(581.879599221457589+t2)))
+		return p / q
+	}
+	if mc > 0.161052 {
+		t := 7.07638962601280827*mc - 1.13966670204861480
+		t2 := t * t
+		p := ((2359.14823394150129 + t2*(2254.30785457761760+t2*35.2259786264917876)) + t*(3983.28520266051676+t2*(492.601686517364701+t2*0.396605124984359783)))
+		q := ((2563.95563932625156 + t2*(2633.23323959119935+t2*52.6711647124832948)) + t*(4450.19076667898892+t2*(622.983787815718489+t2)))
+		return p / q
+	}
+	if mc > 0.083522 {
+		t := 12.8982329420869341*mc - 1.07728621178898491
+		t2 := t * t
+		p := ((2464.65334987833736 + t2*(2541.68516994216007+t2*41.5832527504007778)) + t*(4333.38639187691528+t2*(571.53606797524881+t2*0.465975784547025267)))
+		q := ((2600.66956117247726 + t2*(2823.69445052534842+t2*56.136001230010910)) + t*(4661.64381841490914+t2*(674.25435972414302+t2)))
+		return p / q
+	}
+	if mc > 0.041966 {
+		t := 24.0639137549331023*mc - 1.00986620463952257
+		t2 := t * t
+		p := ((2509.86724450741259 + t2*(2835.27071287535469+t2*48.9701196718008345)) + t*(4631.12336462339975+t2*(659.86172161727281+t2*0.54158304771955794)))
+		q := ((2594.15983397593723 + t2*(3034.20118545214106+t2*60.652838995496991)) + t*(4848.17491604384532+t2*(737.15143838356850+t2)))
+		return p / q
+	}
+	if mc > 0.020313 {
+		t := 46.1829769546944996*mc - 0.938114810880709371
+		t2 := t * t
+		p := ((2480.58307884128017 + t2*(3122.00900554841322+t2*57.541132641218839)) + t*(4845.57861173250699+t2*(757.31633816400643+t2*0.62119950515996627)))
+		q := ((2528.85218300581396 + t2*(3253.86151324157460+t2*66.496093157522450)) + t*(4979.31783250484768+t2*(812.40556572486862+t2)))
+		return p / q
+	}
+	if mc > 0.009408 {
+		t := 91.7010545621274645*mc - 0.862723521320495186
+		t2 := t * t
+		p := ((2365.25385348859592 + t2*(3381.09304915246175+t2*67.442026950538221)) + t*(4939.53925884558687+t2*(862.16657576129841+t2*0.70143698925710129)))
+		q := ((2390.48737882063755 + t2*(3462.34808443022907+t2*73.934680452209164)) + t*(5015.4675579215077+t2*(898.99542983710459+t2)))
+		return p / q
+	}
+	if mc > 0.004136 {
+		t := 189.681335356600910*mc - 0.784522003034901366
+		t2 := t * t
+		p := ((2160.82916040868119 + t2*(3584.53058926175721+t2*78.769178005879162)) + t*(4877.14832623847052+t2*(970.53716686804832+t2*0.77797110431753920)))
+		q := ((2172.70451405048305 + t2*(3630.52345460629336+t2*83.173163222639080)) + t*(4916.35263668839769+t2*(993.36676027886685+t2)))
+		return p / q
+	}
+	if mc > 0 {
+		t := 1 - 106.292517006802721*mc
+		p := mc * math.Log(mc*0.0625) * (6607.46457640413908 + t*(19.0287633783211078-t*0.00625368946932704460)) / (26150.3443630974309 + t*(354.603981274536040+t))
+		q := (26251.5678902584870 + t*(168.788023807915689+t*0.352150236262724288)) / (26065.7912239203873 + t*(353.916840382280456+t))
+		return p + q
+	}
+
+	return 1
+}
+
+// CompleteD computes an associate complete elliptic integral of the 2nd kind, 0≤m≤1. It returns math.NaN() if m is not in [0,1].
+//
+//	D(m) = \int_{0}^{π/2} {\sin^2θ} / {\sqrt{1-m{\sin^2θ}}} dθ
+func CompleteD(m float64) float64 {
+	// Reference:
+	// Toshio Fukushima, Precise and fast computation of complete elliptic integrals
+	// by piecewise minimax rational function approximation,
+	// Journal of Computational and Applied Mathematics, Volume 282, 2015, Pages 71-76.
+	// https://doi.org/10.1016/j.cam.2014.12.038
+	// Original Fortran code available at:
+	// https://www.researchgate.net/publication/295857819_xceitxt_F90_package_of_complete_elliptic_integral_computation
+	if m < 0 || 1 < m || math.IsNaN(m) {
+		return math.NaN()
+	}
+
+	mc := 1 - m
+
+	if mc > 0.599909 {
+		t := 2.49943137936119533*mc - 1.49943137936119533
+		t2 := t * t
+		p := ((1593.39813781813498 + t2*(1058.56241259843217+t2*11.7584241242587571)) + t*(2233.25576544961714+t2*(195.247394601357872+t2*0.101486443490307517)))
+		q := ((1685.47865546030468 + t2*(1604.88100543517015+t2*38.6743012128666717)) + t*(2756.20968383181114+t2*(397.504162950935944+t2)))
+		return p / q
+	}
+	if mc > 0.359180 {
+		t := 4.15404874360795750*mc - 1.49205122772910617
+		t2 := t * t
+		p := ((1967.01442513777287 + t2*(1329.30058268219177+t2*15.0447805948342760)) + t*(2779.87604145516343+t2*(247.475085945854673+t2*0.130547566005491628)))
+		q := ((1749.70634057327467 + t2*(1654.40804288486242+t2*39.1895256017535337)) + t*(2853.92630369567765+t2*(406.925098588378587+t2)))
+		return p / q
+	}
+	if mc > 0.214574 {
+		t := 6.91534237860116454*mc - 1.48385267554596628
+		t2 := t * t
+		p := ((2409.64196912091452 + t2*(1659.30176823041376+t2*19.1942111405094383)) + t*(3436.40744503228691+t2*(312.186468430688790+t2*0.167847673021897479)))
+		q := ((1824.89205701262525 + t2*(1715.38574780156913+t2*39.8798253173462218)) + t*(2971.02216287936566+t2*(418.929791715319490+t2)))
+		return p / q
+	}
+	if mc > 0.127875 {
+		t := 11.5341584101316047*mc - 1.47493050669557896
+		t2 := t * t
+		p := ((2926.81143179637839 + t2*(2056.45624281065334+t2*24.3811986813439843)) + t*(4214.52119721241319+t2*(391.420514384925370+t2*0.215574280659075512)))
+		q := ((1910.33091918583314 + t2*(1787.99942542734799+t2*40.7663012893484449)) + t*(3107.04531802441481+t2*(433.673494280825971+t2)))
+		return p / q
+	}
+	if mc > 0.076007 {
+		t := 19.2797100331611013*mc - 1.46539292049047582
+		t2 := t * t
+		p := ((3520.63614251102960 + t2*(2526.67111759550923+t2*30.7739877519417978)) + t*(5121.2842239226937+t2*(486.926821696342529+t2*0.276315678908126399)))
+		q := ((2003.81997889501324 + t2*(1871.05914195570669+t2*41.8489850490387023)) + t*(3259.09205279874214+t2*(451.007555352632053+t2)))
+		return p / q
+	}
+	if mc > 0.045052 {
+		t := 32.3049588111775157*mc - 1.45540300436116944
+		t2 := t * t
+		p := ((4188.00087087025347 + t2*(3072.05695847158556+t2*38.5070211470790031)) + t*(6156.0080960857764+t2*(599.76666155374012+t2*0.352955925261363680)))
+		q := ((2101.60113938424690 + t2*(1961.76794074710108+t2*43.0997999502743622)) + t*(3421.55151253792527+t2*(470.407158843118117+t2)))
+		return p / q
+	}
+	if mc > 0.026626 {
+		t := 54.2711386084880061*mc - 1.44502333658960165
+		t2 := t * t
+		p := ((4916.74442376570733 + t2*(3688.12811638360551+t2*47.6447145147811350)) + t*(7304.6632479558695+t2*(729.75841970840314+t2*0.448422756936257635)))
+		q := ((2197.49982676612397 + t2*(2055.19657857622715+t2*44.4576261146308645)) + t*(3584.94502590860852+t2*(490.880160668822953+t2)))
+		return p / q
+	}
+	if mc > 0.015689 {
+		t := 91.4327512114839536*mc - 1.43448843375697175
+		t2 := t * t
+		p := ((5688.7542903989517 + t2*(4364.21513060078954+t2*58.159468141567195)) + t*(8542.6096475195826+t2*(875.35992968472914+t2*0.56528145509695951)))
+		q := ((2285.44062680812883 + t2*(2145.80779422696555+t2*45.8427480379028781)) + t*(3739.30422133833258+t2*(511.23253971875808+t2)))
+		return p / q
+	}
+	if mc > 0.009216 {
+		t := 154.487872701992894*mc - 1.42376023482156651
+		t2 := t * t
+		p := ((6475.3392225234969 + t2*(5081.2997108708577+t2*69.910123337464043)) + t*(9829.1138694605662+t2*(1033.32687775311981+t2*0.70526087421186325)))
+		q := ((2357.74885505777295 + t2*(2226.89527217032394+t2*47.1609071069631012)) + t*(3872.32565152553360+t2*(530.03943432061149+t2)))
+		return p / q
+	}
+	if mc > 0 {
+		t := 1 - 108.506944444444444*mc
+		p := -math.Log(mc*0.0625) * (6.2904323649908115e6 + t*(58565.284164780476+t*(131.176674599188545+t*0.0426826410911220304))) / (1.24937550257219890e7 + t*(203580.534005225410+t*(921.17729845011868+t)))
+		q := -(27356.1090344387530 + t*(107.767403612304371-t*0.0827769227048233593)) / (27104.0854889805978 + t*(358.708172147752755+t))
+		return p + q
+	}
+
+	return math.Inf(1)
+}

vendor/gonum.org/v1/gonum/mathext/erf.go

@@ -0,0 +1,91 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import "math"
+
+/*
+Copyright (c) 2012 The Probab Authors. All rights reserved.
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+* Redistributions of source code must retain the above copyright
+notice, this list of conditions and the following disclaimer.
+* Redistributions in binary form must reproduce the above
+copyright notice, this list of conditions and the following disclaimer
+in the documentation and/or other materials provided with the
+distribution.
+* Neither the name of Google Inc. nor the names of its
+contributors may be used to endorse or promote products derived from
+this software without specific prior written permission.
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+*/
+
+// NormalQuantile computes the quantile function (inverse CDF) of the standard
+// normal. NormalQuantile panics if the input p is less than 0 or greater than 1.
+func NormalQuantile(p float64) float64 {
+	switch {
+	case p < 0 || 1 < p:
+		panic("mathext: quantile out of bounds")
+	case p == 1:
+		return math.Inf(1)
+	case p == 0:
+		return math.Inf(-1)
+	}
+	// Compute rational approximation based on the value of p.
+
+	dp := p - 0.5
+	if math.Abs(dp) <= 0.425 {
+		z := 0.180625 - dp*dp
+		z1 := ((((((zQSA[0]*z+zQSA[1])*z+zQSA[2])*z+zQSA[3])*z+zQSA[4])*z+zQSA[5])*z+zQSA[6])*z + zQSA[7]
+		z2 := ((((((zQSB[0]*z+zQSB[1])*z+zQSB[2])*z+zQSB[3])*z+zQSB[4])*z+zQSB[5])*z+zQSB[6])*z + zQSB[7]
+		return dp * z1 / z2
+	}
+
+	if p < 0.5 {
+		r := math.Sqrt(-math.Log(p))
+		if r <= 5.0 {
+			z := r - 1.6
+			z1 := ((((((zQIA[0]*z+zQIA[1])*z+zQIA[2])*z+zQIA[3])*z+zQIA[4])*z+zQIA[5])*z+zQIA[6])*z + zQIA[7]
+			z2 := ((((((zQIB[0]*z+zQIB[1])*z+zQIB[2])*z+zQIB[3])*z+zQIB[4])*z+zQIB[5])*z+zQIB[6])*z + zQIB[7]
+			return -z1 / z2
+		}
+		z := r - 5
+		z1 := ((((((zQTA[0]*z+zQTA[1])*z+zQTA[2])*z+zQTA[3])*z+zQTA[4])*z+zQTA[5])*z+zQTA[6])*z + zQTA[7]
+		z2 := ((((((zQTB[0]*z+zQTB[1])*z+zQTB[2])*z+zQTB[3])*z+zQTB[4])*z+zQTB[5])*z+zQTB[6])*z + zQTB[7]
+		return -z1 / z2
+	}
+	r := math.Sqrt(-math.Log(1 - p))
+	if r <= 5.0 {
+		z := r - 1.6
+		z1 := ((((((zQIA[0]*z+zQIA[1])*z+zQIA[2])*z+zQIA[3])*z+zQIA[4])*z+zQIA[5])*z+zQIA[6])*z + zQIA[7]
+		z2 := ((((((zQIB[0]*z+zQIB[1])*z+zQIB[2])*z+zQIB[3])*z+zQIB[4])*z+zQIB[5])*z+zQIB[6])*z + zQIB[7]
+		return z1 / z2
+	}
+
+	z := r - 5
+	z1 := ((((((zQTA[0]*z+zQTA[1])*z+zQTA[2])*z+zQTA[3])*z+zQTA[4])*z+zQTA[5])*z+zQTA[6])*z + zQTA[7]
+	z2 := ((((((zQTB[0]*z+zQTB[1])*z+zQTB[2])*z+zQTB[3])*z+zQTB[4])*z+zQTB[5])*z+zQTB[6])*z + zQTB[7]
+	return z1 / z2
+}
+
+var (
+	zQSA = [...]float64{2509.0809287301226727, 33430.575583588128105, 67265.770927008700853, 45921.953931549871457, 13731.693765509461125, 1971.5909503065514427, 133.14166789178437745, 3.387132872796366608}
+	zQSB = [...]float64{5226.495278852854561, 28729.085735721942674, 39307.89580009271061, 21213.794301586595867, 5394.1960214247511077, 687.1870074920579083, 42.313330701600911252, 1.0}
+	zQIA = [...]float64{7.7454501427834140764e-4, 0.0227238449892691845833, 0.24178072517745061177, 1.27045825245236838258, 3.64784832476320460504, 5.7694972214606914055, 4.6303378461565452959, 1.42343711074968357734}
+	zQIB = [...]float64{1.05075007164441684324e-9, 5.475938084995344946e-4, 0.0151986665636164571966, 0.14810397642748007459, 0.68976733498510000455, 1.6763848301838038494, 2.05319162663775882187, 1.0}
+	zQTA = [...]float64{2.01033439929228813265e-7, 2.71155556874348757815e-5, 0.0012426609473880784386, 0.026532189526576123093, 0.29656057182850489123, 1.7848265399172913358, 5.4637849111641143699, 6.6579046435011037772}
+	zQTB = [...]float64{2.04426310338993978564e-15, 1.4215117583164458887e-7, 1.8463183175100546818e-5, 7.868691311456132591e-4, 0.0148753612908506148525, 0.13692988092273580531, 0.59983220655588793769, 1.0}
+)

vendor/gonum.org/v1/gonum/mathext/gamma_inc.go

@@ -0,0 +1,58 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import (
+	"gonum.org/v1/gonum/mathext/internal/cephes"
+)
+
+// GammaIncReg computes the regularized incomplete Gamma integral.
+//
+//	GammaIncReg(a,x) = (1/ Γ(a)) \int_0^x e^{-t} t^{a-1} dt
+//
+// The input argument a must be positive and x must be non-negative or GammaIncReg
+// will panic.
+//
+// See http://mathworld.wolfram.com/IncompleteGammaFunction.html
+// or https://en.wikipedia.org/wiki/Incomplete_gamma_function for more detailed
+// information.
+func GammaIncReg(a, x float64) float64 {
+	return cephes.Igam(a, x)
+}
+
+// GammaIncRegComp computes the complemented regularized incomplete Gamma integral.
+//
+//	GammaIncRegComp(a,x) = 1 - GammaIncReg(a,x)
+//	                     = (1/ Γ(a)) \int_x^\infty e^{-t} t^{a-1} dt
+//
+// The input argument a must be positive and x must be non-negative or
+// GammaIncRegComp will panic.
+func GammaIncRegComp(a, x float64) float64 {
+	return cephes.IgamC(a, x)
+}
+
+// GammaIncRegInv computes the inverse of the regularized incomplete Gamma integral. That is,
+// it returns the x such that:
+//
+//	GammaIncReg(a, x) = y
+//
+// The input argument a must be positive and y must be between 0 and 1
+// inclusive or GammaIncRegInv will panic. GammaIncRegInv should return a positive
+// number, but can return NaN if there is a failure to converge.
+func GammaIncRegInv(a, y float64) float64 {
+	return gammaIncRegInv(a, y)
+}
+
+// GammaIncRegCompInv computes the inverse of the complemented regularized incomplete Gamma
+// integral. That is, it returns the x such that:
+//
+//	GammaIncRegComp(a, x) = y
+//
+// The input argument a must be positive and y must be between 0 and 1
+// inclusive or GammaIncRegCompInv will panic. GammaIncRegCompInv should return a
+// positive number, but can return 0 even with non-zero y due to underflow.
+func GammaIncRegCompInv(a, y float64) float64 {
+	return cephes.IgamI(a, y)
+}

vendor/gonum.org/v1/gonum/mathext/gamma_inc_inv.go

@@ -0,0 +1,58 @@
+// Derived from SciPy's special/c_misc/gammaincinv.c
+// https://github.com/scipy/scipy/blob/master/scipy/special/c_misc/gammaincinv.c
+
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/mathext/internal/cephes"
+)
+
+const (
+	allowedATol = 1e-306
+	allowedRTol = 1e-6
+)
+
+func gammaIncReg(x float64, params []float64) float64 {
+	return cephes.Igam(params[0], x) - params[1]
+}
+
+// gammaIncRegInv is the inverse of the regularized incomplete Gamma integral. That is, it
+// returns x such that:
+//
+//	Igam(a, x) = y
+//
+// The input argument a must be positive and y must be between 0 and 1
+// inclusive or gammaIncRegInv will panic. gammaIncRegInv should return a
+// positive number, but can return NaN if there is a failure to converge.
+func gammaIncRegInv(a, y float64) float64 {
+	// For y not small, we just use
+	//  IgamI(a, 1-y)
+	// (inverse of the complemented incomplete Gamma integral). For y small,
+	// however, 1-y is about 1, and we lose digits.
+	if a <= 0 || y <= 0 || y >= 0.25 {
+		return cephes.IgamI(a, 1-y)
+	}
+
+	lo := 0.0
+	flo := -y
+	hi := cephes.IgamI(a, 0.75)
+	fhi := 0.25 - y
+
+	params := []float64{a, y}
+
+	// Also, after we generate a small interval by bisection above, false
+	// position will do a large step from an interval of width ~1e-4 to ~1e-14
+	// in one step (a=10, x=0.05, but similar for other values).
+	result, bestX, _, errEst := falsePosition(lo, hi, flo, fhi, 2*machEp, 2*machEp, 1e-2*a, gammaIncReg, params)
+	if result == fSolveMaxIterations && errEst > allowedATol+allowedRTol*math.Abs(bestX) {
+		bestX = math.NaN()
+	}
+
+	return bestX
+}

vendor/gonum.org/v1/gonum/mathext/internal/amos/doc.go

@@ -0,0 +1,6 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package amos implements functions originally in the Netlib code by Donald Amos.
+package amos // import "gonum.org/v1/gonum/mathext/internal/amos"

vendor/gonum.org/v1/gonum/mathext/internal/amos/staticcheck.conf

@@ -0,0 +1 @@
+checks = []

vendor/gonum.org/v1/gonum/mathext/internal/cephes/cephes.go

@@ -0,0 +1,28 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package cephes
+
+import "math"
+
+/*
+Additional copyright information:
+
+Code in this package is adapted from the Cephes library (http://www.netlib.org/cephes/).
+There is no explicit licence on Netlib, but the author has agreed to a BSD release.
+See https://github.com/deepmind/torch-cephes/blob/master/LICENSE.txt and
+https://lists.debian.org/debian-legal/2004/12/msg00295.html
+*/
+
+const (
+	paramOutOfBounds            = "cephes: parameter out of bounds"
+	errParamFunctionSingularity = "cephes: function singularity"
+)
+
+const (
+	machEp  = 1.0 / (1 << 53)
+	maxLog  = 1024 * math.Ln2
+	minLog  = -1075 * math.Ln2
+	maxIter = 2000
+)

vendor/gonum.org/v1/gonum/mathext/internal/cephes/doc.go

@@ -0,0 +1,6 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package cephes implements functions originally in the Netlib code by Stephen Mosher.
+package cephes // import "gonum.org/v1/gonum/mathext/internal/cephes"

vendor/gonum.org/v1/gonum/mathext/internal/cephes/igam.go

@@ -0,0 +1,320 @@
+// Derived from SciPy's special/cephes/igam.c and special/cephes/igam.h
+// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/igam.c
+// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/igam.h
+// Made freely available by Stephen L. Moshier without support or guarantee.
+
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+// Copyright ©1985, ©1987 by Stephen L. Moshier
+// Portions Copyright ©2016 The Gonum Authors. All rights reserved.
+
+package cephes
+
+import "math"
+
+const (
+	igamDimK       = 25
+	igamDimN       = 25
+	igam           = 1
+	igamC          = 0
+	igamSmall      = 20
+	igamLarge      = 200
+	igamSmallRatio = 0.3
+	igamLargeRatio = 4.5
+)
+
+var igamCoefs = [igamDimK][igamDimN]float64{
+	{-3.3333333333333333e-1, 8.3333333333333333e-2, -1.4814814814814815e-2, 1.1574074074074074e-3, 3.527336860670194e-4, -1.7875514403292181e-4, 3.9192631785224378e-5, -2.1854485106799922e-6, -1.85406221071516e-6, 8.296711340953086e-7, -1.7665952736826079e-7, 6.7078535434014986e-9, 1.0261809784240308e-8, -4.3820360184533532e-9, 9.1476995822367902e-10, -2.551419399494625e-11, -5.8307721325504251e-11, 2.4361948020667416e-11, -5.0276692801141756e-12, 1.1004392031956135e-13, 3.3717632624009854e-13, -1.3923887224181621e-13, 2.8534893807047443e-14, -5.1391118342425726e-16, -1.9752288294349443e-15},
+	{-1.8518518518518519e-3, -3.4722222222222222e-3, 2.6455026455026455e-3, -9.9022633744855967e-4, 2.0576131687242798e-4, -4.0187757201646091e-7, -1.8098550334489978e-5, 7.6491609160811101e-6, -1.6120900894563446e-6, 4.6471278028074343e-9, 1.378633446915721e-7, -5.752545603517705e-8, 1.1951628599778147e-8, -1.7543241719747648e-11, -1.0091543710600413e-9, 4.1627929918425826e-10, -8.5639070264929806e-11, 6.0672151016047586e-14, 7.1624989648114854e-12, -2.9331866437714371e-12, 5.9966963656836887e-13, -2.1671786527323314e-16, -4.9783399723692616e-14, 2.0291628823713425e-14, -4.13125571381061e-15},
+	{4.1335978835978836e-3, -2.6813271604938272e-3, 7.7160493827160494e-4, 2.0093878600823045e-6, -1.0736653226365161e-4, 5.2923448829120125e-5, -1.2760635188618728e-5, 3.4235787340961381e-8, 1.3721957309062933e-6, -6.298992138380055e-7, 1.4280614206064242e-7, -2.0477098421990866e-10, -1.4092529910867521e-8, 6.228974084922022e-9, -1.3670488396617113e-9, 9.4283561590146782e-13, 1.2872252400089318e-10, -5.5645956134363321e-11, 1.1975935546366981e-11, -4.1689782251838635e-15, -1.0940640427884594e-12, 4.6622399463901357e-13, -9.905105763906906e-14, 1.8931876768373515e-17, 8.8592218725911273e-15},
+	{6.4943415637860082e-4, 2.2947209362139918e-4, -4.6918949439525571e-4, 2.6772063206283885e-4, -7.5618016718839764e-5, -2.3965051138672967e-7, 1.1082654115347302e-5, -5.6749528269915966e-6, 1.4230900732435884e-6, -2.7861080291528142e-11, -1.6958404091930277e-7, 8.0994649053880824e-8, -1.9111168485973654e-8, 2.3928620439808118e-12, 2.0620131815488798e-9, -9.4604966618551322e-10, 2.1541049775774908e-10, -1.388823336813903e-14, -2.1894761681963939e-11, 9.7909989511716851e-12, -2.1782191880180962e-12, 6.2088195734079014e-17, 2.126978363279737e-13, -9.3446887915174333e-14, 2.0453671226782849e-14},
+	{-8.618882909167117e-4, 7.8403922172006663e-4, -2.9907248030319018e-4, -1.4638452578843418e-6, 6.6414982154651222e-5, -3.9683650471794347e-5, 1.1375726970678419e-5, 2.5074972262375328e-10, -1.6954149536558306e-6, 8.9075075322053097e-7, -2.2929348340008049e-7, 2.956794137544049e-11, 2.8865829742708784e-8, -1.4189739437803219e-8, 3.4463580499464897e-9, -2.3024517174528067e-13, -3.9409233028046405e-10, 1.8602338968504502e-10, -4.356323005056618e-11, 1.2786001016296231e-15, 4.6792750266579195e-12, -2.1492464706134829e-12, 4.9088156148096522e-13, -6.3385914848915603e-18, -5.0453320690800944e-14},
+	{-3.3679855336635815e-4, -6.9728137583658578e-5, 2.7727532449593921e-4, -1.9932570516188848e-4, 6.7977804779372078e-5, 1.419062920643967e-7, -1.3594048189768693e-5, 8.0184702563342015e-6, -2.2914811765080952e-6, -3.252473551298454e-10, 3.4652846491085265e-7, -1.8447187191171343e-7, 4.8240967037894181e-8, -1.7989466721743515e-14, -6.3061945000135234e-9, 3.1624176287745679e-9, -7.8409242536974293e-10, 5.1926791652540407e-15, 9.3589442423067836e-11, -4.5134262161632782e-11, 1.0799129993116827e-11, -3.661886712685252e-17, -1.210902069055155e-12, 5.6807435849905643e-13, -1.3249659916340829e-13},
+	{5.3130793646399222e-4, -5.9216643735369388e-4, 2.7087820967180448e-4, 7.9023532326603279e-7, -8.1539693675619688e-5, 5.6116827531062497e-5, -1.8329116582843376e-5, -3.0796134506033048e-9, 3.4651553688036091e-6, -2.0291327396058604e-6, 5.7887928631490037e-7, 2.338630673826657e-13, -8.8286007463304835e-8, 4.7435958880408128e-8, -1.2545415020710382e-8, 8.6496488580102925e-14, 1.6846058979264063e-9, -8.5754928235775947e-10, 2.1598224929232125e-10, -7.6132305204761539e-16, -2.6639822008536144e-11, 1.3065700536611057e-11, -3.1799163902367977e-12, 4.7109761213674315e-18, 3.6902800842763467e-13},
+	{3.4436760689237767e-4, 5.1717909082605922e-5, -3.3493161081142236e-4, 2.812695154763237e-4, -1.0976582244684731e-4, -1.2741009095484485e-7, 2.7744451511563644e-5, -1.8263488805711333e-5, 5.7876949497350524e-6, 4.9387589339362704e-10, -1.0595367014026043e-6, 6.1667143761104075e-7, -1.7562973359060462e-7, -1.2974473287015439e-12, 2.695423606288966e-8, -1.4578352908731271e-8, 3.887645959386175e-9, -3.8810022510194121e-17, -5.3279941738772867e-10, 2.7437977643314845e-10, -6.9957960920705679e-11, 2.5899863874868481e-17, 8.8566890996696381e-12, -4.403168815871311e-12, 1.0865561947091654e-12},
+	{-6.5262391859530942e-4, 8.3949872067208728e-4, -4.3829709854172101e-4, -6.969091458420552e-7, 1.6644846642067548e-4, -1.2783517679769219e-4, 4.6299532636913043e-5, 4.5579098679227077e-9, -1.0595271125805195e-5, 6.7833429048651666e-6, -2.1075476666258804e-6, -1.7213731432817145e-11, 3.7735877416110979e-7, -2.1867506700122867e-7, 6.2202288040189269e-8, 6.5977038267330006e-16, -9.5903864974256858e-9, 5.2132144922808078e-9, -1.3991589583935709e-9, 5.382058999060575e-16, 1.9484714275467745e-10, -1.0127287556389682e-10, 2.6077347197254926e-11, -5.0904186999932993e-18, -3.3721464474854592e-12},
+	{-5.9676129019274625e-4, -7.2048954160200106e-5, 6.7823088376673284e-4, -6.4014752602627585e-4, 2.7750107634328704e-4, 1.8197008380465151e-7, -8.4795071170685032e-5, 6.105192082501531e-5, -2.1073920183404862e-5, -8.8585890141255994e-10, 4.5284535953805377e-6, -2.8427815022504408e-6, 8.7082341778646412e-7, 3.6886101871706965e-12, -1.5344695190702061e-7, 8.862466778790695e-8, -2.5184812301826817e-8, -1.0225912098215092e-14, 3.8969470758154777e-9, -2.1267304792235635e-9, 5.7370135528051385e-10, -1.887749850169741e-19, -8.0931538694657866e-11, 4.2382723283449199e-11, -1.1002224534207726e-11},
+	{1.3324454494800656e-3, -1.9144384985654775e-3, 1.1089369134596637e-3, 9.932404122642299e-7, -5.0874501293093199e-4, 4.2735056665392884e-4, -1.6858853767910799e-4, -8.1301893922784998e-9, 4.5284402370562147e-5, -3.127053674781734e-5, 1.044986828530338e-5, 4.8435226265680926e-11, -2.1482565873456258e-6, 1.329369701097492e-6, -4.0295693092101029e-7, -1.7567877666323291e-13, 7.0145043163668257e-8, -4.040787734999483e-8, 1.1474026743371963e-8, 3.9642746853563325e-18, -1.7804938269892714e-9, 9.7480262548731646e-10, -2.6405338676507616e-10, 5.794875163403742e-18, 3.7647749553543836e-11},
+	{1.579727660730835e-3, 1.6251626278391582e-4, -2.0633421035543276e-3, 2.1389686185689098e-3, -1.0108559391263003e-3, -3.9912705529919201e-7, 3.6235025084764691e-4, -2.8143901463712154e-4, 1.0449513336495887e-4, 2.1211418491830297e-9, -2.5779417251947842e-5, 1.7281818956040463e-5, -5.6413773872904282e-6, -1.1024320105776174e-11, 1.1223224418895175e-6, -6.8693396379526735e-7, 2.0653236975414887e-7, 4.6714772409838506e-14, -3.5609886164949055e-8, 2.0470855345905963e-8, -5.8091738633283358e-9, -1.332821287582869e-16, 9.0354604391335133e-10, -4.9598782517330834e-10, 1.3481607129399749e-10},
+	{-4.0725121195140166e-3, 6.4033628338080698e-3, -4.0410161081676618e-3, -2.183732802866233e-6, 2.1740441801254639e-3, -1.9700440518418892e-3, 8.3595469747962458e-4, 1.9445447567109655e-8, -2.5779387120421696e-4, 1.9009987368139304e-4, -6.7696499937438965e-5, -1.4440629666426572e-10, 1.5712512518742269e-5, -1.0304008744776893e-5, 3.304517767401387e-6, 7.9829760242325709e-13, -6.4097794149313004e-7, 3.8894624761300056e-7, -1.1618347644948869e-7, -2.816808630596451e-15, 1.9878012911297093e-8, -1.1407719956357511e-8, 3.2355857064185555e-9, 4.1759468293455945e-20, -5.0423112718105824e-10},
+	{-5.9475779383993003e-3, -5.4016476789260452e-4, 8.7910413550767898e-3, -9.8576315587856125e-3, 5.0134695031021538e-3, 1.2807521786221875e-6, -2.0626019342754683e-3, 1.7109128573523058e-3, -6.7695312714133799e-4, -6.9011545676562133e-9, 1.8855128143995902e-4, -1.3395215663491969e-4, 4.6263183033528039e-5, 4.0034230613321351e-11, -1.0255652921494033e-5, 6.612086372797651e-6, -2.0913022027253008e-6, -2.0951775649603837e-13, 3.9756029041993247e-7, -2.3956211978815887e-7, 7.1182883382145864e-8, 8.925574873053455e-16, -1.2101547235064676e-8, 6.9350618248334386e-9, -1.9661464453856102e-9},
+	{1.7402027787522711e-2, -2.9527880945699121e-2, 2.0045875571402799e-2, 7.0289515966903407e-6, -1.2375421071343148e-2, 1.1976293444235254e-2, -5.4156038466518525e-3, -6.3290893396418616e-8, 1.8855118129005065e-3, -1.473473274825001e-3, 5.5515810097708387e-4, 5.2406834412550662e-10, -1.4357913535784836e-4, 9.9181293224943297e-5, -3.3460834749478311e-5, -3.5755837291098993e-12, 7.1560851960630076e-6, -4.5516802628155526e-6, 1.4236576649271475e-6, 1.8803149082089664e-14, -2.6623403898929211e-7, 1.5950642189595716e-7, -4.7187514673841102e-8, -6.5107872958755177e-17, 7.9795091026746235e-9},
+	{3.0249124160905891e-2, 2.4817436002649977e-3, -4.9939134373457022e-2, 5.9915643009307869e-2, -3.2483207601623391e-2, -5.7212968652103441e-6, 1.5085251778569354e-2, -1.3261324005088445e-2, 5.5515262632426148e-3, 3.0263182257030016e-8, -1.7229548406756723e-3, 1.2893570099929637e-3, -4.6845138348319876e-4, -1.830259937893045e-10, 1.1449739014822654e-4, -7.7378565221244477e-5, 2.5625836246985201e-5, 1.0766165333192814e-12, -5.3246809282422621e-6, 3.349634863064464e-6, -1.0381253128684018e-6, -5.608909920621128e-15, 1.9150821930676591e-7, -1.1418365800203486e-7, 3.3654425209171788e-8},
+	{-9.9051020880159045e-2, 1.7954011706123486e-1, -1.2989606383463778e-1, -3.1478872752284357e-5, 9.0510635276848131e-2, -9.2828824411184397e-2, 4.4412112839877808e-2, 2.7779236316835888e-7, -1.7229543805449697e-2, 1.4182925050891573e-2, -5.6214161633747336e-3, -2.39598509186381e-9, 1.6029634366079908e-3, -1.1606784674435773e-3, 4.1001337768153873e-4, 1.8365800754090661e-11, -9.5844256563655903e-5, 6.3643062337764708e-5, -2.076250624489065e-5, -1.1806020912804483e-13, 4.2131808239120649e-6, -2.6262241337012467e-6, 8.0770620494930662e-7, 6.0125912123632725e-16, -1.4729737374018841e-7},
+	{-1.9994542198219728e-1, -1.5056113040026424e-2, 3.6470239469348489e-1, -4.6435192311733545e-1, 2.6640934719197893e-1, 3.4038266027147191e-5, -1.3784338709329624e-1, 1.276467178337056e-1, -5.6213828755200985e-2, -1.753150885483011e-7, 1.9235592956768113e-2, -1.5088821281095315e-2, 5.7401854451350123e-3, 1.0622382710310225e-9, -1.5335082692563998e-3, 1.0819320643228214e-3, -3.7372510193945659e-4, -6.6170909729031985e-12, 8.4263617380909628e-5, -5.5150706827483479e-5, 1.7769536448348069e-5, 3.8827923210205533e-14, -3.53513697488768e-6, 2.1865832130045269e-6, -6.6812849447625594e-7},
+	{7.2438608504029431e-1, -1.3918010932653375, 1.0654143352413968, 1.876173868950258e-4, -8.2705501176152696e-1, 8.9352433347828414e-1, -4.4971003995291339e-1, -1.6107401567546652e-6, 1.9235590165271091e-1, -1.6597702160042609e-1, 6.8882222681814333e-2, 1.3910091724608687e-8, -2.146911561508663e-2, 1.6228980898865892e-2, -5.9796016172584256e-3, -1.1287469112826745e-10, 1.5167451119784857e-3, -1.0478634293553899e-3, 3.5539072889126421e-4, 8.1704322111801517e-13, -7.7773013442452395e-5, 5.0291413897007722e-5, -1.6035083867000518e-5, 1.2469354315487605e-14, 3.1369106244517615e-6},
+	{1.6668949727276811, 1.165462765994632e-1, -3.3288393225018906, 4.4692325482864037, -2.6977693045875807, -2.600667859891061e-4, 1.5389017615694539, -1.4937962361134612, 6.8881964633233148e-1, 1.3077482004552385e-6, -2.5762963325596288e-1, 2.1097676102125449e-1, -8.3714408359219882e-2, -7.7920428881354753e-9, 2.4267923064833599e-2, -1.7813678334552311e-2, 6.3970330388900056e-3, 4.9430807090480523e-11, -1.5554602758465635e-3, 1.0561196919903214e-3, -3.5277184460472902e-4, 9.3002334645022459e-14, 7.5285855026557172e-5, -4.8186515569156351e-5, 1.5227271505597605e-5},
+	{-6.6188298861372935, 1.3397985455142589e+1, -1.0789350606845146e+1, -1.4352254537875018e-3, 9.2333694596189809, -1.0456552819547769e+1, 5.5105526029033471, 1.2024439690716742e-5, -2.5762961164755816, 2.3207442745387179, -1.0045728797216284, -1.0207833290021914e-7, 3.3975092171169466e-1, -2.6720517450757468e-1, 1.0235252851562706e-1, 8.4329730484871625e-10, -2.7998284958442595e-2, 2.0066274144976813e-2, -7.0554368915086242e-3, 1.9402238183698188e-12, 1.6562888105449611e-3, -1.1082898580743683e-3, 3.654545161310169e-4, -5.1290032026971794e-11, -7.6340103696869031e-5},
+	{-1.7112706061976095e+1, -1.1208044642899116, 3.7131966511885444e+1, -5.2298271025348962e+1, 3.3058589696624618e+1, 2.4791298976200222e-3, -2.061089403411526e+1, 2.088672775145582e+1, -1.0045703956517752e+1, -1.2238783449063012e-5, 4.0770134274221141, -3.473667358470195, 1.4329352617312006, 7.1359914411879712e-8, -4.4797257159115612e-1, 3.4112666080644461e-1, -1.2699786326594923e-1, -2.8953677269081528e-10, 3.3125776278259863e-2, -2.3274087021036101e-2, 8.0399993503648882e-3, -1.177805216235265e-9, -1.8321624891071668e-3, 1.2108282933588665e-3, -3.9479941246822517e-4},
+	{7.389033153567425e+1, -1.5680141270402273e+2, 1.322177542759164e+2, 1.3692876877324546e-2, -1.2366496885920151e+2, 1.4620689391062729e+2, -8.0365587724865346e+1, -1.1259851148881298e-4, 4.0770132196179938e+1, -3.8210340013273034e+1, 1.719522294277362e+1, 9.3519707955168356e-7, -6.2716159907747034, 5.1168999071852637, -2.0319658112299095, -4.9507215582761543e-9, 5.9626397294332597e-1, -4.4220765337238094e-1, 1.6079998700166273e-1, -2.4733786203223402e-8, -4.0307574759979762e-2, 2.7849050747097869e-2, -9.4751858992054221e-3, 6.419922235909132e-6, 2.1250180774699461e-3},
+	{2.1216837098382522e+2, 1.3107863022633868e+1, -4.9698285932871748e+2, 7.3121595266969204e+2, -4.8213821720890847e+2, -2.8817248692894889e-2, 3.2616720302947102e+2, -3.4389340280087117e+2, 1.7195193870816232e+2, 1.4038077378096158e-4, -7.52594195897599e+1, 6.651969984520934e+1, -2.8447519748152462e+1, -7.613702615875391e-7, 9.5402237105304373, -7.5175301113311376, 2.8943997568871961, -4.6612194999538201e-7, -8.0615149598794088e-1, 5.8483006570631029e-1, -2.0845408972964956e-1, 1.4765818959305817e-4, 5.1000433863753019e-2, -3.3066252141883665e-2, 1.5109265210467774e-2},
+	{-9.8959643098322368e+2, 2.1925555360905233e+3, -1.9283586782723356e+3, -1.5925738122215253e-1, 1.9569985945919857e+3, -2.4072514765081556e+3, 1.3756149959336496e+3, 1.2920735237496668e-3, -7.525941715948055e+2, 7.3171668742208716e+2, -3.4137023466220065e+2, -9.9857390260608043e-6, 1.3356313181291573e+2, -1.1276295161252794e+2, 4.6310396098204458e+1, -7.9237387133614756e-6, -1.4510726927018646e+1, 1.1111771248100563e+1, -4.1690817945270892, 3.1008219800117808e-3, 1.1220095449981468, -7.6052379926149916e-1, 3.6262236505085254e-1, 2.216867741940747e-1, 4.8683443692930507e-1},
+}
+
+// Igam computes the incomplete Gamma integral.
+//
+//	Igam(a,x) = (1/ Γ(a)) \int_0^x e^{-t} t^{a-1} dt
+//
+// The input argument a must be positive and x must be non-negative or Igam
+// will panic.
+func Igam(a, x float64) float64 {
+	// The integral is evaluated by either a power series or continued fraction
+	// expansion, depending on the relative values of a and x.
+	// Sources:
+	// [1] "The Digital Library of Mathematical Functions", dlmf.nist.gov
+	// [2] Maddock et. al., "Incomplete Gamma Functions",
+	// http://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html
+
+	// Check zero integration limit first
+	if x == 0 {
+		return 0
+	}
+
+	if x < 0 || a <= 0 {
+		panic(paramOutOfBounds)
+	}
+
+	// Asymptotic regime where a ~ x; see [2].
+	absxmaA := math.Abs(x-a) / a
+	if (igamSmall < a && a < igamLarge && absxmaA < igamSmallRatio) ||
+		(igamLarge < a && absxmaA < igamLargeRatio/math.Sqrt(a)) {
+		return asymptoticSeries(a, x, igam)
+	}
+
+	if x > 1 && x > a {
+		return 1 - IgamC(a, x)
+	}
+
+	return igamSeries(a, x)
+}
+
+// IgamC computes the complemented incomplete Gamma integral.
+//
+//	IgamC(a,x) = 1 - Igam(a,x)
+//	           = (1/ Γ(a)) \int_0^\infty e^{-t} t^{a-1} dt
+//
+// The input argument a must be positive and x must be non-negative or
+// IgamC will panic.
+func IgamC(a, x float64) float64 {
+	// The integral is evaluated by either a power series or continued fraction
+	// expansion, depending on the relative values of a and x.
+	// Sources:
+	// [1] "The Digital Library of Mathematical Functions", dlmf.nist.gov
+	// [2] Maddock et. al., "Incomplete Gamma Functions",
+	// http://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/sf_gamma/igamma.html
+
+	switch {
+	case x < 0, a <= 0:
+		panic(paramOutOfBounds)
+	case x == 0:
+		return 1
+	case math.IsInf(x, 0):
+		return 0
+	}
+
+	// Asymptotic regime where a ~ x; see [2].
+	absxmaA := math.Abs(x-a) / a
+	if (igamSmall < a && a < igamLarge && absxmaA < igamSmallRatio) ||
+		(igamLarge < a && absxmaA < igamLargeRatio/math.Sqrt(a)) {
+		return asymptoticSeries(a, x, igamC)
+	}
+
+	// Everywhere else; see [2].
+	if x > 1.1 {
+		if x < a {
+			return 1 - igamSeries(a, x)
+		}
+		return igamCContinuedFraction(a, x)
+	} else if x <= 0.5 {
+		if -0.4/math.Log(x) < a {
+			return 1 - igamSeries(a, x)
+		}
+		return igamCSeries(a, x)
+	}
+
+	if x*1.1 < a {
+		return 1 - igamSeries(a, x)
+	}
+	return igamCSeries(a, x)
+}
+
+// igamFac computes
+//
+//	x^a * e^{-x} / Γ(a)
+//
+// corrected from (15) and (16) in [2] by replacing
+//
+//	e^{x - a}
+//
+// with
+//
+//	e^{a - x}
+func igamFac(a, x float64) float64 {
+	if math.Abs(a-x) > 0.4*math.Abs(a) {
+		ax := a*math.Log(x) - x - lgam(a)
+		return math.Exp(ax)
+	}
+
+	fac := a + lanczosG - 0.5
+	res := math.Sqrt(fac/math.Exp(1)) / lanczosSumExpgScaled(a)
+
+	if a < 200 && x < 200 {
+		res *= math.Exp(a-x) * math.Pow(x/fac, a)
+	} else {
+		num := x - a - lanczosG + 0.5
+		res *= math.Exp(a*log1pmx(num/fac) + x*(0.5-lanczosG)/fac)
+	}
+
+	return res
+}
+
+// igamCContinuedFraction computes IgamC using DLMF 8.9.2.
+func igamCContinuedFraction(a, x float64) float64 {
+	ax := igamFac(a, x)
+	if ax == 0 {
+		return 0
+	}
+
+	// Continued fraction
+	y := 1 - a
+	z := x + y + 1
+	c := 0.0
+	pkm2 := 1.0
+	qkm2 := x
+	pkm1 := x + 1.0
+	qkm1 := z * x
+	ans := pkm1 / qkm1
+
+	for i := 0; i < maxIter; i++ {
+		c += 1.0
+		y += 1.0
+		z += 2.0
+		yc := y * c
+		pk := pkm1*z - pkm2*yc
+		qk := qkm1*z - qkm2*yc
+		var t float64
+		if qk != 0 {
+			r := pk / qk
+			t = math.Abs((ans - r) / r)
+			ans = r
+		} else {
+			t = 1.0
+		}
+		pkm2 = pkm1
+		pkm1 = pk
+		qkm2 = qkm1
+		qkm1 = qk
+		if math.Abs(pk) > big {
+			pkm2 *= biginv
+			pkm1 *= biginv
+			qkm2 *= biginv
+			qkm1 *= biginv
+		}
+		if t <= machEp {
+			break
+		}
+	}
+
+	return ans * ax
+}
+
+// igamSeries computes Igam using DLMF 8.11.4.
+func igamSeries(a, x float64) float64 {
+	ax := igamFac(a, x)
+	if ax == 0 {
+		return 0
+	}
+
+	// Power series
+	r := a
+	c := 1.0
+	ans := 1.0
+
+	for i := 0; i < maxIter; i++ {
+		r += 1.0
+		c *= x / r
+		ans += c
+		if c <= machEp*ans {
+			break
+		}
+	}
+
+	return ans * ax / a
+}
+
+// igamCSeries computes IgamC using DLMF 8.7.3. This is related to the series
+// in igamSeries but extra care is taken to avoid cancellation.
+func igamCSeries(a, x float64) float64 {
+	fac := 1.0
+	sum := 0.0
+
+	for n := 1; n < maxIter; n++ {
+		fac *= -x / float64(n)
+		term := fac / (a + float64(n))
+		sum += term
+		if math.Abs(term) <= machEp*math.Abs(sum) {
+			break
+		}
+	}
+
+	logx := math.Log(x)
+	term := -expm1(a*logx - lgam1p(a))
+	return term - math.Exp(a*logx-lgam(a))*sum
+}
+
+// asymptoticSeries computes Igam/IgamC using DLMF 8.12.3/8.12.4.
+func asymptoticSeries(a, x float64, fun int) float64 {
+	maxpow := 0
+	lambda := x / a
+	sigma := (x - a) / a
+	absoldterm := math.MaxFloat64
+	etapow := [igamDimN]float64{1}
+	sum := 0.0
+	afac := 1.0
+
+	var sgn float64
+	if fun == igam {
+		sgn = -1
+	} else {
+		sgn = 1
+	}
+
+	var eta float64
+	if lambda > 1 {
+		eta = math.Sqrt(-2 * log1pmx(sigma))
+	} else if lambda < 1 {
+		eta = -math.Sqrt(-2 * log1pmx(sigma))
+	} else {
+		eta = 0
+	}
+	res := 0.5 * math.Erfc(sgn*eta*math.Sqrt(a/2))
+
+	for k := 0; k < igamDimK; k++ {
+		ck := igamCoefs[k][0]
+		for n := 1; n < igamDimN; n++ {
+			if n > maxpow {
+				etapow[n] = eta * etapow[n-1]
+				maxpow++
+			}
+			ckterm := igamCoefs[k][n] * etapow[n]
+			ck += ckterm
+			if math.Abs(ckterm) < machEp*math.Abs(ck) {
+				break
+			}
+		}
+		term := ck * afac
+		absterm := math.Abs(term)
+		if absterm > absoldterm {
+			break
+		}
+		sum += term
+		if absterm < machEp*math.Abs(sum) {
+			break
+		}
+		absoldterm = absterm
+		afac /= a
+	}
+	res += sgn * math.Exp(-0.5*a*eta*eta) * sum / math.Sqrt(2*math.Pi*a)
+
+	return res
+}

vendor/gonum.org/v1/gonum/mathext/internal/cephes/igami.go

@@ -0,0 +1,155 @@
+// Derived from SciPy's special/cephes/igami.c
+// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/igami.c
+// Made freely available by Stephen L. Moshier without support or guarantee.
+
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+// Copyright ©1984, ©1987, ©1995 by Stephen L. Moshier
+// Portions Copyright ©2017 The Gonum Authors. All rights reserved.
+
+package cephes
+
+import "math"
+
+// IgamI computes the inverse of the incomplete Gamma function. That is, it
+// returns the x such that:
+//
+//	IgamC(a, x) = p
+//
+// The input argument a must be positive and p must be between 0 and 1
+// inclusive or IgamI will panic. IgamI should return a positive number, but
+// can return 0 even with non-zero y due to underflow.
+func IgamI(a, p float64) float64 {
+	// Bound the solution
+	x0 := math.MaxFloat64
+	yl := 0.0
+	x1 := 0.0
+	yh := 1.0
+	dithresh := 5.0 * machEp
+
+	if p < 0 || p > 1 || a <= 0 {
+		panic(paramOutOfBounds)
+	}
+
+	if p == 0 {
+		return math.Inf(1)
+	}
+
+	if p == 1 {
+		return 0.0
+	}
+
+	// Starting with the approximate value
+	//  x = a y^3
+	// where
+	//  y = 1 - d - ndtri(p) sqrt(d)
+	// and
+	//  d = 1/9a
+	// the routine performs up to 10 Newton iterations to find the root of
+	//  IgamC(a, x) - p = 0
+	d := 1.0 / (9.0 * a)
+	y := 1.0 - d - Ndtri(p)*math.Sqrt(d)
+	x := a * y * y * y
+
+	lgm := lgam(a)
+
+	for i := 0; i < 10; i++ {
+		if x > x0 || x < x1 {
+			break
+		}
+
+		y = IgamC(a, x)
+
+		if y < yl || y > yh {
+			break
+		}
+
+		if y < p {
+			x0 = x
+			yl = y
+		} else {
+			x1 = x
+			yh = y
+		}
+
+		// Compute the derivative of the function at this point
+		d = (a-1)*math.Log(x) - x - lgm
+		if d < -maxLog {
+			break
+		}
+		d = -math.Exp(d)
+
+		// Compute the step to the next approximation of x
+		d = (y - p) / d
+		if math.Abs(d/x) < machEp {
+			return x
+		}
+		x = x - d
+	}
+
+	d = 0.0625
+	if x0 == math.MaxFloat64 {
+		if x <= 0 {
+			x = 1
+		}
+		for x0 == math.MaxFloat64 {
+			x = (1 + d) * x
+			y = IgamC(a, x)
+			if y < p {
+				x0 = x
+				yl = y
+				break
+			}
+			d = d + d
+		}
+	}
+
+	d = 0.5
+	dir := 0
+	for i := 0; i < 400; i++ {
+		x = x1 + d*(x0-x1)
+		y = IgamC(a, x)
+
+		lgm = (x0 - x1) / (x1 + x0)
+		if math.Abs(lgm) < dithresh {
+			break
+		}
+
+		lgm = (y - p) / p
+		if math.Abs(lgm) < dithresh {
+			break
+		}
+
+		if x <= 0 {
+			break
+		}
+
+		if y >= p {
+			x1 = x
+			yh = y
+			if dir < 0 {
+				dir = 0
+				d = 0.5
+			} else if dir > 1 {
+				d = 0.5*d + 0.5
+			} else {
+				d = (p - yl) / (yh - yl)
+			}
+			dir++
+		} else {
+			x0 = x
+			yl = y
+			if dir > 0 {
+				dir = 0
+				d = 0.5
+			} else if dir < -1 {
+				d = 0.5 * d
+			} else {
+				d = (p - yl) / (yh - yl)
+			}
+			dir--
+		}
+	}
+
+	return x
+}

vendor/gonum.org/v1/gonum/mathext/internal/cephes/incbeta.go

@@ -0,0 +1,312 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+/*
+ * Cephes Math Library, Release 2.3:  March, 1995
+ * Copyright 1984, 1995 by Stephen L. Moshier
+ */
+
+package cephes
+
+import (
+	"math"
+
+	"gonum.org/v1/gonum/mathext/internal/gonum"
+)
+
+const (
+	maxGam = 171.624376956302725
+	big    = 4.503599627370496e15
+	biginv = 2.22044604925031308085e-16
+)
+
+// Incbet computes the regularized incomplete beta function.
+func Incbet(aa, bb, xx float64) float64 {
+	if aa <= 0 || bb <= 0 {
+		panic(paramOutOfBounds)
+	}
+	if xx <= 0 || xx >= 1 {
+		if xx == 0 {
+			return 0
+		}
+		if xx == 1 {
+			return 1
+		}
+		panic(paramOutOfBounds)
+	}
+
+	var flag int
+	if bb*xx <= 1 && xx <= 0.95 {
+		t := pseries(aa, bb, xx)
+		return transformT(t, flag)
+	}
+
+	w := 1 - xx
+
+	// Reverse a and b if x is greater than the mean.
+	var a, b, xc, x float64
+	if xx > aa/(aa+bb) {
+		flag = 1
+		a = bb
+		b = aa
+		xc = xx
+		x = w
+	} else {
+		a = aa
+		b = bb
+		xc = w
+		x = xx
+	}
+
+	if flag == 1 && (b*x) <= 1.0 && x <= 0.95 {
+		t := pseries(a, b, x)
+		return transformT(t, flag)
+	}
+
+	// Choose expansion for better convergence.
+	y := x*(a+b-2.0) - (a - 1.0)
+	if y < 0.0 {
+		w = incbcf(a, b, x)
+	} else {
+		w = incbd(a, b, x) / xc
+	}
+
+	// Multiply w by the factor
+	// x^a * (1-x)^b * Γ(a+b) / (a*Γ(a)*Γ(b))
+	var t float64
+	y = a * math.Log(x)
+	t = b * math.Log(xc)
+	if (a+b) < maxGam && math.Abs(y) < maxLog && math.Abs(t) < maxLog {
+		t = math.Pow(xc, b)
+		t *= math.Pow(x, a)
+		t /= a
+		t *= w
+		t *= 1.0 / gonum.Beta(a, b)
+		return transformT(t, flag)
+	}
+
+	// Resort to logarithms.
+	y += t - gonum.Lbeta(a, b)
+	y += math.Log(w / a)
+	if y < minLog {
+		t = 0.0
+	} else {
+		t = math.Exp(y)
+	}
+
+	return transformT(t, flag)
+}
+
+func transformT(t float64, flag int) float64 {
+	if flag == 1 {
+		if t <= machEp {
+			t = 1.0 - machEp
+		} else {
+			t = 1.0 - t
+		}
+	}
+	return t
+}
+
+// incbcf returns the incomplete beta integral evaluated by a continued fraction
+// expansion.
+func incbcf(a, b, x float64) float64 {
+	var xk, pk, pkm1, pkm2, qk, qkm1, qkm2 float64
+	var k1, k2, k3, k4, k5, k6, k7, k8 float64
+	var r, t, ans, thresh float64
+	var n int
+
+	k1 = a
+	k2 = a + b
+	k3 = a
+	k4 = a + 1.0
+	k5 = 1.0
+	k6 = b - 1.0
+	k7 = k4
+	k8 = a + 2.0
+
+	pkm2 = 0.0
+	qkm2 = 1.0
+	pkm1 = 1.0
+	qkm1 = 1.0
+	ans = 1.0
+	r = 1.0
+	thresh = 3.0 * machEp
+
+	for n = 0; n <= 300; n++ {
+
+		xk = -(x * k1 * k2) / (k3 * k4)
+		pk = pkm1 + pkm2*xk
+		qk = qkm1 + qkm2*xk
+		pkm2 = pkm1
+		pkm1 = pk
+		qkm2 = qkm1
+		qkm1 = qk
+
+		xk = (x * k5 * k6) / (k7 * k8)
+		pk = pkm1 + pkm2*xk
+		qk = qkm1 + qkm2*xk
+		pkm2 = pkm1
+		pkm1 = pk
+		qkm2 = qkm1
+		qkm1 = qk
+
+		if qk != 0 {
+			r = pk / qk
+		}
+		if r != 0 {
+			t = math.Abs((ans - r) / r)
+			ans = r
+		} else {
+			t = 1.0
+		}
+
+		if t < thresh {
+			return ans
+		}
+
+		k1 += 1.0
+		k2 += 1.0
+		k3 += 2.0
+		k4 += 2.0
+		k5 += 1.0
+		k6 -= 1.0
+		k7 += 2.0
+		k8 += 2.0
+
+		if (math.Abs(qk) + math.Abs(pk)) > big {
+			pkm2 *= biginv
+			pkm1 *= biginv
+			qkm2 *= biginv
+			qkm1 *= biginv
+		}
+		if (math.Abs(qk) < biginv) || (math.Abs(pk) < biginv) {
+			pkm2 *= big
+			pkm1 *= big
+			qkm2 *= big
+			qkm1 *= big
+		}
+	}
+
+	return ans
+}
+
+// incbd returns the incomplete beta integral evaluated by a continued fraction
+// expansion.
+func incbd(a, b, x float64) float64 {
+	var xk, pk, pkm1, pkm2, qk, qkm1, qkm2 float64
+	var k1, k2, k3, k4, k5, k6, k7, k8 float64
+	var r, t, ans, z, thresh float64
+	var n int
+
+	k1 = a
+	k2 = b - 1.0
+	k3 = a
+	k4 = a + 1.0
+	k5 = 1.0
+	k6 = a + b
+	k7 = a + 1.0
+	k8 = a + 2.0
+
+	pkm2 = 0.0
+	qkm2 = 1.0
+	pkm1 = 1.0
+	qkm1 = 1.0
+	z = x / (1.0 - x)
+	ans = 1.0
+	r = 1.0
+	thresh = 3.0 * machEp
+	for n = 0; n <= 300; n++ {
+
+		xk = -(z * k1 * k2) / (k3 * k4)
+		pk = pkm1 + pkm2*xk
+		qk = qkm1 + qkm2*xk
+		pkm2 = pkm1
+		pkm1 = pk
+		qkm2 = qkm1
+		qkm1 = qk
+
+		xk = (z * k5 * k6) / (k7 * k8)
+		pk = pkm1 + pkm2*xk
+		qk = qkm1 + qkm2*xk
+		pkm2 = pkm1
+		pkm1 = pk
+		qkm2 = qkm1
+		qkm1 = qk
+
+		if qk != 0 {
+			r = pk / qk
+		}
+		if r != 0 {
+			t = math.Abs((ans - r) / r)
+			ans = r
+		} else {
+			t = 1.0
+		}
+
+		if t < thresh {
+			return ans
+		}
+
+		k1 += 1.0
+		k2 -= 1.0
+		k3 += 2.0
+		k4 += 2.0
+		k5 += 1.0
+		k6 += 1.0
+		k7 += 2.0
+		k8 += 2.0
+
+		if (math.Abs(qk) + math.Abs(pk)) > big {
+			pkm2 *= biginv
+			pkm1 *= biginv
+			qkm2 *= biginv
+			qkm1 *= biginv
+		}
+		if (math.Abs(qk) < biginv) || (math.Abs(pk) < biginv) {
+			pkm2 *= big
+			pkm1 *= big
+			qkm2 *= big
+			qkm1 *= big
+		}
+	}
+	return ans
+}
+
+// pseries returns the incomplete beta integral evaluated by a power series. Use
+// when b*x is small and x not too close to 1.
+func pseries(a, b, x float64) float64 {
+	var s, t, u, v, n, t1, z, ai float64
+	ai = 1.0 / a
+	u = (1.0 - b) * x
+	v = u / (a + 1.0)
+	t1 = v
+	t = u
+	n = 2.0
+	s = 0.0
+	z = machEp * ai
+	for math.Abs(v) > z {
+		u = (n - b) * x / n
+		t *= u
+		v = t / (a + n)
+		s += v
+		n += 1.0
+	}
+	s += t1
+	s += ai
+
+	u = a * math.Log(x)
+	if (a+b) < maxGam && math.Abs(u) < maxLog {
+		t = 1.0 / gonum.Beta(a, b)
+		s = s * t * math.Pow(x, a)
+	} else {
+		t = -gonum.Lbeta(a, b) + u + math.Log(s)
+		if t < minLog {
+			s = 0.0
+		} else {
+			s = math.Exp(t)
+		}
+	}
+	return (s)
+}

vendor/gonum.org/v1/gonum/mathext/internal/cephes/incbi.go

@@ -0,0 +1,247 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+/*
+ * Cephes Math Library Release 2.4:  March,1996
+ * Copyright 1984, 1996 by Stephen L. Moshier
+ */
+
+package cephes
+
+import "math"
+
+// Incbi computes the inverse of the regularized incomplete beta integral.
+func Incbi(aa, bb, yy0 float64) float64 {
+	var a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt float64
+	var i, rflg, dir, nflg int
+
+	if yy0 <= 0 {
+		return (0.0)
+	}
+	if yy0 >= 1.0 {
+		return (1.0)
+	}
+	x0 = 0.0
+	yl = 0.0
+	x1 = 1.0
+	yh = 1.0
+	nflg = 0
+
+	if aa <= 1.0 || bb <= 1.0 {
+		dithresh = 1.0e-6
+		rflg = 0
+		a = aa
+		b = bb
+		y0 = yy0
+		x = a / (a + b)
+		y = Incbet(a, b, x)
+		goto ihalve
+	} else {
+		dithresh = 1.0e-4
+	}
+	// Approximation to inverse function
+	yp = -Ndtri(yy0)
+
+	if yy0 > 0.5 {
+		rflg = 1
+		a = bb
+		b = aa
+		y0 = 1.0 - yy0
+		yp = -yp
+	} else {
+		rflg = 0
+		a = aa
+		b = bb
+		y0 = yy0
+	}
+
+	lgm = (yp*yp - 3.0) / 6.0
+	x = 2.0 / (1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0))
+	d = yp*math.Sqrt(x+lgm)/x - (1.0/(2.0*b-1.0)-1.0/(2.0*a-1.0))*(lgm+5.0/6.0-2.0/(3.0*x))
+	d = 2.0 * d
+	if d < minLog {
+		// mtherr("incbi", UNDERFLOW)
+		x = 0
+		goto done
+	}
+	x = a / (a + b*math.Exp(d))
+	y = Incbet(a, b, x)
+	yp = (y - y0) / y0
+	if math.Abs(yp) < 0.2 {
+		goto newt
+	}
+
+	/* Resort to interval halving if not close enough. */
+ihalve:
+
+	dir = 0
+	di = 0.5
+	for i = 0; i < 100; i++ {
+		if i != 0 {
+			x = x0 + di*(x1-x0)
+			if x == 1.0 {
+				x = 1.0 - machEp
+			}
+			if x == 0.0 {
+				di = 0.5
+				x = x0 + di*(x1-x0)
+				if x == 0.0 {
+					// mtherr("incbi", UNDERFLOW)
+					goto done
+				}
+			}
+			y = Incbet(a, b, x)
+			yp = (x1 - x0) / (x1 + x0)
+			if math.Abs(yp) < dithresh {
+				goto newt
+			}
+			yp = (y - y0) / y0
+			if math.Abs(yp) < dithresh {
+				goto newt
+			}
+		}
+		if y < y0 {
+			x0 = x
+			yl = y
+			if dir < 0 {
+				dir = 0
+				di = 0.5
+			} else if dir > 3 {
+				di = 1.0 - (1.0-di)*(1.0-di)
+			} else if dir > 1 {
+				di = 0.5*di + 0.5
+			} else {
+				di = (y0 - y) / (yh - yl)
+			}
+			dir += 1
+			if x0 > 0.75 {
+				if rflg == 1 {
+					rflg = 0
+					a = aa
+					b = bb
+					y0 = yy0
+				} else {
+					rflg = 1
+					a = bb
+					b = aa
+					y0 = 1.0 - yy0
+				}
+				x = 1.0 - x
+				y = Incbet(a, b, x)
+				x0 = 0.0
+				yl = 0.0
+				x1 = 1.0
+				yh = 1.0
+				goto ihalve
+			}
+		} else {
+			x1 = x
+			if rflg == 1 && x1 < machEp {
+				x = 0.0
+				goto done
+			}
+			yh = y
+			if dir > 0 {
+				dir = 0
+				di = 0.5
+			} else if dir < -3 {
+				di = di * di
+			} else if dir < -1 {
+				di = 0.5 * di
+			} else {
+				di = (y - y0) / (yh - yl)
+			}
+			dir -= 1
+		}
+	}
+	// mtherr("incbi", PLOSS)
+	if x0 >= 1.0 {
+		x = 1.0 - machEp
+		goto done
+	}
+	if x <= 0.0 {
+		// mtherr("incbi", UNDERFLOW)
+		x = 0.0
+		goto done
+	}
+
+newt:
+	if nflg > 0 {
+		goto done
+	}
+	nflg = 1
+	lgm = lgam(a+b) - lgam(a) - lgam(b)
+
+	for i = 0; i < 8; i++ {
+		/* Compute the function at this point. */
+		if i != 0 {
+			y = Incbet(a, b, x)
+		}
+		if y < yl {
+			x = x0
+			y = yl
+		} else if y > yh {
+			x = x1
+			y = yh
+		} else if y < y0 {
+			x0 = x
+			yl = y
+		} else {
+			x1 = x
+			yh = y
+		}
+		if x == 1.0 || x == 0.0 {
+			break
+		}
+		/* Compute the derivative of the function at this point. */
+		d = (a-1.0)*math.Log(x) + (b-1.0)*math.Log(1.0-x) + lgm
+		if d < minLog {
+			goto done
+		}
+		if d > maxLog {
+			break
+		}
+		d = math.Exp(d)
+		/* Compute the step to the next approximation of x. */
+		d = (y - y0) / d
+		xt = x - d
+		if xt <= x0 {
+			y = (x - x0) / (x1 - x0)
+			xt = x0 + 0.5*y*(x-x0)
+			if xt <= 0.0 {
+				break
+			}
+		}
+		if xt >= x1 {
+			y = (x1 - x) / (x1 - x0)
+			xt = x1 - 0.5*y*(x1-x)
+			if xt >= 1.0 {
+				break
+			}
+		}
+		x = xt
+		if math.Abs(d/x) < 128.0*machEp {
+			goto done
+		}
+	}
+	/* Did not converge.  */
+	dithresh = 256.0 * machEp
+	goto ihalve
+
+done:
+
+	if rflg > 0 {
+		if x <= machEp {
+			x = 1.0 - machEp
+		} else {
+			x = 1.0 - x
+		}
+	}
+	return (x)
+}
+
+func lgam(a float64) float64 {
+	lg, _ := math.Lgamma(a)
+	return lg
+}

vendor/gonum.org/v1/gonum/mathext/internal/cephes/lanczos.go

@@ -0,0 +1,153 @@
+// Derived from SciPy's special/cephes/lanczos.c
+// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/lanczos.c
+
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+// Copyright ©2006 John Maddock
+// Portions Copyright ©2003 Boost
+// Portions Copyright ©2016 The Gonum Authors. All rights reserved.
+
+package cephes
+
+// Optimal values for G for each N are taken from
+// http://web.mala.bc.ca/pughg/phdThesis/phdThesis.pdf,
+// as are the theoretical error bounds.
+
+// Constants calculated using the method described by Godfrey
+// http://my.fit.edu/~gabdo/gamma.txt and elaborated by Toth at
+// http://www.rskey.org/gamma.htm using NTL::RR at 1000 bit precision.
+
+var lanczosNum = [...]float64{
+	2.506628274631000270164908177133837338626,
+	210.8242777515793458725097339207133627117,
+	8071.672002365816210638002902272250613822,
+	186056.2653952234950402949897160456992822,
+	2876370.628935372441225409051620849613599,
+	31426415.58540019438061423162831820536287,
+	248874557.8620541565114603864132294232163,
+	1439720407.311721673663223072794912393972,
+	6039542586.35202800506429164430729792107,
+	17921034426.03720969991975575445893111267,
+	35711959237.35566804944018545154716670596,
+	42919803642.64909876895789904700198885093,
+	23531376880.41075968857200767445163675473,
+}
+
+var lanczosDenom = [...]float64{
+	1,
+	66,
+	1925,
+	32670,
+	357423,
+	2637558,
+	13339535,
+	45995730,
+	105258076,
+	150917976,
+	120543840,
+	39916800,
+	0,
+}
+
+var lanczosSumExpgScaledNum = [...]float64{
+	0.006061842346248906525783753964555936883222,
+	0.5098416655656676188125178644804694509993,
+	19.51992788247617482847860966235652136208,
+	449.9445569063168119446858607650988409623,
+	6955.999602515376140356310115515198987526,
+	75999.29304014542649875303443598909137092,
+	601859.6171681098786670226533699352302507,
+	3481712.15498064590882071018964774556468,
+	14605578.08768506808414169982791359218571,
+	43338889.32467613834773723740590533316085,
+	86363131.28813859145546927288977868422342,
+	103794043.1163445451906271053616070238554,
+	56906521.91347156388090791033559122686859,
+}
+
+var lanczosSumExpgScaledDenom = [...]float64{
+	1,
+	66,
+	1925,
+	32670,
+	357423,
+	2637558,
+	13339535,
+	45995730,
+	105258076,
+	150917976,
+	120543840,
+	39916800,
+	0,
+}
+
+var lanczosSumNear1D = [...]float64{
+	0.3394643171893132535170101292240837927725e-9,
+	-0.2499505151487868335680273909354071938387e-8,
+	0.8690926181038057039526127422002498960172e-8,
+	-0.1933117898880828348692541394841204288047e-7,
+	0.3075580174791348492737947340039992829546e-7,
+	-0.2752907702903126466004207345038327818713e-7,
+	-0.1515973019871092388943437623825208095123e-5,
+	0.004785200610085071473880915854204301886437,
+	-0.1993758927614728757314233026257810172008,
+	1.483082862367253753040442933770164111678,
+	-3.327150580651624233553677113928873034916,
+	2.208709979316623790862569924861841433016,
+}
+
+var lanczosSumNear2D = [...]float64{
+	0.1009141566987569892221439918230042368112e-8,
+	-0.7430396708998719707642735577238449585822e-8,
+	0.2583592566524439230844378948704262291927e-7,
+	-0.5746670642147041587497159649318454348117e-7,
+	0.9142922068165324132060550591210267992072e-7,
+	-0.8183698410724358930823737982119474130069e-7,
+	-0.4506604409707170077136555010018549819192e-5,
+	0.01422519127192419234315002746252160965831,
+	-0.5926941084905061794445733628891024027949,
+	4.408830289125943377923077727900630927902,
+	-9.8907772644920670589288081640128194231,
+	6.565936202082889535528455955485877361223,
+}
+
+const lanczosG = 6.024680040776729583740234375
+
+func lanczosSum(x float64) float64 {
+	return ratevl(x,
+		lanczosNum[:],
+		len(lanczosNum)-1,
+		lanczosDenom[:],
+		len(lanczosDenom)-1)
+}
+
+func lanczosSumExpgScaled(x float64) float64 {
+	return ratevl(x,
+		lanczosSumExpgScaledNum[:],
+		len(lanczosSumExpgScaledNum)-1,
+		lanczosSumExpgScaledDenom[:],
+		len(lanczosSumExpgScaledDenom)-1)
+}
+
+func lanczosSumNear1(dx float64) float64 {
+	var result float64
+
+	for i, val := range lanczosSumNear1D {
+		k := float64(i + 1)
+		result += (-val * dx) / (k*dx + k*k)
+	}
+
+	return result
+}
+
+func lanczosSumNear2(dx float64) float64 {
+	var result float64
+	x := dx + 2
+
+	for i, val := range lanczosSumNear2D {
+		k := float64(i + 1)
+		result += (-val * dx) / (x + k*x + k*k - 1)
+	}
+
+	return result
+}

vendor/gonum.org/v1/gonum/mathext/internal/cephes/ndtri.go

@@ -0,0 +1,150 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+/*
+ * Cephes Math Library Release 2.1:  January, 1989
+ * Copyright 1984, 1987, 1989 by Stephen L. Moshier
+ * Direct inquiries to 30 Frost Street, Cambridge, MA 02140
+ */
+
+package cephes
+
+import "math"
+
+// TODO(btracey): There is currently an implementation of this functionality
+// in gonum/stat/distuv. Find out which implementation is better, and rectify
+// by having distuv call this, or moving this implementation into
+// gonum/mathext/internal/gonum.
+
+// math.Sqrt(2*pi)
+const s2pi = 2.50662827463100050242e0
+
+// approximation for 0 <= |y - 0.5| <= 3/8
+var P0 = [5]float64{
+	-5.99633501014107895267e1,
+	9.80010754185999661536e1,
+	-5.66762857469070293439e1,
+	1.39312609387279679503e1,
+	-1.23916583867381258016e0,
+}
+
+var Q0 = [8]float64{
+	/* 1.00000000000000000000E0, */
+	1.95448858338141759834e0,
+	4.67627912898881538453e0,
+	8.63602421390890590575e1,
+	-2.25462687854119370527e2,
+	2.00260212380060660359e2,
+	-8.20372256168333339912e1,
+	1.59056225126211695515e1,
+	-1.18331621121330003142e0,
+}
+
+// Approximation for interval z = math.Sqrt(-2 log y ) between 2 and 8
+// i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
+var P1 = [9]float64{
+	4.05544892305962419923e0,
+	3.15251094599893866154e1,
+	5.71628192246421288162e1,
+	4.40805073893200834700e1,
+	1.46849561928858024014e1,
+	2.18663306850790267539e0,
+	-1.40256079171354495875e-1,
+	-3.50424626827848203418e-2,
+	-8.57456785154685413611e-4,
+}
+
+var Q1 = [8]float64{
+	/*  1.00000000000000000000E0, */
+	1.57799883256466749731e1,
+	4.53907635128879210584e1,
+	4.13172038254672030440e1,
+	1.50425385692907503408e1,
+	2.50464946208309415979e0,
+	-1.42182922854787788574e-1,
+	-3.80806407691578277194e-2,
+	-9.33259480895457427372e-4,
+}
+
+// Approximation for interval z = math.Sqrt(-2 log y ) between 8 and 64
+// i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
+var P2 = [9]float64{
+	3.23774891776946035970e0,
+	6.91522889068984211695e0,
+	3.93881025292474443415e0,
+	1.33303460815807542389e0,
+	2.01485389549179081538e-1,
+	1.23716634817820021358e-2,
+	3.01581553508235416007e-4,
+	2.65806974686737550832e-6,
+	6.23974539184983293730e-9,
+}
+
+var Q2 = [8]float64{
+	/*  1.00000000000000000000E0, */
+	6.02427039364742014255e0,
+	3.67983563856160859403e0,
+	1.37702099489081330271e0,
+	2.16236993594496635890e-1,
+	1.34204006088543189037e-2,
+	3.28014464682127739104e-4,
+	2.89247864745380683936e-6,
+	6.79019408009981274425e-9,
+}
+
+// Ndtri returns the argument, x, for which the area under the
+// Gaussian probability density function (integrated from
+// minus infinity to x) is equal to y.
+func Ndtri(y0 float64) float64 {
+	// For small arguments 0 < y < exp(-2), the program computes
+	// z = math.Sqrt( -2.0 * math.Log(y) );  then the approximation is
+	// x = z - math.Log(z)/z  - (1/z) P(1/z) / Q(1/z).
+	// There are two rational functions P/Q, one for 0 < y < exp(-32)
+	// and the other for y up to exp(-2).  For larger arguments,
+	// w = y - 0.5, and  x/math.Sqrt(2pi) = w + w**3 R(w**2)/S(w**2)).
+	var x, y, z, y2, x0, x1 float64
+	var code int
+
+	if y0 <= 0.0 {
+		if y0 < 0 {
+			panic(paramOutOfBounds)
+		}
+		return math.Inf(-1)
+	}
+	if y0 >= 1.0 {
+		if y0 > 1 {
+			panic(paramOutOfBounds)
+		}
+		return math.Inf(1)
+	}
+	code = 1
+	y = y0
+	if y > (1.0 - 0.13533528323661269189) { /* 0.135... = exp(-2) */
+		y = 1.0 - y
+		code = 0
+	}
+
+	if y > 0.13533528323661269189 {
+		y = y - 0.5
+		y2 = y * y
+		x = y + y*(y2*polevl(y2, P0[:], 4)/p1evl(y2, Q0[:], 8))
+		x = x * s2pi
+		return (x)
+	}
+
+	x = math.Sqrt(-2.0 * math.Log(y))
+	x0 = x - math.Log(x)/x
+
+	z = 1.0 / x
+	if x < 8.0 { /* y > exp(-32) = 1.2664165549e-14 */
+		x1 = z * polevl(z, P1[:], 8) / p1evl(z, Q1[:], 8)
+	} else {
+		x1 = z * polevl(z, P2[:], 8) / p1evl(z, Q2[:], 8)
+	}
+	x = x0 - x1
+	if code != 0 {
+		x = -x
+	}
+	return (x)
+}

vendor/gonum.org/v1/gonum/mathext/internal/cephes/polevl.go

@@ -0,0 +1,84 @@
+// Derived from SciPy's special/cephes/polevl.h
+// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/polevl.h
+// Made freely available by Stephen L. Moshier without support or guarantee.
+
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+// Copyright ©1984, ©1987, ©1988 by Stephen L. Moshier
+// Portions Copyright ©2016 The Gonum Authors. All rights reserved.
+
+package cephes
+
+import "math"
+
+// polevl evaluates a polynomial of degree N
+//
+//	y = c_0 + c_1 x_1 + c_2 x_2^2 ...
+//
+// where the coefficients are stored in reverse order, i.e. coef[0] = c_n and
+// coef[n] = c_0.
+func polevl(x float64, coef []float64, n int) float64 {
+	ans := coef[0]
+	for i := 1; i <= n; i++ {
+		ans = ans*x + coef[i]
+	}
+	return ans
+}
+
+// p1evl is the same as polevl, except c_n is assumed to be 1 and is not included
+// in the slice.
+func p1evl(x float64, coef []float64, n int) float64 {
+	ans := x + coef[0]
+	for i := 1; i <= n-1; i++ {
+		ans = ans*x + coef[i]
+	}
+	return ans
+}
+
+// ratevl evaluates a rational function
+func ratevl(x float64, num []float64, m int, denom []float64, n int) float64 {
+	// Source: Holin et. al., "Polynomial and Rational Function Evaluation",
+	// http://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/roots/rational.html
+	absx := math.Abs(x)
+
+	var dir, idx int
+	var y float64
+	if absx > 1 {
+		// Evaluate as a polynomial in 1/x
+		dir = -1
+		idx = m
+		y = 1 / x
+	} else {
+		dir = 1
+		idx = 0
+		y = x
+	}
+
+	// Evaluate the numerator
+	numAns := num[idx]
+	idx += dir
+	for i := 0; i < m; i++ {
+		numAns = numAns*y + num[idx]
+		idx += dir
+	}
+
+	// Evaluate the denominator
+	if absx > 1 {
+		idx = n
+	} else {
+		idx = 0
+	}
+
+	denomAns := denom[idx]
+	idx += dir
+	for i := 0; i < n; i++ {
+		denomAns = denomAns*y + denom[idx]
+		idx += dir
+	}
+
+	if absx > 1 {
+		pow := float64(n - m)
+		return math.Pow(x, pow) * numAns / denomAns
+	}
+	return numAns / denomAns
+}

vendor/gonum.org/v1/gonum/mathext/internal/cephes/staticcheck.conf

@@ -0,0 +1 @@
+checks = []

vendor/gonum.org/v1/gonum/mathext/internal/cephes/unity.go

@@ -0,0 +1,184 @@
+// Derived from SciPy's special/cephes/unity.c
+// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/unity.c
+// Made freely available by Stephen L. Moshier without support or guarantee.
+
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+// Copyright ©1984, ©1996 by Stephen L. Moshier
+// Portions Copyright ©2016 The Gonum Authors. All rights reserved.
+
+package cephes
+
+import "math"
+
+// Relative error approximations for function arguments near unity.
+//  log1p(x) = log(1+x)
+//  expm1(x) = exp(x) - 1
+//  cosm1(x) = cos(x) - 1
+//  lgam1p(x) = lgam(1+x)
+
+const (
+	invSqrt2 = 1 / math.Sqrt2
+	pi4      = math.Pi / 4
+	euler    = 0.577215664901532860606512090082402431 // Euler constant
+)
+
+// Coefficients for
+//
+//	log(1+x) = x - \frac{x^2}{2} + \frac{x^3 lP(x)}{lQ(x)}
+//
+// for
+//
+//	\frac{1}{\sqrt{2}} <= x < \sqrt{2}
+//
+// Theoretical peak relative error = 2.32e-20
+var lP = [...]float64{
+	4.5270000862445199635215e-5,
+	4.9854102823193375972212e-1,
+	6.5787325942061044846969e0,
+	2.9911919328553073277375e1,
+	6.0949667980987787057556e1,
+	5.7112963590585538103336e1,
+	2.0039553499201281259648e1,
+}
+
+var lQ = [...]float64{
+	1.5062909083469192043167e1,
+	8.3047565967967209469434e1,
+	2.2176239823732856465394e2,
+	3.0909872225312059774938e2,
+	2.1642788614495947685003e2,
+	6.0118660497603843919306e1,
+}
+
+// log1p computes
+//
+//	log(1 + x)
+func log1p(x float64) float64 {
+	z := 1 + x
+	if z < invSqrt2 || z > math.Sqrt2 {
+		return math.Log(z)
+	}
+	z = x * x
+	z = -0.5*z + x*(z*polevl(x, lP[:], 6)/p1evl(x, lQ[:], 6))
+	return x + z
+}
+
+// log1pmx computes
+//
+//	log(1 + x) - x
+func log1pmx(x float64) float64 {
+	if math.Abs(x) < 0.5 {
+		xfac := x
+		res := 0.0
+
+		var term float64
+		for n := 2; n < maxIter; n++ {
+			xfac *= -x
+			term = xfac / float64(n)
+			res += term
+			if math.Abs(term) < machEp*math.Abs(res) {
+				break
+			}
+		}
+		return res
+	}
+	return log1p(x) - x
+}
+
+// Coefficients for
+//
+//	e^x = 1 + \frac{2x eP(x^2)}{eQ(x^2) - eP(x^2)}
+//
+// for
+//
+//	-0.5 <= x <= 0.5
+var eP = [...]float64{
+	1.2617719307481059087798e-4,
+	3.0299440770744196129956e-2,
+	9.9999999999999999991025e-1,
+}
+
+var eQ = [...]float64{
+	3.0019850513866445504159e-6,
+	2.5244834034968410419224e-3,
+	2.2726554820815502876593e-1,
+	2.0000000000000000000897e0,
+}
+
+// expm1 computes
+//
+//	expm1(x) = e^x - 1
+func expm1(x float64) float64 {
+	if math.IsInf(x, 0) {
+		if math.IsNaN(x) || x > 0 {
+			return x
+		}
+		return -1
+	}
+	if x < -0.5 || x > 0.5 {
+		return math.Exp(x) - 1
+	}
+	xx := x * x
+	r := x * polevl(xx, eP[:], 2)
+	r = r / (polevl(xx, eQ[:], 3) - r)
+	return r + r
+}
+
+var coscof = [...]float64{
+	4.7377507964246204691685e-14,
+	-1.1470284843425359765671e-11,
+	2.0876754287081521758361e-9,
+	-2.7557319214999787979814e-7,
+	2.4801587301570552304991e-5,
+	-1.3888888888888872993737e-3,
+	4.1666666666666666609054e-2,
+}
+
+// cosm1 computes
+//
+//	cosm1(x) = cos(x) - 1
+func cosm1(x float64) float64 {
+	if x < -pi4 || x > pi4 {
+		return math.Cos(x) - 1
+	}
+	xx := x * x
+	xx = -0.5*xx + xx*xx*polevl(xx, coscof[:], 6)
+	return xx
+}
+
+// lgam1pTayler computes
+//
+//	lgam(x + 1)
+//
+// around x = 0 using its Taylor series.
+func lgam1pTaylor(x float64) float64 {
+	if x == 0 {
+		return 0
+	}
+	res := -euler * x
+	xfac := -x
+	for n := 2; n < 42; n++ {
+		nf := float64(n)
+		xfac *= -x
+		coeff := Zeta(nf, 1) * xfac / nf
+		res += coeff
+		if math.Abs(coeff) < machEp*math.Abs(res) {
+			break
+		}
+	}
+
+	return res
+}
+
+// lgam1p computes
+//
+//	lgam(x + 1)
+func lgam1p(x float64) float64 {
+	if math.Abs(x) <= 0.5 {
+		return lgam1pTaylor(x)
+	} else if math.Abs(x-1) < 0.5 {
+		return math.Log(x) + lgam1pTaylor(x-1)
+	}
+	return lgam(x + 1)
+}

vendor/gonum.org/v1/gonum/mathext/internal/cephes/zeta.go

@@ -0,0 +1,117 @@
+// Derived from SciPy's special/cephes/zeta.c
+// https://github.com/scipy/scipy/blob/master/scipy/special/cephes/zeta.c
+// Made freely available by Stephen L. Moshier without support or guarantee.
+
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+// Copyright ©1984, ©1987 by Stephen L. Moshier
+// Portions Copyright ©2016 The Gonum Authors. All rights reserved.
+
+package cephes
+
+import "math"
+
+// zetaCoegs are the expansion coefficients for Euler-Maclaurin summation
+// formula:
+//
+//	\frac{(2k)!}{B_{2k}}
+//
+// where
+//
+//	B_{2k}
+//
+// are Bernoulli numbers.
+var zetaCoefs = [...]float64{
+	12.0,
+	-720.0,
+	30240.0,
+	-1209600.0,
+	47900160.0,
+	-1.307674368e12 / 691,
+	7.47242496e10,
+	-1.067062284288e16 / 3617,
+	5.109094217170944e18 / 43867,
+	-8.028576626982912e20 / 174611,
+	1.5511210043330985984e23 / 854513,
+	-1.6938241367317436694528e27 / 236364091,
+}
+
+// Zeta computes the Riemann zeta function of two arguments.
+//
+//	Zeta(x,q) = \sum_{k=0}^{\infty} (k+q)^{-x}
+//
+// Note that Zeta returns +Inf if x is 1 and will panic if x is less than 1,
+// q is either zero or a negative integer, or q is negative and x is not an
+// integer.
+//
+// Note that:
+//
+//	zeta(x,1) = zetac(x) + 1
+func Zeta(x, q float64) float64 {
+	// REFERENCE: Gradshteyn, I. S., and I. M. Ryzhik, Tables of Integrals, Series,
+	// and Products, p. 1073; Academic Press, 1980.
+	if x == 1 {
+		return math.Inf(1)
+	}
+
+	if x < 1 {
+		panic(paramOutOfBounds)
+	}
+
+	if q <= 0 {
+		if q == math.Floor(q) {
+			panic(errParamFunctionSingularity)
+		}
+		if x != math.Floor(x) {
+			panic(paramOutOfBounds) // Because q^-x not defined
+		}
+	}
+
+	// Asymptotic expansion: http://dlmf.nist.gov/25.11#E43
+	if q > 1e8 {
+		return (1/(x-1) + 1/(2*q)) * math.Pow(q, 1-x)
+	}
+
+	// The Euler-Maclaurin summation formula is used to obtain the expansion:
+	//  Zeta(x,q) = \sum_{k=1}^n (k+q)^{-x} + \frac{(n+q)^{1-x}}{x-1} - \frac{1}{2(n+q)^x} + \sum_{j=1}^{\infty} \frac{B_{2j}x(x+1)...(x+2j)}{(2j)! (n+q)^{x+2j+1}}
+	// where
+	//  B_{2j}
+	// are Bernoulli numbers.
+	// Permit negative q but continue sum until n+q > 9. This case should be
+	// handled by a reflection formula. If q<0 and x is an integer, there is a
+	// relation to the polyGamma function.
+	s := math.Pow(q, -x)
+	a := q
+	i := 0
+	b := 0.0
+	for i < 9 || a <= 9 {
+		i++
+		a += 1.0
+		b = math.Pow(a, -x)
+		s += b
+		if math.Abs(b/s) < machEp {
+			return s
+		}
+	}
+
+	w := a
+	s += b * w / (x - 1)
+	s -= 0.5 * b
+	a = 1.0
+	k := 0.0
+	for _, coef := range zetaCoefs {
+		a *= x + k
+		b /= w
+		t := a * b / coef
+		s = s + t
+		t = math.Abs(t / s)
+		if t < machEp {
+			return s
+		}
+		k += 1.0
+		a *= x + k
+		b /= w
+		k += 1.0
+	}
+	return s
+}

vendor/gonum.org/v1/gonum/mathext/internal/gonum/beta.go

@@ -0,0 +1,58 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package gonum
+
+import (
+	"math"
+)
+
+// Beta returns the value of the complete beta function B(a, b). It is defined as
+//
+//	Γ(a)Γ(b) / Γ(a+b)
+//
+// Special cases are:
+//
+//	B(a,b) returns NaN if a or b is Inf
+//	B(a,b) returns NaN if a and b are 0
+//	B(a,b) returns NaN if a or b is NaN
+//	B(a,b) returns NaN if a or b is < 0
+//	B(a,b) returns +Inf if a xor b is 0.
+//
+// See http://mathworld.wolfram.com/BetaFunction.html for more detailed information.
+func Beta(a, b float64) float64 {
+	return math.Exp(Lbeta(a, b))
+}
+
+// Lbeta returns the natural logarithm of the complete beta function B(a,b).
+// Lbeta is defined as:
+//
+//	Ln(Γ(a)Γ(b)/Γ(a+b))
+//
+// Special cases are:
+//
+//	Lbeta(a,b) returns NaN if a or b is Inf
+//	Lbeta(a,b) returns NaN if a and b are 0
+//	Lbeta(a,b) returns NaN if a or b is NaN
+//	Lbeta(a,b) returns NaN if a or b is < 0
+//	Lbeta(a,b) returns +Inf if a xor b is 0.
+func Lbeta(a, b float64) float64 {
+	switch {
+	case math.IsInf(a, +1) || math.IsInf(b, +1):
+		return math.NaN()
+	case a == 0 && b == 0:
+		return math.NaN()
+	case a < 0 || b < 0:
+		return math.NaN()
+	case math.IsNaN(a) || math.IsNaN(b):
+		return math.NaN()
+	case a == 0 || b == 0:
+		return math.Inf(+1)
+	}
+
+	la, _ := math.Lgamma(a)
+	lb, _ := math.Lgamma(b)
+	lab, _ := math.Lgamma(a + b)
+	return la + lb - lab
+}

vendor/gonum.org/v1/gonum/mathext/internal/gonum/doc.go

@@ -0,0 +1,7 @@
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+// Package gonum contains functions implemented by the gonum team.
+// It is here to avoid circular imports and/or double coding of functions.
+package gonum // import "gonum.org/v1/gonum/mathext/internal/gonum"

vendor/gonum.org/v1/gonum/mathext/internal/gonum/gonum.go

@@ -0,0 +1,5 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package gonum

vendor/gonum.org/v1/gonum/mathext/mvgamma.go

@@ -0,0 +1,32 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import "math"
+
+const (
+	logPi = 1.14472988584940017414342735135305871164729481 // http://oeis.org/A053510
+)
+
+// MvLgamma returns the log of the multivariate Gamma function. Dim
+// must be greater than zero, and MvLgamma will return NaN if v < (dim-1)/2.
+//
+// See https://en.wikipedia.org/wiki/Multivariate_gamma_function for more
+// information.
+func MvLgamma(v float64, dim int) float64 {
+	if dim < 1 {
+		panic("mathext: negative dimension")
+	}
+	df := float64(dim)
+	if v < (df-1)*0.5 {
+		return math.NaN()
+	}
+	ans := df * (df - 1) * 0.25 * logPi
+	for i := 1; i <= dim; i++ {
+		lg, _ := math.Lgamma(v + float64(1-i)*0.5)
+		ans += lg
+	}
+	return ans
+}

vendor/gonum.org/v1/gonum/mathext/roots.go

@@ -0,0 +1,181 @@
+// Derived from SciPy's special/c_misc/fsolve.c and special/c_misc/misc.h
+// https://github.com/scipy/scipy/blob/master/scipy/special/c_misc/fsolve.c
+// https://github.com/scipy/scipy/blob/master/scipy/special/c_misc/misc.h
+
+// Copyright ©2017 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import "math"
+
+type objectiveFunc func(float64, []float64) float64
+
+type fSolveResult uint8
+
+const (
+	// An exact solution was found, in which case the first point on the
+	// interval is the value
+	fSolveExact fSolveResult = iota + 1
+	// Interval width is less than the tolerance
+	fSolveConverged
+	// Root-finding didn't converge in a set number of iterations
+	fSolveMaxIterations
+)
+
+const (
+	machEp = 1.0 / (1 << 53)
+)
+
+// falsePosition uses a combination of bisection and false position to find a
+// root of a function within a given interval. This is guaranteed to converge,
+// and always keeps a bounding interval, unlike Newton's method. Inputs are:
+//
+//	x1, x2:   initial bounding interval
+//	f1, f2: value of f() at x1 and x2
+//	absErr, relErr: absolute and relative errors on the bounding interval
+//	bisectTil: if > 0.0, perform bisection until the width of the bounding
+//	           interval is less than this
+//	f, fExtra: function to find root of is f(x, fExtra)
+//
+// Returns:
+//
+//	result: whether an exact root was found, the process converged to a
+//	        bounding interval small than the required error, or the max number
+//	        of iterations was hit
+//	bestX: best root approximation
+//	bestF: function value at bestX
+//	errEst: error estimation
+func falsePosition(x1, x2, f1, f2, absErr, relErr, bisectTil float64, f objectiveFunc, fExtra []float64) (fSolveResult, float64, float64, float64) {
+	// The false position steps are either unmodified, or modified with the
+	// Anderson-Bjorck method as appropriate. Theoretically, this has a "speed of
+	// convergence" of 1.7 (bisection is 1, Newton is 2).
+	// Note that this routine was designed initially to work with gammaincinv, so
+	// it may not be tuned right for other problems. Don't use it blindly.
+
+	if f1*f2 >= 0 {
+		panic("Initial interval is not a bounding interval")
+	}
+
+	const (
+		maxIterations = 100
+		bisectIter    = 4
+		bisectWidth   = 4.0
+	)
+
+	const (
+		bisect = iota + 1
+		falseP
+	)
+
+	var state uint8
+	if bisectTil > 0 {
+		state = bisect
+	} else {
+		state = falseP
+	}
+
+	gamma := 1.0
+
+	w := math.Abs(x2 - x1)
+	lastBisectWidth := w
+
+	var nFalseP int
+	var x3, f3, bestX, bestF float64
+	for i := 0; i < maxIterations; i++ {
+		switch state {
+		case bisect:
+			x3 = 0.5 * (x1 + x2)
+			if x3 == x1 || x3 == x2 {
+				// i.e., x1 and x2 are successive floating-point numbers
+				bestX = x3
+				if x3 == x1 {
+					bestF = f1
+				} else {
+					bestF = f2
+				}
+				return fSolveConverged, bestX, bestF, w
+			}
+
+			f3 = f(x3, fExtra)
+			if f3 == 0 {
+				return fSolveExact, x3, f3, w
+			}
+
+			if f3*f2 < 0 {
+				x1 = x2
+				f1 = f2
+			}
+			x2 = x3
+			f2 = f3
+			w = math.Abs(x2 - x1)
+			lastBisectWidth = w
+			if bisectTil > 0 {
+				if w < bisectTil {
+					bisectTil = -1.0
+					gamma = 1.0
+					nFalseP = 0
+					state = falseP
+				}
+			} else {
+				gamma = 1.0
+				nFalseP = 0
+				state = falseP
+			}
+		case falseP:
+			s12 := (f2 - gamma*f1) / (x2 - x1)
+			x3 = x2 - f2/s12
+			f3 = f(x3, fExtra)
+			if f3 == 0 {
+				return fSolveExact, x3, f3, w
+			}
+
+			nFalseP++
+			if f3*f2 < 0 {
+				gamma = 1.0
+				x1 = x2
+				f1 = f2
+			} else {
+				// Anderson-Bjorck method
+				g := 1.0 - f3/f2
+				if g <= 0 {
+					g = 0.5
+				}
+				gamma *= g
+			}
+			x2 = x3
+			f2 = f3
+			w = math.Abs(x2 - x1)
+
+			// Sanity check. For every 4 false position checks, see if we really are
+			// decreasing the interval by comparing to what bisection would have
+			// achieved (or, rather, a bit more lenient than that -- interval
+			// decreased by 4 instead of by 16, as the fp could be decreasing gamma
+			// for a bit). Note that this should guarantee convergence, as it makes
+			// sure that we always end up decreasing the interval width with a
+			// bisection.
+			if nFalseP > bisectIter {
+				if w*bisectWidth > lastBisectWidth {
+					state = bisect
+				}
+				nFalseP = 0
+				lastBisectWidth = w
+			}
+		}
+
+		tol := absErr + relErr*math.Max(math.Max(math.Abs(x1), math.Abs(x2)), 1.0)
+		if w <= tol {
+			if math.Abs(f1) < math.Abs(f2) {
+				bestX = x1
+				bestF = f1
+			} else {
+				bestX = x2
+				bestF = f2
+			}
+			return fSolveConverged, bestX, bestF, w
+		}
+	}
+
+	return fSolveMaxIterations, x3, f3, w
+}

vendor/gonum.org/v1/gonum/mathext/zeta.go

@@ -0,0 +1,22 @@
+// Copyright ©2016 The Gonum Authors. All rights reserved.
+// Use of this source code is governed by a BSD-style
+// license that can be found in the LICENSE file.
+
+package mathext
+
+import "gonum.org/v1/gonum/mathext/internal/cephes"
+
+// Zeta computes the Riemann zeta function of two arguments.
+//
+//	Zeta(x,q) = \sum_{k=0}^{\infty} (k+q)^{-x}
+//
+// Note that Zeta returns +Inf if x is 1 and will panic if x is less than 1,
+// q is either zero or a negative integer, or q is negative and x is not an
+// integer.
+//
+// See http://mathworld.wolfram.com/HurwitzZetaFunction.html
+// or https://en.wikipedia.org/wiki/Multiple_zeta_function#Two_parameters_case
+// for more detailed information.
+func Zeta(x, q float64) float64 {
+	return cephes.Zeta(x, q)
+}