freepascal.compiler/packages/numlib/spe.pas

{
    $Id$
    This file is part of the Numlib package.
    Copyright (c) 1986-2000 by
     Kees van Ginneken, Wil Kortsmit and Loek van Reij of the
     Computational centre of the Eindhoven University of Technology

    FPC port Code          by Marco van de Voort ([email protected])
             documentation by Michael van Canneyt ([email protected])

    Special functions. (Currently, most of these functions only work for
            ArbFloat=REAL/DOUBLE, not for Extended(at least not at full
            precision, implement the tables in the diverse functions
            for extended)

    See the file COPYING.FPC, included in this distribution,
    for details about the copyright.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 **********************************************************************}
unit spe;
{$I DIRECT.INC}

interface

uses typ;

{  Calculate modified Besselfunction "of the first kind" I0(x) }
function spebi0(x: ArbFloat): ArbFloat;

{  Calculate modified Besselfunction "of the first kind" I1(x) }
function spebi1(x: ArbFloat): ArbFloat;

{  Calculate Besselfunction "of the first kind" J0(x) }
function spebj0(x: ArbFloat): ArbFloat;

{  Calculate Besselfunction "of the first kind" J1(x) }
function spebj1(x: ArbFloat): ArbFloat;

{  Calculate modified Besselfunction "of the second kind" K0(x) }
function spebk0(x: ArbFloat): ArbFloat;

{  Calculate modified Besselfunction "of the second kind" K1(x) }
function spebk1(x: ArbFloat): ArbFloat;

{  Calculate Besselfunction "of the second kind" Y0(x) }
function speby0(x: ArbFloat): ArbFloat;

{  Calculate Besselfunction "of the second kind" Y1(x) }
function speby1(x: ArbFloat): ArbFloat;

{  Entier function, calculates first integer greater or equal than X}
function speent(x: ArbFloat): longint;

{  Errorfunction ( 2/sqrt(pi)* Int(t,0,pi,exp(sqr(t)) )}
function speerf(x: ArbFloat): ArbFloat;

{  Errorfunction's complement ( 2/sqrt(pi)* Int(t,pi,inf,exp(sqr(t)) )}
function speefc(x: ArbFloat): ArbFloat;

{  Function to calculate the Gamma function ( int(t,0,inf,t^(x-1)*exp(-t)) }
function spegam(x: ArbFloat): ArbFloat;

{  Function to calculate the natural logaritm of the Gamma function}
function spelga(x: ArbFloat): ArbFloat;

{  "Calculates" the maximum of two ArbFloat values     }
function spemax(a, b: ArbFloat): ArbFloat;

{  Calculates the functionvalue of a polynomalfunction with n coefficients in a
for variable X }
function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;

{ Calc a^b with a and b real numbers}
function spepow(a, b: ArbFloat): ArbFloat;

{ Returns sign of x (-1 for x<0, 0 for x=0 and 1 for x>0)  }
function spesgn(x: ArbFloat): ArbInt;

{  ArcSin(x) }
function spears(x: ArbFloat): ArbFloat;

{  ArcCos(x) }
function spearc(x: ArbFloat): ArbFloat;

{  Sinh(x) }
function spesih(x: ArbFloat): ArbFloat;

{  Cosh(x) }
function specoh(x: ArbFloat): ArbFloat;

{  Tanh(x) }
function spetah(x: ArbFloat): ArbFloat;

{  ArcSinH(x) }
function speash(x: ArbFloat): ArbFloat;

{  ArcCosh(x) }
function speach(x: ArbFloat): ArbFloat;

{  ArcTanH(x) }
function speath(x: ArbFloat): ArbFloat;

implementation

function spebi0(x: ArbFloat): ArbFloat;

const

     xvsmall = 3.2e-9;
          a1 : array[0..23] of ArbFloat =
               (  3.08508322553671039e-1, -1.86478066609466760e-1,
                  1.57686843969995904e-1, -1.28895621330524993e-1,
                  9.41616340200868389e-2, -6.04316795007737183e-2,
                  3.41505388391452157e-2, -1.71317947935716536e-2,
                  7.70061052263382555e-3, -3.12923286656374358e-3,
                  1.15888319775791686e-3, -3.93934532072526720e-4,
                  1.23682594989692688e-4, -3.60645571444886286e-5,
                  9.81395862769787105e-6, -2.50298975966588680e-6,
                  6.00566861079330132e-7, -1.36042013507151017e-7,
                  2.92096163521178835e-8, -5.94856273204259507e-9,
                  1.13415934215369209e-9, -2.10071360134551962e-10,
                  4.44484446637868974e-11,-7.48150165756234957e-12);
          a2 : array[0..26] of ArbFloat =
               (  1.43431781856850311e-1, -3.71571542566085323e-2,
                  1.44861237337359455e-2, -6.30121694459896307e-3,
                  2.89362046530968701e-3, -1.37638906941232170e-3,
                  6.72508592273773611e-4, -3.35833513200679384e-4,
                  1.70524543267970595e-4, -8.74354291104467762e-5,
                  4.48739019580173804e-5, -2.28278155280668483e-5,
                  1.14032404021741277e-5, -5.54917762110482949e-6,
                  2.61457634142262604e-6, -1.18752840689765504e-6,
                  5.18632519069546106e-7, -2.17653548816447667e-7,
                  8.75291839187305722e-8, -3.34900221934314738e-8,
                  1.24131668344616429e-8, -4.66215489983794905e-9,
                  1.58599776268172290e-9, -3.80370174256271589e-10,
                  1.23188158175419302e-10,-8.46900307934754898e-11,
                  2.45185252963941089e-11);
           a3: array[0..19] of ArbFloat =
               (  4.01071065066847416e-1,  2.18216817211694382e-3,
                  5.59848253337377763e-5,  2.79770701849785597e-6,
                  2.17160501061222148e-7,  2.36884434055843528e-8,
                  3.44345025431425567e-9,  6.47994117793472057e-10,
                  1.56147127476528831e-10, 4.82726630988879388e-11,
                  1.89599322920800794e-11, 1.05863621425699789e-11,
                  8.27719401266046976e-12, 2.82807056475555021e-12,
                 -4.34624739357691085e-12,-4.29417106720584499e-12,
                  4.30812577328136192e-13, 1.44572313799118029e-12,
                  4.73229306831831040e-14,-1.95679809047625728e-13);


var t : ArbFloat;

begin
  t:=abs(x);
  if t <=xvsmall
  then
    spebi0:=1
  else
  if t <= 4
  then
    spebi0 := exp(t)*spepol(t/2-1, a1[0], SizeOf(a1) div SizeOf(ArbFloat) -1)
  else
  if t <= 12
  then
    spebi0:=exp(t)*spepol(t/4-2, a2[0], SizeOf(a2) div SizeOf(ArbFloat) -1)
  else { t > 12}
    spebi0:=(exp(t)/sqrt(t))*
            spepol(24/t-1, a3[0], SizeOf(a3) div SizeOf(ArbFloat) -1)
end; {spebi0}

function spebi1(x: ArbFloat): ArbFloat;


const xvsmall = 3.2e-9;
      a1: array[0..11] of ArbFloat =
      ( 1.19741654963670236e+0, 9.28758890114609554e-1,
        2.68657659522092832e-1, 4.09286371827770484e-2,
        3.84763940423809498e-3, 2.45224314039278904e-4,
        1.12849795779951847e-5, 3.92368710996392755e-7,
        1.06662712314503955e-8, 2.32856921884663846e-10,
        4.17372709788222413e-12,6.24387910353848320e-14);

      a2: array[0..26] of ArbFloat =
      ( 1.34142493292698178e-1, -2.99140923897405570e-2,
        9.76021102528646704e-3, -3.40759647928956354e-3,
        1.17313412855965374e-3, -3.67626180992174570e-4,
        8.47999438119288094e-5,  5.21557319070236939e-6,
       -2.62051678511418163e-5,  2.47493270133518925e-5,
       -1.79026222757948636e-5,  1.13818992442463952e-5,
       -6.63144162982509821e-6,  3.60186151617732531e-6,
       -1.83910206626348772e-6,  8.86951515545183908e-7,
       -4.05456611578551130e-7,  1.76305222240064495e-7,
       -7.28978293484163628e-8,  2.84961041291017650e-8,
       -1.07563514207617768e-8,  4.11321223904934809e-9,
       -1.41575617446629553e-9,  3.38883570696523350e-10,
       -1.10970391104678003e-10, 7.79929176497056645e-11,
       -2.27061376122617856e-11);

       a3: array[0..19] of ArbFloat =
       ( 3.92624494204116555e-1, -6.40545360348237412e-3,
        -9.12475535508497109e-5, -3.82795135453556215e-6,
        -2.72684545741400871e-7, -2.82537120880041703e-8,
        -3.96757162863209348e-9, -7.28107961041827952e-10,
        -1.72060490748583241e-10,-5.23524129533553498e-11,
        -2.02947854602758139e-11,-1.11795516742222899e-11,
        -8.69631766630563635e-12,-3.05957293450420224e-12,
         4.42966462319664333e-12, 4.47735589657057690e-12,
        -3.95353303949377536e-13,-1.48765082315961139e-12,
        -5.77176811730370560e-14, 1.99448557598015488e-13);

var t : ArbFloat;

begin
  t:=abs(x);
  if t <= xvsmall
  then
    spebi1:=x/2
  else
  if t <= 4
  then
    spebi1:=x*spepol(sqr(t)/8-1, a1[0], sizeof(a1) div sizeof(ArbFloat)-1)
  else
  if t <= 12
  then
    spebi1:=
      exp(t)*spepol(t/4-2, a2[0], sizeof(a2) div sizeof(ArbFloat)-1)*spesgn(x)
  else { t > 12}
    spebi1:=
      (exp(t)/sqrt(t))*
      spepol(24/t-1, a3[0], sizeof(a3) div sizeof(ArbFloat)-1)*spesgn(x)
end; {spebi1}

function spebj0(x: ArbFloat): ArbFloat;
const

       xvsmall = 3.2e-9;
          tbpi = 6.36619772367581343e-1;
           a1 : array[0..5] of ArbFloat =
           ( 1.22200000000000000e-17, -1.94383469000000000e-12,
             7.60816359241900000e-8,  -4.60626166206275050e-4,
             1.58067102332097261e-1,  -8.72344235285222129e-3);

            b1 : array[0..5] of ArbFloat =
            ( - 7.58850000000000000e-16, 7.84869631400000000e-11,
              - 1.76194690776215000e-6,  4.81918006946760450e-3,
              - 3.70094993872649779e-1,  1.57727971474890120e-1);

            c1 : array[0..4] of ArbFloat =
            ( 4.12532100000000000e-14, - 2.67925353056000000e-9,
              3.24603288210050800e-5,  - 3.48937694114088852e-2,
              2.65178613203336810e-1);

            d1 : array[0..13] of ArbFloat =
            ( 9.99457275788251954e-1, -5.36367319213004570e-4,
              6.13741608010926000e-6, -2.05274481565160000e-7,
              1.28037614434400000e-8, -1.21211819632000000e-9,
              1.55005642880000000e-10,-2.48827276800000000e-11,
              4.78702080000000000e-12,-1.06365696000000000e-12,
              2.45294080000000000e-13,-6.41843200000000000e-14,
              3.34028800000000000e-14,-1.17964800000000000e-14);

             d2 : array[0..16] of ArbFloat =
             ( -1.55551138795135187e-2,  6.83314909934390000e-5,
               -1.47713883264594000e-6,  7.10621485930000000e-8,
               -5.66871613024000000e-9,  6.43278173280000000e-10,
               -9.47034774400000000e-11, 1.70330918400000000e-11,
               -3.59094272000000000e-12, 8.59855360000000000e-13,
               -2.28807680000000000e-13, 6.95193600000000000e-14,
               -2.27942400000000000e-14, 4.75136000000000000e-15,
               -1.14688000000000000e-15, 2.12992000000000000e-15,
               -9.83040000000000000e-16);

var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
    i, bov                       : ArbInt;

begin
  t:=abs(x);
  if t<=xvsmall
  then
    spebj0:=1
  else
  if t<=8
  then
    begin
      t:=0.03125*sqr(t)-1; t2:=2*t;
      b:=0; c:=0;
      bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
      for i:=0 to bov do
        begin
          a:=t2*c-b+a1[i];
          if i<5
          then
            b:=t2*a-c+b1[i]
          else
            spebj0:=t*a-c+b1[i];
          if i<bov
          then
            c:=t2*b-a+c1[i]
          else
            if i<5
            then
              spebj0:=t*b-a+c1[i]
        end {i}
    end {abs(x)<=8}
  else
    begin
      g:=t-1/(2*tbpi); y:=sqrt(tbpi/t);
      cx:=cos(g)*y; sx:=-sin(g)*y*8/t;
      t:=128/sqr(t)-1;
      spebj0:=cx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
              + sx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
    end {abs(x)>8}
    
end {spebj0};

function spebj1(x: ArbFloat): ArbFloat;
const

       xvsmall = 3.2e-9;
          tbpi = 6.36619772367581343e-1;
      a1 : array[0..5] of ArbFloat =
      ( 2.95000000000000000e-18, -5.77740420000000000e-13,
        2.94970700727800000e-8,  -2.60444389348580680e-4,
        1.77709117239728283e-1,  -1.19180116054121687e+0);

      b1 : array[0..5] of ArbFloat =
      ( -1.95540000000000000e-16, 2.52812366400000000e-11,
        -7.61758780540030000e-7,  3.24027018268385747e-3,
        -6.61443934134543253e-1,  6.48358770605264921e-1);

      c1 : array[0..4] of ArbFloat =
      ( 1.13857200000000000e-14, -9.42421298160000000e-10,
        1.58870192399321300e-5,  -2.91755248061542077e-2,
        1.28799409885767762e+0);

       d1 : array[0..13] of ArbFloat =
       ( 1.00090702627808217e+0,  8.98804941670557880e-4,
        -7.95969469843846000e-6,  2.45367662227560000e-7,
        -1.47085129889600000e-8,  1.36030580128000000e-9,
        -1.71310758400000000e-10, 2.72040729600000000e-11,
        -5.19113984000000000e-12, 1.14622464000000000e-12,
        -2.63372800000000000e-13, 6.86387200000000000e-14,
        -3.54508800000000000e-14, 1.24928000000000000e-14);

       d2 : array[0..15] of ArbFloat =
       ( 4.67768740274489776e-2,  -9.62145882205441600e-5,
         1.82120185123076000e-6,  -8.29196070929200000e-8,
         6.42013250344000000e-9,  -7.15110504800000000e-10,
         1.03950931840000000e-10, -1.85248000000000000e-11,
         3.87554432000000000e-12, -9.23228160000000000e-13,
         2.50224640000000000e-13, -7.43936000000000000e-14,
         1.75718400000000000e-14, -4.83328000000000000e-15,
         5.32480000000000000e-15, -2.29376000000000000e-15);

var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
    i, bov                       : ArbInt;

begin
  t:=abs(x);
  if t<xvsmall
  then
    spebj1:=x/2
  else
  if t<=8
  then
    begin
      t:=0.03125*sqr(t)-1; t2:=2*t;
      b:=0; c:=0;
      bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
      for i:=0 to bov do
        begin
          a:=t2*c-b+a1[i];
          if i<bov
          then
            begin
              b:=t2*a-c+b1[i];
              c:=t2*b-a+c1[i]
            end
          else
            spebj1:=(x/8)*(t*a-c+b1[i])
        end {i}
    end {abs(x)<=8}
  else
    begin
      g:=t-1.5/tbpi; y:=sqrt(tbpi/t)*spesgn(x);
      cx:=cos(g)*y; sx:=-sin(g)*y*8/t;
      t:=128/sqr(t)-1;
      spebj1:=cx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
              + sx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
    end {abs(x)>8}
end {spebj1};

function spebk0(x: ArbFloat): ArbFloat;

const

     egam = 0.57721566490153286;
     xvsmall = 3.2e-9;
     highexp = 745;

      a0: array[0..7] of ArbFloat =
      ( 1.12896092945412762e+0,  1.32976966478338191e-1,
        4.07157485171389048e-3,  5.59702338227915383e-5,
        4.34562671546158210e-7,  2.16382411824721532e-9,
        7.49110736894134794e-12, 1.90674197514561280e-14);

      a1: array[0..8] of ArbFloat =
      ( 2.61841879258687055e-1,  1.52436921799395196e-1,
        6.63513979313943827e-3,  1.09534292632401542e-4,
        9.57878493265929443e-7,  5.19906865800665633e-9,
        1.92405264219706684e-11, 5.16867886946332160e-14,
        1.05407718191360000e-16);

      a2: array[0..22] of ArbFloat =
      ( 9.58210053294896496e-1, -1.42477910128828254e-1,
        3.23582010649653009e-2, -8.27780350351692662e-3,
        2.24709729617770471e-3, -6.32678357460594866e-4,
        1.82652460089342789e-4, -5.37101208898441760e-5,
        1.60185974149720562e-5, -4.83134250336922161e-6,
        1.47055796078231691e-6, -4.51017292375200017e-7,
        1.39217270224614153e-7, -4.32185089841834127e-8,
        1.34790467361340101e-8, -4.20597329258249948e-9,
        1.32069362385968867e-9, -4.33326665618780914e-10,
        1.37999268074442719e-10, -3.19241059198852137e-11,
        9.74410152270679245e-12, -7.83738609108569293e-12,
        2.57466288575820595e-12);

      a3: array[0..22] of ArbFloat =
     ( 6.97761598043851776e-1, -1.08801882084935132e-1,
       2.56253646031960321e-2, -6.74459607940169198e-3,
       1.87292939725962385e-3, -5.37145622971910027e-4,
       1.57451516235860573e-4, -4.68936653814896712e-5,
       1.41376509343622727e-5, -4.30373871727268511e-6,
       1.32052261058932425e-6, -4.07851207862189007e-7,
       1.26672629417567360e-7, -3.95403255713518420e-8,
       1.23923137898346852e-8, -3.88349705250555658e-9,
       1.22424982779432970e-9, -4.03424607871960089e-10,
       1.28905587479980147e-10,-2.97787564633235128e-11,
       9.11109430833001267e-12,-7.39672783987933184e-12,
       2.43538242247537459e-12);
      a4: array[0..16] of ArbFloat =
      ( 1.23688664769425422e+0,  -1.72683652385321641e-2,
       -9.25551464765637133e-4,  -9.02553345187404564e-5,
       -6.31692398333746470e-6,  -7.69177622529272933e-7,
       -4.16044811174114579e-8,  -9.41555321137176073e-9,
        1.75359321273580603e-10, -2.22829582288833265e-10,
        3.49564293256545992e-11, -1.11391758572647639e-11,
        2.85481235167705907e-12, -7.31344482663931904e-13,
        2.06328892562554880e-13, -1.28108310826991616e-13,
        4.43741979886551040e-14);


var t: ArbFloat;

begin
  if x<=0
  then
    RunError(401);
  if x<=xvsmall
  then
    spebk0:=-(ln(x/2)+egam)
  else
  if x<=1
  then
    begin
      t:=2*sqr(x)-1;
      spebk0:=-ln(x)*spepol(t, a0[0], sizeof(a0) div sizeof(ArbFloat) - 1)
              + spepol(t, a1[0], sizeof(a1) div sizeof(ArbFloat) - 1)
    end
  else
  if x<=2
  then
    spebk0:=exp(-x)*spepol(2*x-3, a2[0], sizeof(a2) div sizeof(ArbFloat) - 1)
  else
  if x<=4
  then
    spebk0:=exp(-x)*spepol(x-3, a3[0], sizeof(a3) div sizeof(ArbFloat) - 1)
  else
  if x <= highexp
  then
    spebk0:=exp(-x)*
            spepol(10/(1+x)-1, a4[0], sizeof(a4) div sizeof(ArbFloat) - 1)/sqrt(x)
  else
    spebk0:=0
end; {spebk0}

function spebk1(x: ArbFloat): ArbFloat;

const

   xsmall = 7.9e-10;
  highexp = 745;
   a0: array[0..7] of ArbFloat =
   ( 5.31907865913352762e-1,  3.25725988137110495e-2,
     6.71642805873498653e-4,  6.95300274548206237e-6,
     4.32764823642997753e-8,  1.79784792380155752e-10,
     5.33888268665658944e-13, 1.18964962439910400e-15);

   a1: array[0..7] of ArbFloat =
   ( 3.51825828289325536e-1,  4.50490442966943726e-2,
     1.20333585658219028e-3,  1.44612432533006139e-5,
     9.96686689273781531e-8,  4.46828628435618679e-10,
     1.40917103024514301e-12, 3.29881058019865600e-15);

   a2: array[0..23] of ArbFloat =
   ( 1.24316587355255299e+0, -2.71910714388689413e-1,
     8.20250220860693888e-2, -2.62545818729427417e-2,
     8.57388087067410089e-3, -2.82450787841655951e-3,
     9.34594154387642940e-4, -3.10007681013626626e-4,
     1.02982746700060730e-4, -3.42424912211942134e-5,
     1.13930169202553526e-5, -3.79227698821142908e-6,
     1.26265578331941923e-6, -4.20507152338934956e-7,
     1.40138351985185509e-7, -4.66928912168020101e-8,
     1.54456653909012693e-8, -5.13783508140332214e-9,
     1.82808381381205361e-9, -6.15211416898895086e-10,
     1.28044023949946257e-10, -4.02591066627023831e-11,
     4.27404330568767242e-11, -1.46639291782948454e-11);

   a3: array[0..23] of ArbFloat =
   ( 8.06563480128786903e-1,  -1.60052611291327173e-1,
     4.58591528414023064e-2,  -1.42363136684423646e-2,
     4.55865751206724687e-3,  -1.48185472032688523e-3,
     4.85707174778663652e-4,  -1.59994873621599146e-4,
     5.28712919123131781e-5,  -1.75089594354079944e-5,
     5.80692311842296724e-6,  -1.92794586996432593e-6,
     6.40581814037398274e-7,  -2.12969229346310343e-7,
     7.08723366696569880e-8,  -2.35855618461025265e-8,
     7.79421651144832709e-9,  -2.59039399308009059e-9,
     9.20781685906110546e-10, -3.09667392343245062e-10,
     6.44913423545894175e-11, -2.02680401514735862e-11,
     2.14736751065133220e-11, -7.36478297050421658e-12);

    a4: array[0..16] of ArbFloat =
    ( 1.30387573604230402e+0,   5.44845254318931612e-2,
      4.31639434283445364e-3,   4.29973970898766831e-4,
      4.04720631528495020e-5,   4.32776409784235211e-6,
      4.07563856931843484e-7,   4.86651420008153956e-8,
      3.82717692121438315e-9,   6.77688943857588882e-10,
      6.97075379117731379e-12,  1.72026097285930936e-11,
     -2.60774502020271104e-12,  8.58211523713560576e-13,
     -2.19287104441802752e-13,  1.39321122940600320e-13,
     -4.77850238111580160e-14);

var t: ArbFloat;

begin
  if x<=0
  then
    RunError(402);
  if x<=xsmall
  then
    spebk1:=1/x
  else
  if x<=1
  then
    begin
      t:=2*sqr(x)-1;
      spebk1:=( ln(x)*spepol(t, a0[0], sizeof(a0) div sizeof(ArbFloat) - 1)
              -spepol(t, a1[0], sizeof(a1) div sizeof(ArbFloat) -1) )*x + 1/x
    end
  else
  if x<=2
  then
    spebk1:=exp(-x)*spepol(2*x-3, a2[0], sizeof(a2) div sizeof(ArbFloat) - 1)
  else
  if x<=4
  then
    spebk1:=exp(-x)*spepol(x-3, a3[0], sizeof(a3) div sizeof(ArbFloat) - 1)
  else
  if x <= highexp
  then
    spebk1:=exp(-x)*spepol(10/(1+x)-1, a4[0],
            sizeof(a4) div sizeof(ArbFloat) - 1)/sqrt(x)
  else
    spebk1:=0
end; {spebk1}

function speby0(x: ArbFloat): ArbFloat;

const

      tbpi = 6.36619772367581343e-1;
      egam = 5.77215664901532861e-1;
   xvsmall = 3.2e-9;
   a1 : array[0..5] of ArbFloat =
   ( 3.90000000000000000e-19, -8.74734100000000000e-14,
     5.24879478733000000e-9,  -5.63207914105698700e-5,
     4.71966895957633869e-2,   1.79034314077182663e-1);

   b1 : array[0..5] of ArbFloat =
   ( -2.69800000000000000e-17, 4.02633082000000000e-12,
     -1.44072332740190000e-7,  7.53113593257774230e-4,
     -1.77302012781143582e-1, -2.74474305529745265e-1);

   c1 : array[0..5] of ArbFloat =
   ( 1.64349000000000000e-15, -1.58375525420000000e-10,
     3.20653253765480000e-6,  -7.28796247955207918e-3,
     2.61567346255046637e-1,  -3.31461132032849417e-2);

    d1 : array[0..13] of ArbFloat =
    ( 9.99457275788251954e-1, -5.36367319213004570e-4,
      6.13741608010926000e-6, -2.05274481565160000e-7,
      1.28037614434400000e-8, -1.21211819632000000e-9,
      1.55005642880000000e-10,-2.48827276800000000e-11,
      4.78702080000000000e-12,-1.06365696000000000e-12,
      2.45294080000000000e-13,-6.41843200000000000e-14,
      3.34028800000000000e-14,-1.17964800000000000e-14);

    d2 : array[0..16] of ArbFloat =
    (-1.55551138795135187e-2,  6.83314909934390000e-5,
     -1.47713883264594000e-6,  7.10621485930000000e-8,
     -5.66871613024000000e-9,  6.43278173280000000e-10,
     -9.47034774400000000e-11, 1.70330918400000000e-11,
     -3.59094272000000000e-12, 8.59855360000000000e-13,
     -2.28807680000000000e-13, 6.95193600000000000e-14,
     -2.27942400000000000e-14, 4.75136000000000000e-15,
     -1.14688000000000000e-15, 2.12992000000000000e-15,
     -9.83040000000000000e-16);

var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
    i, bov                       : ArbInt;

begin
  if x<=0
  then
    RunError(403);
  if x<=xvsmall
  then
    speby0:=(ln(x/2)+egam)*tbpi
  else
  if x<=8
  then
    begin
      t:=0.03125*sqr(x)-1; t2:=2*t;
      b:=0; c:=0;
      bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
      for i:=0 to bov do
        begin
          a:=t2*c-b+a1[i];
          b:=t2*a-c+b1[i];
          if i<bov
          then
            c:=t2*b-a+c1[i]
          else
            speby0:=t*b-a+c1[i]+tbpi*spebj0(x)*ln(x)
        end {i}
    end {x<=8}
  else
    begin
      g:=x-1/(2*tbpi); y:=sqrt(tbpi/x);
      cx:=cos(g)*y*8/x; sx:=sin(g)*y;
      t:=128/sqr(x)-1;
      speby0:=sx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
            + cx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
    end {x>8}
end {speby0};

function speby1(x: ArbFloat): ArbFloat;

const
    tbpi = 6.36619772367581343e-1;
    xsmall = 7.9e-10;
   a1 : array[0..5] of ArbFloat =
   (-6.58000000000000000e-18, 1.21143321000000000e-12,
    -5.68844003991900000e-8,  4.40478629867099510e-4,
    -2.26624991556754924e-1, -1.28697384381350001e-1);

   b1 : array[0..5] of ArbFloat =
   ( 4.27730000000000000e-16,-5.17212147300000000e-11,
     1.41662436449235000e-6, -5.13164116106108479e-3,
     6.75615780772187667e-1,  2.03041058859342538e-2);

   c1 : array[0..4] of ArbFloat =
   (-2.44094900000000000e-14, 1.87547032473000000e-9,
    -2.83046401495148000e-5,  4.23191803533369041e-2,
    -7.67296362886645940e-1);
    
   d1 : array[0..13] of ArbFloat =
   ( 1.00090702627808217e+0,  8.98804941670557880e-4,
    -7.95969469843846000e-6,  2.45367662227560000e-7,
    -1.47085129889600000e-8,  1.36030580128000000e-9,
    -1.71310758400000000e-10, 2.72040729600000000e-11,
    -5.19113984000000000e-12, 1.14622464000000000e-12,
    -2.63372800000000000e-13, 6.86387200000000000e-14,
    -3.54508800000000000e-14, 1.24928000000000000e-14);

    d2 : array[0..15] of ArbFloat =
    ( 4.67768740274489776e-2, -9.62145882205441600e-5,
      1.82120185123076000e-6, -8.29196070929200000e-8,
      6.42013250344000000e-9, -7.15110504800000000e-10,
      1.03950931840000000e-10,-1.85248000000000000e-11,
      3.87554432000000000e-12,-9.23228160000000000e-13,
      2.50224640000000000e-13,-7.43936000000000000e-14,
      1.75718400000000000e-14,-4.83328000000000000e-15,
      5.32480000000000000e-15,-2.29376000000000000e-15);

var t, g, y, t2, a, b, c, cx, sx : ArbFloat;
    i, bov                       : ArbInt;

begin
  if x<=0
  then
    RunError(404);
  if x<=xsmall
  then
    speby1:=-tbpi/x
  else
  if x<=8
  then
    begin
      t:=0.03125*sqr(x)-1; t2:=2*t;
      b:=0; c:=0;
      bov:=sizeof(a1) div sizeof(ArbFloat) - 1;
      for i:=0 to bov do
        begin
          a:=t2*c-b+a1[i];
          if i<bov
          then
            begin
              b:=t2*a-c+b1[i];
              c:=t2*b-a+c1[i]
            end
          else
          if bov=3   {single}
          then
            begin
              b:=t2*a-c+b1[i];
              speby1:=(t*b-a+c1[i])*x/8 + spebj1(x)*ln(x)*tbpi - tbpi/x
            end
          else
            speby1:=(t*a-c+b1[i])*x/8 + spebj1(x)*ln(x)*tbpi - tbpi/x
        end {i}
    end {x<=8}
  else
    begin
      g:=x-3/(2*tbpi); y:=sqrt(tbpi/x);
      cx:=cos(g)*y*8/x; sx:=sin(g)*y;
      t:=128/sqr(x)-1;
      speby1:=sx*spepol(t, d1[0], sizeof(d1) div sizeof(ArbFloat) - 1)
            + cx*spepol(t, d2[0], sizeof(d2) div sizeof(ArbFloat) - 1)
    end {x>8}
end {speby1};

function speent(x : ArbFloat): longint;

var tx : longint;

begin
  tx:=trunc(x);
  if x>=0
  then
    speent:=tx
  else
    if x=tx
    then
      speent:=tx
    else
      speent:=tx-1
end; {speent}

function speerf(x : ArbFloat) : ArbFloat;
const

        xup = 6.25;
     sqrtpi = 1.7724538509055160;
     c : array[1..18] of ArbFloat =
     ( 1.9449071068178803e0,  4.20186582324414e-2, -1.86866103976769e-2,
       5.1281061839107e-3,   -1.0683107461726e-3,   1.744737872522e-4,
      -2.15642065714e-5,      1.7282657974e-6,     -2.00479241e-8,
      -1.64782105e-8,         2.0008475e-9,         2.57716e-11,
      -3.06343e-11,           1.9158e-12,           3.703e-13,
      -5.43e-14,             -4.0e-15,              1.2e-15);

     d : array[1..17] of ArbFloat =
     ( 1.4831105640848036e0, -3.010710733865950e-1, 6.89948306898316e-2,
      -1.39162712647222e-2,   2.4207995224335e-3,  -3.658639685849e-4,
       4.86209844323e-5,     -5.7492565580e-6,      6.113243578e-7,
      -5.89910153e-8,         5.2070091e-9,        -4.232976e-10,
       3.18811e-11,          -2.2361e-12,           1.467e-13,
      -9.0e-15,               5.0e-16);

  var t, s, s1, s2, x2: ArbFloat;
         bovc, bovd, j: ArbInt;
begin
  bovc:=sizeof(c) div sizeof(ArbFloat);
  bovd:=sizeof(d) div sizeof(ArbFloat);
  t:=abs(x);
  if t <= 2
  then
    begin
      x2:=sqr(x)-2;
      s1:=d[bovd]; s2:=0; j:=bovd-1;
      s:=x2*s1-s2+d[j];
      while j > 1 do
        begin
          s2:=s1; s1:=s; j:=j-1;
          s:=x2*s1-s2+d[j]
        end;
      speerf:=(s-s2)*x/2
    end
  else
    if t < xup
    then
      begin
        x2:=2-20/(t+3);
        s1:=c[bovc]; s2:=0; j:=bovc-1;
        s:=x2*s1-s2+c[j];
        while j > 1 do
          begin
            s2:=s1; s1:=s; j:=j-1;
            s:=x2*s1-s2+c[j]
          end;
        x2:=((s-s2)/(2*t))*exp(-sqr(x))/sqrtpi;
        speerf:=(1-x2)*spesgn(x)
      end
    else
      speerf:=spesgn(x)
end;  {speerf}

function spemax(a, b: ArbFloat): ArbFloat;
begin
  if a>b
  then
    spemax:=a
  else
    spemax:=b
end; {spemax}

function speefc(x : ArbFloat) : ArbFloat;
const

   xlow = -6.25;
  xhigh = 27.28;
      c : array[0..22] of ArbFloat =
      ( 1.455897212750385e-1, -2.734219314954260e-1,
        2.260080669166197e-1, -1.635718955239687e-1,
        1.026043120322792e-1, -5.480232669380236e-2,
        2.414322397093253e-2, -8.220621168415435e-3,
        1.802962431316418e-3, -2.553523453642242e-5,
       -1.524627476123466e-4,  4.789838226695987e-5,
        3.584014089915968e-6, -6.182369348098529e-6,
        7.478317101785790e-7,  6.575825478226343e-7,
       -1.822565715362025e-7, -7.043998994397452e-8,
        3.026547320064576e-8,  7.532536116142436e-9,
       -4.066088879757269e-9, -5.718639670776992e-10,
        3.328130055126039e-10);

  var t, s : ArbFloat;
begin
  if x <= xlow
  then
    speefc:=2
  else
  if x >= xhigh
  then
    speefc:=0
  else
    begin
      t:=1-7.5/(abs(x)+3.75);
      s:=exp(-x*x)*spepol(t, c[0], sizeof(c) div sizeof(ArbFloat) - 1);
      if x < 0
      then
        speefc:=2-s
      else
        speefc:=s
    end
end {speefc};

function spegam(x: ArbFloat): ArbFloat;
const

    tmax = 170;
    a: array[0..23] of ArbFloat =
    ( 8.86226925452758013e-1,  1.61691987244425092e-2,
      1.03703363422075456e-1, -1.34118505705967765e-2,
      9.04033494028101968e-3, -2.42259538436268176e-3,
      9.15785997288933120e-4, -2.96890121633200000e-4,
      1.00928148823365120e-4, -3.36375833240268800e-5,
      1.12524642975590400e-5, -3.75499034136576000e-6,
      1.25281466396672000e-6, -4.17808776355840000e-7,
      1.39383522590720000e-7, -4.64774927155200000e-8,
      1.53835215257600000e-8, -5.11961333760000000e-9,
      1.82243164160000000e-9, -6.13513953280000000e-10,
      1.27679856640000000e-10,-4.01499750400000000e-11,
      4.26560716800000000e-11,-1.46381209600000000e-11);

var tvsmall, t, g: ArbFloat;
             m, i: ArbInt;
begin
  if sizeof(ArbFloat) = 6
  then
    tvsmall:=2*midget
  else
    tvsmall:=midget;
  t:=abs(x);
  if t > tmax
  then
    RunError(407);
  if t < macheps
  then
    begin
      if t < tvsmall
      then
        RunError(407);
      spegam:=1/x
    end
  else  { abs(x) >= macheps }
    begin
      m:=trunc(x);
      if x > 0
      then
        begin
          t:=x-m; m:=m-1; g:=1;
          if m<0
          then
            g:=g/x
          else
            if m>0
            then
              for i:=1 to m do
                g:=(x-i)*g
        end
      else { x < 0 }
        begin
          t:=x-m+1;
          if t=1
          then
            RunError(407);
          m:=1-m;
          g:=x;
          for i:=1 to m do
            g:=(i+x)*g;
          g:=1/g
        end;
      spegam:=spepol(2*t-1, a[0], sizeof(a) div sizeof(ArbFloat) - 1)*g
    end { abs(x) >= macheps }
end; {spegam}

function spelga(x: ArbFloat): ArbFloat;

const

    xbig = 7.7e7;
    xmax = 2.559e305;
  lnr2pi = 9.18938533204672742e-1;
    a: array[0..23] of ArbFloat =
    ( 8.86226925452758013e-1,  1.61691987244425092e-2,
      1.03703363422075456e-1, -1.34118505705967765e-2,
      9.04033494028101968e-3, -2.42259538436268176e-3,
      9.15785997288933120e-4, -2.96890121633200000e-4,
      1.00928148823365120e-4, -3.36375833240268800e-5,
      1.12524642975590400e-5, -3.75499034136576000e-6,
      1.25281466396672000e-6, -4.17808776355840000e-7,
      1.39383522590720000e-7, -4.64774927155200000e-8,
      1.53835215257600000e-8, -5.11961333760000000e-9,
      1.82243164160000000e-9, -6.13513953280000000e-10,
      1.27679856640000000e-10,-4.01499750400000000e-11,
      4.26560716800000000e-11,-1.46381209600000000e-11);
    b: array[0..5] of ArbFloat =
    ( 8.33271644065786580e-2,  -6.16502049453716986e-6,
      3.89978899876484712e-9,  -6.45101975779653651e-12,
      2.00201927337982364e-14, -9.94561064728159347e-17);


var t, g : ArbFloat;
    m, i : ArbInt;

begin
  if x <= 0
  then
    RunError(408);
  if x <= macheps
  then
    spelga:=-ln(x)
  else
  if x <= 15
  then
    begin
      m:=trunc(x); t:=x-m; m:=m-1; g:=1;
      if m < 0
      then
        g:=g/x
      else
      if m > 0
      then
        for i:=1 to m do
          g:=(x-i)*g;
      spelga:=ln(g*spepol(2*t-1, a[0], sizeof(a) div sizeof(ArbFloat) - 1))
    end
  else { x > 15 }
  if x <= xbig
  then
    spelga:=(x-0.5)*ln(x) - x + lnr2pi
            + spepol(450/sqr(x)-1, b[0], sizeof(b) div sizeof(ArbFloat) - 1)/x
  else { x > xbig }
  if x <= xmax
  then
    spelga:=(x-0.5)*ln(x) - x + lnr2pi
  else  { x > xmax => x*ln(x) > giant }
    RunError(408)
end; {spelga}

function spepol(x: ArbFloat; var a: ArbFloat; n: ArbInt): ArbFloat;
var   pa : ^arfloat0;
       i : ArbInt;
    polx : ArbFloat;
begin
  pa:=@a;
  polx:=0;
  for i:=n downto 0 do
    polx:=polx*x+pa^[i];
  spepol:=polx
end {spepol};

function spepow(a, b: ArbFloat): ArbFloat;

   function PowInt(a: double; n: longint): double;
   var a1 : double;
   begin
     if n=0 then PowInt := 1 else
     begin
        a1 := 1;
        if n<0 then begin a := 1/a; n := -n end;
        while n>0
        do if Odd(n)
           then begin Dec(n); a1 := a1*a end
           else begin n := n div 2; a := sqr(a) end;
        PowInt := a1
     end
   end;

var tb : longint;
    fb : double;
begin

  { (a < 0, b niet geheel) of (a = 0, b <= 0), dan afbreken}
  if (a=0) then if (b<=0) then RunError(400) else begin SpePow :=0; exit end;
  tb := Trunc(b); fb := b-tb;
  if (a<0) and (fb<>0) then RunError(400);

  if a>0
  then if fb=0 then SpePow := PowInt(a, tb)
               else SpePow := PowInt(a, tb)*exp(fb*ln(a))
  else if odd(tb) then Spepow := -PowInt(-a, tb)
                  else SpePow := PowInt(-a, tb)

end; {spepow}

function spesgn(x: ArbFloat): ArbInt;

begin
  if x<0
  then
    spesgn:=-1
  else
    if x=0
    then
      spesgn:=0
    else
      spesgn:=1
end; {spesgn}

function spears(x: ArbFloat): ArbFloat;
const

    pi2 = 1.570796326794897;
    a : array[0..17] of ArbFloat =
    (  1.047197551196598e+0, 5.375149359132719e-2, 7.798902238957732e-3,
       1.519668539582420e-3, 3.408637238430600e-4, 8.302317819598986e-5,
       2.134554822576075e-5, 5.701781046148566e-6, 1.566985123962741e-6,
       4.402076871418002e-7, 1.257811162594110e-7, 3.646577948300129e-8,
       1.081021746966715e-8, 3.212744286269388e-9, 8.515014302985799e-10,
       2.513296398553196e-10, 1.342121568282535e-10, 4.210346761190271e-11);

var    y, u, t, s : ArbFloat;
    uprang        : boolean;
begin
  if abs(x) > 1
  then
    RunError(401);
  u:=sqr(x); uprang:= u > 0.5;
  if uprang
  then
    u:=1-u;
  t:=4*u-1; y:=spepol(t, a[0], sizeof(a) div sizeof(ArbFloat) - 1);
  if uprang
  then
    begin
      s:=pi2-sqrt(u)*y;
      if x < 0
      then
        s:=-s;
      spears:=s
    end
  else
    spears:=x*y
end;  {spears}

function spearc(x: ArbFloat): ArbFloat;
const

    pi2 = 1.570796326794897;
    a : array[0..17] of ArbFloat =
    ( 1.047197551196598e+0,  5.375149359132719e-2,  7.798902238957732e-3,
      1.519668539582420e-3,  3.408637238430600e-4,  8.302317819598986e-5,
      2.134554822576075e-5,  5.701781046148566e-6,  1.566985123962741e-6,
      4.402076871418002e-7,  1.257811162594110e-7,  3.646577948300129e-8,
      1.081021746966715e-8,  3.212744286269388e-9,  8.515014302985799e-10,
      2.513296398553196e-10, 1.342121568282535e-10, 4.210346761190271e-11);

var u, t, y, s    : ArbFloat;
    uprang        : boolean;
begin
  if abs(x)>1.0
  then
    RunError(402);
  u:=sqr(x); uprang:=u>0.5;
  if uprang
  then
    u:=1-u;
  t:=4*u-1; y:=spepol(t, a[0], sizeof(a) div sizeof(ArbFloat) - 1);
  if uprang
  then
    begin
      s:=sqrt(u)*y;
      if x<0
      then
        s:=2*pi2-s;
      spearc:=s
    end
  else
    spearc:=pi2-x*y
end;  {spearc}

function spesih(x: ArbFloat): ArbFloat;
const

    a : array[0..6] of ArbFloat =
    ( 1.085441641272607e+0,  8.757509762437522e-2,  2.158779361257021e-3,
      2.549839945498292e-5,  1.761854853281383e-7,  7.980704288665359e-10,
      2.551377137317034e-12);

var u : ArbFloat;
begin
  if abs(x)<=1.0
  then
    begin
      u:=2*sqr(x)-1;
      spesih:=x*spepol(u, a[0], sizeof(a) div sizeof(ArbFloat) - 1)
    end
  else
  begin
    u:=exp(x); spesih:=(u-1/u)/2
  end
end; {spesih}

function specoh(x: ArbFloat): ArbFloat;
var u: ArbFloat;
begin
  u:=exp(x); specoh:=(u+1/u)/2
end; {specoh}

function spetah(x: ArbFloat): ArbFloat;
const
    xhi = 18.50;
    a : array[0..15] of ArbFloat =
    ( 8.610571715805476e-1, -1.158834489728470e-1,  1.918072383973950e-2,
     -3.225255180728459e-3,  5.433071386922689e-4, -9.154289983175165e-5,
      1.542469328074432e-5, -2.599022539340038e-6,  4.379282308765732e-7,
     -7.378980192173815e-8,  1.243517352745986e-8, -2.095373768837420e-9,
      3.509758916273561e-10,-5.908745181531817e-11, 1.124199312776748e-11,
     -1.907888434471600e-12);

var t, y: ArbFloat;

begin
  t:=abs(x);
  if t <= 1
  then
    begin
      y:=2*sqr(x)-1;
      spetah:=x*spepol(y, a[0], sizeof(a) div sizeof(ArbFloat) - 1)
    end
  else
  if t < xhi
  then
    begin
      y:=exp(2*x); spetah:=(y-1)/(y+1)
    end
  else
    spetah:=spesgn(x)
end; {spetah}

function speash(x: ArbFloat): ArbFloat;
const

    xhi = 1e9;
    c : array[0..18] of ArbFloat =
    (  9.312298594527122e-1,  -5.736663926249348e-2,
       9.004288574881897e-3,  -1.833458667045431e-3,
       4.230023450529706e-4,  -1.050715136470630e-4,
       2.740790473603819e-5,  -7.402952157663977e-6,
       2.052474396638805e-6,  -5.807433412373489e-7,
       1.670117348345774e-7,  -4.863477336087045e-8,
       1.432753532351304e-8,  -4.319978113584910e-9,
       1.299779213740398e-9,  -3.394726871170490e-10,
       1.008344962167889e-10, -5.731943029121004e-11,
       1.810792296549804e-11);


var t : ArbFloat;

begin
  t:=abs(x);
  if t <= 1 then
    speash:=x*spepol(2*sqr(x)-1, c[0], sizeof(c) div sizeof(ArbFloat) - 1)
  else
  if t < xhi then
    speash:=ln(sqrt(sqr(x)+1)+t)*spesgn(x)
  else
    speash:=ln(2*t)*spesgn(x)
end; {speash}

function speach(x: ArbFloat): ArbFloat;
const

    xhi = 1e9;

begin
  if x<1 then
    RunError(405);
  if x=1 then
    speach:=0
  else
  if x<=xhi then
    speach:=ln(x+sqrt(sqr(x)-1))
  else
    speach:=ln(2*x)
end; {speach}

function speath(x: ArbFloat): ArbFloat;
const

    c : array[0..19] of ArbFloat =
    ( 1.098612288668110e+0,  1.173605223326117e-1,  2.309071936165689e-2,
      5.449091889986991e-3,  1.404884102286929e-3,  3.816948426588058e-4,
      1.073604335435426e-4,  3.095027782918129e-5,  9.088050814470148e-6,
      2.706881064641104e-6,  8.155200644023077e-7,  2.479830612463254e-7,
      7.588067811453948e-8,  2.339295963220429e-8,  7.408486568719348e-9,
      2.319454882064018e-9,  5.960921368486746e-10, 1.820410351379402e-10,
      1.184977617320312e-10, 3.856235316559190e-11);

var t, u: ArbFloat;
begin
  t:=abs(x);
  if t >= 1 then
    RunError(406);
  u:=sqr(x);
  if u < 0.5 then
    speath:=x*spepol(4*u-1, c[0], sizeof(c) div sizeof(ArbFloat) - 1)
  else { 0.5 < x*x < 1 }
    speath:=ln((1+t)/(1-t))/2*spesgn(x)
end; {speath}

var exitsave : pointer;

procedure MyExit; Far;
const ErrorS : array[400..408,1..6] of char =
     ('spepow',
      'spebk0',
      'spebk1',
      'speby0',
      'speby1',
      'speach',
      'speath',
      'spegam',
      'spelga');

var ErrFil : text;

begin
     ExitProc := ExitSave;
     Assign(ErrFil, 'CON');
     ReWrite(ErrFil);
     if (ExitCode>=400) AND (ExitCode<=408) then
       begin
         write(ErrFil, 'critical error in ', ErrorS[ExitCode]);
         ExitCode := 201
       end;
     Close(ErrFil)
end;

begin
   ExitSave := ExitProc;
   ExitProc := @MyExit;
end.
{
  $Log$
  Revision 1.1  2000-07-13 06:34:16  michael
  + Initial import

  Revision 1.2  2000/01/25 20:21:42  marco
   * small updates, crlf fix, and RTE 207 problem

  Revision 1.1  2000/01/24 22:08:58  marco
   * initial version


}