''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' http://www.alglib.net/
'Cephes Math Library Release 2.8: June, 2000
'Copyright by Stephen L. Moshier
'
'Contributors:
' * Sergey Bochkanov (ALGLIB project). Translation from C to
' pseudocode.
'
'See subroutines comments for additional copyrights.
'
'>>> SOURCE LICENSE >>>
'This program is free software; you can redistribute it and/or modify
'it under the terms of the GNU General Public License as published by
'the Free Software Foundation (www.fsf.org); either version 2 of the
'License, or (at your option) any later version.
'
'This program is distributed in the hope that it will be useful,
'but WITHOUT ANY WARRANTY; without even the implied warranty of
'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
'GNU General Public License for more details.
'
'A copy of the GNU General Public License is available at
'http://www.fsf.org/licensing/licenses
'
'>>> END OF LICENSE >>>
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Gamma function
'
'Input parameters:
' X - argument
'
'Domain:
' 0 < X < 171.6
' -170 < X < 0, X is not an integer.
'
'Relative error:
' arithmetic domain # trials peak rms
' IEEE -170,-33 20000 2.3e-15 3.3e-16
' IEEE -33, 33 20000 9.4e-16 2.2e-16
' IEEE 33, 171.6 20000 2.3e-15 3.2e-16
'
'Cephes Math Library Release 2.8: June, 2000
'Original copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
'Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function Gamma(ByVal X As Double) As Double
Dim Result As Double
Dim p As Double
Dim PP As Double
Dim q As Double
Dim QQ As Double
Dim z As Double
Dim i As Long
Dim SgnGam As Double
SgnGam = 1#
If q > 33# Then
If X < 0# Then
i = Round(p)
If i Mod 2# = 0# Then
SgnGam = -1#
End If
z = q - p
If z > 0.5 Then
p = p + 1#
z = q - p
End If
z = Pi / (z * GammaStirF(q))
Else
z = GammaStirF(X)
End If
Result = SgnGam * z
Gamma = Result
Exit Function
End If
z = 1#
Do While X >= 3#
X = X - 1#
z = z * X
Loop
Do While X < 0#
If X > -0.000000001 Then
Result = z / ((1# + 0.577215664901533 * X) * X)
Gamma = Result
Exit Function
End If
z = z / X
X = X + 1#
Loop
Do While X < 2#
If X < 0.000000001 Then
Result = z / ((1# + 0.577215664901533 * X) * X)
Gamma = Result
Exit Function
End If
z = z / X
X = X + 1#
Loop
If X = 2# Then
Result = z
Gamma = Result
Exit Function
End If
X = X - 2#
PP = 1.60119522476752E-04
PP = 1.19135147006586E-03 + X * PP
PP = 1.04213797561762E-02 + X * PP
PP = 4.76367800457137E-02 + X * PP
PP = 0.207448227648436 + X * PP
PP = 0.494214826801497 + X * PP
PP = 1# + X * PP
QQ = -2.3158187332412E-05
QQ = 5.39605580493303E-04 + X * QQ
QQ = -4.45641913851797E-03 + X * QQ
QQ = 0.011813978522206 + X * QQ
QQ = 3.58236398605499E-02 + X * QQ
QQ = -0.234591795718243 + X * QQ
QQ = 7.14304917030273E-02 + X * QQ
QQ = 1# + X * QQ
Result = z * PP / QQ
Gamma = Result
Exit Function
Gamma = Result
End Function
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Natural logarithm of gamma function
'
'Input parameters:
' X - argument
'
'Result:
' logarithm of the absolute value of the Gamma(X).
'
'Output parameters:
' SgnGam - sign(Gamma(X))
'
'Domain:
' 0 < X < 2.55e305
' -2.55e305 < X < 0, X is not an integer.
'
'ACCURACY:
'arithmetic domain # trials peak rms
' IEEE 0, 3 28000 5.4e-16 1.1e-16
' IEEE 2.718, 2.556e305 40000 3.5e-16 8.3e-17
'The error criterion was relative when the function magnitude
'was greater than one but absolute when it was less than one.
'
'The following test used the relative error criterion, though
'at certain points the relative error could be much higher than
'indicated.
' IEEE -200, -4 10000 4.8e-16 1.3e-16
'
'Cephes Math Library Release 2.8: June, 2000
'Copyright 1984, 1987, 1989, 1992, 2000 by Stephen L. Moshier
'Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007).
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function LnGamma(ByVal X As Double, ByRef SgnGam As Double) As Double
Dim Result As Double
Dim A As Double
Dim B As Double
Dim C As Double
Dim p As Double
Dim q As Double
Dim u As Double
Dim w As Double
Dim z As Double
Dim i As Long
Dim LogPi As Double
Dim LS2PI As Double
Dim Tmp As Double
SgnGam = 1#
LogPi = 1.1447298858494
LS2PI = 0.918938533204673
If X < -34# Then
q = -X
w = LnGamma(q, Tmp)
i = Round(p)
If i Mod 2# = 0# Then
SgnGam = -1#
Else
SgnGam = 1#
End If
z = q - p
If z > 0.5 Then
p = p + 1#
z = p - q
End If
Result
= LogPi
- Log(z
) - w
LnGamma = Result
Exit Function
End If
If X < 13# Then
z = 1#
p = 0#
u = X
Do While u >= 3#
p = p - 1#
u = X + p
z = z * u
Loop
Do While u < 2#
z = z / u
p = p + 1#
u = X + p
Loop
If z < 0# Then
SgnGam = -1#
z = -z
Else
SgnGam = 1#
End If
If u = 2# Then
LnGamma = Result
Exit Function
End If
p = p - 2#
X = X + p
B = -1378.25152569121
B = -38801.6315134638 + X * B
B = -331612.992738871 + X * B
B = -1162370.97492762 + X * B
B = -1721737.0082084 + X * B
B = -853555.664245765 + X * B
C = 1#
C = -351.815701436523 + X * C
C = -17064.2106651881 + X * C
C = -220528.590553854 + X * C
C = -1139334.44367983 + X * C
C = -2532523.07177583 + X * C
C = -2018891.41433533 + X * C
p = X * B / C
LnGamma = Result
Exit Function
End If
q
= (X
- 0.5) * Log(X
) - X
+ LS2PI
If X > 100000000# Then
Result = q
LnGamma = Result
Exit Function
End If
p = 1# / (X * X)
If X >= 1000# Then
q = q + ((7.93650793650794 * 0.0001 * p - 2.77777777777778 * 0.001) * p + 8.33333333333333E-02) / X
Else
A = 8.11614167470508 * 0.0001
A = -(5.95061904284301 * 0.0001) + p * A
A = 7.93650340457717 * 0.0001 + p * A
A = -(2.777777777301 * 0.001) + p * A
A = 8.33333333333332 * 0.01 + p * A
q = q + A / X
End If
Result = q
LnGamma = Result
End Function
Private Function GammaStirF(ByVal X As Double) As Double
Dim Result As Double
Dim y As Double
Dim w As Double
Dim v As Double
Dim Stir As Double
w = 1# / X
Stir = 7.87311395793094E-04
Stir = -2.29549961613378E-04 + w * Stir
Stir = -2.68132617805781E-03 + w * Stir
Stir = 3.47222221605459E-03 + w * Stir
Stir = 8.33333333333482E-02 + w * Stir
w = 1# + w * Stir
If X > 143.01608 Then
v = X ^ (0.5 * X - 0.25)
y = v * (v / y)
Else
y = X ^ ((X - 0.5) / y)
End If
Result = 2.506628274631 * y * w
GammaStirF = Result
End Function
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' http://www.alglib.net/
'Cephes Math Library Release 2.8: June, 2000
'Copyright by Stephen L. Moshier
'
'Contributors:
' * Sergey Bochkanov (ALGLIB project). Translation from C to
' pseudocode.
'
'See subroutines comments for additional copyrights.
'
'>>> SOURCE LICENSE >>>
'This program is free software; you can redistribute it and/or modify
'it under the terms of the GNU General Public License as published by
'the Free Software Foundation (www.fsf.org); either version 2 of the
'License, or (at your option) any later version.
'
'This program is distributed in the hope that it will be useful,
'but WITHOUT ANY WARRANTY; without even the implied warranty of
'MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
'GNU General Public License for more details.
'
'A copy of the GNU General Public License is available at
'http://www.fsf.org/licensing/licenses
'
'>>> END OF LICENSE >>>
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Routines
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Lower Incomplete gamma integral
'
'The function is defined by
'
' x
' -
' 1 | | -t a-1
' igam(a,x) = ----- | e t dt.
' - | |
' | (a) -
' 0
'
'
'In this implementation both arguments must be positive.
'The integral is evaluated by either a power series or
'continued fraction expansion, depending on the relative
'values of a and x.
'
'ACCURACY:
'
' Relative error:
'arithmetic domain # trials peak rms
' IEEE 0,30 200000 3.6e-14 2.9e-15
' IEEE 0,100 300000 9.9e-14 1.5e-14
'
'Cephes Math Library Release 2.8: June, 2000
'Copyright 1985, 1987, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function IncompleteGamma(ByVal A As Double, ByVal X As Double) As Double
Dim Result As Double
Dim IGammaEpsilon As Double
Dim ans As Double
Dim ax As Double
Dim C As Double
Dim r As Double
Dim Tmp As Double
IGammaEpsilon = 0.000000000000001
If X <= 0# Or A <= 0# Then
Result = 0#
IncompleteGamma = Result
Exit Function
End If
If X > 1# And X > A Then
Result = 1# - IncompleteGammaC(A, X)
IncompleteGamma = Result
Exit Function
End If
ax
= A
* Log(X
) - X
- LnGamma
(A, Tmp
) If ax < -709.782712893384 Then
Result = 0#
IncompleteGamma = Result
Exit Function
End If
r = A
C = 1#
ans = 1#
Do
r = r + 1#
C = C * X / r
ans = ans + C
Loop Until C / ans <= IGammaEpsilon
Result = ans * ax / A
IncompleteGamma = Result
End Function
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Upper / Complemented incomplete gamma integral
'
'The function is defined by
'
'
' igamc(a,x) = 1 - igam(a,x)
'
' inf.
' -
' 1 | | -t a-1
' = ----- | e t dt.
' - | |
' | (a) -
' x
'
'
'In this implementation both arguments must be positive.
'The integral is evaluated by either a power series or
'continued fraction expansion, depending on the relative
'values of a and x.
'
'ACCURACY:
'
'Tested at random a, x.
' a x Relative error:
'arithmetic domain domain # trials peak rms
' IEEE 0.5,100 0,100 200000 1.9e-14 1.7e-15
' IEEE 0.01,0.5 0,100 200000 1.4e-13 1.6e-15
'
'Cephes Math Library Release 2.8: June, 2000
'Copyright 1985, 1987, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function IncompleteGammaC(ByVal A As Double, _
ByVal X As Double) As Double
Dim Result As Double
Dim IGammaEpsilon As Double
Dim IGammaBigNumber As Double
Dim IGammaBigNumberInv As Double
Dim ans As Double
Dim ax As Double
Dim C As Double
Dim yc As Double
Dim r As Double
Dim t As Double
Dim y As Double
Dim z As Double
Dim pk As Double
Dim pkm1 As Double
Dim pkm2 As Double
Dim qk As Double
Dim qkm1 As Double
Dim qkm2 As Double
Dim Tmp As Double
IGammaEpsilon = 0.000000000000001
IGammaBigNumber = 4.5035996273705E+15
IGammaBigNumberInv = 2.22044604925031 * 1E-16
If X <= 0# Or A <= 0# Then
Result = 1#
IncompleteGammaC = Result
Exit Function
End If
If X < 1# Or X < A Then
Result = 1# - IncompleteGamma(A, X)
IncompleteGammaC = Result
Exit Function
End If
ax
= A
* Log(X
) - X
- LnGamma
(A, Tmp
) If ax < -709.782712893384 Then
Result = 0#
IncompleteGammaC = Result
Exit Function
End If
y = 1# - A
z = X + y + 1#
C = 0#
pkm2 = 1#
qkm2 = X
pkm1 = X + 1#
qkm1 = z * X
ans = pkm1 / qkm1
Do
C = C + 1#
y = y + 1#
z = z + 2#
yc = y * C
pk = pkm1 * z - pkm2 * yc
qk = qkm1 * z - qkm2 * yc
If qk <> 0# Then
r = pk / qk
ans = r
Else
t = 1#
End If
pkm2 = pkm1
pkm1 = pk
qkm2 = qkm1
qkm1 = qk
If Abs(pk
) > IGammaBigNumber
Then pkm2 = pkm2 * IGammaBigNumberInv
pkm1 = pkm1 * IGammaBigNumberInv
qkm2 = qkm2 * IGammaBigNumberInv
qkm1 = qkm1 * IGammaBigNumberInv
End If
Loop Until t <= IGammaEpsilon
Result = ans * ax
IncompleteGammaC = Result
End Function
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'Inverse of complemented imcomplete gamma integral
'
'Given p, the function finds x such that
'
' igamc( a, x ) = p.
'
'Starting with the approximate value
'
' 3
' x = a t
'
' where
'
' t = 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.
'
'ACCURACY:
'
'Tested at random a, p in the intervals indicated.
'
' a p Relative error:
'arithmetic domain domain # trials peak rms
' IEEE 0.5,100 0,0.5 100000 1.0e-14 1.7e-15
' IEEE 0.01,0.5 0,0.5 100000 9.0e-14 3.4e-15
' IEEE 0.5,10000 0,0.5 20000 2.3e-13 3.8e-14
'
'Cephes Math Library Release 2.8: June, 2000
'Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
'
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Public Function InvIncompleteGammaC(ByVal A As Double, _
ByVal y0 As Double) As Double
Dim Result As Double
Dim IGammaEpsilon As Double
Dim IInvGammaBigNumber As Double
Dim x0 As Double
Dim x1 As Double
Dim X As Double
Dim yl As Double
Dim yh As Double
Dim y As Double
Dim d As Double
Dim lgm As Double
Dim dithresh As Double
Dim i As Long
Dim dir As Long
Dim Tmp As Double
IGammaEpsilon = 0.000000000000001
IInvGammaBigNumber = 4.5035996273705E+15
x0 = IInvGammaBigNumber
yl = 0#
x1 = 0#
yh = 1#
dithresh = 5# * IGammaEpsilon
d = 1# / (9# * A)
y = 1# - d - InvNormalDistribution(y0) * Sqr(d)
X = A * y * y * y
lgm = LnGamma(A, Tmp)
i = 0#
Do While i < 10#
If X > x0 Or X < x1 Then
d = 0.0625
Exit Do
End If
y = IncompleteGammaC(A, X)
If y < yl Or y > yh Then
d = 0.0625
Exit Do
End If
If y < y0 Then
x0 = X
yl = y
Else
x1 = X
yh = y
End If
d
= (A
- 1#
) * Log(X
) - X
- lgm
If d < -709.782712893384 Then
d = 0.0625
Exit Do
End If
d = (y - y0) / d
If Abs(d
/ X
) < IGammaEpsilon
Then Result = X
InvIncompleteGammaC = Result
Exit Function
End If
X = X - d
i = i + 1#
Loop
If x0 = IInvGammaBigNumber Then
If X <= 0# Then
X = 1#
End If
Do While x0 = IInvGammaBigNumber
X = (1# + d) * X
y = IncompleteGammaC(A, X)
If y < y0 Then
x0 = X
yl = y
Exit Do
End If
d = d + d
Loop
End If
d = 0.5
dir = 0#
i = 0#
Do While i < 400#
X = x1 + d * (x0 - x1)
y = IncompleteGammaC(A, X)
lgm = (x0 - x1) / (x1 + x0)
If Abs(lgm
) < dithresh
Then Exit Do
End If
lgm = (y - y0) / y0
If Abs(lgm
) < dithresh
Then Exit Do
End If
If X <= 0# Then
Exit Do
End If
If y >= y0 Then
x1 = X
yh = y
If dir < 0# Then
dir = 0#
d = 0.5
Else
If dir > 1# Then
d = 0.5 * d + 0.5
Else
d = (y0 - yl) / (yh - yl)
End If
End If
dir = dir + 1#
Else
x0 = X
yl = y
If dir > 0# Then
dir = 0#
d = 0.5
Else
If dir < -1# Then
d = 0.5 * d
Else
d = (y0 - yl) / (yh - yl)
End If
End If
dir = dir - 1#
End If
i = i + 1#
Loop
Result = X
InvIncompleteGammaC = Result
End Function
