Page 1 of 2 12 LastLast
Results 1 to 10 of 15

Thread: Gamma function

  1. #1
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152

    Gamma function

    The factorial function (!) only works for integers, it is discrete.

    The Gamma function works for integers (n),

    gamma(n) = (n-1)!

    but also for values between the integers, it is continuous,

    gamma(n + 0.32457), etc.

    Here is a C program which has an approximation of the Gamma function, from,

    http://en.wikipedia.org/wiki/Gamma_function (the equation just above the heading, "History")

    If you run the program, you see that,

    gamma(2) = 1!
    gamma(2) < gamma(2.5) < gamma(3)
    gamma(3) = 2!
    gamma(3) < gamma(3.5) < gamma(4)
    gamma(4) = 3!
    gamma(4) < gamma(4.5) < gamma(5)
    gamma(5) = 4!
    etc.

    #include <stdio.h>
    #include <math.h>
    
    double gamma(double x)
    // rough approximation of Gamma function
    {
    static double pi;
    pi = 4 * atan(1);
    static double e;
    e = exp(1);
    int i;
    double sum = 0;
    static double c[5];
    c[0] = 1;
    c[1] = 1.0/12;
    c[2] = 1.0/288;
    c[3] = -139.0/51840;
    c[4] = -571.0/2488320;
    for(i=0;i<5;i++) sum += c[i]/pow(x, i);
    return sum *= pow(x,x-0.5) * pow(e,-x) * sqrt(2*pi);
    }
    
    int main()
    {
    char c;
    double i;
    for(i = 0.5; i<=40; i += 0.5) printf("%04.1f %060.10f\n", i, gamma(i));
    c = getchar();
    return 0;
    }
    
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  2. #2
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    I guess it's pretty stupid to write the gamma function, when C already includes it (tgamma).

    Remember that, for integers (n),

    gamma(n) = (n-1)!.

    (Below, you can see that as n grows, so does the error.)

    (Also notice that, the output simulates the appearance of rain falling from thunderstorm clouds, when viewed from a distance. (I don't know if that is a mathematical property.))

    ' code ------------------------------------------------------------------------------------------------
    
    #include <stdio.h>
    #include <math.h>
    
    int main()
    {
    char c;
    double i;
    for(i = 0.5; i<=40; i += 0.5) printf(".1f 0.10f\n", i, tgamma(i));
    c = getchar();
    return 0;
    }
    
    ' output ----------------------------------------------------------------------------------------------
    
    00.5 0000000000000000000000000000000000000000000000001.7724538509
    01.0 0000000000000000000000000000000000000000000000001.0000000000
    01.5 0000000000000000000000000000000000000000000000000.8862269255
    02.0 0000000000000000000000000000000000000000000000001.0000000000
    02.5 0000000000000000000000000000000000000000000000001.3293403882
    03.0 0000000000000000000000000000000000000000000000002.0000000000
    03.5 0000000000000000000000000000000000000000000000003.3233509704
    04.0 0000000000000000000000000000000000000000000000006.0000000000
    04.5 0000000000000000000000000000000000000000000000011.6317283966
    05.0 0000000000000000000000000000000000000000000000024.0000000000
    05.5 0000000000000000000000000000000000000000000000052.3427777846
    06.0 0000000000000000000000000000000000000000000000120.0000000000
    06.5 0000000000000000000000000000000000000000000000287.8852778150
    07.0 0000000000000000000000000000000000000000000000720.0000000000
    07.5 0000000000000000000000000000000000000000000001871.2543057978
    08.0 0000000000000000000000000000000000000000000005040.0000000000
    08.5 0000000000000000000000000000000000000000000014034.4072934834
    09.0 0000000000000000000000000000000000000000000040320.0000000000
    09.5 0000000000000000000000000000000000000000000119292.4619946090
    10.0 0000000000000000000000000000000000000000000362880.0000000001
    10.5 0000000000000000000000000000000000000000001133278.3889487856
    11.0 0000000000000000000000000000000000000000003628800.0000000009
    11.5 0000000000000000000000000000000000000000011899423.0839622486
    12.0 0000000000000000000000000000000000000000039916800.0000000075
    12.5 0000000000000000000000000000000000000000136843365.4655658574
    13.0 0000000000000000000000000000000000000000479001600.0000001077
    13.5 0000000000000000000000000000000000000001710542068.3195733000
    14.0 0000000000000000000000000000000000000006227020800.0000007451
    14.5 0000000000000000000000000000000000000023092317922.3142378032
    15.0 0000000000000000000000000000000000000087178291200.0000104308
    15.5 0000000000000000000000000000000000000334838609873.5564574599
    16.0 0000000000000000000000000000000000001307674368000.0003223500
    16.5 0000000000000000000000000000000000005189998453040.1263327700
    17.0 0000000000000000000000000000000000020922789888000.0051576600
    17.5 0000000000000000000000000000000000085634974475162.0572060300
    18.0 0000000000000000000000000000000000355687428096000.0585764600
    18.5 0000000000000000000000000000000001498612053315336.0709547900
    19.0 0000000000000000000000000000000006402373705728001.1475086200
    19.5 0000000000000000000000000000000027724322986333718.2905000000
    20.0 0000000000000000000000000000000121645100408832033.2812000000
    20.5 0000000000000000000000000000000540624298233507433.9061000000
    21.0 0000000000000000000000000000002432902008176640607.4166000000
    21.5 0000000000000000000000000000011082798113786906003.9520000000
    22.0 0000000000000000000000000000051090942171709448099.1363000000
    22.5 0000000000000000000000000000238280159446418474544.0000000000
    23.0 0000000000000000000000000001124000727777607971802.0000000000
    23.5 0000000000000000000000000005361303587544413749128.0000000000
    24.0 0000000000000000000000000025852016738884984515607.0000000000
    24.5 0000000000000000000000000125990634307293761521577.0000000000
    25.0 0000000000000000000000000620448401733239552410000.0000000000
    25.5 0000000000000000000000003086770540528697529220000.0000000000
    26.0 0000000000000000000000015511210043330988264640000.0000000000
    26.5 0000000000000000000000078712648783481796272090000.0000000000
    27.0 0000000000000000000000403291461126605793833730000.0000000000
    27.5 0000000000000000000002085885192762267589569090000.0000000000
    28.0 0000000000000000000010888869450418354972400000000.0000000000
    28.5 0000000000000000000057361842800962331239100000000.0000000000
    29.0 0000000000000000000304888344611713895574200000000.0000000000
    29.5 0000000000000000001634812519827426644042100000000.0000000000
    30.0 0000000000000000008841761993739703670144000000000.0000000000
    30.5 0000000000000000048226969334909066557884200000000.0000000000
    31.0 0000000000000000265252859812191200035000000000000.0000000000
    31.5 0000000000000001470922564714727341197000000000000.0000000000
    32.0 0000000000000008222838654177925782278000000000000.0000000000
    32.5 0000000000000046334060788513915613293000000000000.0000000000
    33.0 0000000000000263130836933693625032901000000000000.0000000000
    33.5 0000000000001505856975626701932920000000000000000.0000000000
    34.0 0000000000008683317618811891588840000000000000000.0000000000
    34.5 0000000000050446208683494507567950000000000000000.0000000000
    35.0 0000000000295232799039604142308230000000000000000.0000000000
    35.5 0000000001740394199580549821257590000000000000000.0000000000
    36.0 0000000010333147966386148254900000000000000000000.0000000000
    36.5 0000000061783994085110651212700000000000000000000.0000000000
    37.0 0000000371993326789906364865600000000000000000000.0000000000
    37.5 0000002255115784106522798538200000000000000000000.0000000000
    38.0 0000013763753091226549819111800000000000000000000.0000000000
    38.5 0000084566841903994709253311100000000000000000000.0000000000
    39.0 0000523022617466599549516000000000000000000000000.0000000000
    39.5 0003255823413303802954033000000000000000000000000.0000000000
    40.0 0020397882081197472289204000000000000000000000000.0000000000
    
    Last edited by danbaron; 01-08-2011 at 06:32.
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  3. #3
    actually it's interesting to find new ways to approximate the gamma function, when I get around to it I will post some implementations.

  4. #4
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    http://en.wikipedia.org/wiki/Gamma_function

    http://en.wikipedia.org/wiki/Euler%E2%80%93Mascheroni_constant

    I can say this.

    Before stumbling onto the realization that C now includes the Gamma function, I tried to implement it using the infinite product definition of Euler and Weierstrass, using the Euler-Mascheroni constant.

    In C, I had real problems getting the constant to converge, and with the accuracy of the partial product.

    My experience really made me wonder how the intrinsic C Gamma function is so accurate.

    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  5. #5
    the Lanczos approximation is a popular method but there are other interesting formulas, my idea is to round up the newer methods and see how they perform.

  6. #6
    here'a Lanczos implementation
    Uses "console"
    Function gamma(ByVal y As Ext) As Ext
        Dim As Ext Pi    = 3.1415926535897932385
        Dim As Ext sq2pi = 2.50662827463100050241577 'sqrt(2Pi)
        Dim As Ext g     = 607/128 ' best resu'ts when 4<=g<=5
        Dim As Ext t, w, gam, x = y-1
        Dim As Integer i, cg = 14
        Dim c(15) As Ext
        
        '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  'thiBasic arrays start at 1 ?
          c( 2) =       57.156235665862923517
          c( 3) =      -59.597960355475491248
          c( 4) =       14.136097974741747174
          c( 5) =       -0.49191381609762019978
          c( 6) =        0.33994649984811888699/10000 
          c( 7) =        0.46523628927048575665/10000
          c( 8) =       -0.98374475304879564677/10000
          c( 9) =        0.15808870322491248884/1000
          c(10) =       -0.21026444172410488319/1000
          c(11) =        0.21743961811521264320/1000
          c(12) =       -0.16431810653676389022/1000
          c(13) =        0.84418223983852743293/10000
          c(14) =       -0.26190838401581408670/10000
          c(15) =        0.36899182659531622704/100000
          
        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
    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; 07-08-2011 at 04:46.

  7. #7
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    Wow, it works good, Johan.

    I didn't know anything about the Lanczos implementation.

    I had heard of him before.

    http://www.gap-system.org/~history/Biographies/Lanczos.html

    100! = 93326215443944152681699238856266700490715968264381621468592963895217599993229915608941463976156518286253697920827223758251185210916864000000000000000000000000

    From your implementation,

    Gamma 101 = 9.332621544394413E+157.

    So, the first 15 digits are correct.

    (I'll look at the implementation of the Gamma function, in "Numerical Recipes in C", now.)

    Dan



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

  8. #8
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    For 100!, the function from, "Numerical Recipes in C", is only correct for 9 digits, so the Lanczos implementation beats it.

    ' code --------------------------------------------------------------------------------------
    
    #include <stdio.h>
    #include <math.h>
    
    double gammaln(double xx)
    // From "Numerical Recipes in C", 2nd ed., p.214.
    // Returns ln(Gamma(xx)), for xx > 0.
    {
    int j;
    double x,y,tmp,ser;
    static double cof[6];
    cof[0] =  76.18009172947146;
    cof[1] = -86.50532032941677;
    cof[2] =  24.01409824083091;
    cof[3] =  -1.231739572450155;
    cof[4] =   0.1208650973866179e-2;
    cof[5] =  -0.5395239384953e-5;
    
    y=x=xx;
    tmp=x+5.5;
    tmp -= (x+0.5)*log(tmp);
    ser=1.000000000190015;
    for(j=0;j<=5;j++) ser += cof[j]/++y;
    return -tmp+log(2.5066282746310005*ser/x);
    }
    
    int main()
    {
    char c;
    double i;
    i = 101;
    printf("%05.1f %25.20e\n", i, exp(gammaln(i)));
    c = getchar();
    return 0;
    }
    
    ' output ------------------------------------------------------------------------------------
    
    101.0 9.33262154537299415097e+157
    
    "You can't cheat an honest man. Never give a sucker an even break, or smarten up a chump." - W.C.Fields

  9. #9
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    On the other hand, "Super-Mathematica" will give the exact result for Gamma(1001) (maybe it calculates 1000!).

    http://www.wolframalpha.com/input/?i=Gamma[1001]

    Result:
    402387260077093773543702433923003985719374864210714632543799910429938512398629020592044208486969404800479988610197196058631666872994808558901323829669944590997424504087073759918823627727188732519779505950995276120874975462497043601418278094646496291056393887437886487337119181045825783647849977012476632889835955735432513185323958463075557409114262417474349347553428646576611667797396668820291207379143853719588249808126867838374559731746136085379534524221586593201928090878297308431392844403281231558611036976801357304216168747609675871348312025478589320767169132448426236131412508780208000261683151027341827977704784635868170164365024153691398281264810213092761244896359928705114964975419909342221566832572080821333186116811553615836546984046708975602900950537616475847728421889679646244945160765353408198901385442487984959953319101723355556602139450399736280750137837615307127761926849034352625200015888535147331611702103968175921510907788019393178114194545257223865541461062892187960223838971476088506276862967146674697562911234082439208160153780889893964518263243671616762179168909779911903754031274622289988005195444414282012187361745992642956581746628302955570299024324153181617210465832036786906117260158783520751516284225540265170483304226143974286933061690897968482590125458327168226458066526769958652682272807075781391858178889652208164348344825993266043367660176999612831860788386150279465955131156552036093988180612138558600301435694527224206344631797460594682573103790084024432438465657245014402821885252470935190620929023136493273497565513958720559654228749774011413346962715422845862377387538230483865688976461927383814900140767310446640259899490222221765904339901886018566526485061799702356193897017860040811889729918311021171229845901641921068884387121855646124960798722908519296819372388642614839657382291123125024186649353143970137428531926649875337218940694281434118520158014123344828015051399694290153483077644569099073152433278288269864602789864321139083506217095002597389863554277196742822248757586765752344220207573630569498825087968928162753848863396909959826280956121450994871701244516461260379029309120889086942028510640182154399457156805941872748998094254742173582401063677404595741785160829230135358081840096996372524230560855903700624271243416909004153690105933983835777939410970027753472000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

    Check result, from Racket.

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

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

  10. #10
    thinBasic MVPs danbaron's Avatar
    Join Date
    Jan 2010
    Location
    California
    Posts
    1,378
    Rep Power
    152
    For Gamma(1001 + 1/1000000),

    http://www.wolframalpha.com/input/?i=Gamma[1001.000001]

    this is all I could get from "Super-Mathematica".

    4.0239003988057584997902173976163592011969374179001511239151002642392471030979284415068711496796131858912114129894951019307516229727006305338006965443711395900922510707344833143982161291576196145549524252116586498402388174031102033111615654658920197684590636911143385941614831004326546290077112736324575126945175550119098139642713208122214518177314109490782649292284077778713228049034887764970619206691682137254901666889850542038221053151508464693805173317329282324843144943943599436514279514538607552755619053756475634096657267703608254014755211290420692809243460322686928748132961796450943572995373926927227559737480726681980983218643224702815506486387652208117712019519946798982740183845603220534925666831744612502114908382710346769016257450862275257405860597569351553374047531992774021108553153132441252879047714632798432033625895109963398433028869133545373648578084478742621362263016745157981846648908850249353140681346060239970206896405377407332686209471584545211937434617377658459807959756083415727145958537520435373456645804904729989329964598572790374526220294313751473179336564051622005259346833431900944629570839237229481513254432209634074253387839700361771632471703417078441912584157269029731250759827936047129505797040335623306987799596245979371220717534340058326409719830787501939605783323783893179283319665132503930707413454484359916812366340570334427980972071481375201194203796140336701067494124741436407947480809292277983651740542286852851769396310225165686126061205580214600926014228591882855121520450907878434749341566595088164908352656635137221802845387984029365296643760097823289428298591031940194745031433675760122709191816281215841944892670582736814231986110261890538652395205193148570945253573720796741616911796790452097412733314334762831996654388583958842108704860923901748565161209141596165235353480433962920389547091119284361184464231936462484553514420363877766052084394699696456390950010286890909749030319272522477519054298031886134749961363250122571370745518417715230196582522727560938697643549511273538629747079520406413534104121981790847012314953152367591394890854899650620559354790016382878566029849228881525569110021829561650569024779437891761222842924279083000460739467674535783721585102675448762616885235579170233920356944428551215896798274345929092227557674918376075932337861270314616782583093758957562631739284704672146519207109744977914313663591015850751595443459206763794303939770317473799857171590222662301548722568155715675746100133944439923277239113204532536037954523472602398569689151629072755297077462874071783240172319640630147138674416129455126541858280762275787308341634326062407520397616731640922794849171635359556882358169551264574095302587695705568308648876430372219809533500063843080428965101561654802023673373158391310104669887045308869837699763815934065102121965101508796092056855661749281613039002338281146882969916854109020038123529481436584940630288994678359267874045978361539072677171117303849094719953611327095943318932698734386917593119024538031854596120942550157270174457336010032617790660722732015345076544097871929876417549941543216004206034777057170912796699287536720531688865073593195894207758319489176583240052580932558196302682862695669489650826090073859613164147950379619784492273985221444482402596660180723828594850186630616315085145571389650601927654671155352412964610239592652562298130511431240532574515358034040225214190832152608962048874360481181163117706057960456785858703174930508635662497576032855328818899495348962384480754728207869918433136816825060953289916424884475186026464979296732458581572805992968768306412122770943206407590057950197962062011156536819841696174893179842028679782235569275374134554732020123806556821231345176... × 10^2567

    (And, my experience with "Super-Mathematica" is that it usually gives the answer for almost any input, NOW.)

    (The answers for Gamma(1001) and Gamma(1001.000001) diverge after only 4 digits. It seems strange. I guess it must be correct.)

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

Page 1 of 2 12 LastLast

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
  •