yet another gamma function implementation

Post your FreeBASIC source, examples, tips and tricks here. Please don’t post code without including an explanation.
Post Reply
srvaldez
Posts: 3415
Joined: Sep 25, 2005 21:54

yet another gamma function implementation

Post by 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
 
Last edited by srvaldez on May 21, 2024 0:50, edited 1 time in total.
jdebord
Posts: 551
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact:

Re: yet another gamma function implementation

Post by 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.
Post Reply