Gamma function: Difference between revisions
Added solution for Free Pascal (Lazarus)
(Gamma function en BASIC256) |
(Added solution for Free Pascal (Lazarus)) |
||
Line 2,629:
===Stirling approximation===
<lang parigp>Stirling(x)=x--;sqrt(2*Pi*x)*(x/exp(1))^x</lang>
=={{header|Pascal}}==
A console application in Free Pascal, created with the Lazarus IDE.
Based on the algorithm for ln(Gamma(x)) (x > 0) in Press et al., Numerical Recipes, 3rd edition, pp. 256-7. For x >= 1/2, we simply take the exponential of their value; for x < 1/2 we calculate Gamma(1 - x) and use the reflection formula
Gamma(x)*Gamma(1 - x) = pi/sin(pi*x).
<lang pascal>
program GammaTest;
{$mode objfpc}{$H+}
uses SysUtils;
function Gamma( x : extended) : extended;
const COF : array [0..14] of extended =
( 0.999999999999997092, // may as well include this in the array
57.1562356658629235,
-59.5979603554754912,
14.1360979747417471,
-0.491913816097620199,
0.339946499848118887e-4,
0.465236289270485756e-4,
-0.983744753048795646e-4,
0.158088703224912494e-3,
-0.210264441724104883e-3,
0.217439618115212643e-3,
-0.164318106536763890e-3,
0.844182239838527433e-4,
-0.261908384015814087e-4,
0.368991826595316234e-5);
const
K = 2.5066282746310005;
PI_OVER_K = PI / K;
var
j : integer;
tmp, w, ser : extended;
reflect : boolean;
begin
reflect := (x < 0.5);
if reflect then w := 1.0 - x else w := x;
tmp := w + 5.2421875;
tmp := (w + 0.5)*Ln(tmp) - tmp;
ser := COF[0];
for j := 1 to 14 do ser := ser + COF[j]/(w + j);
try
if reflect then
result := PI_OVER_K * w * Exp(-tmp) / (Sin(PI*x) * ser)
else
result := K * Exp(tmp) * ser / w;
except
raise SysUtils.Exception.CreateFmt(
'Gamma(%g) is undefined or out of floating-point range', [x]);
end;
end;
// Main routine for testing the Gamma function
var
x, k : extended;
begin
WriteLn( 'Is it seamless at x = 1/2 ?');
x := 0.49999999999999;
WriteLn( SysUtils.Format( 'Gamma(%g) = %g', [x, Gamma(x)]));
x := 0.50000000000001;
WriteLn( SysUtils.Format( 'Gamma(%g) = %g', [x, Gamma(x)]));
WriteLn( 'Test a few values:');
WriteLn( SysUtils.Format( 'Gamma(1) = %g', [Gamma(1)]));
WriteLn( SysUtils.Format( 'Gamma(2) = %g', [Gamma(2)]));
WriteLn( SysUtils.Format( 'Gamma(3) = %g', [Gamma(3)]));
WriteLn( SysUtils.Format( 'Gamma(10) = %g', [Gamma(10)]));
WriteLn( SysUtils.Format( 'Gamma(101) = %g', [Gamma(101)]));
WriteLn( ' 100! = 9.33262154439442E157');
WriteLn( SysUtils.Format( 'Gamma(1/2) = %g', [Gamma(0.5)]));
WriteLn( SysUtils.Format( 'Sqrt(pi) = %g', [Sqrt(PI)]));
WriteLn( SysUtils.Format( 'Gamma(-7/2) = %g', [Gamma(-3.5)]));
(*
Note here a bug or misfeature in Lazarus (doesn't occur in Delphi):
Putting (16.0/105.0)*Sqrt(PI) does not give the required precision.
We have to explicitly define the integers as extended floating-point.
*)
k := extended(16.0)/extended(105.0);
WriteLn( SysUtils.Format( ' (16/105)Sqrt(pi) = %g', [k*Sqrt(PI)]));
WriteLn( SysUtils.Format( 'Gamma(-9/2) = %g', [Gamma(-4.5)]));
k := extended(32.0)/extended(945.0);
WriteLn( SysUtils.Format( '-(32/945)Sqrt(pi) = %g', [-k*Sqrt(PI)]));
end.
</lang>
{{out}}
<pre>
Is it seamless at x = 1/2 ?
Gamma(0.49999999999999) = 1.77245385090555
Gamma(0.50000000000001) = 1.77245385090548
Test a few values:
Gamma(1) = 1
Gamma(2) = 1
Gamma(3) = 2
Gamma(10) = 362880
Gamma(101) = 9.33262154439452E157
100! = 9.33262154439442E157
Gamma(1/2) = 1.77245385090552
Sqrt(pi) = 1.77245385090552
Gamma(-7/2) = 0.270088205852269
(16/105)Sqrt(pi) = 0.270088205852269
Gamma(-9/2) = -0.0600196013005043
-(32/945)Sqrt(pi) = -0.0600196013005042
</pre>
=={{header|Perl}}==
| |||