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
```

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
```