unit2.pas

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

unit Unit2;

interface
uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs;
Function GAMMLN(xx:real):real;
Function FACTRL(N:integer):REAL;
procedure DDPOLY(C:array of real; NC:integer; X:real;
                                  var PD:array of real; ND:integer);
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 FACTRL(N:integer):REAL;
var
  GLNTOP:integer;
  GLA:ARRAY[1..33] OF real;
  J:INTEGER;
BEGIN
  GLNTOP:=0;
  GLA[1]:= 1;
  If N < 0 Then  ShowMessage('警告：现在N<0，请修改程序');
  If N = 0 Then
    FACTRL:= GLA[N+1]
  else If N <= 32 Then
    BEGIN
      For J:= GLNTOP + 1 To N do
        GLA[J + 1]:= J * GLA[J];
      GLNTOP:= N;
      FACTRL:= GLA[N + 1];
    END
  else If N>32 then
    FACTRL:= Exp(GAMMLN(N + 1));
END;

procedure DDPOLY(C:array of real; NC:integer; X:real;
                                  var PD:array of real; ND:integer);
var
    J,I,NND:integer;  CONST1:real;
begin
    PD[1]:=C[NC];
    For J:=2 To ND do
        PD[J]:=0; 
    For I:=NC - 1 downTo 1 do
       begin
        If NC + 1 - I < ND Then
           NND:=NC + 1 - I
        Else
           NND:=ND;
        For J:=NND DownTo 2 do
              begin
            PD[J]:=PD[J] * X + PD[J - 1];
        end;
        PD[1]:=PD[1] * X + C[I];
    end;
    CONST1:=2 ;
    For I:=3 To ND do
       begin
        PD[I]:=CONST1 * PD[I];
        CONST1:=CONST1 * I;
    end;
end;
end.