## yet another gamma function implementation

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

### yet another gamma function implementation

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: 550
Joined: May 27, 2005 6:20
Location: Limoges, France
Contact: