1 files changed, 96 insertions, 0 deletions
diff --git a/vendor/github.com/aclements/go-moremath/mathx/gamma.go b/vendor/github.com/aclements/go-moremath/mathx/gamma.go
new file mode 100644
index 0000000..d11096e
--- /dev/null
+++ b/vendor/github.com/aclements/go-moremath/mathx/gamma.go
@@ -0,0 +1,96 @@
+// 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")
+}