'********************************************* '* Program to demonstrate the Gamma Function * '* * '* Basic Version By J-P Moreau * '* (www.jpmoreau.fr) * '* ----------------------------------------- * '* Reference: * '* "Numerical Recipes, By W.H. Press, B.P. * '* Flannery, S.A. Teukolsky and T. Vetter- * '* ling, Cambridge University Press, 1986" * '* [BIBLI 08]. * '* ----------------------------------------- * '* SAMPLE RUN: * '* * '* X Gamma(X) * '* ------------------------- * '* 0.5000 1.772453851 * '* 1.0000 1.000000000 * '* 1.5000 0.886226925 * '* 2.0000 1.000000000 * '* 2.5000 1.329340388 * '* 3.0000 2.000000000 * '* 3.5000 3.323350970 * '* 4.0000 5.999999999 * '* 4.5000 11.631728395 * '* 5.0000 23.999999996 * '* * '********************************************* defint i-n defdbl a-h,o-z cls print print " X Gamma(X) " print " -------------------------" x=0# for i=1 to 10 x=x+0.5# 'call Gamma function gosub 1000 print using " ##.#### ####.########"; x; y next i print end '******************************************* '* FUNCTION GAMMA(X) * '* --------------------------------------- * '* Returns the value of Gamma(x) in double * '* precision as EXP(LN(GAMMA(X))) for X>0. * '* (input: x output: y) * '******************************************* 1000 'Function Gamma(x) dim cof(6) cof(1)=76.18009173# : cof(2)=-86.50532033# cof(3)=24.01409822# : cof(4)=-1.231739516# cof(5)=0.00120858003# : cof(6)=-0.00000536382# stp=2.50662827465# half=0.5# : one=1# : fpf=5.5# xx=x-one tmp=xx+fpf tmp=(xx+half)*LOG(tmp)-tmp ser=one for j=1 to 6 xx=xx+one ser=ser+cof(j)/xx next j y = EXP(tmp+LOG(stp*ser)) return 'end of file gamma.bas