Gamma function

[jack]

For 100!, the function from, "Numerical Recipes in C", is only correct for 9 digits, so the Lanczos implementation beats it.

Dan, both are Lanczos approximations, the one I posted uses more terms, but there's something wrong with thinBasic,
it's supposed to use extended precision in it's calculations but

Val("1e+4932")

gives

1.7014118346046923E+38

[jack]

here's the Stieltjes' approximation.

uses "console"

Function gamma ( x As Ext ) As Ext
'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 Ext cf, logpi2half, xt, gam, pp 
    Dim As Integer i
    
    logpi2half = 0.91893853320467274178032973640562
    ' COEFFICIENTS FOR CONTINUED FRACTION
    Dim As Ext 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

Dim As Ext y, x = 1
While x>0
    Console_Write "x "
    x=Console_ReadLine()
    If x=0 Then Exit While
    y = gamma (x)
    Console_WriteLine "Gamma    "&LTrim$(Str$(x))&"    = "&Format$(y,17)
Wend
Console_WriteLine "All done. Press any key to finish"
Console_WaitKey

[danbaron]

x 101

Gamma 101 = 9.3326215443944156E+157
x 1001

Gamma 1001 = 4.023872600770938E+2567
x

Welcome to DrRacket, version 5.1 [3m].
Language: racket.
> (fact 100)
93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000
> (fact 1000)
402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
>

(Your script needs a, 'uses "console"'.)

Wow, that is really good, Johan.

You know a lot of math.

I just was looking at continued fractions, for the first time the other night.

http://en.wikipedia.org/wiki/Continued_fraction

Considering Stieltjes, all I know is that integrals are named after him.

http://en.wikipedia.org/wiki/Riemann–Stieltjes_integral

http://en.wikipedia.org/wiki/Lebesgue–Stieltjes_integration

I have a book, by Daniel Stroock, "A Concise Introduction to the Theory of Integration", 2nd ed., Birkhauser. But, it seems to be beyond my current capacity/motivation.

Dan

-----------------------------------------------------------------

I agree with you about Val("1e+4932").

But, try this.

Uses "console"

Dim e As Ext

e = Val("9.9999999999999999999999999999999999e+4931")

PrintL e

Console_ReadLine

[Petr Schreiber]

Hi Jack,

I can confirm that, interesting. It seems there is some kind of overflow(?) taking place.
VAL("1e4931") is still okay.

Anyway - if you need to specify 1e4932, you can do it like:

Dim myNumber As Ext
myNumber = 10 ^ 4932

... and the assigned value will be correct.


Petr

P.S. One wish - if you find some odd behavior, please do not hesitate to report the problem via Support in the top forum bar. It is easier to track and fix problems then. Thanks!

Dan, both are Lanczos approximations, the one I posted uses more terms, but there's something wrong with thinBasic,
it's supposed to use extended precision in it's calculations but

Val("1e+4932")

gives

[jack]

Dan, added the missing uses "console" guess I was in a hurry.

P.S. One wish - if you find some odd behavior, please do not hesitate to report the problem via Support in the top forum bar. It is easier to track and fix problems then. Thanks!

agree Petr.