specfunc.pas

Delphi math processing compononets and sources. Release.

{
@abstract(EBK&NVS Pascal-Delphi Math Library: Real Special functions)
@author(Nikolai V. Shokhirev <[e-mail]> <[e-mail]>)
@author(Eugene B. Krissinel <[e-mail]> <[e-mail]>)
@created(09.09.1991)
@lastmod(04.04.2002)
This is a temporary publication (reduced variant), will be updated later
check: http://www.shokhirev.com/nikolai/programs/samplecode.html
㎞ikolai V. Shokhirev, 2002
}
unit SpecFunc;

interface

uses
  MathTypes;

const
  ln10     = 2.3025850929940456840179915;
  Gamma_Eu = 0.577215664901532860606512;

{ Basic functions for the Real Calculations }
function iSign(x: IntType)      : IntType;
function Sign(x: RealType)      : RealType;
function Sign2(a, b: RealType)  : RealType;
function Combinations(m, n: IntType): IntType;
function Fac(m: IntType)        : IntType;
function power10(x: RealType)   : RealType;
function lg(x: RealType)        : RealType;
function amin2(x1, x2: RealType): RealType;
function amax2(x1, x2: RealType): RealType;
function sqroot1(x, y: RealType): RealType; { safe sqrt(x + y) }
function sqroot2(x, y: RealType): RealType; { safe sqrt(x*x+y*y) }

{ Hyperbolic Functions }
function ArcSin(x: RealType): RealType;
function ArcCos(X: RealType): RealType;
function exp1(X: RealType)  : RealType;
function ch(x: RealType)    : RealType;
function sh(x: RealType)    : RealType;
function th    (x: RealType): RealType;
function arch  (x: RealType): RealType;
function arsh  (x: RealType): RealType;
function arth  (x: RealType): RealType;

{ Special functions for the Real Calculations }
function er(x: Realtype)        : RealType; { er(x) = exp(x**2)*erfc(x)    }
function erf(x: Realtype)       : RealType; { error function: 1 = erfc+erf }
function erfc(x: Realtype)      : RealType; { complementary error function }
function I0_(x: RealType)       : RealType; { Modified Bessel function     }
function I1_(x: RealType)       : RealType; { Modified Bessel function     }
function K0_(x: RealType)       : RealType; { Modified Bessel function     }
function K1_(x: RealType)       : RealType; { Modified Bessel function     }
function fk0_(x: RealType)      : RealType; { fk0(x)=K0(x)+ln(x/2)*I0(x)   }
function fk1_(x: RealType)      : RealType; { fk1(x)=K1(x)-ln(x/2)*I1(x)   }
function gamma_(x: Realtype)    : RealType; { gamma-function               }
function gamma_in(a,x: Realtype): RealType; { incomplete gamma-function    }
function E1_(x: Realtype)       : RealType; { integral exp-function        }
function I_nu(nu, x: RealType)  : RealType; { Modified Bessel function     }

(*)
{ Undocumented Functions }
function Gamma(x: RealType): RealType;
function  PolyGamma(n: IntType; z: RealType): RealType;
function  GammaX(a,x: RealType): RealType;
function  BJn   (n,x: RealType): RealType;
function  BIn   (n,x: RealType): RealType;
function  BKn   (n,x: RealType): RealType;
function  Sn(n: IntType; z: RealType): RealType;
function  Un(n: IntType; z,a: RealType): RealType;
(*)
{ ===================================================================  }

implementation

{ Basic functions for the Real Calculations }

function  iSign (x: IntType): IntType;
begin
  if x > 0 then result := 1 else
    if x = 0 then result := 0 else
      result := -1;
end;

function sign(x: RealType): RealType;
begin
if x > 0.0 then sign := 1.0
           else if x < 0.0 then sign := -1.0
                           else sign := 0.0;
end;

function Sign2(a, b: RealType): RealType;
begin
  if b < 0.0 then result := - abs(a)
             else result :=   abs(a) ;
end;

function Combinations(m, n: IntType): IntType;
var
  i, C   : IntType;
begin
  C := 1;
  if m > 0  then
    for i := 1 to m  do
      C := (C*(n-i+1)) div i;
  Combinations := C;
end;

function Fac(m: IntType): IntType;
var
  i, z: IntType;
begin
  z := 1;
  if m > 1  then
    for i := 2 to m  do
      z := z*i;
  Fac := z;
end;

function power10(x: RealType): RealType;
begin
if x > -MaxExp then power10 := exp1(x*ln(10.0))
               else power10 := 1.0;
end;

function lg(x: RealType): RealType;
begin
  lg := ln(x)/ln(10.0);
end;

function amin2(x1,x2: RealType): RealType;
begin
if x1 < x2 then amin2 := x1
           else amin2 := x2;
end;

function amax2(x1,x2: RealType): RealType;
begin
if x1 > x2 then amax2 := x1
           else amax2 := x2;
end;

function  ArcSin (x: RealType): RealType;
var
  R,S,y    : RealType;
  l        : IntType;
begin
  if  x>=1.0    then    ArcSin := Pi/2.0
  else if x<=-1.0 then  ArcSin := -Pi/2.0
  else
    begin
      if abs(x)<0.71  then  y := x
                      else  y := sqrt(1.0-sqr(x));
      S := 0.0;
      R := y;
      l := 1;
      repeat
        S := S+R;
        R := R*sqr(y*l)/((l+1) * (l+2));
        l := l+2;
      until (abs(R)<1.0e-35) or
            (abs(R)<=1.0e-16*abs(S));
      S := S+R;
      if abs(x-y)>1.0e-12  then
        begin
          if x>0.0  then   S := Pi/2.0-S
                    else   S := -Pi/2.0+S;
        end;
      ArcSin := S;
    end;
end;

function  ArcCos(X: RealType): RealType;
begin
  ArcCos := Pi*0.5 - ArcSin(X);
end;

function  exp1(X: RealType): RealType;{ exp with underflow protection}
begin
  if X = 0.0 then exp1 := 1.0
             else if X > -ExpArg then  exp1 := exp(X)
                                 else  exp1 := 0.0;
end;


function sqroot1(x, y: RealType) : RealType;
{ safe sqrt: sqroot = sqrt(x + y) }
begin
  if x > y then
    result := sqrt(x)*sqrt(1.0 + y/x)
  else
    result := sqrt(y)*sqrt(1.0 + x/y);
end;

function sqroot2(x, y: RealType) : RealType;
{ safe sqrt: sqroot = sqrt(x*x+y*y) }
begin
  if abs(x) > abs(y) then
    result := abs(x)*sqrt(1.0+sqr(y/x))
  else
    result := abs(y)*sqrt(1.0+sqr(x/y));
end;

{ Hyperbolic Functions }

function ch(x: RealType): RealType;
begin
  ch := 0.5*(exp1(x)+exp1(-x));
end;

function sh(x: RealType): RealType;
begin
  sh := 0.5*(exp1(x)-exp1(-x));
end;

function th(x: RealType): RealType;
begin
  th := (1.0-exp1(-2.0*x))/(1.0+exp1(-2.0*x));
end;

function arch(x: RealType): RealType;
begin
  arch := ln(x + sqrt(sqr(x)-1.0) );
end;

function arsh(x: RealType): RealType;
begin
  arsh := ln(x + sqrt(sqr(x)+1.0) );
end;

function arth(x: RealType): RealType;
begin
  arth := 0.5 * ln((1.0+x)/(1.0-x));
end;

{ Special functions for the Real Calculations }

function er(x: Realtype): RealType;
const
  epsilon = 1.0E-14;
  sp = 1.7724538509055160;{ sqrt(pi) }
  a1 = 0.4613135;    b1 = 0.1901635;
  a2 = 0.09999216;   b2 = 1.7844927;
  a3 = 0.002883894;  b3 = 5.5253437;
var
  s, t, m, xx, xx2           : RealType;
begin
  xx := x*x;
  if abs(x) < 3.0
  then {  series expansion of er(x)   }
  begin
    m := 1.0; t   := 1.0;
    s := 0.0; xx2 := xx*2.0;
    repeat
      s := s + t; m := m + 2.0;
      t := t*xx2; t := t/m;
    until (abs(t) < epsilon);
    er := exp1(xx) - s*x*2.0/sp;
  end
  else{  asymptotics from Abramowitz }
    if x > 0.0 then
      er := x*(a1/(b1+xx)+a2/(b2+xx)+a3/(b3+xx))
    else
      er := exp1(xx)*2.0 +
            x*(a1/(b1+xx)+a2/(b2+xx)+a3/(b3+xx));
end;{ of er }

function erf(x: Realtype): RealType;
const
  epsilon = 1.0E-14;
  sp = 1.7724538509055160;{ sqrt(pi) }
  a1 = 0.4613135;    b1 = 0.1901635;
  a2 = 0.09999216;   b2 = 1.7844927;
  a3 = 0.002883894;  b3 = 5.5253437;
var
  s, t, m, xx, xx2           : RealType;
begin
  xx := x*x;
  if abs(x) < 3.0
  then {  series expansion of erf(x)   }
  begin
    m := 1.0; t   := 1.0;
    s := 0.0; xx2 := xx*2.0;
    repeat
      s := s + t;  m := m + 2.0;  t := t*xx2/m;
    until (abs(t)<epsilon);
    erf := exp1(-xx)*s*x*2.0/sp;
  end
  else{  asimptotics from Abramowitz }
    erf := sign(x) - exp1(-xx)*x*(a1/(b1+xx)+a2/(b2+xx)+a3/(b3+xx));
end;

function erfc(x: Realtype): RealType;
begin
  erfc := 1.0 - erf(x);
end;

function I0_(x: RealType): RealType;
var
  r               : RealType;
begin
  if x < 3.75 then
  begin
    r  := sqr(x/3.75);
    I0_:= 1.0+r*(3.5156229+r*(3.0899424+r*(1.2067492
             +r*(0.2659732+r*(0.0360768+r*0.0045813)))));
  end else
  begin
    r  := 3.75/x;
    I0_:=exp1(x)*(0.39894228+r*(0.01328592+r*(0.00225319
              +r*(-0.00157565+r*(0.00916281+r*(-0.02057706
              +r*(0.02635537+r*(-0.01647633+r*0.00392377))))))))/sqrt(x);
  end;
end;

function I1_(x: RealType): RealType;
var
  r               : RealType;
begin
  if x < 3.75 then
  begin
    r  := sqr(x/3.75);
    I1_:= (0.5+r*(0.87890594+r*(0.51498869+r*(0.15084934
              +r*(0.02658733+r*(0.00301532+r*0.00032411))))))*x;
  end else
  begin
    r  := 3.75/x;
    I1_:=exp1(x)*(0.39894228+r*(-0.03988024+r*(-0.00362018
              +r*(0.00163801+r*(-0.01031555+r*(0.02282967
              +r*(-0.02895312+r*(0.01787654-r*0.00420059))))))))/sqrt(x);
  end;
end;

function fk0_(x: RealType): RealType;
var
  r                : RealType;
begin
    r   := sqr(x/2.0);
    fk0_:= -0.57721566+r*(0.42278420+r*(0.23069756+r*(0.03488590
        +r*(0.00262698+r*(0.0001075+r*0.0000074)))));
end;

function K0_(x: RealType): RealType;
var
  r               : RealType;
begin
  if x < 2.0
  then
    K0_:= -ln(x/2.0)*I0_(x) + fk0_(x)
  else
  begin
    r := 2.0/x;
    K0_:=exp1(-x)*(1.25331414+r*(-0.07832358
              +r*(0.02189568+r*(-0.01062446+r*(0.00587872
              +r*(-0.00251540+r*0.00053208))))))/sqrt(x);
  end;
end;

function fk1_(x: RealType): RealType;
var
  r: RealType;
begin
  r  := sqr(x/2.0);
  fk1_:= (1.0+r*(0.15443144-r*(0.67278579+r*(0.18156897
            +r*(0.01919402+r*(0.00110404+r*0.00004686))))))/x;
end;

function K1_(x: RealType): RealType;
var
  r: RealType;
begin
  if x < 2.0
  then
    K1_:= ln(x/2.0)*I1_(x)+fk1_(x)
  else
  begin
    r  := 2.0/x;
    K1_:=exp1(-x)*(1.25331414+r*(0.23498619
               +r*(-0.03655620+r*(0.01504268+r*(-0.00780353
               +r*(0.00325614-r*0.00068245))))))/sqrt(x)
  end;
end;

function gamma_(x: Realtype) : RealType;
{ gamma-function  error < 1 e-6    }
const{  for expansion from Abramowitz }
  b1 = -0.577191652;    b2 = 0.988205891;
  b3 = -0.897056937;    b4 = 0.918206857;
  b5 = -0.756704078;    b6 = 0.482199394;
  b7 = -0.193527818;    b8 = 0.035868343;
var
  t               : RealType;
begin
  t := 1.0;
  while x >= 2.0 do begin x := x - 1.0; t := t*x; end;
  while x < 1.0  do begin t := t/x; x := x + 1.0; end;
  x := x - 1.0;
  gamma_:=t*((((((((b8*x+b7)*x+b6)*x+b5)*x+b4)*x+b3)*x+b2)*x+b1)*x+1.0);
end;

function gamma_in(a, x :Realtype): RealType;
const
  epsilon = 1.0E-8;
  Xmax    = 10.0;
var
  s, t, m: RealType;
begin
  if x < Xmax then {  series expansion  }
  begin
    m := 0.0;  t := exp1(a*ln(x))/a;  s := t;
    repeat
      m := m + 1.0;  t := -t*x*a/(a+1.0)/m;
      s := s + t;    a := a + 1.0;
    until (abs(t)<epsilon);
    gamma_in := s;
  end
  else{  asimptotics from Abramowitz }
  begin
    m := a - 1.0;  t := exp1(m*ln(x)-x);  s := t;
    repeat
      t := t*m/x; s := s + t;  m := m - 1.0;
    until (abs(t)<epsilon) or (abs(a) > x);
    gamma_in := gamma_(a) - s;
  end;
end;

function E1_(x: Realtype)    : RealType;
const
  epsilon = 1.0E-9;
  Xmax    = 10.0;
  gamma   = 0.57721566490153286;{ Euler's constant }
var
  s, t, m           : RealType;
begin
  if x < Xmax then {  series expansion  }
  begin
    m := 1.0;  t := x;  s := t - gamma - ln(x);
    repeat
      t := -t*x*m/sqr(m+1.0);
      s := s + t;    m := m + 1.0;
    until (abs(t)<epsilon);
    E1_ := s;
  end
  else{  asymptotics from Abramowitz }
  begin
    m := 1.0;  t := exp1(-x)/x;  s := t;
    repeat
      t := -t*m/x; s := s + t;  m := m + 1.0;
    until (abs(t)<epsilon) or (m > x);
    E1_ := s;
  end;
end;

function I_nu(nu, x: RealType): RealType;
const
  epsilon = 1.0E-7;
  Xmax    = 10.0;
  s2pi    = 2.506628274631;{ sqrt(2*pi) }
var
  s, t, m, f, n, x8: RealType;
begin
  if x < Xmax then {  series expansion  }
  begin
    x  := x/2.0;     t := exp1(nu*ln(x));  x := sqr(x);
    nu := nu + 1.0;  t := t/gamma_(nu);    m := 1.0;
    s  := t;
    repeat
      t  := t*x/m/nu;  s := s + t;
      nu := nu +1.0;   m := m + 1.0;
    until (abs(t)<epsilon);
    I_nu := s;
  end
  else{  asimptotics from Abramowitz }
  begin