gamma function using Oxygen

[jack]

here are two implementations of the gamma function, the first uses the Lanczos approximation and it’s a bit slower than the second which uses the Stieltjes approximation, but the first function also handles negative arguments.
please point out any mistakes or simplifications that could be done.



Uses "Console" , "oxygen"
Dim i , j, t As Long
Dim scr As String
Dim As Long pGamma
Dim As Long pFinish
Dim x, y As Double
scr = "Function gamma(ByVal y As Double) As Double link #pGamma
Dim As double Pi = 3.1415926535897932385
Dim As double sq2pi = 2.50662827463100050241577 'sqrt(2Pi)
Dim As double g = 607/128 ' best resu'ts when 4<=g<=5
Dim As double t, w, gam, x = y-1
Dim As long i, cg = 14
Dim c(15) As double

'Lanczos approximation for the complex plane
'calculated using vpa digits(256)
'the best set of coeffs was selected from
'a solution space of g=0 to 32 with 1 to 32 terms
'these coeffs really give superb performance
'of 15 sig. digits for 0<=real(z)<=171
'coeffs should sum to about g*g/2+23/24

'http://www.numericana.com/answer/info/godfrey.htm

c( 1) = 0.99999999999999709182
c( 2) = 57.156235665862923517
c( 3) = -59.597960355475491248
c( 4) = 14.136097974741747174
c( 5) = -0.49191381609762019978
c( 6) = 0.000033994649984811888699
c( 7) = 0.000046523628927048575665
c( 8) = -0.000098374475304879564677
c( 9) = 0.00015808870322491248884
c(10) = -0.00021026444172410488319
c(11) = 0.00021743961811521264320
c(12) = -0.00016431810653676389022
c(13) = 0.000084418223983852743293
c(14) = -0.000026190838401581408670
c(15) = 0.0000036899182659531622704

If ( x < 0 ) Then
x=-x
If Frac(x)=0 Then
gam=10^4932
Else
t = c(1)
For i=1 To cg
t = t + c(i+1)/(x+i)
Next
w = x + g + 0.5
gam=w^(x+0.5)*Exp(-w)*sq2pi*t
gam = Pi*x/(gam*Sin(Pi*x))
End If
Else
t = c(1)
For i=1 To cg
t = t + c(i+1)/(x+i)
Next
w = x + g + 0.5
gam=w^(x+0.5)*Exp(-w)*sq2pi*t
End If
Function = gam
End Function
sub finish() link #pFinish
terminate
end sub
"
O2_Asm scr
If O2_Error Then
MsgBox 0,O2_Error
Stop
Else
O2_Exec
End If
Declare Function Gamma(ByVal Double ) As Double At pGamma
Declare Sub Finish() At pFinish
t=GetTickCount
For i=1 To 10000
For j=1 To 170
x=gamma(j)
'Console_WriteLine x
Next
Next
t=GetTickCount-t

Console_WriteLine t
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
Finish()


and the second.




Uses "Console" , "oxygen"
Dim i , j, t As Long
Dim scr As String
Dim As Long pGamma
Dim As Long pFinish
Dim x, y As Double
scr = "Function gamma (byval x As Double ) As Double link #pGamma
'Stieltjes' approximation to the Gamma function
'x!=exp(log(2*Pi)/2-x+(x+1/2)*log(x)+cf(x))
'cf(x) = c1
' ---------------------
' c2
' x + ----------------
' c3
' x + ------------
' c4
' x + --------
' x + '
' '
' '
Dim As Double cf, logpi2half, xt, gam, pp
Dim As long i

logpi2half = 0.91893853320467274178032973640562
' COEFFICIENTS FOR CONTINUED FRACTION
Dim As Double C(14)
c( 1) = 0.0833333333333333333333333333333333
c( 2) = 0.0333333333333333333333333333333333
c( 3) = 0.252380952380952380952380952380952
c( 4) = 0.525606469002695417789757412398922
c( 5) = 1.01152306812684171174737212473062
c( 6) = 1.51747364915328739842849151949548
c( 7) = 2.26948897420495996090915067220988
c( 8) = 3.00991738325939817007314073420772
c( 9) = 4.02688719234390122616887595318144
c(10) = 5.00276808075403005168850241227619
c(11) = 6.28391137081578218007266315495165
c(12) = 7.49591912238403392975235470826751
c(13) = 9.04066023436772669953113936026048
c(14) = 10.4893036545094822771883713045926
' CONTINUED FRACTION
xt=x+8
CF = C(9)+Xt
For i = 8 To 2 Step -1
cf=xt+c(i)/cf
Next
cf=c(1)/cf

' POCHHAMMER POLYNOMIAL
PP = X * (X+1) * (X+2) * (X+3) * (X+4) * (X+5) * (X+6) * (X+7) * (X+8)
gam = Exp(logpi2half-xt+(xt+.5)*Log(xt)+cf)
Return gam/pp
End Function
Sub finish() link #pFinish
terminate
End Sub
"
O2_Asm scr
If O2_Error Then
MsgBox 0,O2_Error
Stop
Else
O2_Exec
End If
Declare Function Gamma(ByVal Double ) As Double At pGamma
Declare Sub Finish() At pFinish
t=GetTickCount
For i=1 To 10000
For j=1 To 170
x=gamma(j)
'Console_WriteLine x
Next
Next
t=GetTickCount-t

Console_WriteLine t
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey
Finish()

[RobbeK]

Thanks Jack !!

I'll study them and report back.

My idea is to use the Striling on complex numbers -- (ok, it gives poor results for low numbers -- but this can be tuned by expanding the series if needed , it also holds for complex numbers if R(z) > 0 (perfect for the interesting parts of the complex Zeta).
Strangely (? maybe), the old formula got a boost last decade and interesting methods have been found.

I attached something, but the name is wrong -- it has been found after the Windschitl formula (more recent) -- should be lightning fast and this one should be correct upto 8 digits.
Worth trying and suitable for computer graphics.
For the Zeta there's one integration left then ... (and a pole at z=1 -- but that's ok , the function changes form divergent to convergent there ).

best , thanks for your help .. Rob