gamma(z)
Func
© Γ(Z) by Stirling's formula ©M.Dave1
If real(floor(z))=z:Return when(z>0,(z-1)!,undef)
Local n,x,y
when(real(z)<0,1-z,z)→x
10-min(floor(abs(x)),10)→n
If n≠0:x+n→x
approx(ℯ^(⁻x)*x^(x-0.5)*√(2*Œ)*polyEval({0.00083949872087209,⁻5.1717909082606⁻005,⁻0.000
59216643735369,6.9728137583659⁻005,0.00078403922172007,⁻0.0002294720936214,⁻0.0026813271
604938,0.0034722222222222,0.083333333333333,1},1/x)*when(n=0,1,‹(1/(x-k),k,1,n)))→y
approx(when(real(z)<0,Œ/(sin(z*when(sin(Œ)=0,Œ,180)))/y,y))
EndFunc
This function works for all real and complex arguments. The accuracy is near full machine precision,
except for very large negative arguments.
Since the gamma function is similar to the factorial function, the result overflows the calculator floating
point range for relatively small arguments. For example, gamma(z) returns infinity for z>450. This
limitation can be overcome by using a function that returns the natural log of the gamma function,
instead of the gamma function itself. The log-gamma function is:
lngamma(z)
Func
© lnΓ(Z) by asymptotic series ©M.Dave1
If fPart(z)=0 and z<1:Return undef
Local n,x,y
when(real(z)<0,1-z,z)→x
10-min(floor(abs(x)),10)→n
If n≠0:x+n→x
approx((x-0.5)*ln(x)-x+0.5*ln(2*Œ)+polyEval({⁻0.0005952380952381,0.00079365079365079,⁻0.0
027777777777778,0.083333333333333},x^(⁻2))/x+when(n=0,0,ln(‹(1/(x-k),k,1,n))))→y
approx(when(real(z)<0,ln(Œ/(sin(z*when(sin(Œ)=0,Œ,180))))-y,y))
EndFunc
The program author, Martin Daveluy, has these additional comments:
These two series use asymptotic series combined with the recurrence formula ( gamma(Z+1) =
Z*gamma(Z) ) for Z<10 to keep full precision and the reflection formula ( gamma(Z)*gamma(1-Z)
= pi/(sin(pi*Z)) ) to extend domain to the entire complex plane. Note that the Gamma Stirling's
fomula is obtained by this LnGamma formula. The Stirling's coefficients are obtained by collecting
X power of the Maclaurin series for e^X ( 1+X+(X^2)/2!+... with X= LnGamma_asymptotic_series )
to reach higher precision.
With a gamma function, it is possible to write a factorial function for non-integer and complex
arguments:
factrl(xx)
func
©(x)factorial, complex & non-integer arguments
©1jan00/dburkett@infinet.com
when(imag(xx)=0 and fpart(xx)=0 and xx≤450,xx!,ℯ^(math\lngamma(xx+1)))
Endfunc
If the input argument is a real integer less than 450, the result is found with the built-in factorial function
is used, otherwise the log-gamma function is used. Note that lngamma() is installed in the math\ folder.
(credit to Martin Daveluy)
6 - 4
To Next Page
Loading...
Need help?
General
