Procedure.d Gamma(x.d) ; AKJ  01-May-10
; Complete Gamma function for x>0 and x<2 (approx)
; Extended outside this range via: Gamma(x+1) = x*Gamma(x)
; Based on http://rosettacode.org/wiki/Gamma_function  [Ada]
Protected Dim A.d(28)
A(0) = 1.0
A(1) = 0.5772156649015328606
A(2) =-0.6558780715202538811
A(3) =-0.0420026350340952355
A(4) = 0.1665386113822914895
A(5) =-0.0421977345555443368 ; was ...33675
A(6) =-0.0096219715278769736
A(7) = 0.0072189432466630995
A(8) =-0.0011651675918590651
A(9) =-0.0002152416741149510
A(10) = 0.0001280502823881162
A(11) =-0.0000201348547807882
A(12) =-0.0000012504934821427
A(13) = 0.0000011330272319817
A(14) =-0.0000002056338416978
A(15) = 0.0000000061160951045
A(16) = 0.0000000050020076445
A(17) =-0.0000000011812745705
A(18) = 0.0000000001043426712
A(19) = 0.0000000000077822634
A(20) =-0.0000000000036968056
A(21) = 0.0000000000005100370
A(22) =-0.0000000000000205833
A(23) =-0.0000000000000053481
A(24) = 0.0000000000000012268
A(25) =-0.0000000000000001181
A(26) = 0.0000000000000000012
A(27) = 0.0000000000000000014
A(28) =-0.0000000000000000002
;A(29) = 0.00000000000000000002
Protected y.d, Prod.d, Sum.d, N
If x<=0.0: ProcedureReturn 0.0: EndIf ; Error
y = x-1.0: Prod = 1.0
While y>=1.0
  Prod*y: y-1.0 ; Recurse using Gamma(x+1) = x*Gamma(x)
Wend
Sum= A(28)
For N = 27 To 0 Step -1
  Sum*y+A(N)
Next N
If Sum=0.0: ProcedureReturn Infinity(): EndIf
ProcedureReturn Prod / Sum
EndProcedure

Procedure.d GammLn(x.d) ; AKJ  01-May-10
; Returns Ln(Gamma()) or 0 on error
; Based on Numerical Recipes gamma.h
Protected j, tmp.d, y.d, ser.d
Protected Dim cof.d(13)
cof(0) = 57.1562356658629235
cof(1) = -59.5979603554754912
cof(2) = 14.1360979747417471
cof(3) = -0.491913816097620199
cof(4) = 0.339946499848118887e-4
cof(5) = 0.465236289270485756e-4
cof(6) = -0.983744753048795646e-4
cof(7) = 0.158088703224912494e-3
cof(8) = -0.210264441724104883e-3
cof(9) = 0.217439618115212643e-3
cof(10) = -0.164318106536763890e-3
cof(11) = 0.844182239838527433e-4
cof(12) = -0.261908384015814087e-4
cof(13) = 0.368991826595316234e-5
If x<=0: ProcedureReturn 0: EndIf ; Bad argument
y = x
tmp = x+5.2421875
tmp = (x+0.5)*Log(tmp)-tmp
ser = 0.999999999999997092
For j=0 To 13
  y + 1: ser + cof(j)/y
Next j
ProcedureReturn tmp+Log(2.5066282746310005*ser/x)
EndProcedure

Procedure Factorial(x) ; AKJ  01-May-10
  ProcedureReturn Gamma(x+1)
EndProcedure