'*********************************************
'* 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