module mpi_ari;
{============================================================================}
{===  mpi_ari  ======   Multi-Precision Interval Arithmetic   ===============}
{====================   =     =         =        ===          ===============}
{============================================================================}
{                                                                            }
{    File    :  MPI_ARI.P                                                    }
{                                                                            }
{=   Version :  1.1                                                          }
{                                                                            }
{=   Date    :  10/07/91  ( first implementation )                           }
{=              03/21/95  ( last change )                                    }
{                                                                            }
{    Language:  PASCAL-XSC, Vers. 3.0 or later                               }
{                                                                            }
{    Authors :  W. Kraemer, W. Hofschuster, R. Hammer                        }
{                                                                            }
{    Site    :  Institute for Applied Mathematics (IAM),                     }
{               University of Karlsruhe,                                     }
{               76128 Karlsruhe, Germany.                                    }
{                                                                            }
{    Tel.    :  +49 721 608-6013 (W. Kraemer)                                }
{    Fax     :  +49 721 69 52 83 (IAM)                                       }
{----------------------------------------------------------------------------}
{                                                                            }
{    Short Description:                                                      }
{    ------------------                                                      }
{      Module for multi-precision interval arithmetic. This module is based  }
{      on module MP_ARI for real multi-precision arithmetic.                 }
{                                                                            }
{    Literature:                                                             }
{    -----------                                                             }
{      [1] W. Kraemer: Eine portable Langzahl- und Langzahlintervallarithme- }
{                      tik mit Anwendungen. ZAMM 73, 1992.                   }
{                                                                            }
{      [2] W. Kraemer: Die Berechnung von Standardfunktionen in Rechenanla-  }
{                      gen. In Chatterji, S., Fuchssteiner, B., Kulisch, U., }
{                      Liedl, R., Purkert, W. (Eds):  Jahrbuch Ueberblicke   }
{                      Mathematik 1992, Vieweg, Braunschweig, 1992.          }
{                                                                            }
{      [3] W. Kraemer: Multiple-Precision Computations with Result Verifica- }
{                      tion. In Adams, E., Kulisch, U. (Eds.): Scientific    }
{                      Computing with Automatic Result Verification.         }
{                      Academic Press, San Diego, 1993.                      }
{                                                                            }
{      [4] R. Hammer : Multi-Precision Arithmetic in PASCAL-XSC, Implementa- }
{                      tion and Apllication. Special Issue of the SCAN-93    }
{                      and the MMSC'93, Math. and Comp. for Simulation, to   }
{                      appear 1994.                                          }
{============================================================================}
use
  iostd,    { For exception handling          }
  i_ari;    { Interval arithmetic             }

use global
  mp_ari;   { Real multi-precision arithmetic }


{============================================================================}
{==========   Datatypes and Variables   =====================================}
{============================================================================}
{ Multi-precision interval data type: }
global type  mpinterval  =  global record   inf,  sup: mpreal;  end;

var
  one  : mpreal;                         { mpreal and mpinterval }
  ione : mpinterval;                     { constants for 1       }


{============================================================================}
{==========   Memory Management   ===========================================}
{============================================================================}
global procedure mpinit( var x : mpinterval );     { Initialize a mpinterval }
begin                                              { variable                }
  mpinit(x.inf);
  mpinit(x.sup);
end;

global procedure mpvlcp ( var x : mpinterval );    { Get a local copy of a   }
begin                                              { mpinterval variable     }
  mpvlcp(x.inf);
  mpvlcp(x.sup);
end;

global procedure mpfree ( var x : mpinterval );    { Free a mpinterval       }
begin                                              { variable                }
  mpfree(x.inf);
  mpfree(x.sup);
end;


{============================================================================}
{==========   I/O - Routines   ==============================================}
{============================================================================}
global procedure read( var t : text; var x : mpinterval );   { Input without }
begin                                                        { brackets      }
  read(t,x.inf:-1);  { round downwards }                     {---------------}
  read(t,x.sup:+1);  { round upwards   }
end;

global procedure write( var t : text; x : mpinterval );
begin
  mpvlcp(x);
  writeln(t,x.inf:0:0:-1);
  write  (t,x.sup:0:0:+1);
  mpfree(x);
end;


{============================================================================}
{==========   Exception Handling   ==========================================}
{============================================================================}
const                     { Error codes }
  div_by_zero  =  1;      {-------------}
  sqrt_range   = 10;
  ln_range     = 20;
  loga_range   = 21;
  log2_range   = 22;
  log10_range  = 23;
  power_range  = 30;
  power_notdef = 31;
  arcsin_range = 40;
  arccos_range = 41;
  sin_range    = 50;
  cos_range    = 51;
  tan_range    = 52;
  cot_range    = 53;
  coth_range   = 60;
  arcosh_range = 70;
  artanh_range = 71;
  arcoth_range = 72;

procedure traceback( var t : text );
  external f_back ;

procedure exit( ReturnCode : integer );
  external a_exit;

procedure _math_message( ErrNo : integer; var Arg1, Arg2 : mpinterval );
var
  Msg : text;
begin
  rewrite(Msg,stderr);
  write(Msg,'MPI_ARI.P: !!! ERROR !!! ');

  case (ErrNo) of
    div_by_zero  : writeln(Msg,'Division by 0 in  MPI / MPI !');
    sqrt_range   : writeln(Msg,'Range error in sqrt(MPI)!');
    ln_range     : writeln(Msg,'Range error in ln(MPI)!');
    loga_range   : writeln(Msg,'Range error in loga(MPI)!');
    log2_range   : writeln(Msg,'Range error in log2(MPI)!');
    log10_range  : writeln(Msg,'Range error in log10(MPI)!');
    power_range  : writeln(Msg,'Range error in power(MPI,MPI)!');
    power_notdef : writeln(Msg,'0^0 not defined for power(MPI,MPI)!');
    arcsin_range : writeln(Msg,'Range error in arcsin(MPI)!');
    arccos_range : writeln(Msg,'Range error in arccos(MPI)!');
    sin_range    : writeln(Msg,'Range error in sin(MPI)!');
    cos_range    : writeln(Msg,'Range error in cos(MPI)!');
    tan_range    : writeln(Msg,'Range error in tan(MPI)!');
    cot_range    : writeln(Msg,'Range error in tan(MPI)!');
    coth_range   : writeln(Msg,'Range error in coth(MPI)!');
    arcosh_range : writeln(Msg,'Range error in arcosh(MPI)!');
    artanh_range : writeln(Msg,'Range error in artanh(MPI)!');
    arcoth_range : writeln(Msg,'Range error in arcoth(MPI)!');
    else         : writeln(Msg,'Unknown error!');
  end ;

  if (ErrNo in [div_by_zero,power_range,power_notdef]) then
    begin                              { Error message for two arguments }
      writeln(Msg,'First argument:');  writeln(Msg,Arg1);
      writeln(Msg,'Second argument:'); writeln(Msg,Arg2);
    end
  else                                 { Second argument is a dummy }
    begin
      writeln(Msg,'Argument: '); writeln(Msg,Arg1);
    end;

  traceback(Msg);
  exit(-1);
end;


{============================================================================}
{==========   Special Assignment Operators  =================================}
{============================================================================}
{= A special assignment operator is used to initialize the result variable of }
{= a function (InitFunc) or operator (InitOp) or to mark the result variable  }
{= to be temporary (TempFunc, TempOp). 'InitFunc' and 'InitOp' are synonyms.  }
{= 'TempFunc' and 'TempOp' are also synonyms. For compatibility to older      }
{= versions of the multi-precision modules there also exists an assignment    }
{= operator of type boolean. The value 'true' is equivalent to 'InitFunc' and }
{= 'InitOp', the value 'false' is equivalent to 'TempFunc' and 'TempOp'. For  }
{= definition of type 'mpmode' see module mp_ari.                             }
{=---------------------------------------------------------------------------}
global operator := ( var x : mpinterval; mode : mpmode );
begin
  case mode of
    InitFunc, InitOp : begin
                         mpinit(x.inf);
                         mpinit(x.sup);
                       end;
    TempFunc, TempOp : begin
                         mptemp(x.inf);
                         mptemp(x.sup);
                       end;
    else:  { No action }
  end;
end;

global operator := ( var x : mpinterval; mode : boolean );
begin
  if (mode) then
    begin
      mpinit(x.inf);
      mpinit(x.sup);
    end
  else
    begin
      mptemp(x.inf);
      mptemp(x.sup);
    end;
end;


{============================================================================}
{==========   Assignment Operators   ========================================}
{============================================================================}
global operator := ( var x : mpinterval;          { mpinterval := mpinterval }
                         y : mpinterval );        {--------------------------}
begin
  mpvlcp(y);
  x.inf := y.inf;
  x.sup := y.sup;
  mpfree(y);
end;

global operator := ( var x : mpinterval;            { mpinterval := interval }
                         y : interval );            {------------------------}
begin
  x.inf := y.inf;
  x.sup := y.sup;
end;

global operator := ( var x : mpinterval;              { mpinterval := mpreal }
                         y : mpreal );                {----------------------}
begin
  mpvlcp(y);
  x.inf := y;
  x.sup := y;
  mpfree(y);
end;

global operator := ( var x : mpinterval;                { mpinterval := real }
                         r : real );                    {--------------------}
begin
  x.inf := _mpreal(r,-1);
  x.sup := _mpreal(r,+1);
end;

global operator := ( var x : interval;              { interval := mpinterval }
                         y : mpinterval );          {------------------------}
begin
  mpvlcp(y);
  x.inf := _real(y.inf,-1);
  x.sup := _real(y.sup,+1);
  mpfree(y);
end;


{============================================================================}
{==========   Functions for Type Coercions   ================================}
{============================================================================}
global function _mpinterval( x : mpreal ) : mpinterval;    { mpreal-coercion }
begin                                                      {-----------------}
  mpvlcp(x);
  _mpinterval := InitFunc;
  _mpinterval := x;
  _mpinterval := TempFunc;
  mpfree(x);
end;

global function _mpinterval( x : interval ) : mpinterval;  { interval-coerc. }
begin                                                      {-----------------}
  _mpinterval := InitFunc;
  _mpinterval.inf := x.inf;
  _mpinterval.sup := x.sup;
  _mpinterval := TempFunc;
end;

global function _mpinterval( r : real ) : mpinterval;        { real-coercion }
begin                                                        {---------------}
  _mpinterval := InitFunc;
  _mpinterval := r;
  _mpinterval := TempFunc;
end;


{============================================================================}
{==========   Logical Operators: mpinterval <--> mpinterval   ===============}
{============================================================================}
global operator = ( x, y : mpinterval ) eq : boolean;             { Equality }
begin                                                             {----------}
  mpvlcp(x); mpvlcp(y);
  eq := (x.inf = y.inf) and (x.sup = y.sup);
  mpfree(x); mpfree(y);
end;

global operator <> ( x, y : mpinterval ) neq : boolean;         { Unequality }
begin                                                           {------------}
  mpvlcp(x); mpvlcp(y);
  neq := (x.inf <> y.inf) or (x.sup <> y.sup);
  mpfree(x); mpfree(y);
end;

global operator <= ( x, y : mpinterval ) sub : boolean;             { Subset }
begin                                                               {--------}
  mpvlcp(x); mpvlcp(y);
  sub := (y.inf <= x.inf) and (x.sup <= y.sup);
  mpfree(x); mpfree(y);
end;

global operator < ( x, y : mpinterval ) psub : boolean;      { Proper subset }
begin                                                        {---------------}
  mpvlcp(x); mpvlcp(y);
  psub := (x <= y) and (x <> y);
  mpfree(x); mpfree(y);
end;

global operator >= ( x, y : mpinterval ) sup : boolean;           { Superset }
begin                                                             {----------}
  mpvlcp(x); mpvlcp(y);
  sup := (x.inf <= y.inf) and (y.sup <= x.sup);
  mpfree(x); mpfree(y);
end;

global operator > ( x, y : mpinterval ) psup : boolean;    { Proper Superset }
begin                                                      {-----------------}
  mpvlcp(x); mpvlcp(y);
  psup := (x >= y) and (x <> y);
  mpfree(x); mpfree(y);
end;

global operator in ( x, y : mpinterval ) cont : boolean; { Contained-in test }
begin                                                    {-------------------}
  mpvlcp(x); mpvlcp(y);
  cont := (y.inf < x.inf) and (x.sup < y.sup);
  mpfree(x); mpfree(y);
end;

global operator >< ( x, y : mpinterval ) disj : boolean;    { Disjointedness }
begin                                                       {----------------}
  mpvlcp(x); mpvlcp(y);
  disj := (x.sup < y.inf) or (y.sup < x.inf);
  mpfree(x); mpfree(y);
end;

global function point( x : mpinterval ) : boolean;     { Point interval test }
begin                                                  {---------------------}
  mpvlcp(x);
  point := (x.inf = x.sup);
  mpfree(x);
end;


{============================================================================}
{==========   Logical Operators: mpreal     <--> mpinterval   ===============}
{==========                      mpinterval <--> mpreal       ===============}
{============================================================================}
global operator = ( x : mpinterval; a : mpreal ) eq : boolean;
begin
  mpvlcp(x); mpvlcp(a);
  eq := (x.inf = a) and (x.sup = a);
  mpfree(a); mpfree(x);
end;

global operator = ( a : mpreal; x : mpinterval ) eq : boolean;
begin
  mpvlcp(x); mpvlcp(a);
  eq := (a = x.inf) and (a = x.sup);
  mpfree(a); mpfree(x);
end;

global operator in ( a : mpreal; x : mpinterval ) cont : boolean;
begin
  mpvlcp(x); mpvlcp(a);
  cont := (x.inf <= a) and (a <= x.sup);
  mpfree(a); mpfree(x);
end;


{============================================================================}
{==========   Logical Operators: real       <--> mpinterval   ===============}
{==========                      mpinterval <--> real         ===============}
{============================================================================}
global operator = ( x : mpinterval; r : real ) eq : boolean;
begin
  mpvlcp(x);
  eq := (x.inf = r) and (x.sup = r);
  mpfree(x);
end;

global operator = ( r : real; x : mpinterval ) eq : boolean;
begin
  mpvlcp(x);
  eq := (r = x.inf) and (r = x.sup);
  mpfree(x);
end;

global operator in ( r : real; x : mpinterval ) cont : boolean;
begin
  mpvlcp(x);
  cont := (x.inf <= r) and (r <= x.sup);
  mpfree(x);
end;


{============================================================================}
{==========   Logical Operators: interval   <--> mpinterval   ===============}
{==========                      mpinterval <--> interval     ===============}
{============================================================================}
global operator = ( x : mpinterval; a : interval ) eq : boolean;
begin
  mpvlcp(x);
  eq := (x.inf = a.inf) and (x.sup = a.sup);
  mpfree(x);
end;

global operator = ( a : interval; x : mpinterval ) eq : boolean;
begin
  mpvlcp(x);
  eq := (a.inf = x.inf) and (a.sup = x.sup);
  mpfree(x);
end;

global operator in ( a : interval; x : mpinterval ) cont : boolean;
begin
  mpvlcp(x);
  cont := (x.inf < _mpreal(a.inf,-1)) and (_mpreal(a.sup,+1) < x.sup);
  mpfree(x);
end;


{============================================================================}
{==========   Arithmetic Operators: mpinterval <--> mpinterval   ============}
{============================================================================}
global operator + ( x : mpinterval ) mplus : mpinterval;
begin
  mpvlcp(x);
  mplus := InitOp;
  mplus := x;
  mplus := TempOp;
  mpfree(x);
end;

global operator - ( x : mpinterval ) mminus : mpinterval;
begin
  mpvlcp(x);
  mminus := InitOp;
  mminus.inf := -x.sup;
  mminus.sup := -x.inf;
  mminus := TempOp;
  mpfree(x);
end;

global operator + ( x, y : mpinterval ) add : mpinterval;
begin
  mpvlcp(x); mpvlcp(y);
  add := InitOp;
  add.inf := x.inf +< y.inf;
  add.sup := x.sup +> y.sup;
  add := TempOp;
  mpfree(x); mpfree(y);
end;

global operator - ( x, y : mpinterval ) sub : mpinterval;
begin
  mpvlcp(x); mpvlcp(y);
  sub := InitOp;
  sub.inf := x.inf -< y.sup;
  sub.sup := x.sup -> y.inf;
  sub := TempOp;
  mpfree(x); mpfree(y);
end;

global operator * ( x, y : mpinterval ) mul : mpinterval;
{----------------------------------------------------------------------------}
{ For the multiplication of two intervals [z] = [x] * [y] the following      }
{ rules are applied:                                                         }
{                                                                            }
{   x.inf  x.sup  y.inf  y.sup  |       z.inf         |       z.sup          }
{   ----------------------------+---------------------+--------------------  }
{     +      +      +      +    |   x.inf * y.inf     |   x.sup * y.sup      }
{     +      +      -      +    |   x.sup * y.inf     |   x.sup * y.sup      }
{     +      +      -      -    |   x.sup * y.inf     |   x.inf * y.sup      }
{                               |                     |                      }
{     -      -      +      +    |   x.inf * y.sup     |   x.sup * y.inf      }
{     -      -      -      +    |   x.inf * y.sup     |   x.inf * y.inf      }
{     -      -      -      -    |   x.sup * y.sup     |   x.inf * y.inf      }
{                               |                     |                      }
{     -      +      +      +    |   x.inf * y.sup     |   x.sup * y.sup      }
{     -      +      -      -    |   x.sup * y.inf     |   x.inf * y.inf      }
{                               |                     |                      }
{     -      +      -      +    | min(x.inf * y.sup,  | max(x.inf * y.inf,   }
{                               |     x.sup * y.inf ) |     x.sup * y.sup )  }
{----------------------------------------------------------------------------}
begin
   mpvlcp(x);  mpvlcp(y);
   mul := InitOp;

   if (sign(x.inf) >= 0) then
     begin
       if (sign(y.inf) >= 0) then mul.inf := x.inf *< y.inf
                             else mul.inf := x.sup *< y.inf;
       if (sign(y.sup) >= 0) then mul.sup := x.sup *> y.sup
                             else mul.sup := x.inf *> y.sup;
     end
   else if (sign(x.sup) <= 0) then
     begin
       if (sign(y.sup) >= 0) then mul.inf := x.inf *< y.sup
                             else mul.inf := x.sup *< y.sup;
       if (sign(y.inf) >= 0) then mul.sup := x.sup *> y.inf
                             else mul.sup := x.inf *> y.inf;
     end
   else if (sign(y.inf) >= 0) then
     begin
       mul.inf := x.inf *< y.sup;
       mul.sup := x.sup *> y.sup;
     end
   else if (sign(y.sup) <= 0) then
     begin
       mul.inf := x.sup *< y.inf;
       mul.sup := x.inf *> y.inf;
     end
   else
     begin
       mul.inf := min( x.inf *< y.sup, x.sup *< y.inf);
       mul.sup := max( x.inf *> y.inf, x.sup *> y.sup);
     end;

   mul := TempOp;
   mpfree(x); mpfree(y);
end;

global operator / ( x, y : mpinterval ) frac : mpinterval;
{----------------------------------------------------------------------------}
{ Suppose that 0 is not contained in [y]. For an interval division it holds  }
{ [z] = [x] / [y] = [x] * [1/y.sup,1/y.inf]. Applying the rules for multi-   }
{ plication the following rules can be derived:                              }
{                                                                            }
{        y.inf, y.sup   x.inf  x.sup  |     z.inf     |     z.sup            }
{       ------------------------------+---------------+---------------       }
{             +           +      +    | x.inf / y.sup | x.sup / y.inf        }
{             +           -      +    | x.inf / y.inf | x.sup / y.inf        }
{             +           -      -    | x.inf / y.inf | x.sup / y.sup        }
{                                     |               |                      }
{             -           +      +    | x.sup / y.sup | x.inf / y.inf        }
{             -           -      +    | x.sup / y.sup | x.inf / y.sup        }
{             -           -      -    | x.sup / y.inf | x.inf / y.sup        }
{----------------------------------------------------------------------------}
begin
  mpvlcp(x); mpvlcp(y);
  frac := InitOp;

  if (sign(y.inf) > 0) then
    begin
      if (sign(x.inf) >= 0) then frac.inf := x.inf /< y.sup
                            else frac.inf := x.inf /< y.inf;
      if (sign(x.sup) >= 0) then frac.sup := x.sup /> y.inf
                            else frac.sup := x.sup /> y.sup;
    end
  else if (sign(y.sup) < 0) then
    begin
      if (sign(x.sup) >= 0) then frac.inf := x.sup /< y.sup
                            else frac.inf := x.sup /< y.inf;
      if (sign(x.inf) >= 0) then frac.sup := x.inf /> y.inf
                            else frac.sup := x.inf /> y.sup;
    end
  else
    _math_message(div_by_zero,x,y);

  frac := TempOp;
  mpfree(x); mpfree(y);
end;


{============================================================================}
{==========   Arithmetic Operators: mpreal     <--> mpinterval   ============}
{==========                         mpinterval <--> mpreal       ============}
{============================================================================}
global operator + ( a : mpreal; x : mpinterval ) add : mpinterval;
begin
  mpvlcp(x); mpvlcp(a);
  add := InitOp;
  add := _mpinterval(a) + x;
  add := TempOp;
  mpfree(a); mpfree(x);
end;

global operator + ( x : mpinterval ; a : mpreal ) add : mpinterval;
begin
  mpvlcp(x); mpvlcp(a);
  add := InitOp;
  add := x + _mpinterval(a);
  add := TempOp;
  mpfree(a); mpfree(x);
end;

global operator - ( a : mpreal; x : mpinterval ) sub : mpinterval;
begin
  mpvlcp(x); mpvlcp(a);
  sub := InitOp;
  sub := _mpinterval(a) - x;
  sub := TempOp;
  mpfree(a); mpfree(x);
end;

global operator - ( x : mpinterval ; a : mpreal ) sub : mpinterval;
begin
  mpvlcp(x); mpvlcp(a);
  sub := InitOp;
  sub := x - _mpinterval(a);
  sub := TempOp;
  mpfree(a); mpfree(x);
end;

global operator * ( a : mpreal; x : mpinterval ) mul : mpinterval;
begin
  mpvlcp(x); mpvlcp(a);
  mul := InitOp;
  mul := _mpinterval(a) * x;
  mul := TempOp;
  mpfree(a); mpfree(x);
end;

global operator * ( x : mpinterval ; a : mpreal ) mul : mpinterval;
begin
  mpvlcp(x); mpvlcp(a);
  mul := InitOp;
  mul := x * _mpinterval(a);
  mul := TempOp;
  mpfree(a); mpfree(x);
end;

global operator / ( a : mpreal; x : mpinterval ) frac : mpinterval;
begin
  mpvlcp(x); mpvlcp(a);
  frac := InitOp;
  frac := _mpinterval(a) / x;
  frac := TempOp;
  mpfree(a); mpfree(x);
end;

global operator / ( x : mpinterval ; a : mpreal ) frac : mpinterval;
begin
  mpvlcp(x); mpvlcp(a);
  frac := InitOp;
  frac := x / _mpinterval(a);
  frac := TempOp;
  mpfree(a); mpfree(x);
end;


{============================================================================}
{==========   Arithmetic Operators: real       <--> mpinterval   ============}
{==========                         mpinterval <--> real         ============}
{============================================================================}
global operator + ( r : real; x : mpinterval ) add : mpinterval;
begin
  mpvlcp(x);
  add := InitOp;
  add := _mpinterval(r) + x;
  add := TempOp;
  mpfree(x);
end;

global operator + ( x : mpinterval; r : real ) add : mpinterval;
begin
  mpvlcp(x);
  add := InitOp;
  add := x + _mpinterval(r);
  add := TempOp;
  mpfree(x);
end;

global operator - ( r : real; x : mpinterval ) sub : mpinterval;
begin
  mpvlcp(x);
  sub := InitOp;
  sub := _mpinterval(r) - x;
  sub := TempOp;
  mpfree(x);
end;

global operator - ( x : mpinterval; r : real ) sub : mpinterval;
begin
  mpvlcp(x);
  sub := InitOp;
  sub := x - _mpinterval(r);
  sub := TempOp;
  mpfree(x);
end;

global operator * ( r : real; x : mpinterval ) mul : mpinterval;
begin
  mpvlcp(x);
  mul := InitOp;
  mul := _mpinterval(r) * x;
  mul := TempOp;
  mpfree(x);
end;

global operator * ( x : mpinterval; r : real ) mul : mpinterval;
begin
  mpvlcp(x);
  mul := InitOp;
  mul := x * _mpinterval(r);
  mul := TempOp;
  mpfree(x);
end;

global operator / ( r : real; x : mpinterval ) frac : mpinterval;
begin
  mpvlcp(x);
  frac := InitOp;
  frac := _mpinterval(r) / x;
  frac := TempOp;
  mpfree(x);
end;

global operator / ( x : mpinterval; r : real ) frac : mpinterval;
begin
  mpvlcp(x);
  frac := InitOp;
  frac := x / _mpinterval(r);
  frac := TempOp;
  mpfree(x);
end;


{============================================================================}
{==========   Arithmetic Operators: interval   <--> mpinterval   ============}
{==========                         mpinterval <--> interval     ============}
{============================================================================}
global operator + ( a : interval; x : mpinterval ) add : mpinterval;
begin
  mpvlcp(x);
  add := InitOp;
  add := _mpinterval(a) + x;
  add := TempOp;
  mpfree(x);
end;

global operator + ( x : mpinterval ; a : interval ) add : mpinterval;
begin
  mpvlcp(x);
  add := InitOp;
  add := x + _mpinterval(a);
  add := TempOp;
  mpfree(x);
end;

global operator - ( a : interval; x : mpinterval ) sub : mpinterval;
begin
  mpvlcp(x);
  sub := InitOp;
  sub := _mpinterval(a) - x;
  sub := TempOp;
  mpfree(x);
end;

global operator - ( x : mpinterval ; a : interval ) sub : mpinterval;
begin
  mpvlcp(x);
  sub := InitOp;
  sub := x - _mpinterval(a);
  sub := TempOp;
  mpfree(x);
end;

global operator * ( a : interval; x : mpinterval ) mul : mpinterval;
begin
  mpvlcp(x);
  mul := InitOp;
  mul := _mpinterval(a) * x;
  mul := TempOp;
  mpfree(x);
end;

global operator * ( x : mpinterval ; a : interval ) mul : mpinterval;
begin
  mpvlcp(x);
  mul := InitOp;
  mul := x * _mpinterval(a);
  mul := TempOp;
  mpfree(x);
end;

global operator / ( a : interval; x : mpinterval ) frac : mpinterval;
begin
  mpvlcp(x);
  frac := InitOp;
  frac := _mpinterval(a) / x;
  frac := TempOp;
  mpfree(x);
end;

global operator / ( x : mpinterval ; a : interval ) frac : mpinterval;
begin
  mpvlcp(x);
  frac := InitOp;
  frac := x / _mpinterval(a);
  frac := TempOp;
  mpfree(x);
end;


{============================================================================}
{==========   Lattice Operators   ===========================================}
{============================================================================}
global operator +* ( x, y : mpinterval ) hull : mpinterval;    { Convex hull }
begin                                                          {-------------}
  mpvlcp(x); mpvlcp(y);
  hull := InitFunc;
  if (x.inf <= y.inf) then hull.inf := x.inf
                      else hull.inf := y.inf;
  if (x.sup >= y.sup) then hull.sup := x.sup
                      else hull.sup := y.sup;
  hull := TempFunc;
  mpfree(x); mpfree(y);
end;

global operator ** ( x, y : mpinterval ) intsec : mpinterval; { Intersection }
begin                                                         {--------------}
  mpvlcp(x); mpvlcp(y);
  intsec := InitFunc;
  if (x.inf >= y.inf) then intsec.inf:= x.inf
                      else intsec.inf:= y.inf;
  if (x.sup <= y.sup) then intsec.sup:= x.sup
                      else intsec.sup:= y.sup;
  intsec := TempFunc;
  mpfree(x); mpfree(y);
end;


{============================================================================}
{==========   Access and Transfer Functions   ===============================}
{============================================================================}
global function inf( x : mpinterval ) : mpreal;
begin
  mpvlcp(x);
  inf := InitFunc;
  inf := x.inf;
  inf := TempFunc;
  mpfree(x);
end;

global function sup( x : mpinterval ) : mpreal;
begin
  mpvlcp(x);
  sup := InitFunc;
  sup := x.sup;
  sup := TempFunc;
  mpfree(x);
end;

global function intval( lb, ub : mpreal) : mpinterval;
begin
  mpvlcp(lb); mpvlcp(ub);
  intval := InitFunc;
  intval.inf := lb;
  intval.sup := ub;
  intval := TempFunc;
  mpfree(lb); mpfree(ub);
end;


{============================================================================}
{==========   Utilities   ===================================================}
{============================================================================}
global procedure iround( value: mpreal; var lb, ub : mpreal );
  external b_irnd;     { Creates a lower and an upper bound of type 'mpreal' }
                       { according to the number of ulps indicated by the    }
                       { representation of the mpreal value.                 }

function is_integer( x : mpinterval ) : boolean;      { Check for an integer }
begin                                                 { point interval       }
  mpvlcp(x);                                          {----------------------}
  is_integer := point(x) and (x.inf = trunc_mp(x.inf));
  mpfree(x);
end;

{ Function increasing() is used to compute an enclosure of a monotonically   }
{ increasing function f on the interval [x]. For [z] := f([x]) it holds      }
{ z.inf := inf( f(x.inf) ) and z.sup := sup( f(x.sup) ).                     }
{----------------------------------------------------------------------------}
function increasing( var x : mpinterval;
                     function f(y : mpreal) : mpreal ) : mpinterval;
var
  lb, ub : mpreal;
begin
  mpinit(lb); mpinit(ub);
  increasing := InitFunc;

  iround(f(x.inf),lb,ub);
  increasing.inf := lb;

  if point(x) then
    increasing.sup := ub
  else
    begin
      iround(f(x.sup),lb,ub);
      increasing.sup := ub;
    end;

  increasing := TempFunc;
  mpfree(lb); mpfree(ub);
end;

{ Function decreasing() is used to compute an enclosure of a monotonically   }
{ decreasing function f on the interval [x]. For [z] := f([x]) it holds      }
{ z.inf := inf( f(x.sup) ) and z.sup := sup( f(x.inf) ).                     }
{----------------------------------------------------------------------------}
function decreasing( var x : mpinterval;
                     function f(y : mpreal) : mpreal ) : mpinterval;
var
  lb, ub : mpreal;
begin
  mpinit(lb); mpinit(ub);
  decreasing := InitFunc;

  iround(f(x.sup),lb,ub);
  decreasing.inf := lb;

  if point(x) then
    decreasing.sup := ub
  else
    begin
      iround(f(x.inf),lb,ub);
      decreasing.sup := ub;
    end;

  decreasing := TempFunc;
  mpfree(lb); mpfree(ub);
end;


{============================================================================}
{==========   Standard Functions   ==========================================}
{============================================================================}

{----------------------------------------------------------------------------}
{ Local functions for computing an interval enclosure for point arguments    }
{----------------------------------------------------------------------------}
function pt_loga( x, a : mpreal ) : mpinterval;      { Logarithm, base a for }
var                                                  { point arguments       }
  lb, ub         : mpreal;                           {-----------------------}
  x1, x2, z1, z2 : mpreal;
  exact          : boolean;
  j              : integer;
begin
  mpvlcp(x);  mpvlcp(a);  mpinit(lb); mpinit(ub);
  mpinit(x1); mpinit(x2); mpinit(z1); mpinit(z2);
  pt_loga := InitFunc;

  { Check if x = a^j with an integer j }
  {------------------------------------}
  if (x = 1) then
    begin
      j := 0;
      exact := true;
    end
  else if (a < 1) then
    if (x < 1) then
        if (a < x) then
          exact := false
        else
          begin
            z1 := a;
            z2 := a;
            j  := 1;
            while (z1 > x) and (z1 = z2) do
            begin
              z2 := z2 *> a;
              z1 := z1 *< a;
              j  := succ(j);
            end;
            exact := (z1 = z2) and (z1 = x);
          end
    else  { a < 1 and x > 1 }
      begin
        x1 := 1 /< x;
        x2 := 1 /> x;
        z1 := a;
        z2 := a;
        j := -1;
        if (x1 = x2) then
          begin
            while (z1 > x1) and (z1 = z2) do
            begin
              z1 := z1 *< a;
              z2 := z2 *> a;
              j  := pred(j);
            end;
            exact := (z1 = z2) and (z1 = x1);
          end
        else
          exact := false;
      end
  else { a > 1 }
    if (x < 1) then
      begin
        z1 := 1 /< a;
        z2 := 1 /> a;
        j  := -1;
        while (z1 > x) and (z1 = z2) do
        begin
          z1 := z1 /< a;
          z2 := z2 /> a;
          j  := pred(j);
        end;
        exact := (z1 = z2) and (z1 = x);
      end
    else if (a > x) then
      exact := false
    else
      begin { a < x and  a, x > 1 }
        x1 := x;
        x2 := x;
        j  := 1;
        while (a < x1) and (x1 = x2) do
        begin
          x1 := x1 /< a;
          x2 := x2 /> a;
          j  := succ(j);
        end;
        exact := (x1 = x2) and (x1 = a);
      end;

  if exact then { x = a^j }
    begin
      pt_loga.inf := j;
      pt_loga.sup := j;
    end
  else
    begin
      iround(loga(x,a),lb,ub);
      pt_loga.inf := lb;
      pt_loga.sup := ub;
    end;

  pt_loga := TempFunc;
  mpfree(x);  mpfree(a);  mpfree(lb); mpfree(ub);
  mpfree(x1); mpfree(x2); mpfree(z1); mpfree(z2);
end;  {pt_loga}

function pt_power( x, y : mpreal) : mpinterval;    { x^y for point arguments }
var                                                {-------------------------}
  lb, ub : mpreal;
begin
  mpvlcp(x); mpvlcp(y); mpinit(lb); mpinit(ub);
  iround(power(x,y),lb,ub);
  pt_power := InitFunc;
  pt_power.inf := lb;
  pt_power.sup := ub;
  pt_power:= TempFunc;
  mpfree(x); mpfree(y); mpfree(lb); mpfree(ub);
end;

{----------------------------------------------------------------------------}
{ Arg(x,y) computes an enclosure for the angle of the point p = (x,y) to the }
{ x-axis. For the origin p = (0,0) the interval [-Pi,+Pi] is returned.       }
{----------------------------------------------------------------------------}
global function arctan( x : mpinterval ) : mpinterval;
  forward;

function pt_arg( x, y : mpreal ) : mpinterval;       { argument function for }
var                                                  { point parameters      }
  long_pi : mpinterval;                              {-----------------------}
begin
  mpvlcp(x); mpvlcp(y); mpinit(long_pi);
  pt_arg := InitFunc;

  if (x = 0) and (y = 0) then   { Return [-Pi,+Pi] for the origin p = (0,0) }
    begin
      long_pi := 4*arctan(ione);
      pt_arg.inf:= -long_pi.sup;
      pt_arg.sup:=  long_pi.sup;
    end

  else if (x > 0) then          { p lies in the right halve plane ...     }
    if (y = 0) then             { ... and is part of the pos. real axis   }
      pt_arg := 0
    else                        { .. and lies in quadrant I or IV         }
      pt_arg := arctan(_mpinterval(y) / _mpinterval(x))

  else if (x = 0) then          { p lies not in the right halve plane ... }
    if (y > 0) then             { ... and lies on the positive imag. axis }
      pt_arg :=  2*long_pi
    else                        { ... and lies on the negative imag. axis }
      pt_arg := -2*long_pi
  else {x < 0}                  { p in the left halve plane ...           }
    if (y = 0) then
      pt_arg := long_pi
    else if (y >= 0) then       { ... and lies in quadrant II             }
      pt_arg := arctan(_mpinterval(y) / _mpinterval(x)) + long_pi
    else                        { ... and lies in quadrant III            }
      pt_arg := arctan(_mpinterval(y) / _mpinterval(x)) - long_pi;

  pt_arg := TempFunc;
  mpfree(x); mpfree(y); mpfree(long_pi);
end; {pt_arg}


{----------------------------------------------------------------------------}
{ Special version of the above arg(.,.) function for a point p = (x,y) in    }
{ quadrant III. In contrast to the above implementation the angle is not     }
{ normalized between -Pi and +Pi. The function returns an enclosure of an    }
{ angle > Pi.                                                                }
{----------------------------------------------------------------------------}
function special_arg( x, y : mpreal ) : mpinterval;
var
  long_pi : mpinterval;
begin
  mpvlcp(x); mpvlcp(y); mpinit(long_pi);
  special_arg:= InitFunc;

  if (x = 0) then
    special_arg := 1.5 * long_pi     { p lies on the negative imaginary axis }
  else
    special_arg := arctan(_mpinterval(y) / _mpinterval(x)) + long_pi;

  special_arg := TempFunc;
  mpfree(x); mpfree(y); mpfree(long_pi);
end; {special_arg}


{-------------------}
{ String conversion }
{-------------------}
global function mpival( s : string ) : mpinterval;
begin
  mpival := InitFunc;
  mpival.inf := mpval(s,-1);
  mpival.sup := mpval(s,+1);
  mpival := TempFunc;
end;


{----------------------------------------}
{ Midpoint, diameter, and blow functions }
{----------------------------------------}
global function mid( x : mpinterval ) : mpreal;          { Interval midpoint }
begin                                                    {-------------------}
  mpvlcp(x);
  mid := InitFunc;
  mid := x.inf + 0.5*(x.sup-x.inf);
  mid := TempFunc;
  mpfree(x);
end; {mid}

global function diam ( x : mpinterval ) : real;     { Real interval diameter }
begin                                               {------------------------}
  mpvlcp(x);
  diam := x.sup -> x.inf;
  mpfree(x);
end; {diam}

global function blow( x : mpinterval ) : mpinterval;           { Blow 2 ULPs }
begin                                                          {-------------}
  mpvlcp(x);
  blow := InitFunc;
  blow.inf := pred(pred(x.inf));
  blow.sup := succ(succ(x.sup));
  blow := TempFunc;
  mpfree(x);
end; {blow}


{----------------------------------------}
{ Absolute value, square and square root }
{----------------------------------------}
global function abs( x : mpinterval ) : mpinterval;         { Absolute value }
begin                                                       {----------------}
  mpvlcp(x);
  abs := InitFunc;

  if (sign(x.inf) >= 0) then abs := x             { Interval positiv }
  else if (sign(x.sup) <= 0) then abs := -x       { Interval negativ }
  else
    begin                                         { 0 in interval    }
      abs.inf := 0;
      if (-x.inf > x.sup) then abs.sup := -x.inf
                          else abs.sup :=  x.sup;
    end;
  abs := TempFunc;
  mpfree(x);
end; {abs}

global function sqr( x : mpinterval ) : mpinterval;        { Square function }
begin                                                      {-----------------}
  mpvlcp(x);
  sqr := InitFunc;

  if (sign(x.inf) >= 0) then
    begin                           { Interval positiv }
      sqr.inf := x.inf *< x.inf;
      sqr.sup := x.sup *> x.sup;
    end
  else if (sign(x.sup) <= 0) then
    begin                           { Interval negativ }
      sqr.inf := x.sup *< x.sup;
      sqr.sup := x.inf *> x.inf;
    end
  else
    begin                           { 0 in interval    }
      sqr.inf := 0;
      if (-x.inf > x.sup) then
        sqr.sup := x.inf *> x.inf
      else
        sqr.sup := x.sup *> x.sup;
    end;

  sqr := TempFunc;
  mpfree(x);
end; {sqr}

global function sqrt( x : mpinterval ) : mpinterval;           { Square root }
begin                                                          {-------------}
  mpvlcp(x);
  if (x.inf < 0) then _math_message(sqrt_range,x,x);
  sqrt := InitFunc;
  sqrt := increasing(x,sqrt);
  sqrt := TempFunc;
  mpfree(x);
end; {sqrt}


{-----------------------}
{ Logarithmic functions }
{-----------------------}
global function ln( x : mpinterval ) : mpinterval;       { Natural logarithm }
begin                                                    {-------------------}
  mpvlcp(x);
  if (x.inf <= 0) then _math_message(ln_range,x,x);
  ln := InitFunc;
  ln := increasing(x,ln);
  ln := TempFunc;
  mpfree(x);
end; {ln}

global function log2( x : mpinterval ) : mpinterval;     { Logarithm, base 2 }
begin                                                    {-------------------}
  mpvlcp(x);
  if (x.inf <= 0) then _math_message(log2_range,x,x);
  log2 := InitFunc;
  log2 := increasing(x,log2);
  log2 := TempFunc;
  mpfree(x);
end; {log2}

global function log10( x : mpinterval ): mpinterval;    { Logarithm, base 10 }
begin                                                   {--------------------}
  mpvlcp(x);
  if (x.inf <= 0) then _math_message(log10_range,x,x);
  log10 := InitFunc;
  log10 := increasing(x,log10);
  log10 := TempFunc;
  mpfree(x);
end; {log10}

global function loga( x, a : mpinterval ) : mpinterval;  { Logarithm, base a }
begin                                                    {-------------------}
  mpvlcp(x); mpvlcp(a);
  loga := InitFunc;

  if (x.inf <= 0) or (a.inf <= 0) or (1 in a) then
    _math_message(loga_range,x,x)
  else
    if point(x) and point(a) then
      loga := pt_loga(x.inf,a.inf)
    else if (a.sup < 1) then
      if (x.sup <= 1) then
        begin
          loga.inf := inf( pt_loga(x.inf,a.sup) );
          loga.sup := sup( pt_loga(x.sup,a.inf) );
        end
      else if (x.inf >= 1) then
        begin
          loga.inf := inf( pt_loga(x.sup,a.sup) );
          loga.sup := sup( pt_loga(x.inf,a.inf) );
        end
      else
        begin
          loga.inf := inf( pt_loga(x.sup,a.sup) );
          loga.sup := sup( pt_loga(x.sup,a.inf) );
        end
    else if (x.sup <= 1) then
      begin
        loga.inf := inf( pt_loga(x.inf,a.inf) );
        loga.sup := sup( pt_loga(x.sup,a.sup) );
      end
    else if (x.inf >= 1) then
      begin
        loga.inf := inf( pt_loga(x.sup,a.inf) );
        loga.sup := sup( pt_loga(x.inf,a.sup) );
      end
    else
      begin
        loga.inf := inf( pt_loga(x.inf,a.inf) );
        loga.sup := sup( pt_loga(x.inf,a.sup) );
      end;

  loga := TempFunc;
  mpfree(x); mpfree(a);
end; {loga}


{-----------------------}
{ Exponential functions }
{-----------------------}
global function exp( x : mpinterval ) : mpinterval;   { Exponential function }
begin                                                 {----------------------}
  mpvlcp(x);
  exp := InitFunc;
  exp := increasing(x,exp);
  exp := TempFunc;
  mpfree(x);
end; {exp}

global function exp2( x : mpinterval ) : mpinterval;    { Power fct., base 2 }
begin                                                   {--------------------}
  mpvlcp(x);
  exp2 := InitFunc;
  exp2 := increasing(x,exp2);
  exp2 := TempFunc;
  mpfree(x);
end; {exp2}

global function exp10( x : mpinterval ) : mpinterval;  { Power fct., base 10 }
begin                                                  {---------------------}
  mpvlcp(x);
  exp10 := InitFunc;
  exp10 := increasing(x,exp10);
  exp10 := TempFunc;
  mpfree(x);
end; {exp10}

{----------------------------------------------------------------------------}
{ For computing the power function [z] = [x]^[y] for interval arguments [x]  }
{ and [y] the following cases are to be distinguished:                       }
{                                                                            }
{ 1. [y] = 0:                                                                }
{      0 in [x]  ==>  Error: contains 0^0, not defined!                      }
{      else      ==>  [z] := 1                                               }
{                                                                            }
{ 2. x.inf < 0:                                                              }
{      y is an integer point interval:                                       }
{        y <= -1:                                                            }
{          x.sup >= 0:         ==> Error: contains 0^0, not defined!         }
{          x.sup <  0:                                                       }
{            y is odd:         ==> [z] := [ x.sup^y , x.inf^y ]              }
{            y is even:        ==> [z] := [ x.inf^y , x.sup^y ]              }
{        y >= +1:                                                            }
{          y is odd:           ==> [z] := [ x.inf^y , x.sup^y ]              }
{          y is even:                                                        }
{            x.sup < 0:        ==> [z] := [ x.sup^y , x.inf^y ]              }
{            x.sup >= 0:       ==> [z] := [ 0, max(|x.inf|,x.sup)^y ]        }
{      y is not an integer point interval:                                   }
{                              ==> Error: Illegal range!                     }
{                                                                            }
{ 3. x.inf = 0:                                                              }
{      y.inf <= 0:             ==> Error: contains 0^0, not defined!         }
{      y.inf > 0:                                                            }
{        x.sup <= 1:           ==> [z] := [ 0, x.sup^y.inf ]                 }
{        x.sup > 1:            ==> [z] := [ 0, x.sup^y.sup ]                 }
{                                                                            }
{ 4. x.inf > 0:                                                              }
{      x.sup <= 1:                                                           }
{        y.sup <= 0:         ==> [z] := [ x.sup^y.sup , x.inf^y.inf ]        }
{        y.inf >= 0:         ==> [z] := [ x.inf^y.sup , x.sup^y.inf ]        }
{        y.inf < 0 < y.sup:  ==> [z] := [ x.inf^y.sup , x.inf^y.inf ]        }
{      x.inf >= 1:                                                           }
{        y.sup <= 0:         ==> [z] := [ x.sup^y.inf , x.inf^y.sup ]        }
{        y.inf >= 0:         ==> [z] := [ x.inf^y.inf , x.sup^y.sup ]        }
{        y.inf < 0 < y.sup:  ==> [z] := [ x.sup^y.inf , x.sup^y.sup ]        }
{      x.inf < 0 < x.sup:                                                    }
{        y.sup <= 0:         ==> [z] := [ x.sup^y.inf , x.inf^y.inf ]        }
{        y.inf >= 0:         ==> [z] := [ x.inf^y.sup , x.sup^y.sup ]        }
{        y.inf < 0 < y.sup:  ==> [z] := [ min(x.inf^y.sup,x.sup^y.inf),      }
{                                         max(x.inf^y.inf,x.sup^y.sup) ]     }
{                                                                            }
{----------------------------------------------------------------------------}
global function power( x, y : mpinterval ) : mpinterval;
begin
  mpvlcp(x); mpvlcp(y);
  power := InitFunc;

  if point(y) and (y.inf = 0) then       { Zero exponent inteval }
    if (0 in x) then                     {-----------------------}
      _math_message(power_notdef,x,y)    { 0^0 not defined! }
    else
      power := 1                         { [x]^0 := 1 }

  else if (x.inf < 0) then               { [x] contains negative values }
                                         {------------------------------}
    if is_integer(y) then                { y = [y] is an integer point  }
                                         { interval with y <> 0 ...     }
      if (y.inf <= -1) then              { ... and y is negative }
        if (x.sup >= 0) then
          _math_message(power_range,x,y)
        else {x.sup < 0}
          if odd_mp(y.inf) then   { y is an odd integer  }
            begin
              power.inf := inf(pt_power(x.sup,y.inf));
              power.sup := sup(pt_power(x.inf,y.inf));
            end
          else                    { y is an even integer }
            begin
              power.inf := inf(pt_power(x.inf,y.inf));
              power.sup := sup(pt_power(x.sup,y.sup));
            end

      else { y >= 1 }                    { ... and y is positive }

        if odd_mp(y.inf) then     { y is an odd integer }
          begin
            power.inf := inf(pt_power(x.inf,y.inf));
            power.sup := sup(pt_power(x.sup,y.inf));
          end
        else                      { y is an even integer }
          if (x.sup < 0) then
            begin
              power.inf := inf(pt_power(x.sup,y.inf));
              power.sup := sup(pt_power(x.inf,y.inf));
            end
          else {x.sup >= 0}
            begin
              power.inf := 0;
              x.inf := -x.inf;    { x.inf := |x.inf| }

              if (x.inf > x.sup) then
                power.sup := sup(pt_power(x.inf,y.inf))
              else
                power.sup := sup(pt_power(x.sup,y.inf))
            end

    else                                 { [x] contains neg. values and y is }
      _math_message(power_range,x,y)     { not an integer point interval     }

  else if (x.inf = 0) then               { Lower bound of [x] is zero ... }
                                         {--------------------------------}
    if (y.inf <= 0) then
     _math_message(power_notdef,x,y)     { 0^0 not defined! }
    else {y.inf > 0}
      if (x.sup <= 1) then               { ... and [x] lies in [0,1] }
        begin
          power.inf := 0;
          power.sup := sup(pt_power(x.sup,y.inf));
        end
      else {1 < x.sup}                   { ... and 1 in the interior of [x] }
        begin
          power.inf := 0;
          power.sup := sup(pt_power(x.sup,y.sup));
        end

  else {0 < x.inf}                       { [x] is positive ... }
                                         {---------------------}
      if (x.sup <= 1) then               { ... and [x] lies in (0,1]  }
        if (y.sup <= 0) then      { negative exponents }
          begin
            power.inf := inf(pt_power(x.sup,y.sup));
            power.sup := sup(pt_power(x.inf,y.inf));
          end
        else if (0 <= y.inf) then { positive exponents }
          begin
            power.inf := inf(pt_power(x.inf,y.sup));
            power.sup := sup(pt_power(x.sup,y.inf));
          end
        else {y.inf < 0 < y.sup}  { pos. and neg. exponents }
          begin
            power.inf := inf(pt_power(x.inf,y.sup));
            power.sup := sup(pt_power(x.inf,y.inf));
          end
      else if (1 <= x.inf) then          { ... and [x] lies right of 1 }
        if (y.sup <= 0) then      { negative exponents }
          begin
            power.inf := inf(pt_power(x.sup,y.inf));
            power.sup := sup(pt_power(x.inf,y.sup));
          end
        else if (0 <= y.inf) then { positive exponents }
          begin
            power.inf := inf(pt_power(x.inf,y.inf));
            power.sup := sup(pt_power(x.sup,y.sup));
          end
        else {y.inf < 0 < y.sup}  { pos. and neg. exponents }
          begin
            power.inf := inf(pt_power(x.sup,y.inf));
            power.sup := sup(pt_power(x.sup,y.sup));
          end

      else {x.inf < 1 < x.sup}           { ... and 1 in the interior of [x] }
        if (y.sup <= 0) then      { negative exponents }
          begin
            power.inf := inf(pt_power(x.sup,y.inf));
            power.sup := sup(pt_power(x.inf,y.inf));
          end
        else if (0 <= y.inf) then { positive exponents }
          begin
            power.inf := inf(pt_power(x.inf,y.sup));
            power.sup := sup(pt_power(x.sup,y.sup));
          end
        else {y.inf < 0 < y.sup}  { pos. and neg. exponents }
          begin
            power.inf := min( inf(pt_power(x.inf,y.sup)),
                              inf(pt_power(x.sup,y.inf)) );
            power.sup := max( sup(pt_power(x.inf,y.inf)),
                              inf(pt_power(x.sup,y.sup)) );
          end;

  power := TempFunc;
  mpfree(x);  mpfree(y);
end; {power}

{----------------------------------------------------------------------------}
{ Arg(x,y) computes an enclosure for the angle of every point (x0,y0) with   }
{ x0 in x and y0 in y. The set p = ((x0,y0) : x0 in x, y0 in y) corresponds  }
{ to a rectangle in the complex plane.                                       }
{----------------------------------------------------------------------------}
global function arg( x, y : mpinterval ) : mpinterval;
var
  long_pi : mpinterval;
begin
  mpvlcp(x); mpvlcp(y); mpinit(long_pi);
  arg := InitFunc;

  if (point(x) and point(y)) then       { p is a point }
    arg := pt_arg(x.inf,y.inf )

  else if (y = 0) then                  { p lies on the real axis ... }
    if (x.inf > 0) then
      arg:= 0                           { ... and is positive }
    else if (x.sup < 0) then
      arg:= 4*arctan(ione)              { ... and is negative }
    else
      begin
        long_pi := 4*arctan(ione);      { ... and contains 0  }
        arg.inf := -long_pi.sup;
        arg.sup :=  long_pi.sup;
      end

  else if (x = 0) then                  { p lies on the imaginary axis ... }
    if (y.inf > 0) then
      arg:= 2*arctan(ione)              { ... and is positive }
    else if (y.sup < 0) then
      arg:= -2*arctan(ione)             { ... and is negative }
    else
      begin
        long_pi := 4*arctan(ione);      { ... and contains 0  }
        arg.inf := -long_pi.sup;
        arg.sup :=  long_pi.sup;
      end

  { p is a rectangle in the plane ... }
  {-----------------------------------}
  else if (x.inf >= 0) then             { consider right halve plane ... }
    if (y.inf >= 0) then
      begin                             { ... with p in quadrant I  }
        arg.inf := inf( pt_arg(x.sup,y.inf) );
        arg.sup := sup( pt_arg(x.inf,y.sup) );
      end
    else if (y.sup <= 0) then
      begin                             { ... with p in quadrant IV }
        arg.inf := inf( pt_arg(x.inf,y.inf) );
        arg.sup := sup( pt_arg(x.sup,y.sup) );
      end
    else
      begin                             { ... with p in quadrants I and IV }
        arg.inf := inf( pt_arg(x.inf,y.inf) );
        arg.sup := sup( pt_arg(x.inf,y.sup) );
      end

  else if (x.sup <= 0) then             { consider left halve plane ... }
    if (y.inf >= 0) then
      begin                             { ... with p in quadrant II  }
        arg.inf := inf( pt_arg(x.sup,y.sup) );
        arg.sup := sup( pt_arg(x.inf,y.inf) );
      end
    else if (y.sup <= 0) then
      begin                             { ... with p in quadrant III }
        arg.inf := inf( pt_arg(x.inf,y.sup) );
        arg.sup := sup( pt_arg(x.sup,y.inf) );
      end
    else
      begin                             { ... and p in quadrants II and III }
        arg.inf := inf( pt_arg(x.sup,y.sup) );
        arg.sup := sup( special_arg(x.sup,y.inf) );
      end

  else {x.inf < 0 < x.sup }             { consider both halve planes ... }
    if (y.inf >= 0) then
      begin                             { ... with p in quadrants I and II }
        arg.inf := inf( pt_arg(x.sup,y.inf) );
        arg.sup := sup( pt_arg(x.inf,y.inf) );
      end
    else if (y.sup <= 0) then
      begin                             { ... with p in quadrants III and IV }
        arg.inf := inf( pt_arg(x.inf,y.sup) );
        arg.sup := sup( pt_arg(x.sup,y.sup) );
      end
    else
      begin                             { ... with p in all quadrants }
        long_pi := 4*arctan(ione);
        arg.inf := -long_pi.sup;
        arg.sup :=  long_pi.sup;
      end;

  arg := TempFunc;
  mpfree(x); mpfree(y); mpfree(long_pi);
end; {arg}


{---------------------------------}
{ Inverse trigonometric functions }
{---------------------------------}
global function arcsin( x : mpinterval ) : mpinterval;            { Arc sine }
begin                                                             {----------}
  mpvlcp(x);
  if (x.inf < -1) or (x.sup > 1) then _math_message(arcsin_range,x,x);
  arcsin := InitFunc;
  arcsin := increasing(x,arcsin);
  arcsin := TempFunc;
  mpfree(x);
end; {arcsin}

global function arccos( x : mpinterval ) : mpinterval;          { arc cosine }
begin                                                           {------------}
  mpvlcp(x);
  if (x.inf < -1) or (x.sup > 1) then _math_message(arccos_range,x,x);
  arccos := InitFunc;
  arccos := decreasing(x,arccos);
  arccos := TempFunc;
  mpfree(x);
end; {arccos}

global function arctan( x : mpinterval ) : mpinterval;         { Arc tangent }
begin                                                          {-------------}
  mpvlcp(x);
  arctan := InitFunc;
  arctan := increasing(x,arctan);
  arctan := TempFunc;
  mpfree(x);
end; {arctan}


global function arccot( x : mpinterval ) : mpinterval;       { arc cotangent }
begin                                                        {---------------}
  mpvlcp(x);
  arccot := InitFunc;
  arccot := decreasing(x,arccot);
  arccot := TempFunc;
  mpfree(x);
end; {arccot}

global function arctan2( x, y : mpinterval ) : mpinterval;     { arctan(x/y) }
begin                                                          {-------------}
  mpvlcp(x); mpvlcp(y);
  arctan2 := InitFunc;
  arctan2 := arg(y,x);
  arctan2 := TempFunc;
  mpfree(x); mpfree(y);
end; {artan2}


{-------------------------}
{ Trigonometric functions }
{-------------------------}
global function cos( x : mpinterval ) : mpinterval;                 { Cosine }
var                                                                 {--------}
  even_pi, odd_pi                     : boolean;
  oldprec                             : integer;
  long_pi, lb, ub, lb1, ub1, k1, k2   : mpreal;
begin
  mpvlcp(x);   mpinit(long_pi);
  mpinit(lb);  mpinit(ub);
  mpinit(lb1); mpinit(ub1);
  mpinit(k1);  mpinit(k2);
  cos := InitFunc;

  if point(x) then
    iround(cos(x.inf),lb,ub)
  else
    begin
      oldprec := getprec; setprec(oldprec+2);
      long_pi := 4*arctan(one);
      if (diam(x) >= 2*long_pi) then
        begin
          lb := -1;
          ub := +1;
        end
      else
        begin
          { Check first if x contains an even and/or odd multiple of Pi }
          {-------------------------------------------------------------}
          k1 := trunc_mp(x.inf /< long_pi);
          k2 := trunc_mp(x.sup /> long_pi);

          if (0 in x) then
            begin
              even_pi := true;
              if (k1 <= -1) then
                odd_pi := true
              else if ((k1 = 0) and (k2 = 0)) then
                odd_pi := false
              else
                odd_pi := true;
            end
          else  {0 not in x}
            begin
              if (x.sup < 0) then
              begin
                lb :=  k1;     { swap sign and values }
                k1 := -k2;     { of k1 and k2         }
                k2 := -lb;
              end;

              if odd_mp(k1) then         { k1 is odd }
                if (k2 = k1) then        {-----------}
                  begin
                    even_pi := false;
                    odd_pi  := false;
                  end
                else if (k2-1 = k1) then
                  begin
                    even_pi := true;
                    odd_pi  := false;
                  end
                else
                  begin
                    even_pi := true;
                    odd_pi  := true;
                  end
              else                       { k1 is even }
                if (k2 = k1) then        {------------}
                  begin
                    even_pi := false;
                    odd_pi  := false;
                  end
                else if (k2-1 = k1) then
                  begin
                    even_pi := false;
                    odd_pi  := true;
                  end
                else
                  begin
                    even_pi := true;
                    odd_pi  := true;
                  end;
            end; { 0 not in x}

          { Now we known if x contains an even and/or add multiple of Pi }
          {--------------------------------------------------------------}
          if (even_pi and odd_pi) then
            begin
              lb := -1;
              ub := +1;
            end
          else
            begin
              iround(cos(x.inf),lb,ub);
              iround(cos(x.sup),lb1,ub1);
              if (even_pi and not odd_pi) then
                begin
                  if (lb > lb1) then lb := lb1;
                  ub := 1;
                end
              else if (not even_pi and odd_pi) then
                begin
                  if (ub < ub1) then ub := ub1;
                  lb := -1;
                end
              else
                begin
                  if (lb > lb1) then lb := lb1;
                  if (ub < ub1) then ub := ub1;
                end;
            end;
        end;

      setprec(oldprec);
    end;

  cos.inf := lb;
  cos.sup := ub;

  cos := TempFunc;
  mpfree(x);   mpfree(long_pi);
  mpfree(lb);  mpfree(ub);
  mpfree(lb1); mpfree(ub1);
  mpfree(k1);  mpfree(k2);
end;

global function sin( x : mpinterval ) : mpinterval;                   { Sine }
var                                                                   {------}
  lb, ub  : mpreal;
  half_pi : mpinterval;
begin
  mpvlcp(x); mpinit(lb); mpinit(ub); mpinit(half_pi);
  sin := InitFunc;

  if point(x) then
    begin
      iround(sin(x.inf),lb,ub);
      sin.inf := lb;
      sin.sup := ub;
    end
  else
    begin
      half_pi := 2*arctan(ione);
      sin := cos(x - half_pi);
    end;

  sin := TempFunc;
  mpfree(x); mpfree(lb); mpfree(ub); mpfree(half_pi);
end;

global function tan( x : mpinterval ) : mpinterval;                { Tangent }
var                                                                {---------}
  k1, k2, oldprec        : integer;
  lb, ub, half_pi, dummy : mpreal;
begin
  mpvlcp(x); mpinit(lb); mpinit(ub); mpinit(half_pi); mpinit(dummy);
  tan := InitFunc;

  if point(x) then
    iround(tan(x.inf),lb,ub)
  else
    begin
      oldprec := getprec; setprec(oldprec+2);
      half_pi := 2*arctan(one);
      if ( (diam(x) >= 2 * half_pi)    or
           (x.sup >  maxint*half_pi)   or
           (x.inf < -maxint*half_pi) ) then _math_message(tan_range,x,x);

      k1 := trunc(x.inf /< half_pi);
      k2 := trunc(x.sup /> half_pi);
      if (x.inf < 0) then k1 := k1 - 1;
      if (x.sup < 0) then k2 := k2 - 1;

      if (not odd(k1)) and (k1 <> k2) then _math_message(tan_range,x,x);

      iround(tan(x.inf),lb,dummy);
      iround(tan(x.sup),dummy,ub);

      setprec(oldprec);
    end;

  tan.inf:= lb;
  tan.sup:= ub;

  tan := TempFunc;
  mpfree(x); mpfree(lb); mpfree(ub); mpfree(half_pi); mpfree(dummy);
end;

global function cot( x : mpinterval ) : mpinterval;              { Cotangent }
var                                                              {-----------}
  k1, k2, oldprec        : integer;
  lb, ub, long_pi, dummy : mpreal;
begin
  mpvlcp(x); mpinit(lb); mpinit(ub); mpinit(long_pi); mpinit(dummy);
  cot := InitFunc;

  if point(x) then
    iround(cot(x.inf), lb, ub)
  else
    begin
      oldprec := getprec; setprec(oldprec+2);
      long_pi := 4*arctan(one);
      if ( (diam(x) >= long_pi)        or
           (x.sup >  maxint*long_pi)   or
           (x.inf < -maxint*long_pi) ) then _math_message(cot_range,x,x);

      k1 := trunc(x.inf /< long_pi);
      k2 := trunc(x.sup /> long_pi);
      if (x.inf < 0) then k1 := k1 - 1;
      if (x.sup < 0) then k2 := k2 - 1;

      if (k1 <> k2) then _math_message(cot_range,x,x);

      iround(cot(x.inf),dummy,ub);
      iround(cot(x.sup),lb,dummy);

      setprec(oldprec);
    end;

  cot.inf:= lb;
  cot.sup:= ub;

  cot := TempFunc;
  mpfree(x); mpfree(lb); mpfree(ub); mpfree(long_pi); mpfree(dummy);
end;

{----------------------}
{ Hyperbolic functions }
{----------------------}
global function sinh( x : mpinterval ) : mpinterval;       { Hyperbolic sine }
begin                                                      {-----------------}
  mpvlcp(x);
  sinh := InitFunc;
  sinh := increasing(x,sinh);
  sinh := TempFunc;
  mpfree(x);
end; {sinh}

global function cosh( x : mpinterval ) : mpinterval;     { Hyperbolic cosine }
var                                                      {-------------------}
  lb, ub : mpreal;
begin
  mpvlcp(x); mpinit(lb); mpinit(ub);
  cosh := InitFunc;

  if point(x) then
    begin
      iround(cosh(x.inf),lb,ub);
      cosh.inf := lb;
      cosh.sup := ub;
    end
  else if (x.inf >= 0) then
    begin
      iround(cosh(x.inf),lb,ub);
      cosh.inf := lb;
      iround(cosh(x.sup),lb,ub);
      cosh.sup := ub;
    end
  else if (x.sup <= 0) then
    begin
      iround(cosh(x.sup),lb,ub);
      cosh.inf := lb;
      iround(cosh(x.inf),lb,ub);
      cosh.sup := ub;
    end
  else {x.inf < 0 < x.sup}
    begin
      cosh.inf := 1;
      if (-x.inf > x.sup) then
        iround(cosh(x.inf),lb,ub)
      else
        iround(cosh(x.sup),lb,ub);
      cosh.sup := ub;
    end;

  cosh := TempFunc;
  mpfree(x); mpfree(lb); mpfree(ub);
end; {cosh}

global function tanh( x : mpinterval ) : mpinterval;    { Hyperbolic tangent }
begin                                                   {--------------------}
  mpvlcp(x);
  tanh := InitFunc;
  tanh := increasing(x,tanh);
  tanh := TempFunc;
  mpfree(x);
end; {tanh}

global function coth( x: mpinterval ) : mpinterval;   { Hyperbolic cotangent }
var                                                   {----------------------}
  lb, ub : mpreal;
begin
  mpvlcp(x); mpinit(lb); mpinit(ub);
  coth := InitFunc;

  if point(x) and (x.inf <> 0) then
    begin
      iround(coth(x.inf),lb,ub);
      coth.inf := lb;
      coth.sup := ub;
    end
  else if (x.sup < 0) then
    begin
      iround(coth(x.sup),lb,ub);
      coth.inf := lb;
      iround(coth(x.inf),lb,ub);
      coth.sup := ub;
    end
  else if (x.inf > 0) then
    begin
      iround(coth(x.sup),lb,ub);
      coth.inf := lb;
      iround(coth(x.inf),lb,ub);
      coth.sup := ub;
    end
  else  {0 in x}
    _math_message(coth_range,x,x);

  coth := TempFunc;
  mpfree(x); mpfree(lb); mpfree(ub);
end; {coth}


{------------------------------}
{ Inverse hyperbolic functions }
{------------------------------}
global function arsinh( x : mpinterval ) : mpinterval;      { Inv. hyp. sine }
begin                                                       {----------------}
  mpvlcp(x);
  arsinh := InitFunc;
  arsinh := increasing(x,arsinh);
  arsinh := TempFunc;
  mpvlcp(x);
end; {arsinh}

global function arcosh( x : mpinterval ) : mpinterval;    { Inv. hyp. cosine }
var                                                       {------------------}
  lb, ub : mpreal;
begin
  mpvlcp(x); mpinit(lb); mpinit(ub);
  if (x.inf < 1) then _math_message(arcosh_range,x,x);
  arcosh := InitFunc;
  arcosh := increasing(x,arcosh);
  arcosh := TempFunc;
  mpvlcp(x);
end; {arcosh}

global function artanh( x : mpinterval ) : mpinterval;   { Inv. hyp. tangent }
begin                                                    {-------------------}
  mpvlcp(x);
  if (x.inf <= -1) or (x.sup >= 1) then _math_message(artanh_range,x,x);
  artanh := InitFunc;
  artanh := increasing(x,artanh);
  artanh := TempFunc;
  mpvlcp(x);
end; {artanh}

global function arcoth( x : mpinterval ) : mpinterval; { Inv. hyp. cotangent }
begin                                                  {---------------------}
  mpvlcp(x);
  if not ((x.sup < -1) or (x.inf > 1)) then _math_message(arcoth_range,x,x);
  arcoth := InitFunc;
  arcoth := decreasing(x,arcoth);
  arcoth := TempFunc;
  mpvlcp(x);
end; {arcoth}


{============================================================================}
{==========   Module Initalization Part   ===================================}
{============================================================================}
begin
  mpinit(one); mpinit(ione);
  one  := 1;
  ione := 1;
end.
{==============  Last line of module mpi_ari: 2052  =========================}