unit2.pas

用于开发税务票据管理的软件

unit Unit2;

interface
Function GAMMLN(xx:real):real;
Function T(Y:real):real;
procedure IAP(X, A:real; NMAX:integer; VAR F:array of real);
procedure BESIAN(X, A:real; NM, IH:integer; VAR F:array of real);

implementation
Function GAMMLN(xx:real):real;
CONST
  STP=2.50662827465; HALF=0.5; ONE=1.0; FPF=5.5;
var
  x,tmp,ser:double;
  j:integer;
  cof:array[1..6] of double;
BEGIN
  COF[1]:=76.18009173;      COF[2]:=-86.50532033;
  COF[3]:=24.01409822;      COF[4]:=-1.231739516;
  COF[5]:= 0.120858003e-2;  COF[6]:=-0.536382e-5;
  X:=XX-ONE;
  TMP:=X+FPF;
  TMP:=(X+HALF)*Ln(TMP)-TMP;
  SER:=ONE;
  For J:=1 To 6 do
    BEGIN
      X:=X+ONE;
      SER:=SER+COF[J]/X
    END;
  GAMMLN:=TMP+Ln(STP*SER);
end;

Function T(Y:real):real;
var
    AAA,P,Z:REAL;
begin
    If Y <= 10 Then
    begin
        AAA:=0.000057941 * Y - 0.00176148 * Y + 0.0208645;
        T:=((AAA * Y - 0.129013) * Y + 0.85777) * Y + 1.0125;
    end
    Else
    begin
        Z:=Ln(Y) - 0.775;
        P:=(0.775 - Ln(Z)) / (1 + Z);
        T:=Y / (P + 1) / Z;
    end;
end;

procedure IAP(X, A:real; NMAX:integer; var F:array of real);
label 1,2,3;
const
    MMAX = 10;    E = 4;
var
    FA,RR:array[0..MMAX] of real;
    EPS,SUM,R,S,AL,AM,D1:real;
    N,NU,M:integer;
begin
    EPS:= 0.5*EXP(-E * Ln(10));
    For N:=0 To NMAX do  FA[N]:=0;
    Sum:=Exp(GAMMLN(1 + A));
    Sum:=Exp(A * Ln(X / 2)) * Exp(X) / Sum;
    D1:=2.3026 * E + 1.3863;
    If NMAX > 0 Then
        R:=T(0.5 * D1 / NMAX) * NMAX
    Else
        R:=0;
    If X < D1 Then
        S:=T(0.73576 * (D1 - X) / X) * 1.3591 * X
    Else
        S:=1.3591 * X;
    If R <= S Then
        NU:=1 + Trunc(S)
    Else
        NU:=1 + Trunc(R);
1:  N:=0;
    AL:=1;
2:  N:=N + 1;
    AL:=AL * (N + 2 * A) / (N + 1);
    If N < NU Then GoTo 2;
    R:=0;
    S:=0;
3:  R:=1 / (2 * (A + N) / X + R);
    AL:=AL * (N + 1) / (N + 2 * A);
    AM:=2 * (N + A) * AL;
    S:=R * (AM + S);
    If N <= NMAX Then RR[N - 1]:=R;
    N:=N - 1;
    If N >= 1 Then GoTo 3;
    F[0]:=Sum / (1 + S);
    For N:=0 To NMAX - 1 do
        F[N + 1]:=RR[N] * F[N];
    For N:=0 To NMAX do
    begin
        If Abs((F[N] - FA[N]) / F[N]) > EPS Then
        begin
            For M:=0 To NMAX do
                FA[M]:=F[M];
            NU:=NU + 5;
            GoTo 1;
        end;
    end;
end;

procedure BESIAN(X, A:real; NM, IH:integer; var F:array of real);
var
    I:integer;
begin
    If IH = 1 Then  IAP(X, A, NM, F);
    If IH = -1 Then
    begin
        IAP(X, A, 1, F);
        F[1]:=2 * A * F[0] / X + F[1];
        For I:=1 To NM - 1 do
            F[I + 1]:=2 * (A - I) * F[I] / X + F[I - 1];
    end;
end;

end.