?? fspec.pas
字號:
{ **********************************************************************
* Unit FSPEC.PAS *
* Version 1.2d *
* (c) J. Debord, May 2003 *
**********************************************************************
Special functions and probability distributions
**********************************************************************
Error handling: The global variable MathError returns the error code
from the last function evaluation. It must be checked immediately
after a function call:
Y := f(X); (* f is one of the functions of the library *)
if MathError = FN_OK then ...
The possible error codes are the following:
---------------------------------------------
Error code Value Significance
---------------------------------------------
FN_OK 0 No error
FN_DOMAIN -1 Argument domain error
FN_SING -2 Function singularity
FN_OVERFLOW -3 Overflow range error
FN_UNDERFLOW -4 Underflow range error
FN_TLOSS -5 Total loss of precision
FN_PLOSS -6 Partial loss of precision
---------------------------------------------
**********************************************************************
Most functions are translated C code from Cephes math library
by S. Moshier (http://www.moshier.net)
********************************************************************** }
unit fspec;
interface
uses
fmath;
{ ----------------------------------------------------------------------
Table of factorials
---------------------------------------------------------------------- }
const
NFACT = 33;
var
FactArray : array[0..NFACT] of Float =
(1.0,
1.0,
2.0,
6.0,
24.0,
120.0,
720.0,
5040.0,
40320.0,
362880.0,
3628800.0,
39916800.0,
479001600.0,
6227020800.0,
87178291200.0,
1307674368000.0,
20922789888000.0,
355687428096000.0,
6402373705728000.0,
121645100408832000.0,
2432902008176640000.0,
51090942171709440000.0,
1124000727777607680000.0,
25852016738884976640000.0,
620448401733239439360000.0,
15511210043330985984000000.0,
403291461126605635584000000.0,
10888869450418352160768000000.0,
304888344611713860501504000000.0,
8841761993739701954543616000000.0,
265252859812191058636308480000000.0,
8222838654177922817725562880000000.0,
263130836933693530167218012160000000.0,
8683317618811886495518194401280000000.0);
{ ----------------------------------------------------------------------
Special functions
---------------------------------------------------------------------- }
function Fact(N : Integer) : Float; { Factorial }
function Binomial(N, K : Integer) : Float; { Binomial coef. C(N,K) }
function Gamma(X : Float) : Float; { Gamma function }
function SgnGamma(X : Float) : Integer; { Sign of Gamma function }
function IGamma(A, X : Float) : Float; { Incomplete Gamma function }
function JGamma(A, X : Float) : Float; { Complement of IGamma }
function Beta(X, Y : Float) : Float; { Beta function }
function IBeta(A, B, X : Float) : Float; { Incomplete Beta function }
function Erf(X : Float) : Float; { Error function }
function Erfc(X : Float) : Float; { Complement of Erf }
function LnGamma(X : Float) : Float; overload; { Log(|Gamma(X)|) }
function LnGamma(Z : Complex) : Complex; overload;
{ ----------------------------------------------------------------------
Binomial distribution with probability P and number of repetitions N
---------------------------------------------------------------------- }
function PBinom(N : Integer; P : Float; K : Integer) : Float; { Prob(X = K) }
function FBinom(N : Integer; P : Float; K : Integer) : Float; { Prob(X <= K) }
{ ----------------------------------------------------------------------
Poisson distribution with mean Mu
---------------------------------------------------------------------- }
function PPoisson(Mu : Float; K : Integer) : Float; { Prob(X = K) }
function FPoisson(Mu : Float; K : Integer) : Float; { Prob(X <= K) }
{ ----------------------------------------------------------------------
Standard normal distribution
---------------------------------------------------------------------- }
function DNorm(X : Float) : Float; { Density of standard normal }
function FNorm(X : Float) : Float; { Prob(U <= X) }
function PNorm(X : Float) : Float; { Prob(|U| >= |X|) }
function InvNorm(P : Float) : Float; { Inverse of FNorm : returns X
such that Prob(U <= X) = P}
{ ----------------------------------------------------------------------
Student distribution with Nu d.o.f.
---------------------------------------------------------------------- }
function DStudent(Nu : Integer; X : Float) : Float; { Density of t }
function FStudent(Nu : Integer; X : Float) : Float; { Prob(t <= X) }
function PStudent(Nu : Integer; X : Float) : Float; { Prob(|t| >= |X|) }
{ ----------------------------------------------------------------------
Khi-2 distribution with Nu d.o.f.
---------------------------------------------------------------------- }
function DKhi2(Nu : Integer; X : Float) : Float; { Density of Khi2 }
function FKhi2(Nu : Integer; X : Float) : Float; { Prob(Khi2 <= X) }
function PKhi2(Nu : Integer; X : Float) : Float; { Prob(Khi2 >= X) }
{ ----------------------------------------------------------------------
Fisher-Snedecor distribution with Nu1 and Nu2 d.o.f.
---------------------------------------------------------------------- }
function DSnedecor(Nu1, Nu2 : Integer; X : Float) : Float; { Density of F }
function FSnedecor(Nu1, Nu2 : Integer; X : Float) : Float; { Prob(F <= X) }
function PSnedecor(Nu1, Nu2 : Integer; X : Float) : Float; { Prob(F >= X) }
{ ----------------------------------------------------------------------
Exponential distribution
---------------------------------------------------------------------- }
function DExpo(A, X : Float) : Float; { Density of exponential distrib. }
function FExpo(A, X : Float) : Float; { Prob( <= X) }
{ ----------------------------------------------------------------------
Beta distribution
---------------------------------------------------------------------- }
function DBeta(A, B, X : Float) : Float; { Density of Beta distribution }
function FBeta(A, B, X : Float) : Float; { Prob( <= X) }
{ ----------------------------------------------------------------------
Gamma distribution
---------------------------------------------------------------------- }
function DGamma(A, B, X : Float) : Float; { Density of Gamma distribution }
function FGamma(A, B, X : Float) : Float; { Prob( <= X) }
{ ********************************************************************** }
implementation
{ ----------------------------------------------------------------------
Error handling function
---------------------------------------------------------------------- }
function DefaultVal(ErrCode : Integer; Default : Float) : Float;
{ Sets the global variable MathError and the function default value
according to the error code }
begin
MathError := ErrCode;
DefaultVal := Default;
end;
{ ----------------------------------------------------------------------
Special functions
---------------------------------------------------------------------- }
const { Used by IGamma and IBeta }
BIG = 9.223372036854775808E18;
BIGINV = 1.084202172485504434007E-19;
type
TabCoef = array[0..9] of Float;
function PolEvl(var X : Float; Coef : TabCoef; N : Integer) : Float;
{ ----------------------------------------------------------------------
Evaluates polynomial of degree N:
2 N
y = C + C x + C x +...+ C x
0 1 2 N
Coefficients are stored in reverse order:
Coef[0] = C , ..., Coef[N] = C
N 0
The function P1Evl() assumes that Coef[N] = 1.0 and is
omitted from the array. Its calling arguments are
otherwise the same as PolEvl().
---------------------------------------------------------------------- }
var
Ans : Float;
I : Integer;
begin
Ans := Coef[0];
for I := 1 to N do
Ans := Ans * X + Coef[I];
PolEvl := Ans;
end;
function P1Evl(var X : Float; Coef : TabCoef; N : Integer) : Float;
{ ----------------------------------------------------------------------
Evaluate polynomial when coefficient of X is 1.0.
Otherwise same as PolEvl.
---------------------------------------------------------------------- }
var
Ans : Float;
I : Integer;
begin
Ans := X + Coef[0];
for I := 1 to N - 1 do
Ans := Ans * X + Coef[I];
P1Evl := Ans;
end;
function SgnGamma(X : Float) : Integer;
begin
if X > 0.0 then
SgnGamma := 1
else if Odd(Trunc(Abs(X))) then
SgnGamma := 1
else
SgnGamma := - 1;
end;
function Stirf(X : Float) : Float;
{ Stirling's formula for the gamma function
Gamma(x) = Sqrt(2*Pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
where P(x) is a polynomial }
const
STIR : TabCoef = (
7.147391378143610789273E-4,
- 2.363848809501759061727E-5,
- 5.950237554056330156018E-4,
6.989332260623193171870E-5,
7.840334842744753003862E-4,
- 2.294719747873185405699E-4,
- 2.681327161876304418288E-3,
3.472222222230075327854E-3,
8.333333333333331800504E-2,
0);
var
W, P : Float;
begin
W := 1.0 / X;
if X > 1024.0 then
begin
P := 6.97281375836585777429E-5 * W + 7.84039221720066627474E-4;
P := P * W - 2.29472093621399176955E-4;
P := P * W - 2.68132716049382716049E-3;
P := P * W + 3.47222222222222222222E-3;
P := P * W + 8.33333333333333333333E-2;
end
else
P := PolEvl(W, STIR, 8);
Stirf := SQRT2PI * Exp((X - 0.5) * Ln(X) - X) * (1.0 + W * P);
end;
function GamSmall(X1, Z : Float) : Float;
{ Gamma function for small values of the argument }
const
S : TabCoef = (
- 1.193945051381510095614E-3,
7.220599478036909672331E-3,
- 9.622023360406271645744E-3,
- 4.219773360705915470089E-2,
1.665386113720805206758E-1,
- 4.200263503403344054473E-2,
- 6.558780715202540684668E-1,
5.772156649015328608253E-1,
1.000000000000000000000E0,
0);
SN : TabCoef = (
1.133374167243894382010E-3,
7.220837261893170325704E-3,
9.621911155035976733706E-3,
- 4.219773343731191721664E-2,
- 1.665386113944413519335E-1,
- 4.200263503402112910504E-2,
6.558780715202536547116E-1,
5.772156649015328608727E-1,
- 1.000000000000000000000E0,
0);
var
P : Float;
begin
if X1 = 0.0 then
begin
GamSmall := DefaultVal(FN_SING, MAXNUM);
Exit;
end;
if X1 < 0.0 then
begin
X1 := - X1;
P := PolEvl(X1, SN, 8);
end
else
P := PolEvl(X1, S, 8);
GamSmall := Z / (X1 * P);
end;
function StirfL(X : Float) : Float;
{ Approximate Ln(Gamma) by Stirling's formula, for X >= 13 }
const
P : TabCoef = (
4.885026142432270781165E-3,
- 1.880801938119376907179E-3,
8.412723297322498080632E-4,
- 5.952345851765688514613E-4,
7.936507795855070755671E-4,
- 2.777777777750349603440E-3,
8.333333333333331447505E-2,
0, 0, 0);
var
Q, W : Float;
begin
Q := Ln(X) * (X - 0.5) - X;
Q := Q + LNSQRT2PI;
if X > 1.0E+10 then
StirfL := Q
else
begin
W := 1.0 / Sqr(X);
StirfL := Q + PolEvl(W, P, 6) / X;
end;
end;
function Gamma(X : Float) : Float;
const
P : TabCoef = (
4.212760487471622013093E-5,
4.542931960608009155600E-4,
4.092666828394035500949E-3,
2.385363243461108252554E-2,
1.113062816019361559013E-1,
3.629515436640239168939E-1,
8.378004301573126728826E-1,
1.000000000000000000009E0,
0, 0);
Q : TabCoef = (
- 1.397148517476170440917E-5,
2.346584059160635244282E-4,
- 1.237799246653152231188E-3,
- 7.955933682494738320586E-4,
2.773706565840072979165E-2,
- 4.633887671244534213831E-2,
- 2.243510905670329164562E-1,
4.150160950588455434583E-1,
9.999999999999999999908E-1,
0);
var
SgnGam, N : Integer;
A, X1, Z : Float;
begin
MathError := FN_OK;
SgnGam := SgnGamma(X);
if (X = 0.0) or ((X < 0.0) and (Frac(X) = 0.0)) then
begin
Gamma := DefaultVal(FN_SING, SgnGam * MAXNUM);
Exit;
end;
if X > MAXGAM then
begin
Gamma := DefaultVal(FN_OVERFLOW, MAXNUM);
Exit;
end;
A := Abs(X);
if A > 13.0 then
begin
if X < 0.0 then
begin
N := Trunc(A);
Z := A - N;
if Z > 0.5 then
begin
N := N + 1;
Z := A - N;
end;
Z := Abs(A * Sin(PI * Z)) * Stirf(A);
if Z <= PI / MAXNUM then
begin
Gamma := DefaultVal(FN_OVERFLOW, SgnGam * MAXNUM);
Exit;
end;
Z := PI / Z;
end
else
Z := Stirf(X);
Gamma := SgnGam * Z;
end
else
begin
Z := 1.0;
X1 := X;
while X1 >= 3.0 do
begin
X1 := X1 - 1.0;
Z := Z * X1;
end;
while X1 < - 0.03125 do
begin
Z := Z / X1;
X1 := X1 + 1.0;
end;
if X1 <= 0.03125 then
Gamma := GamSmall(X1, Z)
else
begin
while X1 < 2.0 do
begin
Z := Z / X1;
X1 := X1 + 1.0;
end;
if (X1 = 2.0) or (X1 = 3.0) then
Gamma := Z
else
begin
X1 := X1 - 2.0;
Gamma := Z * PolEvl(X1, P, 7) / PolEvl(X1, Q, 8);
end;
end;
end;
end;
function LnGamma(X : Float) : Float; overload;
const
P : TabCoef = (
- 2.163690827643812857640E3,
- 8.723871522843511459790E4,
- 1.104326814691464261197E6,
- 6.111225012005214299996E6,
- 1.625568062543700591014E7,
- 2.003937418103815175475E7,
- 8.875666783650703802159E6,
0, 0, 0);
Q : TabCoef = (
- 5.139481484435370143617E2,
- 3.403570840534304670537E4,
- 6.227441164066219501697E5,
- 4.814940379411882186630E6,
- 1.785433287045078156959E7,
- 3.138646407656182662088E7,
- 2.099336717757895876142E7,
0, 0, 0);
var
N, SgnGam : Integer;
A, X1, Z : Float;
begin
MathError := FN_OK;
SgnGam := SgnGamma(X);
if (X = 0.0) or ((X < 0.0) and (Frac(X) = 0.0)) then
begin
LnGamma := DefaultVal(FN_SING, MAXNUM);
Exit;
end;
if X > MAXLGM then
begin
LnGamma := DefaultVal(FN_OVERFLOW, MAXNUM);
Exit;
end;
A := Abs(X);
if A > 34.0 then
begin
if X < 0.0 then
begin
N := Trunc(A);
Z := A - N;
if Z > 0.5 then
begin
N := N + 1;
Z := N - A;
end;
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -