1. ## 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;
}
```

2. 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
```

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. 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.

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. 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

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 "
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
```

7. 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

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

8. 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
```

9. 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
>

10. 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.)

Page 1 of 2 12 Last

#### Posting Permissions

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