module erf_help;


use x_real;

global const MIN_EXP  = -1022;
             MIN_EXP_ = -1074;

var kl_norm_ :real;    { kleinste pos. normalisierte IEEE-Zahl }

global var FW_exakt : boolean;

global type feld_1  = global dynamic array[*] of real;
global type feld_2  = global dynamic array[*,*] of real;
global type feld_1i = global dynamic array[*] of integer;

var e_z2 : feld_1[0..27];             { e^(-z^2) - Konstanten  }

global function EXPO_0( x : real ) : integer;
{**************************************************************}
{*   EXPO_0(0.0) = 0 im Gegensatz zu expo(0.0) = -2147483647  *}
{**************************************************************}
begin
  EXPO_0 := expo( x );
  if SIGN(x) = 0 then EXPO_0 := 0;
end;

priority M_2HOCH = *;
{**************************************************************************}
{*  y := a M_2HOCH k  liefert das Produkt  a * 2^k;                       *}
{*  Im Fall  |a*2^k|<2^-1022  wird das Ergebnis zur naechsten Rasterzahl  *}
{*  gerundet, d.h. Verhalten wie bei der normalen Multiplikation mit 2^k  *}
{**************************************************************************}
global operator M_2HOCH(a : real; k : integer) res : real;
label 1;
var m       : real;
    ex      : integer;
    zuklein : boolean;
begin
     zuklein := k < MIN_EXP_;
     if zuklein then k := k - MIN_EXP_;
  1: ex := EXPO_0(a) + k;
     if ex < (MIN_EXP+1) then a := a * COMP(0.5,k+1)
                         else a := COMP(MANT(a),ex);
     if zuklein then 
     begin 
       k := MIN_EXP_;
       zuklein := false;
       goto 1;
     end;
     res := a;
end;

GLOBAL FUNCTION rounduod(y:real; alpha:real):real;

 {************************************************************************}
 {*   y: Maschinenergebnis des exakten Wertes yo<>y mit der Fehler-      *}
 {*   schranke epso>0;          y=yo(1+eps),  mit: |eps| <= epso>0;      *}
 {*   alpha:=1.001*>epso,  d.h. alpha>0 mu� vorher berechnet werden!     *}
 {*   ROUNDUOD(y,(+/-)alpha) erm�glicht die Berechnung einer Ein-        *}
 {*   schlie�ung des exakten Wertes yo OHNE MULTIPLIKATION:              *}   
 {*                           yU <= yo <= yO;                            *}
 {*   yO:=ROUNDUOD(y,+alpha);  liefert eine Oberschranke:     yo<=yO;    *}
 {*   yU:=ROUNDUOD(y,-alpha);  liefert eine Unterschranke:    yU<=yo;    *}   
 {*   F�r y=0 wird die Einschlie�ung     yU=-2^(-1073)<yo<+2^(-1073)=yO  *}  
 {*   angegeben f�r den Fall, da� y=0 durch Underflow berechnet wurde.   *}
 {************************************************************************}
var ex,ex1 : integer;
    lr     : real; 
begin  ex:=EXPO_0(y);  y:=MANT(y);
  if y=0 then
     if alpha>0 then rounduod:=2*minreal    { minreal = 2^(-1074) }
                else rounduod:=-2*minreal
         else
     begin  { y<>0 }
       if alpha>0 then y:=y+>alpha
                  else y:=y+<alpha;
       ex1:=EXPO_0(y);
       if ex1<>0 then
       begin ex:=ex+ex1;  y:=MANT(y)  end;
         if ex < min_exp+1 then                                   
         begin
           lr := COMP(0.5,ex+1);
           if alpha > 0 then rounduod := y *> lr
                        else rounduod := y *< lr; 
         end        else  rounduod:=COMP(y,ex);
     end
end;

GLOBAL FUNCTION round_norm(y : real; alpha : real) : real;

 {************************************************************************}
 {*   y: Maschinenergebnis des exakten Wertes yo<>y mit der Fehler-      *}
 {*   schranke epso>0;          y=yo(1+eps),  mit: |eps| <= epso>0;      *}
 {*   alpha:=1.001*>epso,  d.h. alpha>0 mu� vorher berechnet werden!     *}
 {*   ROUND_NORM(y,(+/-)alpha)   erm�glicht die Berechnung einer Ein-    *}
 {*   schlie�ung des exakten Wertes yo OHNE MULTIPLIKATION:              *}   
 {*                           yU <= yo <= yO;                            *}
 {*   yO:=ROUND_NORM(y,+alpha);  liefert eine Oberschranke:     yo<=yO;  *}
 {*   yU:=ROUND_NORM(y,-alpha);  liefert eine Unterschranke:    yU<=yo;  *}
 {*   F�r y=0 wird die Einschlie�ung     yU=-2^(-1021)<yo<+2^(-1021)=yO  *}  
 {*   angegeben f�r den Fall, da� y auf Null gesetzt wurde und yo durch  *}
 {*   normalisierte Schranken   yU, yO   einzuschliessen ist.            *}
 {************************************************************************}
var ex,ex1 : integer;
    lr     : real; 
begin  ex:=EXPO_0(y);  y:=MANT(y);
  if y=0 then
     if alpha>0 then round_norm:=2*kl_norm_    { kl_norm_ = 2^(-1022) }
                else round_norm:=-2*kl_norm_
         else
     begin  { y<>0 }
       if alpha>0 then y:=y+>alpha
                  else y:=y+<alpha;
       ex1:=EXPO_0(y);
       if ex1<>0 then
       begin ex:=ex+ex1;  y:=MANT(y)  end;
         if ex < min_exp+1 then                                   
         begin
           lr := COMP(0.5,ex+1);
           if alpha > 0 then round_norm := y *> lr
                        else round_norm := y *< lr; 
         end        else  round_norm:=COMP(y,ex);
     end
end; { round_norm }

GLOBAL FUNCTION horner(a: feld_1; x: real) : real;

 {************************************************************************}
 {*  Polynom:      P(x) = a[0]*x^0 + a[1]*x^1 + .. + a[N]*x^N            *}
 {*  Auswertung des Polynoms nach dem HORNER-Schema:   horner := P(x)    *}
 {*  Die Indizierung der a[k] ist durch das Feld  a : feld_1  vorgegeben *}
 {*  und muss nicht mit der Indizierung des obigen  P(x) uebereinstimmen *}
 {************************************************************************}

var k : integer;
    y : real;
begin { horner }
  y := a[ubound(a)];
  For k := ubound(a)-1 downto lbound(a) do
    y      := y * x + a[k];
    horner := y
end;  { horner }

GLOBAL FUNCTION INT_NR(b: feld_1; x: real) : integer;

 {************************************************************************}
 {* Ein Intervall  [A,B]  wird durch das Feld    b: feld_1   unterteilt, *}
 {* mit:  A = b[lbound(b)]  und  B = b[ubound(b)];                       *}
 {* [c,d] sei ein Teilintervall von [A,B];                               *}
 {* Ist x enthalten im halboffenen Intervall  [c,d[  , mit:  c=b[k] , so *}
 {* gilt:   INT_NR = k,  d.h. die Intervall-Nr. ist der Index des linken *}
 {* Intervallrandpunktes.                                                *}
 {* Fuer x >= B = b[ubound(b)] gilt:  INT_NR = ubound(b)                 *}
 {* Fuer x <  A = b[lbound(b)] gilt:  INT_NR = lbound(b) - 1             *}
 {************************************************************************}

var i,j,k : integer;
begin { INT_NR }
  i := lbound(b);
  j := ubound(b) - 1;
  repeat
    k := (i+j) div 2;
    if x < b[k] then j := k-1 else i := k+1
  until i > j;
  INT_NR := j;
end;  { INT_NR }

global FUNCTION INT_NR(b: feld_1i; x: integer) : integer;

 {************************************************************************}
 {* Ein Intervall  [A,B]  wird durch das Feld    b: feld_1i  unterteilt, *}
 {* mit:  A = b[lbound(b)]  und  B = b[ubound(b)];                       *}
 {* [c,d] sei ein Teilintervall von [A,B];                               *}
 {* Ist x enthalten im halboffenen Intervall  [c,d[  , mit:  c=b[k] , so *}
 {* gilt:   INT_NR = k,  d.h. die Intervall-Nr. ist der Index des linken *}
 {* Intervallrandpunktes.                                                *}
 {* Fuer x >= B = b[ubound(b)] gilt:  INT_NR = ubound(b)                 *}
 {* Fuer x <  A = b[lbound(b)] gilt:  INT_NR = lbound(b) - 1             *}
 {************************************************************************}

var i,j,k : integer;
begin { INT_NR }
  i := lbound(b);
  j := ubound(b) - 1;
  repeat
    k := (i+j) div 2;
    if x < b[k] then j := k-1 else i := k+1
  until i > j;
  INT_NR := j;
end;  { INT_NR }

global function EXP_X2( x : real ) : real;
 {************************************************************************}
 {*  EXP_X2(x) = e^(-x^2);                                               *}
 {*  rel. Fehler = 8.3243E-16        bei maximalgenauen Grundoperationen *}
 {*  rel. Fehler = 1.0823E-15        bei hochgenauen    Grundoperationen *}
 {*  |x| <= 26.615717;                                                   *}
 {*  |x| >  26.615717  --->   rel. Fehler wird ueberschritten; aus Lauf- *} 
 {*  zeitgruenden erfolgt jedoch keine Fehlermeldung !!!                 *}
 {*                                                                      *}
 {*  Zu benutzen ist die schnelle Exponentialfunktion(Tabellenverfahren) *}
 {*  mit der Fehlerschranke :            eps(exp) = 2.3580E-16           *}
 {*                                                                      *}
 {*  Blomquist,  23.08.95                                                *}
 {************************************************************************}
var z   : integer;
    m,t : real;
begin
  if SIGN(x) < 0 then x := -x;
  z := trunc(x);  { z: ganzzahliger Anteil von x   } 
  m := x - z;     { m: zugehoeriger Nachkomma-Wert }
  if m > 0.5 then
  begin
    z := z + 1;
    m := m - 1;   { m <= 0.5;           x = z + m; }
  end;
  t := e_z2[z] * exp( -(z+z)*m ) * exp( -m*m );
  if z = 27 then t := t M_2HOCH -64;
  EXP_X2 := t;
end;  { EXP_X2 }


var messagepipe : boolean;  { FALSE by default, TRUE if special }
                            { messagehandling has been done     }

global function INTEXTD( function f(x : real): real; i : interval;
                             up : boolean; al : real ) : interval;
var c  : interval;                 
    xr : real;
begin    
  if not up then
    begin xr:=i.INF; i.INF:=i.SUP; i.SUP:=xr end;
  IF messagepipe THEN xr:=-f(i.INF)
                 ELSE xr:=f(i.INF);  
  if FW_exakt then c.INF := xr               
              else c.INF := ROUNDUOD(xr,-al);
  if i.INF<>i.SUP then xr:=f(i.SUP);            
  if FW_exakt then c.SUP := xr
              else c.SUP := ROUNDUOD(xr,+al);
  INTEXTD:=c            
end;  

global function INTEXTD_NORM( function f(x : real): real; i : interval;
                             up : boolean; al : real ) : interval;
var c  : interval;                 
    xr : real;
begin    
  if not up then
    begin xr:=i.INF; i.INF:=i.SUP; i.SUP:=xr end;
  IF messagepipe THEN xr:=-f(i.INF)
                 ELSE xr:=f(i.INF);  
  if FW_exakt then c.INF := xr               
              else c.INF := ROUND_NORM(xr,-al);
  if i.INF<>i.SUP then xr:=f(i.SUP);            
  if FW_exakt then c.SUP := xr
              else c.SUP := ROUND_NORM(xr,+al);
  INTEXTD_NORM := c            
end;  

begin
  messagepipe := FALSE;

  { kl_norm_ : kleinste pos. normalisierte IEEE-Zahl }
  kl_norm_ := (>2.225073858507201383e-308);

  e_z2[0]   := 1;                         {  e^(-0^2)  }
  e_z2[1]   := 0.367879441171442321;      {  e^(-1^2)  }
  e_z2[2]   := 1.83156388887341802e-2;    {  e^(-2^2)  }
  e_z2[3]   := 1.23409804086679549e-4;    {     usw.   }
  e_z2[4]   := 1.12535174719259114e-7;
  e_z2[5]   := 1.38879438649640205e-11;
  e_z2[6]   := 2.319522830243569388e-16;
  e_z2[7]   := 5.24288566336346393e-22;
  e_z2[8]   := 1.60381089054863785e-28;
  e_z2[9]   := 6.63967719958073440e-36;
  e_z2[10]  := 3.72007597602083596e-44;
  e_z2[11]  := 2.82077008846013540e-53;
  e_z2[12]  := 2.89464031164830028e-63;
  e_z2[13]  := 4.02006021574335524e-74;
  e_z2[14]  := 7.55581901971196035e-86;
  e_z2[15]  := 1.92194772782384906e-98;
  e_z2[16]  := 6.61626105670948526e-112;
  e_z2[17]  := 3.08244069694909841e-126;
  e_z2[18]  := 1.94351485004929273e-141;
  e_z2[19]  := 1.65841047768114514e-157;
  e_z2[20]  := 1.91516959671400569e-174;
  e_z2[21]  := 2.99318445226019273e-192;
  e_z2[22]  := 6.33097733621059136e-211;
  e_z2[23]  := 1.81225402579399232e-230;
  e_z2[24]  := 7.02066779850473471e-251;
  e_z2[25]  := 3.68085585480180060e-272;
  e_z2[26]  := 2.61174176128405547e-294;  {  e^(-26^2)      }
  e_z2[27]  := 4.62639185846956420e-298;  {  e^(-27^2)*2^64 }
end.