aboutsummaryrefslogtreecommitdiff
path: root/vendor/github.com/aclements/go-moremath/mathx/gamma.go
diff options
context:
space:
mode:
Diffstat (limited to 'vendor/github.com/aclements/go-moremath/mathx/gamma.go')
-rw-r--r--vendor/github.com/aclements/go-moremath/mathx/gamma.go96
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")
-}