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