Results 1 to 2 of 2

Thread: gamma function using Oxygen

  1. #1

    gamma function using Oxygen

    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()
    
    Last edited by jack; 25-12-2013 at 03:02.

  2. #2
    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
    Attached Images Attached Images

Similar Threads

  1. Gamma function
    By danbaron in forum Math: all about
    Replies: 14
    Last Post: 08-08-2011, 03:51
  2. COM and Oxygen
    By Charles Pegge in forum Programs
    Replies: 21
    Last Post: 08-08-2010, 09:56
  3. UI + TBGL + Oxygen II
    By Charles Pegge in forum Programs
    Replies: 1
    Last Post: 18-03-2010, 14:55
  4. my first try using TBGL and OXYGEN
    By D.J.Peters in forum Shout Box Area
    Replies: 9
    Last Post: 05-03-2010, 05:54
  5. Oxygen module
    By ErosOlmi in forum Legacy
    Replies: 4
    Last Post: 02-12-2008, 20:13

Members who have read this thread: 0

There are no members to list at the moment.

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts
  •