+- QB64 Phoenix Edition (https://qb64phoenix.com/forum)
+-- Forum: QB64 Rising (https://qb64phoenix.com/forum/forumdisplay.php?fid=1)
+--- Forum: Code and Stuff (https://qb64phoenix.com/forum/forumdisplay.php?fid=3)
+---- Forum: Programs (https://qb64phoenix.com/forum/forumdisplay.php?fid=7)
+---- Thread: a very accurate gamma function (/showthread.php?tid=2707)
a very accurate gamma function - Jack - 05-19-2024
this implementation is actually very accurate and well behaved, it's from https://www.researchgate.net/publication/302892889_Chebychev_interpolations_of_the_Gamma_and_Polygamma_Functions_and_their_analytical_properties
it took me some time to figure out how to use the tables because they are presented using the convolution operator which was new for me https://en.wikipedia.org/wiki/Convolution and https://youtu.be/KuXjwB4LzSA?t=346, I first converted the table to a power series and then used Maple to simplify the expression
Code: (Select All)
_Title "gamma function test"
Dim As Long i
For i = 1 To 20
Print i, gamma(i)
Next
Function gamma# (x_in As Double)
Static As Double d(20), f, sum, x, z
x = x_in
d(1) = 0.9999999999999999
d(2) = 0.08333333333338228
d(3) = 0.003472222216158396
d(4) = -0.002681326870868177
d(5) = -0.0002294790668608988
d(6) = 0.0007841331256329749
d(7) = 6.903992060449035E-05
d(8) = -0.0005907032612269776
d(9) = -2.877475047743023E-05
d(10) = 0.0005566293593820988
d(11) = 0.001799738210158344
d(12) = -0.008767670094590723
d(13) = 0.01817828637250317
d(14) = -0.02452583787937907
d(15) = 0.02361068245082701
d(16) = -0.01654210549755366
d(17) = 0.008304315532029655
d(18) = -0.00284326571576103
d(19) = 0.0005961678245858015
d(20) = -5.783378931872318E-05
z = 2.506628274631001 * x ^ (x - .5) * Exp(-x)
f = x * x: f = f * f: f = f * f * x: f = f * f * x 'f = x^19
sum=(d(20)*z+(d(19)*z+(d(18)*z+(d(17)*z+(d(16)*z+(d(15)*z+(d(14)*z _
+(d(13)*z+(d(12)*z+(d(11)*z+(d(10)*z+(d(9)*z+(d(8)*z+(d(7)*z+(d(6)* _
z+(d(5)*z+(d(4)*z+(d(3)*z+(x*z*d(1)+z*d(2))*x)*x)*x)*x)*x)*x)*x)*x) _
*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)/f
gamma = sum
End FunctionRE: a very accurate gamma function - Jack - 05-19-2024
and here's the gamma reciprocal
Code: (Select All)
_Title "gamma_inv function test"
Dim As Long i
For i = 1 To 20
Print i, gamma_inv(i)
Next
Function gamma_inv# (x_in As Double)
Static As Double d(20), f, sum, x, z
x = x_in
d(1) = 1
d(2) = -0.08333333333339432
d(3) = 0.003472222229927532
d(4) = 0.002681326782003515
d(5) = -0.0002294625761279813
d(6) = -0.0007841775716605252
d(7) = 7.093133582704529E-05
d(8) = 0.0005865361360167114
d(9) = -5.038235546703003E-05
d(10) = -0.0006597591049770531
d(11) = -0.001339314419184738
d(12) = 0.008102207604020538
d(13) = -0.01770764460311748
d(14) = 0.02454011794984818
d(15) = -0.02403014074759222
d(16) = 0.01704015119564821
d(17) = -0.008632563394507144
d(18) = 0.002976937308522899
d(19) = -0.000627846535400408
d(20) = 6.120299540340159E-05
z = 0.3989422804014327 * x ^ (.5 - x) * Exp(x)
f = x * x: f = f * f: f = f * f * x: f = f * f * x 'f = x^19
sum=(d(20)*z+(d(19)*z+(d(18)*z+(d(17)*z+(d(16)*z+(d(15)*z+(d(14)*z _
+(d(13)*z+(d(12)*z+(d(11)*z+(d(10)*z+(d(9)*z+(d(8)*z+(d(7)*z+(d(6)* _
z+(d(5)*z+(d(4)*z+(d(3)*z+(x*z*d(1)+z*d(2))*x)*x)*x)*x)*x)*x)*x)*x) _
*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)/f
gamma_inv = sum
End Function