Page 2 of 2 FirstFirst 12
Results 11 to 15 of 15

Thread: Gamma function

  1. #11
    Quote Originally Posted by danbaron View Post
    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

  2. #12
    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
    
    Last edited by jack; 08-08-2011 at 02:07.

  3. #13
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    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
    

    Last edited by danbaron; 07-08-2011 at 11:20.
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  4. #14
    Super Moderator Petr Schreiber's Avatar
    Join Date
    Aug 2005
    Location
    Brno - Czech Republic
    Posts
    7,128
    Rep Power
    732
    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!

    Quote Originally Posted by jack View Post
    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
    Learn 3D graphics with ThinBASIC, learn TBGL!
    Windows 10 64bit - Intel Core i5-3350P @ 3.1GHz - 16 GB RAM - NVIDIA GeForce GTX 1050 Ti 4GB

  5. #15
    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.

Page 2 of 2 FirstFirst 12

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
  •