yet another gamma function implementation

[srvaldez]

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

function gamma(byval x as double) as double
    static as double a(20), f, sum, z

    a( 1) =  0.9999999999999999
    a( 2) =  0.08333333333338228
    a( 3) =  0.003472222216158396
    a( 4) = -0.002681326870868177
    a( 5) = -0.0002294790668608988
    a( 6) =  0.0007841331256329749
    a( 7) =  6.903992060449035e-05
    a( 8) = -0.0005907032612269776
    a( 9) = -2.877475047743023e-05
    a(10) =  0.0005566293593820988
    a(11) =  0.001799738210158344
    a(12) = -0.008767670094590723
    a(13) =  0.01817828637250317
    a(14) = -0.02452583787937907
    a(15) =  0.02361068245082701
    a(16) = -0.01654210549755366
    a(17) =  0.008304315532029655
    a(18) = -0.00284326571576103
    a(19) =  0.0005961678245858015
    a(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=z*(a(20)+(a(19)+(a(18)+(a(17)+(a(16)+(a(15)+(a(14)+(a(13)+ _
    (a(12)+(a(11)+(a(10)+(a(9)+(a(8)+(a(7)+(a(6)+(a(5)+(a(4)+(a(3)+ _
    (x*a(1)+a(2))*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)/f

    return sum
end function

dim as long i

for i=1 to 20
    print i, gamma(i)
next
Sleep

output

Code: Select all

 1             1
 2             1
 3             2
 4             6
 5             24
 6             120
 7             720
 8             5040
 9             40320
 10            362880
 11            3628800
 12            39916800
 13            479001600
 14            6227020800
 15            87178291200
 16            1307674368000
 17            20922789888000
 18            355687428096000
 19            6.402373705728e+015
 20            1.21645100408832e+017
 

lngamma

Code: Select all

Function lngamma(Byval x As Double) As Double
    Static As Double c(20), f, sum

    c( 1) =  0.9189385332046727
    c( 2) =  0.08333333333338824
    c( 3) = -6.875523375600131e-12
    c( 4) = -0.002777777444245219
    c( 5) = -8.229108205635117e-09
    c( 6) =  0.0007937666109114311
    c( 7) = -9.417521476898452e-07
    c( 8) = -0.0005917271548255839
    c( 9) =  1.10417689045145e-05
    c(10) =  0.0006092585812576799
    c(11) =  0.001574019576696141
    c(12) = -0.008450831265113272
    c(13) =  0.01796959267198273
    c(14) = -0.02456297768317419
    c(15) =  0.0238450482919148
    c(16) = -0.01680609543961501
    c(17) =  0.008475041030242971
    c(18) = -0.002912109145921757
    c(19) =  0.000612384508536461
    c(20) = -5.955145748126546e-05

    f=x*x : f=f*f : f=f*f*x : f=f*f*x 'f = x^19
    sum=(c(20)+(c(19)+(c(18)+(c(17)+(c(16)+(c(15)+(c(14)+(c(13)+(c(12) _
    +(c(11)+(c(10)+(c(9)+(c(8)+(c(7)+(c(6)+(c(5)+(c(4)+(c(3)+(c(2) _
    +(c(1)-.5*Log(x)+(Log(x)-1)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)*x)/f

    Return sum
End Function

[jdebord]

Thank you for the reference. I have downloaded it.

For special functions I use mainly the MPFR library (thanks to your wonderful FreeBASIC binding), but it is always interesting to have new approximations.