(*==========================================================================*)
(*                                                                          *)
(*  PASCAL-XSC - MODUL  MPS_ARIL                                    990530  *)
(*                                                                          *)
(*  (c) Markus Neher    Institute of Applied Mathematics                    *)
(*                      D-76128 Karlsruhe University, Germany               *)
(*                      e-mail: markus.neher@math.uni-karlsruhe.de          *)
(*                                                                          *)
(*  MULTIPLE PRECISION STAGGERED CORRECTION ARITHMETIC                      *)
(*  FOR LINEAR ODES WITH ANALYTIC COEFFICIENT FUNCTIONS                     *)
(*                                                                          *)
(*==========================================================================*)

MODULE mps_aril;
                                                                              
USE i_ari , mv_ari , mvi_ari, mp_ari , mpi_ari , service;
 
(*--------------------------------------------------------------------------*)

GLOBAL CONST STAGPREC = 40;

GLOBAL TYPE STAGARRAY = GLOBAL ARRAY[1..STAGPREC] OF REAL;

            MULPREC = GLOBAL RECORD
               PREC : INTEGER;
               STAG : STAGARRAY;
               INT  : INTERVAL; END;

            MPVECTOR  = GLOBAL DYNAMIC ARRAY[*] OF MULPREC;
            MPMATRIX  = GLOBAL DYNAMIC ARRAY[*] OF MPVECTOR;


(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR := ( VAR VRX : STAGARRAY; K : INTEGER );

VAR I : INTEGER;

BEGIN
   VRX[1] := K;
   FOR I := 2 TO STAGPREC DO
       VRX[I] := 0;
END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR := ( VAR VRX : STAGARRAY; RY : REAL );

VAR I : INTEGER;

BEGIN
   VRX[1] := RY;
   FOR I := 2 TO STAGPREC DO
       VRX[I] := 0;
END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR := ( VAR MPX : MULPREC; K : INTEGER );

BEGIN
   MPX.PREC := 1;
   MPX.STAG := K;
   MPX.INT  := 0;
END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR := ( VAR MPX : MULPREC; K : REAL );

BEGIN
   MPX.PREC := 1;
   MPX.STAG := K;
   MPX.INT  := 0;
END;

(*--------------------------------------------------------------------------*)

GLOBAL FUNCTION MPS( K : INTEGER ) : MULPREC;

BEGIN
   MPS := K;
END;

(*--------------------------------------------------------------------------*)

GLOBAL FUNCTION MPS( K : REAL ) : MULPREC;

BEGIN
   MPS := K;
END;

(*--------------------------------------------------------------------------*)

GLOBAL FUNCTION ROUND( VAR MPX : MULPREC; PREC : INTEGER ) : MULPREC;

(*  RUNDUNG AUF GEWšNSCHTE GENAUIGKEIT  *)

VAR I   : INTEGER;
    MPZ : MULPREC;

BEGIN
   IF MPX.PREC <= PREC THEN
      MPZ := MPX
   ELSE
      BEGIN
      MPZ.PREC := PREC;
      FOR I := 1 TO PREC DO
          MPZ.STAG[I] := MPX.STAG[I];
      FOR I := PREC + 1 TO STAGPREC DO
          MPZ.STAG[I] := 0;
      MPZ.INT.INF := #<( FOR I := PREC + 1 TO MPX.PREC SUM( MPX.STAG[I] )
                         + MPX.INT.INF );
      MPZ.INT.SUP := #>( FOR I := PREC + 1 TO MPX.PREC SUM( MPX.STAG[I] )
                         + MPX.INT.SUP );
   END;

   ROUND := MPZ;

END;

(*--------------------------------------------------------------------------*)

GLOBAL FUNCTION CONVERT( VAR DPL, DPU : DOTPRECISION; PREC : INTEGER )
                       : MULPREC;

(*  SPEICHERT DEN INHALT DES DOTPRECISION-INTERVALS  [DPL,DPU]  *)
(*  ALS STAGGERED-CORRECTION-FORMAT DER LAENGE  PREC            *)

VAR I       : INTEGER;
    ZL , ZU : REAL;
    MPZ     : MULPREC;

BEGIN

   MPZ.STAG := 0.0;

   I:= 0;

   REPEAT

      I := I + 1;
      ZL:= #*( DPL );
      ZU:= #*( DPU );

      IF ( ZL = ZU ) OR ( ZL = PRED( ZU ) ) THEN BEGIN
         MPZ.STAG[I]:= ZL;
         DPL:= #( DPL - ZL );
         DPU:= #( DPU - ZL );
      END;

   UNTIL ( I = PREC ) OR ( ZL < PRED( ZU ) ) OR ( I = STAGPREC );

   WHILE ( ( I > 1 ) AND ( MPZ.STAG[I] = 0 ) ) DO
         I := I - 1;

   MPZ.PREC := I;
   MPZ.INT.INF := #<( DPL );
   MPZ.INT.SUP := #>( DPU );

   CONVERT := MPZ;

END;

(*--------------------------------------------------------------------------*)

GLOBAL FUNCTION CONVERT( VAR MPIX : MPINTERVAL; PREC : INTEGER ) : MULPREC;

(*  SPEICHERT DEN INHALT DES MP-INTERVALS  MPIX      *)
(*  ALS STAGGERED-CORRECTION-FORMAT DER LAENGE PREC  *)

VAR I       : INTEGER;
    ZL , ZU : REAL;
    MPZ     : MULPREC;

BEGIN

   MPZ.STAG := 0.0;

   I:= 0;

   REPEAT

      I := I + 1;
      ZL:= MPIX.INF;
      ZU:= MPIX.SUP;

      IF ( ZL = ZU ) OR ( ZL = PRED( ZU ) ) THEN BEGIN
         MPZ.STAG[I]:= ZL;
         MPIX := MPIX - ZL;
      END;

   UNTIL ( I = PREC ) OR ( ZL < PRED( ZU ) ) OR ( I = STAGPREC );

   WHILE ( ( I > 1 ) AND ( MPZ.STAG[I] = 0 ) ) DO
         I := I - 1;

   MPZ.PREC := I;
   MPZ.INT := MPIX;

   CONVERT := MPZ;

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR + ( VAR MPX , MPY : MULPREC ) RES: MULPREC;

(*  ADDITION  MPX + MPY                   *)
(*  GENAUIGKEIT DES ERGEBNISSES:          *)
(*  MAXIMALE GENAUIGKEIT VON MPX UND MPY  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPX.PREC SUM( MPX.STAG[I] )
             + FOR I := 1 TO MPY.PREC SUM( MPY.STAG[I] )
             + MPX.INT.INF + MPY.INT.INF );

   DPU := #( DPL + MPX.INT.SUP - MPX.INT.INF + MPY.INT.SUP - MPY.INT.INF );

   RES := CONVERT( DPL , DPU , MAX( MPX.PREC , MPY.PREC ) );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR - ( VAR MPX , MPY : MULPREC ) RES: MULPREC;

(*  SUBTRAKTION  MPX - MPY                *)
(*  GENAUIGKEIT DES ERGEBNISSES:          *)
(*  MAXIMALE GENAUIGKEIT VON MPX UND MPY  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPX.PREC SUM( MPX.STAG[I] )
             - FOR I := 1 TO MPY.PREC SUM( MPY.STAG[I] )
             + MPX.INT.INF - MPY.INT.SUP );

   DPU := #( DPL + MPX.INT.SUP - MPX.INT.INF + MPY.INT.SUP - MPY.INT.INF );

   RES := CONVERT( DPL , DPU , MAX( MPX.PREC , MPY.PREC ) );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR - ( VAR MPX : MULPREC ) RES: MULPREC;

(*  - MPX MONADISCH                       *)
(*  GENAUIGKEIT DES ERGEBNISSES:          *)
(*  GENAUIGKEIT VON MPX                   *)

VAR I         : INTEGER;

BEGIN
   
   RES.PREC := MPX.PREC;
   FOR I := 1 TO MPX.PREC DO
	RES.STAG[I] := -MPX.STAG[I];
   RES.INT.INF := -MPX.INT.SUP;
   RES.INT.SUP := -MPX.INT.INF;
   
END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *< ( VAR RX : INTEGER; VAR MPY : MULPREC )
                     RES: MULPREC;

(*  PRODUKT  RX * MPY                      *)
(*  GENAUIGKEIT DES ERGEBNISSES:           *)
(*  GENAUIGKEIT VON MPY WIRD NICHT ERH™HT  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPY.PREC SUM( RX * MPY.STAG[I] ) );
   DPU := DPL;

   IF RX > 0 THEN BEGIN
      DPL := #( DPL + RX * MPY.INT.INF);
      DPU := #( DPU + RX * MPY.INT.SUP);
      END
   ELSE BEGIN
      DPL := #( DPL + RX * MPY.INT.SUP);
      DPU := #( DPU + RX * MPY.INT.INF);
   END;

   RES := CONVERT( DPL , DPU , MPY.PREC );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *> ( VAR RX : INTEGER; VAR MPY : MULPREC )
                     RES: MULPREC;

(*  PRODUKT  RX * MPY                *)
(*  GENAUIGKEIT DES ERGEBNISSES:     *)
(*  GENAUIGKEIT VON MPY WIRD ERH™HT  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPY.PREC SUM( RX * MPY.STAG[I] ) );
   DPU := DPL;

   IF RX > 0 THEN BEGIN
      DPL := #( DPL + RX * MPY.INT.INF);
      DPU := #( DPU + RX * MPY.INT.SUP);
      END
   ELSE BEGIN
      DPL := #( DPL + RX * MPY.INT.SUP);
      DPU := #( DPU + RX * MPY.INT.INF);
   END;

   RES := CONVERT( DPL , DPU , MPY.PREC + 1 );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *< ( VAR MPY : MULPREC; VAR RX : INTEGER )
                     RES: MULPREC;

(*  PRODUKT  RX * MPY                      *)
(*  GENAUIGKEIT DES ERGEBNISSES:           *)
(*  GENAUIGKEIT VON MPY WIRD NICHT ERH™HT  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPY.PREC SUM( RX * MPY.STAG[I] ) );
   DPU := DPL;

   IF RX > 0 THEN BEGIN
      DPL := #( DPL + RX * MPY.INT.INF);
      DPU := #( DPU + RX * MPY.INT.SUP);
      END
   ELSE BEGIN
      DPL := #( DPL + RX * MPY.INT.SUP);
      DPU := #( DPU + RX * MPY.INT.INF);
   END;

   RES := CONVERT( DPL , DPU , MPY.PREC );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *> ( VAR MPY : MULPREC; VAR RX : INTEGER )
                     RES: MULPREC;

(*  PRODUKT  RX * MPY                *)
(*  GENAUIGKEIT DES ERGEBNISSES:     *)
(*  GENAUIGKEIT VON MPY WIRD ERH™HT  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPY.PREC SUM( RX * MPY.STAG[I] ) );
   DPU := DPL;

   IF RX > 0 THEN BEGIN
      DPL := #( DPL + RX * MPY.INT.INF);
      DPU := #( DPU + RX * MPY.INT.SUP);
      END
   ELSE BEGIN
      DPL := #( DPL + RX * MPY.INT.SUP);
      DPU := #( DPU + RX * MPY.INT.INF);
   END;

   RES := CONVERT( DPL , DPU , MPY.PREC + 1 );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *< ( VAR RX : REAL; VAR MPY : MULPREC )
                    RES: MULPREC;

(*  PRODUKT  RX * MPY                      *)
(*  GENAUIGKEIT DES ERGEBNISSES:           *)
(*  GENAUIGKEIT VON MPY WIRD NICHT ERH™HT  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPY.PREC SUM( RX * MPY.STAG[I] ) );
   DPU := DPL;

   IF RX > 0 THEN BEGIN
      DPL := #( DPL + RX * MPY.INT.INF);
      DPU := #( DPU + RX * MPY.INT.SUP);
      END
   ELSE BEGIN
      DPL := #( DPL + RX * MPY.INT.SUP);
      DPU := #( DPU + RX * MPY.INT.INF);
   END;

   RES := CONVERT( DPL , DPU , MPY.PREC );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *> ( VAR RX : REAL; VAR MPY : MULPREC )
                    RES: MULPREC;

(*  PRODUKT  RX * MPY                *)
(*  GENAUIGKEIT DES ERGEBNISSES:     *)
(*  GENAUIGKEIT VON MPY WIRD ERH™HT  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPY.PREC SUM( RX * MPY.STAG[I] ) );
   DPU := DPL;

   IF RX > 0 THEN BEGIN
      DPL := #( DPL + RX * MPY.INT.INF);
      DPU := #( DPU + RX * MPY.INT.SUP);
      END
   ELSE BEGIN
      DPL := #( DPL + RX * MPY.INT.SUP);
      DPU := #( DPU + RX * MPY.INT.INF);
   END;

   RES := CONVERT( DPL , DPU , MPY.PREC + 1 );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *< ( VAR MPY : MULPREC; VAR RX : REAL )
                    RES: MULPREC;

(*  PRODUKT  RX * MPY                      *)
(*  GENAUIGKEIT DES ERGEBNISSES:           *)
(*  GENAUIGKEIT VON MPY WIRD NICHT ERH™HT  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPY.PREC SUM( RX * MPY.STAG[I] ) );
   DPU := DPL;

   IF RX > 0 THEN BEGIN
      DPL := #( DPL + RX * MPY.INT.INF );
      DPU := #( DPU + RX * MPY.INT.SUP );
      END
   ELSE BEGIN
      DPL := #( DPL + RX * MPY.INT.SUP );
      DPU := #( DPU + RX * MPY.INT.INF );
   END;

   RES := CONVERT( DPL , DPU , MPY.PREC );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *> ( VAR MPY : MULPREC; VAR RX : REAL )
                    RES: MULPREC;

(*  PRODUKT  RX * MPY                *)
(*  GENAUIGKEIT DES ERGEBNISSES:     *)
(*  GENAUIGKEIT VON MPY WIRD ERH™HT  *)

VAR I         : INTEGER;
    DPL , DPU : DOTPRECISION;

BEGIN

   DPL := #( FOR I := 1 TO MPY.PREC SUM( RX * MPY.STAG[I] ) );
   DPU := DPL;

   IF RX > 0 THEN BEGIN
      DPL := #( DPL + RX * MPY.INT.INF );
      DPU := #( DPU + RX * MPY.INT.SUP );
      END
   ELSE BEGIN
      DPL := #( DPL + RX * MPY.INT.SUP );
      DPU := #( DPU + RX * MPY.INT.INF );
   END;

   RES := CONVERT( DPL , DPU , MPY.PREC + 1 );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *< ( VAR MPX , MPY : MULPREC ) RES: MULPREC;

(*  PRODUKT  MPX * MPY                    *)
(*  GENAUIGKEIT DES ERGEBNISSES:          *)
(*  MAXIMALE GENAUIGKEIT VON MPX UND MPY  *)

VAR I , J     : INTEGER;
    IX        : INTERVAL;
    DPL , DPU : DOTPRECISION;

PROCEDURE ADD_INT_TIMES_MP( VAR IX : INTERVAL; VAR MPY : MULPREC;
                            VAR DPL , DPU : DOTPRECISION );
(*  ADDIERT DAS PRODUKT VON INTERVALL IX UND STAGGERED-ANTEIL VON MPY  *)
(*  ZU DPL, DPU                                                        *)

VAR J : INTEGER;

BEGIN
   FOR J := 1 TO MPY.PREC DO
       IF MPY.STAG[J] >= 0 THEN BEGIN
          DPL := #( DPL + IX.INF * MPY.STAG[J] );
          DPU := #( DPU + IX.SUP * MPY.STAG[J] );
          END
       ELSE BEGIN
          DPL := #( DPL + IX.SUP * MPY.STAG[J] );
          DPU := #( DPU + IX.INF * MPY.STAG[J] );
       END;
END;

BEGIN (*  OPERATOR MPX * MPY  *)

   (*  PRODUKT DER REELLEN ANTEILE  *)

   DPL := #( FOR I := 1 TO MPX.PREC
                 SUM( FOR J := 1 TO MPY.PREC
                      SUM( MPX.STAG[I] * MPY.STAG[J] ) ) );
   DPU := DPL;

   (*  PRODUKT VON INTERVALL- MIT REELLEM ANTEIL  *)

   (*  MPX.INT * MPY.STAG  *)
   ADD_INT_TIMES_MP( MPX.INT , MPY , DPL , DPU );
   (*  MPY.INT * MPX.STAG  *)
   ADD_INT_TIMES_MP( MPY.INT , MPX , DPL , DPU );

   (*  PRODUKT DER INTERVALL-ANTEILE  *)
   IX := MPX.INT * MPY.INT;
   DPL := #( DPL + IX.INF );
   DPU := #( DPU + IX.SUP );

   RES := CONVERT( DPL , DPU , MAX( MPX.PREC , MPY.PREC ) );

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR *> ( VAR MPX , MPY : MULPREC ) RES: MULPREC;

(*  PRODUKT  MPX * MPY  MIT ERH™HTER GENAUIGKEIT  *)
(*  GENAUIGKEIT DES ERGEBNISSES:                  *)
(*  GENAUIGKEIT VON MPX PLUS GENAUIGKEIT VON MPY  *)

VAR I , J     : INTEGER;
    IX        : INTERVAL;
    DPL , DPU : DOTPRECISION;

PROCEDURE ADD_INT_TIMES_MP( VAR IX : INTERVAL; VAR MPY : MULPREC;
                          VAR DPL , DPU : DOTPRECISION );
(*  ADDIERT DAS PRODUKT VON INTERVALL IX UND STAGGERED-ANTEIL VON MPY  *)
(*  ZU DPL, DPU                                                        *)

VAR J : INTEGER;

BEGIN
   FOR J := 1 TO MPY.PREC DO
       IF MPY.STAG[J] >= 0 THEN BEGIN
          DPL := #( DPL + IX.INF * MPY.STAG[J] );
          DPU := #( DPU + IX.SUP * MPY.STAG[J] );
          END
       ELSE BEGIN
          DPL := #( DPL + IX.SUP * MPY.STAG[J] );
          DPU := #( DPU + IX.INF * MPY.STAG[J] );
       END;
END;

BEGIN (*  OPERATOR MPX ** MPY  *)

   (*  PRODUKT DER REELLEN ANTEILE  *)

   DPL := #( FOR I := 1 TO MPX.PREC
                 SUM( FOR J := 1 TO MPY.PREC
                      SUM( MPX.STAG[I] * MPY.STAG[J] ) ) );
   DPU := DPL;

   (*  PRODUKT VON INTERVALL- MIT REELLEM ANTEIL  *)

   (*  MPX.INT * MPY.STAG  *)
   ADD_INT_TIMES_MP( MPX.INT , MPY , DPL , DPU );
   (*  MPY.INT * MPX.STAG  *)
   ADD_INT_TIMES_MP( MPY.INT , MPX , DPL , DPU );

   (*  PRODUKT DER INTERVALL-ANTEILE  *)
   IX := MPX.INT * MPY.INT;
   DPL := #( DPL + IX.INF );
   DPU := #( DPU + IX.SUP );

   RES := CONVERT( DPL , DPU , MPX.PREC + MPY.PREC );

END;

(*--------------------------------------------------------------------------*)

PROCEDURE ADD_INT_TIMES_MP( VAR IX : INTERVAL; VAR MPY : MULPREC;
                          VAR DPL , DPU : DOTPRECISION );
(*  ADDIERT DAS PRODUKT VON INTERVALL IX UND STAGGERED-ANTEIL VON MPY  *)
(*  ZU DPL, DPU                                                        *)

VAR J : INTEGER;

BEGIN
   FOR J := 1 TO MPY.PREC DO
       IF MPY.STAG[J] >= 0 THEN BEGIN
          DPL := #( DPL + IX.INF * MPY.STAG[J] );
          DPU := #( DPU + IX.SUP * MPY.STAG[J] );
          END
       ELSE BEGIN
          DPL := #( DPL + IX.SUP * MPY.STAG[J] );
          DPU := #( DPU + IX.INF * MPY.STAG[J] );
       END;
END;

(*--------------------------------------------------------------------------*)

GLOBAL PROCEDURE ADD_MP_TIMES_MP( VAR DPX, DPY : DOTPRECISION;
                                  VAR MPX , MPY : MULPREC );

(*  [DPX,DPY] := [DPX,DPY] + MPX * MPY  *)

VAR I , J     : INTEGER;
    IX        : INTERVAL;
    DPL , DPU : DOTPRECISION;

BEGIN (*  PROCEDURE ADD_MP_TIMES_MP  *)

   (*  PRODUKT DER REELLEN ANTEILE  *)

   DPL := #( FOR I := 1 TO MPX.PREC
                 SUM( FOR J := 1 TO MPY.PREC
                      SUM( MPX.STAG[I] * MPY.STAG[J] ) ) );
   DPU := DPL;

   (*  PRODUKT VON INTERVALL- MIT REELLEM ANTEIL  *)

   (*  MPX.INT * MPY.STAG  *)
   ADD_INT_TIMES_MP( MPX.INT , MPY , DPL , DPU );
   (*  MPY.INT * MPX.STAG  *)
   ADD_INT_TIMES_MP( MPY.INT , MPX , DPL , DPU );

   (*  PRODUKT DER INTERVALL-ANTEILE  *)
   IX := MPX.INT * MPY.INT;
   DPL := #( DPL + IX.INF );
   DPU := #( DPU + IX.SUP );

   DPX := #(DPX + DPL);
   DPY := #(DPY + DPU);

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR ** ( VAR MPX: MULPREC; POT: INTEGER )
                     RES : MULPREC;

(*  EXAKTE LANGZAHLPOTENZ  MPX ** POT  *)

VAR I   : INTEGER;
    ERR : REAL;
    MPZ : MULPREC;

BEGIN

   IF ( POT < 0 ) THEN BEGIN
      WRITELN( 'Fehler bei der Berechnung einer Langzahlpotenz.' );
      WRITELN( 'Modul MPS_ARI, Operator **.' );
      ERR := 1.0 / 0;
   END;

   MPZ.STAG := 0.0;

   IF POT = 0 THEN BEGIN
      MPZ.PREC := 1;
      MPZ.STAG[1] := 1;
      MPZ.INT := 0;
      END
   ELSE BEGIN
       MPZ := MPX;
       FOR I := 1 TO POT - 1 DO
           MPZ := MPZ *> MPX;  (* EXAKTE MULTIPLIKATION *)
   END;

   RES := MPZ;

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR / ( VAR MPX: MULPREC; RY : REAL; )
                    RES : MULPREC;

(*  QUOTIENT  MPX / RY                   *)
(*  GENAUIGKEIT DES ERGEBNISSES:         *)
(*  GENAUIGKEIT VON MPX BLEIBT ERHALTEN  *)

VAR I    : INTEGER;
    IX   : INTERVAL;
    MPY  : MULPREC;
    DPX  : DOTPRECISION;

BEGIN
   DPX := #( FOR I := 1 TO MPX.PREC SUM( MPX.STAG[I] ) );

   MPY.PREC := MPX.PREC;

   MPY.STAG[1] := #*( DPX ) / RY;
   FOR I := 2 TO MPX.PREC DO BEGIN
       DPX := #( DPX - RY * MPY.STAG[I-1] );
       MPY.STAG[I] := #*( DPX ) / RY;
   END;
   DPX := #( DPX - RY * MPY.STAG[MPY.PREC] );
   
   MPY.INT := INTVAL( #<( DPX + MPX.INT.INF ) , #>( DPX + MPX.INT.SUP ) ) / RY;

   RES := MPY;

END;

(*--------------------------------------------------------------------------*)

GLOBAL OPERATOR / ( VAR MPX, MPY : MULPREC ) RES : MULPREC;

(*  QUOTIENT  MPX / MPY                  *)
(*  GENAUIGKEIT DES ERGEBNISSES:         *)
(*  MAXIMALE GENAUIGKEIT VON MPX UND MPY *)

VAR I,J,PREC       : INTEGER;
    MPZ            : MULPREC;
    DPX, DPL, DPU  : DOTPRECISION;
    RX             : REAL;
    INT1, INT2     : INTERVAL;
	

BEGIN
   
   PREC := MAX(MPX.PREC,MPY.PREC);
   
   DPL := #( FOR I := 1 TO MPX.PREC SUM( MPX.STAG[I] ) );
   RX := #*( FOR I := 1 TO MPY.PREC SUM( MPY.STAG[I] )
             + 0.5*MPY.INT.INF + 0.5*MPY.INT.SUP );
   
   MPZ.PREC := PREC;
   MPZ.STAG[1] := #*(DPL) / RX;
   
   FOR I := 2 TO PREC DO 
   BEGIN
       DPL := #(DPL - FOR J := 1 TO MPY.PREC SUM (MPY.STAG[J]*MPZ.STAG[I-1]) );
       MPZ.STAG[I] := #*( DPL ) / RX;
   END;
   
   DPL := #(DPL - FOR J := 1 TO MPY.PREC SUM (MPY.STAG[J]*MPZ.STAG[PREC]));
   DPU := DPL;
   
   ADD_INT_TIMES_MP(-MPY.INT, MPZ, DPL, DPU);
   
   DPX := #( FOR I := 1 TO MPY.PREC SUM ( MPY.STAG[I] ) );

   INT1 := INTVAL( #<(DPL + MPX.INT.INF), #>(DPU + MPX.INT.SUP) );
   INT2 := INTVAL( #<(DPX + MPY.INT.INF), #>(DPX + MPY.INT.SUP) );

   MPZ.INT := INT1 / INT2;

   RES := MPZ;

END;

(*--------------------------------------------------------------------------*)

GLOBAL FUNCTION ABS_SUP( MPX : MULPREC ) : REAL;

(*  SCHŽTZT ABS( MPX ) NACH OBEN AB  *)
(*  ERSTE MANTISSENSTELLEN DES STAG-ANTEILS WERDEN  *)
(*  MIT 1+1E-15 MULTIPLIZIERT, DIE WEITEREN         *)
(*  MANTISSENSTELLEN SIND DANN VERNACHLŽSSIGBAR     *)

BEGIN
   ABS_SUP := 1.000000000000001 * ABS( MPX.STAG[1] ) + SUP( ABS( MPX.INT ) );
END;

(*--------------------------------------------------------------------------*)

GLOBAL PROCEDURE MPSWRITE( NAME : STRING; VAR MPX : MULPREC );

VAR I : INTEGER;

BEGIN
   WRITELN( NAME );
   FOR I := 1 TO MPX.PREC DO
       WRITELN( MPX.STAG[I] );
   IWRITELN( MPX.INT );
END;

(*--------------------------------------------------------------------------*)
END.
(*==========================================================================*)
(*                                                                          *)
(*  PASCAL-XSC - MODUL  MPS_ARI                                             *)
(*                                                                          *)
(*  MULTIPLE PRECISION STAGGERED CORRECTION ARITHMETIC                      *)
(*                                                                          *)
(*==========================================================================*)