diff options
Diffstat (limited to 'vendor/github.com/aclements/go-moremath/mathx/gamma.go')
| -rw-r--r-- | vendor/github.com/aclements/go-moremath/mathx/gamma.go | 96 |
1 files changed, 0 insertions, 96 deletions
diff --git a/vendor/github.com/aclements/go-moremath/mathx/gamma.go b/vendor/github.com/aclements/go-moremath/mathx/gamma.go deleted file mode 100644 index d11096e..0000000 --- a/vendor/github.com/aclements/go-moremath/mathx/gamma.go +++ /dev/null @@ -1,96 +0,0 @@ -// Copyright 2015 The Go 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 mathx - -import "math" - -// GammaInc returns the value of the incomplete gamma function (also -// known as the regularized gamma function): -// -// P(a, x) = 1 / Γ(a) * ∫₀ˣ exp(-t) t**(a-1) dt -func GammaInc(a, x float64) float64 { - // Based on Numerical Recipes in C, section 6.2. - - if a <= 0 || x < 0 || math.IsNaN(a) || math.IsNaN(x) { - return math.NaN() - } - - if x < a+1 { - // Use the series representation, which converges more - // rapidly in this range. - return gammaIncSeries(a, x) - } else { - // Use the continued fraction representation. - return 1 - gammaIncCF(a, x) - } -} - -// GammaIncComp returns the complement of the incomplete gamma -// function 1 - GammaInc(a, x). This is more numerically stable for -// values near 0. -func GammaIncComp(a, x float64) float64 { - if a <= 0 || x < 0 || math.IsNaN(a) || math.IsNaN(x) { - return math.NaN() - } - - if x < a+1 { - return 1 - gammaIncSeries(a, x) - } else { - return gammaIncCF(a, x) - } -} - -func gammaIncSeries(a, x float64) float64 { - const maxIterations = 200 - const epsilon = 3e-14 - - if x == 0 { - return 0 - } - - ap := a - del := 1 / a - sum := del - for n := 0; n < maxIterations; n++ { - ap++ - del *= x / ap - sum += del - if math.Abs(del) < math.Abs(sum)*epsilon { - return sum * math.Exp(-x+a*math.Log(x)-lgamma(a)) - } - } - panic("a too large; failed to converge") -} - -func gammaIncCF(a, x float64) float64 { - const maxIterations = 200 - const epsilon = 3e-14 - - raiseZero := func(z float64) float64 { - if math.Abs(z) < math.SmallestNonzeroFloat64 { - return math.SmallestNonzeroFloat64 - } - return z - } - - b := x + 1 - a - c := math.MaxFloat64 - d := 1 / b - h := d - - for i := 1; i <= maxIterations; i++ { - an := -float64(i) * (float64(i) - a) - b += 2 - d = raiseZero(an*d + b) - c = raiseZero(b + an/c) - d = 1 / d - del := d * c - h *= del - if math.Abs(del-1) < epsilon { - return math.Exp(-x+a*math.Log(x)-lgamma(a)) * h - } - } - panic("a too large; failed to converge") -} |
