specfunc.pas

Delphi math processing compononets and sources. Release.

    nu := 4.0*sqr(nu);  m := 1.0;  n := 1.0;
    x8 := 8.0*x;        t := 1.0;  s := t;
    repeat
      f := (nu-sqr(m))/n/x8;  t := -t*f;
      s := s + t;
      m := m + 2.0;  n := n + 1.0;
    until (abs(t)<epsilon) or (abs(f) > 1.0);
    I_nu := s*exp1(x)/s2pi/sqrt(x);
  end;
end;

{ Undocumented Functions }

function Gamma(x: RealType): RealType;
{    x <> 0.0, -IntType   }
const
  Ck: array [2..26] of RealType =
  ( 0.5772156649015329, -0.6558780715202538, -0.0420026350340952,  
      0.1665386113822915, -0.0421977345555443, -0.0096219715278770,  
      0.0072189432466630, -0.0011651675918591, -0.0002152416741149,  
      0.0001280502823882, -0.0000201348547807, -0.0000012504934821, 
      0.0000011330272320, -0.0000002056338417,  0.0000000061160950,  
      0.0000000050020075, -0.0000000011812746,  0.0000000001043427,  
      0.0000000000077823, -0.0000000000036968,  0.0000000000005100, 
     -0.0000000000000206, -0.0000000000000054,  0.0000000000000014,  
      0.0000000000000001 );
(*)
const
  a1 = 1.0/12.0;          a2 = 1.0/288.0;
  a3 = 139.0/51840.0;     a4 = 571.0/2488320.0;
(*)
var
  B,D,F,BB,T    : RealType;
  i             : byte;
begin
  B := abs(x);
  D := 1.0;
  while B>1.0  do
  begin
    D := D*B;  B := B-1.0;
  end;
  F  := 1.0;
  BB := B;
  T  := BB;
  for i:=2 to 26  do
    begin
      F := F + Ck[i]*T;     T := T*BB;
    end;
  if x>0.0  then  Gamma := D/x/F
            else  Gamma := F*Pi/sin(Pi*x)/D;
(*)
  F := exp1(-B+ln(B)*(B-0.5))*sqrt(2.0*Pi) *
  (1.0 + (a1 + (a2 - (a3 + a4/B)/B)/B)/B );
  if x>0.0  then  Gamma := F
            else  Gamma := Pi/sin(Pi*x)/x;
(*)
end;

function  PolyGamma(n: IntType; z: RealType): RealType;
{                    n+1                 n+1
   PolyGamma(n,z) = d   ln(Gamma(z)) / dz   ,  z > 10, 0<=n<=3  }      
const
  BernN = 32;
  Bernoulli: array [1..BernN] of RealType =
   ( 1.0,                    -0.5,                      1.0/6.0,
    -1.0/30.0,                1.0/42.0,                -1.0/30.0,
     5.0/66.0,               -691.0/2730.0,             7.0/6.0,
    -3617.0/510.0,            43867.0/798.0,           -174611.0/330.0,
     854513.0/138.0,         -236364091.0/2730.0,       8553103.0/6.0,
    -23749461029.0/870.0,     8615841276005.0/14322.0, -7709321041217.0/510.0,
     2577687858367.0/6.0,    -1.371165521e+13,          4.883323190e+14,
    -1.929657934e+16,         8.416930476e+17,         -4.033807185e+19,
     2.115074864e+21,        -1.208662652e+23,          7.500866746e+24,
    -5.038778101e+26,         3.652877648e+28,         -2.849876930e+30,
     2.386542750e+32,        -2.139994926e+34 );

 function PGam_Assc(n: IntType; z: RealType): RealType;
 var
   S,A,A1,zz,Item,N1fac,Nom,Denom  : RealType;
   i                               : IntType;
 begin   
   zz  := sqr(z);      S     := 0.0;
   Nom := n;           N1fac := 1.0;
   for i:=1 to n-1  do
     N1fac := N1fac*i;    {  (n-1)!  }
   Item := N1fac;      A     := MaxReal;
   i    := 3;          Denom := 2.0;
   repeat
     A1 := A;
     if n>0  then
     begin
       Item := Item*Nom/Denom*(Nom+1.0)/(Denom-1.0)/zz;
       Nom  := Nom+2.0;
       A    := Bernoulli[i]*Item;
     end   else
     begin
       Item := Item/zz;
       A    := Bernoulli[i]*Item/Denom;
     end;
     Denom := Denom+2.0;
     inc(i);
     S := S+A;
   until (S=S+A) or (abs(A)>abs(A1)) or (i>BernN);
     
   if n=0  then  PGam_Assc := ln(z) - 0.5/z - S
   else
   begin
     S := S + N1fac*(1.0+n/(2.0*z));
     for i:=1 to n  do
       S := S/z;
     if Odd(n)  then  PGam_Assc := S
                else  PGam_Assc := -S;
   end;
end;
 
const
  AssLim = 30.0;
var
  S,Nfac,A,zz,PG,Sz,Cz : RealType;
  i,j,k                : IntType;
begin
  zz := abs(z); 
  if z<0.0  then  zz := z+1.0;
  if zz>=AssLim  then  PG := PGam_Assc(n,zz )
  else
  begin
    k := trunc(AssLim-zz) + 1;
    S := 0.0;
    for i:=k-1 downto 0  do
      begin
        A := 1.0;
        for j:=1 to n+1  do
          A := A/(zz+i);
        S := S+A;
      end;
    Nfac := 1.0;
    if n>0  then
      for i:=1 to n  do
        Nfac := Nfac*(-i);
    PG := PGam_Assc(n,k+zz) - Nfac*S;
  end;
  if z>=0.0  then
    PolyGamma := PG
  else     
  begin
    Sz := sin(Pi*zz);     Cz := cos(Pi*zz);
    S  := Pi/Sz;
    case  n  of
      0: PolyGamma := PG + Cz*S;
      1: PolyGamma := -PG - sqr(S);
      2: PolyGamma := PG + 2.0*Cz*S*sqr(S);
      3: PolyGamma := -PG - 2.0*Pi*S*sqr(S)*(1.0+2.0*sqr(Cz));
       else
         PolyGamma := PG;  {  wrong, of course  }
    end;
  end;
end;
 
{  ----------------------------------------------------------  }

function GammaX(a,x: RealType): RealType;
  
 function GammaAbsX(x: RealType): RealType;
 var
   n,S,T,dS : RealType;
 begin
   if x=0.0  then   GammaAbsX := 0.0
   else if x<=12.0  then
     begin
       T := 1.0;        dS := T/a;
       S := 0.0;        n  := 1.0;
       repeat
         S  := S+dS;
         T  := T/n;
         T  := T*(-x);
         dS := T/(a+n);
         n  := n+1.0;
       until ((S>0.0) and (S+dS=S)) or (n>64000.0);
       if a<=3.0  then  GammaAbsX := exp1(a*ln(x))*S
                  else  GammaAbsX := exp1(a*ln(x)+ln(S));
     end
   else
     begin
       T := 1.0;    S := 0.0;
       n := 1.0;
       repeat
         S  := S + T;
         dS := T;
         T  := T/x;
         T  := T*(a-n);
         n  := n+1.0;
       until (S+T=S) or ((dS<0.0) and (abs(dS)<abs(T)));
       GammaAbsX := Gamma(a) - exp1((a-1.0)*ln(x)-x+ln(S));
     end;
 end;

begin
  if x>=0.0  then  GammaX := GammaAbsX(x )
             else  GammaX := -GammaAbsX(-x );
end;

{  ----------------------------------------------------------  }

function IJn(n,x,W: RealType): RealType;
var
  F,R,i  : RealType;
begin
  if  x=0.0  then
    begin
      if n=0.0  then  IJn := 1.0
                else  IJn := 0.0;
      Exit;
    end;
  if n<>0.0  then   F := 1.0/Gamma(n)/n
             else   F := 1.0;
  R := F;
  i := 0.0;
  repeat
    i := i+1.0;
    R := R*W;
    R := R/i/(n+i);
    F := F+R;
  until (F-R=F) or (i>64000.0);
  if  W<=0.0  then   IJn := exp1(n*ln(x/2.0))*F
              else   IJn := exp1(n*ln(x/2.0)+ln(F));
end;

function BJn(n,x: RealType): RealType;
const
  Xmax = 23.0;
var
  k                        : byte;
  s,f,nu,x8,Sn,Cs,Ck,m,p   : RealType;
begin
  if x<=Xmax  then
    BJn := IJn(n,x,-sqr(x)/4.0 )
  else
    begin
      s  := x-Pi*(n+0.5)/2.0;
      Sn := sin(s);         Cs := cos(s);
      Ck := 1.0;            s  := 0.0;
      nu := 4.0*sqr(n);     x8 := 8.0*x;
      k  := 0;              m  := 1.0;
      p  := 1.0;
      repeat
        case  k  of
          0: s := s + Ck*Cs;
          1: s := s - Ck*Sn;
          2: s := s - Ck*Cs;
          3: s := s + Ck*Sn;
        end;
        k  := (k+1) mod 4;
        f  := (nu-sqr(m))/p/x8;
        Ck := Ck*f;
        m  := m + 2.0;
        p  := p + 1.0;
      until ((s+Ck=s) and (p>=n)) or (f<=-1.0);
      BJn := sqrt(2.0/(Pi*x))*s;
    end;
end;

function BIn(n,x: RealType): RealType;
const
  Xmax    = 35.0;
  s2pi    = 2.506628274631;  { sqrt(2*pi) }
var
  s,t,m,f,nu,x8,p    : RealType;
begin
  if x<Xmax  then   BIn := IJn(n,x,+sqr(x)/4.0 )
  else
    begin
      nu := 4.0*sqr(n);   m := 1.0;  p := 1.0;
      x8 := 8.0*x;        t := 1.0;  s := t;
      repeat
        f := (nu-sqr(m))/p/x8;
        t := -t*f;
        s := s + t;
        m := m + 2.0;
        p := p + 1.0;
      until (s-t=s) or (f<=-1.0);
      BIn := s*exp1(x)/s2pi/sqrt(x);
    end;
end;

function BKn(n,x: RealType): RealType;
const
  Dmin = 1.0e-4;
var
  dn: RealType;
 function BK(nu: RealType): RealType;
 begin
   BK := Pi/2.0*(BIn(-nu,x)-BIn(nu,x))/sin(nu*Pi);
 end;
begin
  dn := frac(n);
  if abs(dn)>=Dmin  then BKn := BK(n)
  else
    BKn := BK(round(n)+Dmin)*((Dmin+dn)/(2.0*Dmin)) +
           BK(round(n)-Dmin)*((Dmin-dn)/(2.0*Dmin));
end;
 
{  ----------------------------------------------------------  }

const
  Z0   = 20.0;
  Nmax = 20;
  S_Coef_Done: boolean = false;
var
  S_Coef     : array [1..Nmax+1] of RealType;

function  Sn( n: IntType; z: RealType): RealType;

var
  i          : IntType;
  lnZ,lniZ,S : RealType;
  l          : IntType;

 function  S_n(x: RealType; n: IntType): RealType;
 var
   S,Item,k,A  : RealType;
   i           : IntType;
 begin
   if n=0  then
   begin
     S_n := exp1(-x)-1.0;
     Exit;
   end;
   S := -x;   Item := S;    k := 1.0;
   repeat
     k    := k+1.0;
     Item := Item/k;
     Item := Item*(-x);
     A    := 1.0;
     for i:=1 to n  do  A := A*k;
     A    := Item/A;
     S := S+A;
   until  S=S+A;
   S_n := S;
 end;

 procedure  InitSn;
 var
   lnZ0,Ci,lnjZ0: RealType;
   i,j, l          : IntType;
 begin
   lnZ0 := ln(Z0);
   for i:=0 to Nmax  do
     S_Coef[i+1] := S_n(Z0,i);
   for i:=0 to Nmax  do
   begin
     Ci    := -S_Coef[i+1];
     l     := 1;
     lnjZ0 := 1.0;
     for j:=1 to i  do
     begin
       l     := l*j;
       lnjZ0 := lnjZ0*lnZ0;
       Ci    := Ci - S_Coef[i-j+1]*lnjZ0/l;
     end;
     S_Coef[i+1] := Ci;
   end;
 end;

begin
  if  (n=0) or (z<=Z0) or (n>Nmax)  then   Sn := S_n(z,n)
  else
    begin
      if not S_Coef_Done  then
        begin
          InitSn;  S_Coef_Done := true;
        end;
      l    := 1;
      lnZ  := ln(Z);
      lniZ := 1.0;
      S    := -S_Coef[n+1];
      for i:=1 to n  do
        begin
          l    := l*i;
          lniZ := lniZ*lnZ;
          S    := S-S_Coef[n-i+1]*lniZ/l;
        end;
      Sn := S;
    end;
end;
 
{  ----------------------------------------------------------  }

function  Un( n: IntType; z,a: RealType): RealType;

 function  U_n(x,a: RealType; n: IntType): RealType;
 {    Direct Sum calculations   }
 var
   S,Item,k,B  : RealType;
   i           : IntType;
 begin
   if n=0  then
     begin
       U_n := exp1(-x)-1.0;
       Exit;
     end;
   Item := 1.0;   k := 0.0;   S := 0.0; //  B := MaxReal;
   repeat                                              
     k    := k+1.0;
     Item := Item/k;
     Item := Item*(-x);
     B    := Item;
     for i:=1 to n  do  B := B/(k+a);
     S := S+B;
   until  S=S+B;
   U_n := S;
 end; 

 function  Un_Ass(n: IntType; z,a: RealType): RealType;
 var
   S,G1,PG,G2,P,S1,S2,A1,A2,MachLimit : RealType;
   ZF                                 : boolean;
   k                                  : IntType;
 begin
   {    Asymptotics for large  z,   n=2,3     }
   S := ln(z);
   G1 := exp1(-a*S)*Gamma(a);
   PG := PolyGamma(0,a)-S;
   if n=2  then  G1 := G1*PG
           else  G1 := G1*(sqr(PG)+PolyGamma(1,a));
           
  {    questionable exectness off  GammaX    }
   S := 0.0;    k  := 1;     A1 := MaxReal;
   P := 1.0;    S1 := 0.0;   S2 := 0.0;
   MachLimit := 32*sqrt(MachEps);
   ZF := false;
   repeat
     A2 := A1;
     PG := a-k;
     if abs(PG)>MachLimit  then
     begin
       P  := P*PG;  S1 := S1+1.0/PG;
       S2 := S2+1.0/sqr(PG);
     end else
       ZF := true;
     P := P/z;
     if n=2  then
     begin
       if ZF  then  A1 := P
              else  A1 := P*S1;
     end else
       if k>1  then
       begin
         if ZF  then  A1 := P*S1
                else  A1 := P*(sqr(S1)-S2);
         S := S + A1;
       end;
     inc(k);
   until (S+A1=S) or ((A2<0.0) and (abs(A2)<abs(A1))) or (k>512);
   G2 := exp1(-z)/z*S;                         
  (*   else  G2 := 0.0     *)
   if n=2  then  Un_Ass := -(G1-G2+1.0/sqr(a))
           else  Un_Ass := (G1-G2-2.0/sqr(a)/a)/2.0;
 end;
 
const
  ZLim = 20.0;
    
begin
  if (a=0.0) or (n=0)  then   Un := Sn(n,z)
  else if (z<=ZLim) or (n>3)  then  Un := U_n(z,a,n)
  else if n=1  then              
    {    exact formulae for arbitrary z>0.0   }
    Un := exp1(-a*ln(z))*GammaX(a,z) - 1.0/a
  else Un := Un_Ass(n,z,a);
end;


end.