module erf_mod;

{***************************************************************}
{*                                                             *}
{*  ERF(x) :  eps(ERF) = 2.7153E-15 bei hochgenauer Arithmetik *}
{*            eps(ERF) = 1.5643E-15 bei maximalgen. Arithmetik *}
{*                              |x| in [1.97193E-308,+maxreal] *}
{*                                                             *}
{*  ERFC(x):  eps(ERFC)= 6.0761E-15 bei hochgenauer Arithmetik *}
{*            eps(ERFC)= 3.4063E-15 bei maximalgen. Arithmetik *}
{*                                    x in [-maxreal, 26.5432] *}
{*                                                             *}
{*  Benutzt wird die Standard-Exponentialfunktion mit:         *}
{*                   eps(exp) = 2.3580E-16                     *}
{***************************************************************}

use erf_help, iostd;

const al_erf  = 1.5659E-15;   {* = 1.001*>eps(erf) bei maximal *}
                              {*  genauer Arithmetik           *}
      al_erfc = 3.2985E-15;   {* =1.001*>eps(erfc) bei maximal *}
                              {*  genauer Arithmetik           *}
{***************************************************************}
{*  Die folgenden Konstanten sind bei hochgenauer  Arithmetik  *}
{*                                               zu benutzen:  *}
{*    al_erf  = 2.7181E-15       = 1.001*>eps(erf) bei hochge- *}
{*                                           nauer Arithmetik  *}
{*    al_erfc = 5.8599E-15       = 1.001*>eps(erf) bei hochge- *}
{*                                           nauer Arithmetik  *}
{***************************************************************}

var erf0_P4, erf0_Q4, erfc3_P4, erfc3_Q4 : feld_1[0..4];
    erfc1_P5, erfc2_P5, erfc_ber         : feld_1[0..5];
    erfc1_Q6, erfc2_Q6                   : feld_1[0..6];
    erf_ber                              : feld_1[0..3];

    {*****************  Fehler-Kommentare:   *******************}
    comment : array[1..1] of string[45];

procedure ERROR_SPECF_1( NAME : string; ERROR_NR : integer;
                         x : real );
const c = '@@@@@@@   ';
var out : TEXT;
begin
  rewrite(out,'CON');
  writeln(out);
  writeln(out,c,'error in ',NAME);
  writeln(out,c,'argument = ',x:23:0);
  writeln(out,c,comment[ERROR_NR]);
  writeln;
  EXIT(0)
end;  { ERROR_SPECF_1 }

{********************************************************************}
{*                      Fehlerfunktionen                            *}
{********************************************************************}

function ERF0( x : real) : real;
{ erf(x) :   x in [0,0.65]                }
{ x < 1.97193*10^-308  --> Fehlermeldung! }
var y, x2   : real;
    bereich : integer;
begin { ERF0 }
  bereich := INT_NR(erf_ber,x);
  case bereich of
  0 : begin
        if SIGN(x)=0 then ERF0 := 0 
                     else ERROR_SPECF_1('ERF',1,x);
      end;
  1 : ERF0 := erf0_P4[0] * x;
  2 : begin
        x2 := x * x;
        y  := x * HORNER(erf0_P4,x2);
        x2 := HORNER(erf0_Q4,x2);
        ERF0 := y / x2;
      end;
  end; { case }
end;  { ERF0 }

function ERF0000( x : real) : real;
{ erf(x) :   x in [0,0.65]              }
{ x < 1.97193*10^-308  --> erf(x) = 0   }
var y, x2   : real;
    bereich : integer;
begin { ERF0 }
  bereich := INT_NR(erf_ber,x);
  case bereich of
  0 : ERF0000 := 0;
  1 : ERF0000 := erf0_P4[0] * x;
  2 : begin
        x2 := x * x;
        y  := x * HORNER(erf0_P4,x2);
        x2 := HORNER(erf0_Q4,x2);
        ERF0000 := y / x2;
      end;
  end; { case }
end;  { ERF0000 }

function ERFC_B1( x : real) : real;
{ erfc(x) :  x in B1=[0.65,2.2] }
var P5, Q6 : real;
begin { ERFC_B1 }
  P5 := HORNER(erfc1_P5,x);
  Q6 := HORNER(erfc1_Q6,x);
  ERFC_B1 := exp_x2(x) * P5 / Q6;
end;  { ERFC_B1 }

function ERFC_B2( x : real) : real;
{ erfc(x) :  x in B2=[2.2,6] }
var P5, Q6 : real;
begin { ERFC_B2 }
  P5 := HORNER(erfc2_P5,x);
  Q6 := HORNER(erfc2_Q6,x);
  ERFC_B2 := exp_x2(x) * P5 / Q6;
end;  { ERFC_B2 }

function ERFC_B3( x : real) : real;
{ erfc(x) :  x in B3=[6,26.5432] }
var P4, Q4,x2 : real;
begin { ERFC_B3 }
  x2 := sqr(1/x);
  P4 := HORNER(erfc3_P4,x2);
  Q4 := HORNER(erfc3_Q4,x2);
  ERFC_B3 := ( exp_x2(x) * P4 ) / ( x * Q4 );
end;  { ERFC_B3 }

function ERFC_GR0( x : real ) : real;
{ erfc(x) :  x >= 0;                    }
{ Fuer x > 26.5432 gilt:   ERFC_GR0 = 0 }
var bereich : integer;
begin { ERFC_GR0 }
  bereich := INT_NR(erfc_ber,x);
  case bereich of
  0 : ERFC_GR0 := 1 - ERF0000(x);
  1 : ERFC_GR0 := ERFC_B1(x);
  2 : ERFC_GR0 := ERFC_B2(x);
  3 : ERFC_GR0 := ERFC_B3(x);
  4 : ERFC_GR0 := 0;
  end; { case }
end;  { ERFC_GR0 }

function ERF_GR000( x : real ) : real;
{ erf(x) fuer x >= 0                    }
{ x < 1.97193*10^-308  -->  erf(x) = 0  }
begin { ERF_GR000 }
  if x >= erfc_ber[1] then ERF_GR000 := 1 - ERFC_GR0(x)
                      else ERF_GR000 := ERF0000(x) 
end;  { ERF_GR000 }

function ERFC000I( x : real ) : real;
{ erfc(x) zur Berechnung der Intervallfunktion }
{ Fuer x > 26.5432 wird ERFC000I auf 0 gesetzt }
begin { ERFC000I }
  if SIGN(x) < 0 then ERFC000I := 1 + ERF_GR000(-x)
                 else ERFC000I := ERFC_GR0(x);
  FW_exakt := SIGN(x) = 0;
end;  { ERFC000I }

function ERF_GR0( x : real ) : real;
{ erf(x) fuer x >= 0                      }
{ x < 1.97193*10^-308  -->  Fehlermeldung }
begin { ERF_GR0 }
  if x >= erfc_ber[1] then ERF_GR0 := 1 - ERFC_GR0(x)
                      else ERF_GR0 := ERF0(x) 
end;  { ERF_GR0 }

function ERF000I( x : real ) : real;
{ erf(x) fuer alle x aus dem IEEE double-Format }
{ Relativer Fehler :  eps(ERF) = 1.6236*10^-15  }
{ x < 1.97193*10^-308  -->  erf(x) = 0          }
begin { ERF000I }
  if SIGN(x) < 0 then ERF000I := -ERF_GR000(-x)
                 else ERF000I := ERF_GR000(x);
  FW_exakt := SIGN(x) = 0;
end;  { ERF000I }   


global function ERF( x : real ) : real;
{ erf(x) fuer alle x aus dem IEEE double-Format }
{ Relativer Fehler :  eps(ERF) = 1.6236*10^-15  }
{ x < 1.97193*10^-308  -->  Fehlermeldung       }
begin { ERF }
  if SIGN(x) < 0 then ERF := -ERF_GR0(-x)
                 else ERF := ERF_GR0(x);
end;  { ERF }  


global function ERFC( x : real ) : real;
{ erfc(x) fuer alle x aus dem IEEE double-Format           }
{ Bei Underflow (x > 26.5432) erfolgt eine Fehlermeldung,  }
{ da sonst  eps(ERFC) = 3.7323*10^-15 nicht realisiert ist }
begin { ERFC }
  if x >= erfc_ber[4] then  ERROR_SPECF_1('ERFC',1,x)
                      else
  if SIGN(x) < 0 then ERFC := 1 + ERF_GR000(-x)
                 else ERFC := ERFC_GR0(x)
end;  { ERFC }


global function erf( x : interval ) : interval;
var c : interval;
begin
  c := INTEXTD_NORM(ERF000I,x,TRUE,al_erf);
  if c.sup > +1 then c.sup := +1;
  if c.inf < -1 then c.inf := -1;
  if (SIGN(x.inf) > 0) and (SIGN(c.inf) < 0) then c.inf := 0; 
  if (SIGN(x.sup) < 0) and (SIGN(c.sup) > 0) then c.sup := 0; 
  erf := c;
end;

global function erfc( x : interval) : interval;
var c : interval;
begin
  c := INTEXTD_NORM(ERFC000I,x,FALSE,al_erfc);
  if c.sup > 2 then c.sup := 2;
  if SIGN(c.inf) < 0 then c.inf := 0;
  if (SIGN(x.inf) > 0) and (c.sup > 1) then c.sup := 1;
  if (SIGN(x.sup) < 0) and (c.inf < 1) then c.inf := 1;
  erfc := c;
end;

begin
  erf0_P4[0] := 1.12837916709551256; 
  erf0_P4[1] := 1.35894887627277916e-1; 
  erf0_P4[2] := 4.03259488531795274e-2; 
  erf0_P4[3] := 1.20339380863079457e-3; 
  erf0_P4[4] := 6.49254556481904354e-5; 

  erf0_Q4[0] := 1; 
  erf0_Q4[1] := 4.53767041780002545e-1; 
  erf0_Q4[2] := 8.69936222615385890e-2; 
  erf0_Q4[3] := 8.49717371168693357e-3; 
  erf0_Q4[4] := 3.64915280629351082e-4; 

  erfc1_P5[0] := 9.99999992049799098e-1; 
  erfc1_P5[1] := 1.33154163936765307; 
  erfc1_P5[2] := 8.78115804155881782e-1; 
  erfc1_P5[3] := 3.31899559578213215e-1; 
  erfc1_P5[4] := 7.14193832506776067e-2; 
  erfc1_P5[5] := 7.06940843763253131e-3; 
  
  erfc1_Q6[0] := 1; 
  erfc1_Q6[1] := 2.45992070144245533; 
  erfc1_Q6[2] := 2.65383972869775752; 
  erfc1_Q6[3] := 1.61876655543871376; 
  erfc1_Q6[4] := 5.94651311286481502e-1; 
  erfc1_Q6[5] := 1.26579413030177940e-1; 
  erfc1_Q6[6] := 1.25304936549413393e-2; 
  
  erfc2_P5[0] := 9.99921140009714409e-1; 
  erfc2_P5[1] := 1.62356584489366647; 
  erfc2_P5[2] := 1.26739901455873222; 
  erfc2_P5[3] := 5.81528574177741135e-1; 
  erfc2_P5[4] := 1.57289620742838702e-1; 
  erfc2_P5[5] := 2.25716982919217555e-2; 
  
  erfc2_Q6[0] := 1; 
  erfc2_Q6[1] := 2.75143870676376208; 
  erfc2_Q6[2] := 3.37367334657284535; 
  erfc2_Q6[3] := 2.38574194785344389; 
  erfc2_Q6[4] := 1.05074004614827206; 
  erfc2_Q6[5] := 2.78788439273628983e-1; 
  erfc2_Q6[6] := 4.00072964526861362e-2; 
  
  erfc3_P4[0] := 5.64189583547756078e-1; 
  erfc3_P4[1] := 8.80253746105525775; 
  erfc3_P4[2] := 3.84683103716117320e+1; 
  erfc3_P4[3] := 4.77209965874436377e+1; 
  erfc3_P4[4] := 8.08040729052301677; 
  
  erfc3_Q4[0] := 1; 
  erfc3_Q4[1] := 1.61020914205869003e+1; 
  erfc3_Q4[2] := 7.54843505665954743e+1; 
  erfc3_Q4[3] := 1.12123870801026015e+2; 
  erfc3_Q4[4] := 3.73997570145040850e+1; 
  
  erf_ber[0] := 0;
  erf_ber[1] := (>1.97193e-308);
  erf_ber[2] := (>1e-10);
  erf_ber[3] := (>0.65);

  erfc_ber[0] := 0;
  erfc_ber[1] := (>0.65);
  erfc_ber[2] := (>2.2);
  erfc_ber[3] := 6;
  erfc_ber[4] := (>26.5432);
  erfc_ber[5] := 30;

  comment[1] := 'Underflow --> rel. error exceeded';
end.