83 lines
2.3 KiB
Plaintext
83 lines
2.3 KiB
Plaintext
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
|