a very accurate gamma function

[Jack]

this implementation is actually very accurate and well behaved, it's from https://www.researchgate.net/publication...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 Function

[Jack]

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