mirror of
https://github.com/gonum/gonum.git
synced 2025-10-15 03:30:39 +08:00
Merge pull request #24 from zeroviscosity/master
Adding Incomplete Gamma and Zeta
This commit is contained in:
Brendan Tracey
50
gamma_inc.go
Normal file
50
gamma_inc.go
Normal file
@@ -0,0 +1,50 @@
|
|||||||
|
// 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 (
|
||||||
|
"github.com/gonum/mathext/internal/cephes"
|
||||||
|
)
|
||||||
|
|
||||||
|
// GammaInc computes the incomplete Gamma integral.
|
||||||
|
// GammaInc(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 GammaInc
|
||||||
|
// will panic.
|
||||||
|
//
|
||||||
|
// See http://mathworld.wolfram.com/IncompleteGammaFunction.html
|
||||||
|
// or https://en.wikipedia.org/wiki/Incomplete_gamma_function for more detailed
|
||||||
|
// information.
|
||||||
|
func GammaInc(a, x float64) float64 {
|
||||||
|
return cephes.Igam(a, x)
|
||||||
|
}
|
||||||
|
|
||||||
|
// GammaIncComp computes the complemented incomplete Gamma integral.
|
||||||
|
// GammaIncComp(a,x) = 1 - GammaInc(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
|
||||||
|
// GammaIncComp will panic.
|
||||||
|
func GammaIncComp(a, x float64) float64 {
|
||||||
|
return cephes.IgamC(a, x)
|
||||||
|
}
|
||||||
|
|
||||||
|
// GammaIncInv computes the inverse of the incomplete Gamma integral. That is,
|
||||||
|
// it returns the x such that:
|
||||||
|
// GammaInc(a, x) = y
|
||||||
|
// The input argument a must be positive and y must be between 0 and 1
|
||||||
|
// inclusive or GammaIncInv will panic. GammaIncInv should return a positive
|
||||||
|
// number, but can return NaN if there is a failure to converge.
|
||||||
|
func GammaIncInv(a, y float64) float64 {
|
||||||
|
return gammaIncInv(a, y)
|
||||||
|
}
|
||||||
|
|
||||||
|
// GammaIncCompInv computes the inverse of the complemented incomplete Gamma
|
||||||
|
// integral. That is, it returns the x such that:
|
||||||
|
// GammaIncComp(a, x) = y
|
||||||
|
// The input argument a must be positive and y must be between 0 and 1
|
||||||
|
// inclusive or GammaIncCompInv will panic. GammaIncCompInv should return a
|
||||||
|
// positive number, but can return 0 even with non-zero y due to underflow.
|
||||||
|
func GammaIncCompInv(a, y float64) float64 {
|
||||||
|
return cephes.IgamI(a, y)
|
||||||
|
}
|
||||||
56
gamma_inc_inv.go
Normal file
56
gamma_inc_inv.go
Normal file
@@ -0,0 +1,56 @@
|
|||||||
|
// 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"
|
||||||
|
|
||||||
|
"github.com/gonum/mathext/internal/cephes"
|
||||||
|
)
|
||||||
|
|
||||||
|
const (
|
||||||
|
allowedATol = 1e-306
|
||||||
|
allowedRTol = 1e-6
|
||||||
|
)
|
||||||
|
|
||||||
|
func gammaInc(x float64, params []float64) float64 {
|
||||||
|
return cephes.Igam(params[0], x) - params[1]
|
||||||
|
}
|
||||||
|
|
||||||
|
// gammaIncInv is the inverse of the 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 gammaIncInv will panic. gammaIncInv should return a
|
||||||
|
// positive number, but can return NaN if there is a failure to converge.
|
||||||
|
func gammaIncInv(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 similiar for other values).
|
||||||
|
result, bestX, _, errEst := falsePosition(lo, hi, flo, fhi, 2*machEp, 2*machEp, 1e-2*a, gammaInc, params)
|
||||||
|
if result == fSolveMaxIterations && errEst > allowedATol+allowedRTol*math.Abs(bestX) {
|
||||||
|
bestX = math.NaN()
|
||||||
|
}
|
||||||
|
|
||||||
|
return bestX
|
||||||
|
}
|
||||||
138
gamma_inc_test.go
Normal file
138
gamma_inc_test.go
Normal file
@@ -0,0 +1,138 @@
|
|||||||
|
// 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"
|
||||||
|
"testing"
|
||||||
|
)
|
||||||
|
|
||||||
|
func TestGammaInc(t *testing.T) {
|
||||||
|
for i, test := range []struct {
|
||||||
|
a, x, want float64
|
||||||
|
}{
|
||||||
|
// Results computed using scipy.special.gamminc
|
||||||
|
{0, 0, 0},
|
||||||
|
{0.0001, 1, 0.99997805936186279},
|
||||||
|
{0.001, 0.005, 0.99528424172333985},
|
||||||
|
{0.01, 10, 0.99999995718295021},
|
||||||
|
{0.1, 10, 0.99999944520142825},
|
||||||
|
{0.25, 0.75, 0.89993651328449831},
|
||||||
|
{0.5, 0.5, 0.68268949213708596},
|
||||||
|
{0.5, 2, 0.95449973610364147},
|
||||||
|
{0.75, 2.5, 0.95053039734695643},
|
||||||
|
{1, 0.5, 0.39346934028736652},
|
||||||
|
{1, 1, 0.63212055882855778},
|
||||||
|
{1.5, 0.75, 0.31772966966378746},
|
||||||
|
{2.5, 1, 0.15085496391539038},
|
||||||
|
{3, 0.05, 2.0067493624397931e-05},
|
||||||
|
{3, 20, 0.99999954448504946},
|
||||||
|
{5, 50, 1},
|
||||||
|
{7, 10, 0.86985857911751696},
|
||||||
|
{10, 0.9, 4.2519575433351128e-08},
|
||||||
|
{10, 5, 0.031828057306204811},
|
||||||
|
{25, 10, 4.6949381426799868e-05},
|
||||||
|
} {
|
||||||
|
if got := GammaInc(test.a, test.x); math.Abs(got-test.want) > 1e-10 {
|
||||||
|
t.Errorf("test %d GammaInc(%g, %g) failed: got %g want %g", i, test.a, test.x, got, test.want)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestGammaIncComp(t *testing.T) {
|
||||||
|
for i, test := range []struct {
|
||||||
|
a, x, want float64
|
||||||
|
}{
|
||||||
|
// Results computed using scipy.special.gammincc
|
||||||
|
{0.00001, 0.075, 2.0866541002417804e-05},
|
||||||
|
{0.0001, 1, 2.1940638138146658e-05},
|
||||||
|
{0.001, 0.005, 0.0047157582766601536},
|
||||||
|
{0.01, 0.9, 0.0026263432520514662},
|
||||||
|
{0.25, 0.75, 0.10006348671550169},
|
||||||
|
{0.5, 0.5, 0.31731050786291404},
|
||||||
|
{0.75, 0.25, 0.65343980284081038},
|
||||||
|
{0.9, 0.01, 0.98359881081593148},
|
||||||
|
{1, 0, 1},
|
||||||
|
{1, 0.075, 0.92774348632855297},
|
||||||
|
{1, 1, 0.36787944117144233},
|
||||||
|
{1, 10, 4.5399929762484861e-05},
|
||||||
|
{1, math.Inf(1), 0},
|
||||||
|
{3, 20, 4.5551495055892125e-07},
|
||||||
|
{5, 10, 0.029252688076961127},
|
||||||
|
{10, 3, 0.99889751186988451},
|
||||||
|
{50, 25, 0.99999304669475242},
|
||||||
|
{100, 10, 1},
|
||||||
|
{500, 500, 0.49405285382921321},
|
||||||
|
{500, 550, 0.014614408126291296},
|
||||||
|
} {
|
||||||
|
if got := GammaIncComp(test.a, test.x); math.Abs(got-test.want) > 1e-10 {
|
||||||
|
t.Errorf("test %d GammaIncComp(%g, %g) failed: got %g want %g", i, test.a, test.x, got, test.want)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestGammaIncInv(t *testing.T) {
|
||||||
|
for i, test := range []struct {
|
||||||
|
a, x, want float64
|
||||||
|
}{
|
||||||
|
// Results computed using scipy.special.gammincinv
|
||||||
|
{0.001, 0.99, 2.4259428385570885e-05},
|
||||||
|
{0.01, 0.99, 0.26505255025157959},
|
||||||
|
{0.1, 0.5, 0.00059339110446022798},
|
||||||
|
{0.2, 0.8, 0.26354363204872067},
|
||||||
|
{0.25, 0.5, 0.043673802352873381},
|
||||||
|
{0.5, 0.25, 0.050765522133810789},
|
||||||
|
{0.5, 0.5, 0.22746821155978625},
|
||||||
|
{0.75, 0.25, 0.15340752707472377},
|
||||||
|
{1, 0, 0},
|
||||||
|
{1, 0.075, 0.077961541469711862},
|
||||||
|
{1, 1, math.Inf(1)},
|
||||||
|
{2.5, 0.99, 7.5431362346944937},
|
||||||
|
{10, 0.5, 9.6687146147141299},
|
||||||
|
{25, 0.01, 14.853341349420646},
|
||||||
|
{25, 0.99, 38.076945624506337},
|
||||||
|
{50, 0.75, 54.570620535040511},
|
||||||
|
{100, 0.25, 93.08583383712174},
|
||||||
|
{1000, 0.01, 927.90815979664251},
|
||||||
|
{1000, 0.99, 1075.0328320864389},
|
||||||
|
{10000, 0.5, 9999.6666686420485},
|
||||||
|
} {
|
||||||
|
if got := GammaIncInv(test.a, test.x); math.Abs(got-test.want) > 1e-10 {
|
||||||
|
t.Errorf("test %d GammaIncInv(%g, %g) failed: got %g want %g", i, test.a, test.x, got, test.want)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
func TestGammaIncCompInv(t *testing.T) {
|
||||||
|
for i, test := range []struct {
|
||||||
|
a, x, want float64
|
||||||
|
}{
|
||||||
|
// Results computed using scipy.special.gamminccinv
|
||||||
|
{0.001, 0.01, 2.4259428385570885e-05},
|
||||||
|
{0.01, 0.01, 0.26505255025158292},
|
||||||
|
{0.03, 0.4, 2.316980536227699e-08},
|
||||||
|
{0.1, 0.5, 0.00059339110446022798},
|
||||||
|
{0.1, 0.75, 5.7917132949696076e-07},
|
||||||
|
{0.25, 0.25, 0.26062600197823282},
|
||||||
|
{0.5, 0.1, 1.3527717270477047},
|
||||||
|
{0.5, 0.5, 0.22746821155978625},
|
||||||
|
{0.75, 0.25, 1.0340914067758025},
|
||||||
|
{1, 0, math.Inf(1)},
|
||||||
|
{1, 0.5, 0.69314718055994529},
|
||||||
|
{1, 1, 0},
|
||||||
|
{3, 0.75, 1.727299417860519},
|
||||||
|
{25, 0.4, 25.945791937289371},
|
||||||
|
{25, 0.7, 22.156653488661991},
|
||||||
|
{10, 0.5, 9.6687146147141299},
|
||||||
|
{100, 0.25, 106.5510925269767},
|
||||||
|
{1000, 0.01, 1075.0328320864389},
|
||||||
|
{1000, 0.99, 927.90815979664251},
|
||||||
|
{10000, 0.5, 9999.6666686420485},
|
||||||
|
} {
|
||||||
|
if got := GammaIncCompInv(test.a, test.x); math.Abs(got-test.want) > 1e-10 {
|
||||||
|
t.Errorf("test %d GammaIncCompInv(%g, %g) failed: got %g want %g", i, test.a, test.x, got, test.want)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
@@ -2,9 +2,11 @@
|
|||||||
// Use of this source code is governed by a BSD-style
|
// Use of this source code is governed by a BSD-style
|
||||||
// license that can be found in the LICENSE file.
|
// license that can be found in the LICENSE file.
|
||||||
|
|
||||||
// package cephes implements functions originally in the Netlib code by Stephen Mosher
|
// Package cephes implements functions originally in the Netlib code by Stephen Mosher.
|
||||||
package cephes
|
package cephes
|
||||||
|
|
||||||
|
import "math"
|
||||||
|
|
||||||
/*
|
/*
|
||||||
Additional copyright information:
|
Additional copyright information:
|
||||||
|
|
||||||
@@ -15,11 +17,13 @@ https://lists.debian.org/debian-legal/2004/12/msg00295.html
|
|||||||
*/
|
*/
|
||||||
|
|
||||||
var (
|
var (
|
||||||
badParamOutOfBounds = "cephes: parameter out of bounds"
|
badParamOutOfBounds = "cephes: parameter out of bounds"
|
||||||
|
badParamFunctionSingularity = "cephes: function singularity"
|
||||||
)
|
)
|
||||||
|
|
||||||
const (
|
const (
|
||||||
machEp = 1.11022302462515654042e-16 // 2^-53
|
machEp = 1.0 / (1 << 53)
|
||||||
maxLog = 7.09782712893383996732e2 // log(2^127)
|
maxLog = 1024 * math.Ln2
|
||||||
minLog = -7.451332191019412076235e2 // log(2^-128)
|
minLog = -1075 * math.Ln2
|
||||||
|
maxIter = 2000
|
||||||
)
|
)
|
||||||
|
|||||||
311
internal/cephes/igam.go
Normal file
311
internal/cephes/igam.go
Normal file
@@ -0,0 +1,311 @@
|
|||||||
|
// 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{
|
||||||
|
[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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{-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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{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},
|
||||||
|
[igamDimN]float64{-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(badParamOutOfBounds)
|
||||||
|
}
|
||||||
|
|
||||||
|
// 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(badParamOutOfBounds)
|
||||||
|
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
|
||||||
|
}
|
||||||
153
internal/cephes/igami.go
Normal file
153
internal/cephes/igami.go
Normal file
@@ -0,0 +1,153 @@
|
|||||||
|
// 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(badParamOutOfBounds)
|
||||||
|
}
|
||||||
|
|
||||||
|
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
|
||||||
|
}
|
||||||
153
internal/cephes/lanczos.go
Normal file
153
internal/cephes/lanczos.go
Normal file
@@ -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
|
||||||
|
}
|
||||||
@@ -1,15 +1,16 @@
|
|||||||
// Copyright ©2016 The gonum Authors. All rights reserved.
|
// 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
|
// Use of this source code is governed by a BSD-style
|
||||||
// license that can be found in the LICENSE file.
|
// 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.
|
||||||
* Cephes Math Library Release 2.1: December, 1988
|
|
||||||
* Copyright 1984, 1987, 1988 by Stephen L. Moshier
|
|
||||||
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
|
|
||||||
*/
|
|
||||||
|
|
||||||
package cephes
|
package cephes
|
||||||
|
|
||||||
|
import "math"
|
||||||
|
|
||||||
// polevl evaluates a polynomial of degree N
|
// polevl evaluates a polynomial of degree N
|
||||||
// y = c_0 + c_1 x_1 + c_2 x_2^2 ...
|
// 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
|
// where the coefficients are stored in reverse order, i.e. coef[0] = c_n and
|
||||||
@@ -31,3 +32,51 @@ func p1evl(x float64, coef []float64, n int) float64 {
|
|||||||
}
|
}
|
||||||
return ans
|
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
|
||||||
|
}
|
||||||
|
|||||||
170
internal/cephes/unity.go
Normal file
170
internal/cephes/unity.go
Normal file
@@ -0,0 +1,170 @@
|
|||||||
|
// 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)
|
||||||
|
}
|
||||||
110
internal/cephes/zeta.go
Normal file
110
internal/cephes/zeta.go
Normal file
@@ -0,0 +1,110 @@
|
|||||||
|
// 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(badParamOutOfBounds)
|
||||||
|
}
|
||||||
|
|
||||||
|
if q <= 0 {
|
||||||
|
if q == math.Floor(q) {
|
||||||
|
panic(badParamFunctionSingularity)
|
||||||
|
}
|
||||||
|
if x != math.Floor(x) {
|
||||||
|
panic(badParamOutOfBounds) // 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
|
||||||
|
}
|
||||||
178
roots.go
Normal file
178
roots.go
Normal file
@@ -0,0 +1,178 @@
|
|||||||
|
// 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 appropiate. 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
|
||||||
|
}
|
||||||
20
zeta.go
Normal file
20
zeta.go
Normal file
@@ -0,0 +1,20 @@
|
|||||||
|
// 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 "github.com/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)
|
||||||
|
}
|
||||||
42
zeta_test.go
Normal file
42
zeta_test.go
Normal file
@@ -0,0 +1,42 @@
|
|||||||
|
// 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"
|
||||||
|
"testing"
|
||||||
|
)
|
||||||
|
|
||||||
|
func TestZeta(t *testing.T) {
|
||||||
|
for i, test := range []struct {
|
||||||
|
x, q, want float64
|
||||||
|
}{
|
||||||
|
// Results computed using scipy.special.zeta
|
||||||
|
{1, 1, math.Inf(1)},
|
||||||
|
{1.00001, 0.5, 100001.96352290553},
|
||||||
|
{1.0001, 25, 9996.8017690244506},
|
||||||
|
{1.001, 1, 1000.5772884760117},
|
||||||
|
{1.01, 10, 97.773405639173305},
|
||||||
|
{1.5, 2, 1.6123753486854886},
|
||||||
|
{1.5, 20, 0.45287361712938717},
|
||||||
|
{2, -0.7, 14.28618087263834},
|
||||||
|
{2.5, 0.5, 6.2471106345688137},
|
||||||
|
{5, 2.5, 0.013073166646113805},
|
||||||
|
{7.5, 5, 7.9463377443314306e-06},
|
||||||
|
{10, -0.5, 2048.0174503557578},
|
||||||
|
{10, 0.5, 1024.0174503557578},
|
||||||
|
{10, 7.5, 2.5578265694201971e-9},
|
||||||
|
{12, 2.5, 1.7089167198843551e-5},
|
||||||
|
{17, 0.5, 131072.00101513157},
|
||||||
|
{20, -2.5, 2097152.0006014798},
|
||||||
|
{20, 0.75, 315.3368689825316},
|
||||||
|
{25, 0.25, 1125899906842624.0},
|
||||||
|
{30, 1, 1.0000000009313275},
|
||||||
|
} {
|
||||||
|
if got := Zeta(test.x, test.q); math.Abs(got-test.want) > 1e-10 {
|
||||||
|
t.Errorf("test %d Zeta(%g, %g) failed: got %g want %g", i, test.x, test.q, got, test.want)
|
||||||
|
}
|
||||||
|
}
|
||||||
|
}
|
||||||
Reference in New Issue
Block a user