unit2.pas

《Delphi常用数值算法集》的配书源码

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
Function GAMMLN(xx:real):real;
Function BETACF( A, B, X:real):real;
Function BETAI(A,B,X:real):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 BETACF( A, B, X:real):real;
label 1;
const
  ITMAX=100;    EPS=0.0000003;
var
  TEM,QAP,QAM,QAB,EM,D:REAL; BZ,BP,BPP,BM,AZ,AAP:REAL;
  AM,AOLD,AP:REAL;
  M:INTEGER;
begin
  AM:=1;
  BM:=1;
  AZ:=1;
  QAB:=A + B;
  QAP:=A + 1;
  QAM:=A - 1;
  BZ:=1 - QAB * X / QAP;
  For M:=1 To ITMAX DO
  begin
    EM:=M;
    TEM:=EM + EM;
    D:=EM * (B - M) * X / ((QAM + TEM) * (A + TEM));
    AP:=AZ + D * AM;
    BP:=BZ + D * BM;
    D:=-(A + EM) * (QAB + EM) * X / ((A + TEM) * (QAP + TEM));
    AAP:=AP + D * AZ;
    BPP:=BP + D * BZ;
    AOLD:=AZ;
    AM:=AP / BPP;
    BM:=BP / BPP;
    AZ:=AAP / BPP;
    BZ:=1;
    If Abs(AZ - AOLD) < EPS * Abs(AZ) Then GoTo 1;
  end;
  ShowMessage('A or B too big, or ITMAX too small');
  exit;
1: BETACF:=AZ;
end;

Function BETAI(A,B,X:real):real;
var
  AAA,BT:real;
begin
  If (X < 0) Or (X > 1) Then ShowMessage('bad argument X in BETAI');
  If (X=0) Or (X=1) Then
    BT:=0
  Else
    begin
      AAA:=GAMMLN(A + B) - GAMMLN(A) - GAMMLN(B);
      BT:=Exp(AAA + A * Ln(X) + B * Ln(1 - X));
    end;
  If X < (A + 1) / (A + B + 2) Then
    BETAI:=BT * BETACF(A, B, X) / A
  Else
    BETAI:=1 - BT * BETACF(B, A, 1 - X) / B;
end;

end.