(*==========================================================================*)
(*                                                                          *)
(*  PASCAL-XSC - MODUL  MPS_TAYL                                    990530  *)
(*                                                                          *)
(*  MODUL ZUR BERECHNUNG VON TAYLORKOEFFIZIENTEN                            *)
(*  MIT MEHRFACHER GENAUIGKEIT                                              *)
(*                                                                          *)
(*==========================================================================*)

MODULE mps_tayl;

USE iostd , i_ari, mp_ari, mpi_ari, mps_aril;

(*  GLOBAL TYPE MPTAYLOR = MPVECTOR;       *)
(*                                         *)
(*  FšHRT ZU IRREPARABLEN COMPILERFEHLERN  *)
(*                                         *)

VAR TAYPREC : INTEGER;
GLOBAL PROCEDURE SET_TAYPREC( I : INTEGER );
BEGIN
   TAYPREC := I;
END;

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

OPERATOR := ( VAR MP : MULPREC; MPI : MPINTERVAL );

VAR I : INTEGER;

BEGIN
   MPVLCP( MPI );
   MP.PREC := TAYPREC;
   FOR I := 1 TO MP.PREC DO
   BEGIN
	MP.STAG[I] := MID( MPI );
	MPI := MPI - MP.STAG[I];
   END;
   MP.INT := MPI;
   MPFREE( MPI );
END;

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

OPERATOR := ( VAR MPI : MPINTERVAL; MP : MULPREC );

VAR I : INTEGER;

BEGIN
   MPI := MP.INT;
   FOR I := 1 TO MP.PREC DO
   MPI := MPI + MP.STAG[I];
END;

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

GLOBAL FUNCTION PLUS( MPT1, MPT2: MPVECTOR; K: INTEGER ): MULPREC;

BEGIN
   PLUS := MPT1[K] + MPT2[K];
END;

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

GLOBAL FUNCTION MINUS(MPT1, MPT2: MPVECTOR; K: INTEGER): MULPREC;

BEGIN
   MINUS := MPT1[K] - MPT2[K];
END;

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

GLOBAL FUNCTION TIMES(MPT1, MPT2: MPVECTOR; K: INTEGER): MULPREC;

VAR I  : INTEGER;
    MP : MULPREC;

BEGIN
   MP.PREC := MPT1[0].PREC;
   MP.STAG := 0;
   MP.INT  := 0;
   FOR I := 0 TO K DO
	MP := MP + MPT1[I] *< MPT2[K-I];
   TIMES := MP;
END;

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

GLOBAL FUNCTION DIVIDE(MPT1, MPT2, MPTDIV: MPVECTOR; K: INTEGER): MULPREC;

VAR I  : INTEGER;
    MP : MULPREC;

BEGIN
   MP := MPT1[K];
   FOR I := 1 TO K DO
	MP := MP - MPTDIV[K-I] *< MPT2[I];
   DIVIDE := MP / MPT2[0]; 
END;

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

GLOBAL FUNCTION SQR( MPT: MPVECTOR; K: INTEGER ): MULPREC;

VAR I  : INTEGER;
    MP : MULPREC;

BEGIN
   IF K = 0 THEN
	SQR := MPT[0] *< MPT[0]
   ELSE
   BEGIN
	IF ODD(K) THEN
	BEGIN
	   MP.PREC := MPT[0].PREC;
	   MP.STAG := 0;
	   MP.INT  := 0;
	   FOR I := 0 TO ((K-1) DIV 2) DO
		MP := MP + MPT[I] *< MPT[K-I];
	   SQR := 2 *< MP
	END
	ELSE
	BEGIN
	   MP.PREC := MPT[0].PREC;
	   MP.STAG := 0;
	   MP.INT  := 0;
	   FOR I := 0 TO ((K-2) DIV 2) DO
		MP := MP + MPT[I] *< MPT[K-I];
	   SQR := 2 *< MP + MPT[K DIV 2] *< MPT[K DIV 2];
	END;
   END;	
END;

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

GLOBAL FUNCTION SQRT( MPT, MPTSQ: MPVECTOR; K: INTEGER ): MULPREC;

VAR I   : INTEGER;
    MP  : MULPREC;
    MPI : MPINTERVAL;

BEGIN
   MPINIT( MPI );
   IF K = 0 THEN
   BEGIN
      MPI := MPT[0];
      SQRT := SQRT( MPI );
   END
   ELSE
   BEGIN
	IF ODD( K ) THEN
	BEGIN
	   MP.PREC := MPTSQ[0].PREC;
	   MP.STAG := 0;
	   MP.INT  := 0;
	   FOR I := 1 TO ((K-1) DIV 2) DO
		MP := MP + MPTSQ[I] *< MPTSQ[K-I];
	   SQRT := ( MPT[K] - 2 *< MP ) / ( 2 *< MPTSQ[0] )	
	END
	ELSE
	BEGIN
	   MP.PREC := MPT[0].PREC;
	   MP.STAG := 0;
	   MP.INT  := 0;
	   FOR I := 1 TO ((K-2) DIV 2) DO
		MP := MP + MPTSQ[I] *< MPTSQ[K-I];
	   SQRT := ( MPT[K] - 2 *< MP - MPTSQ[K DIV 2] *< MPTSQ[K DIV 2])
	           / ( 2 *< MPTSQ[0] );
	END;
   END;
   MPFREE(MPI);
END;

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

GLOBAL FUNCTION EXP( VAR MPTEXP : MPVECTOR; A , X_0 : REAL; K : INTEGER ) : MULPREC;

(*  K-TER TAYLORKOEFFIZIENT VON EXP( A*X ) AN DER STELLE X_0 WIRD BERECHNET *)
(*  MPTEXP ENTHŽLT ERSTEN BIS (K-1)-TEN TAYLORKOEFFIZIENTEN VOM EXP( A*X )  *)
(*  GENAUE AUSWERTUNG MIT LANGZAHL-STANDARDFUNKTIONEN DES MODULS MPI_ARI    *)
(*  DIE GEWšNSCHTE GENAUIGKEIT WIRD IM HAUPTPROGRAMM FESTGELEGT             *)

VAR PREC                 : INTEGER;
    MPIA , MPIX , MPIRES : MPINTERVAL;
    RES                  : MULPREC;

BEGIN
   PREC := MPTEXP[K].PREC;
   IF K = 0 THEN BEGIN
      MPINIT( MPIA );
      MPINIT( MPIX );
      MPINIT( MPIRES );
      MPIA := A;
      MPIX := X_0;
      MPIRES := EXP( MPIA * MPIX );
      RES := CONVERT( MPIRES , PREC );
      MPFREE( MPIA );
      MPFREE( MPIX );
      MPFREE( MPIRES );
      END
   ELSE BEGIN
      RES := ( MPTEXP[K-1] *> A );
      IF RES.PREC < PREC THEN
         RES.PREC := PREC;
      RES := RES / K;
      RES := ROUND( RES , PREC );
   END;
   EXP := RES;
END;

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

GLOBAL FUNCTION EXP( VAR MPTEXP , MPTF : MPVECTOR; K : INTEGER ) : MULPREC;

(*  K-TER TAYLORKOEFFIZIENT VON EXP( F(X) ) WIRD BERECHNET                  *)
(*  MPTF ENTHŽLT ERSTEN BIS K-TEN TAYLORKOEFFIZIENTEN VON F(X)              *)
(*  MPTEXP ENTHŽLT ERSTEN BIS (K-1)-TEN TAYLORKOEFFIZIENTEN VON EXP( F(X) ) *)

VAR I , PREC      : INTEGER;
    RES           : MULPREC;
    MPIX , MPIRES : MPINTERVAL;
   
BEGIN
   PREC := MPTEXP[K].PREC;
   IF K = 0 THEN BEGIN
      MPINIT( MPIX );
      MPINIT( MPIRES );
      MPIX := MPTF[0];
      MPIRES := EXP( MPIX );
      RES := CONVERT( MPIRES , PREC );
      MPFREE( MPIX );
      MPFREE( MPIRES );
      END
   ELSE BEGIN
      RES.PREC := PREC;
      RES.STAG := 0;
      RES.INT  := 0;
      FOR I := 0 TO K - 1 DO
          RES := ROUND( RES + ( K - I ) *> MPTEXP[I] *> MPTF[K-I] , PREC );
      IF RES.PREC < PREC THEN
         RES.PREC := PREC;
      RES := RES / K;
      RES := ROUND( RES , PREC );
   END;
   EXP := RES;
END;

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

GLOBAL FUNCTION SIN( VAR MPTSIN : MPVECTOR; A , X_0 : REAL; K : INTEGER ) : MULPREC;

(*  K-TER TAYLORKOEFFIZIENT VON SIN( A*X ) AN DER STELLE X_0 WIRD BERECHNET *)
(*  MPTSIN ENTHŽLT ERSTEN BIS (K-1)-TEN TAYLORKOEFFIZIENTEN VON SIN( A*X )  *)
(*  GENAUE AUSWERTUNG MIT LANGZAHL-STANDARDFUNKTIONEN DES MODULS MPI_ARI    *)
(*  DIE GEWšNSCHTE GENAUIGKEIT WIRD IM HAUPTPROGRAMM FESTGELEGT             *)

VAR PREC                 : INTEGER;
    MPIA , MPIX , MPIRES : MPINTERVAL;
    RES                  : MULPREC;

BEGIN
   PREC := MPTSIN[K].PREC;
   IF K < 2 THEN BEGIN
      MPINIT( MPIA );
      MPINIT( MPIX );
      MPINIT( MPIRES );
      MPIA := A;
      MPIX := X_0;
      IF K = 0 THEN
         MPIRES := SIN( MPIA * MPIX )
      ELSE
         MPIRES := A * COS( MPIA * MPIX );
      RES := CONVERT( MPIRES , PREC );
      MPFREE( MPIA );
      MPFREE( MPIX );
      MPFREE( MPIRES );
      END
   ELSE BEGIN
      RES := - ( MPTSIN[K-2] *> A ) *> A;
      IF RES.PREC < PREC THEN
         RES.PREC := PREC;
      RES := RES / ( K * ( K - 1 ) );
      RES := ROUND( RES , PREC );
   END;
   SIN := RES;
END;

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

GLOBAL FUNCTION SIN( VAR MPTSIN , MPTF : MPVECTOR; K : INTEGER) : MULPREC;

(*  K-TER TAYLORKOEFFIZIENT VON SIN( F(X) ) WIRD BERECHNET                  *)
(*  MPTF ENTHŽLT ERSTEN BIS K-TEN TAYLORKOEFFIZIENTEN VON F(X)              *)
(*  MPTSIN ENTHŽLT ERSTEN BIS (K-1)-TEN TAYLORKOEFFIZIENTEN VON SIN( F(X) ) *)

VAR I , J     : INTEGER;
    MP1 , MP2 : MULPREC;
    MPI       : MPINTERVAL;
    
BEGIN
   MPINIT( MPI );
   IF K < 2 THEN
   BEGIN
      IF K = 0 THEN
      BEGIN
         MPI := MPTF[0];
         SIN := SIN( MPI )
      END
      ELSE
      BEGIN
         MPI := MPTF[0];
         MP1 := COS( MPI );
         SIN := MP1 *< MPTF[1]
	END
   END
   ELSE
   BEGIN
	MP1.PREC := MPTSIN[0].PREC;
	MP1.STAG := 0;
	MP1.INT  := 0;
      FOR J := 1 TO K-1 DO
      BEGIN
	   MP2.PREC := MPTSIN[0].PREC;
	   MP2.STAG := 0;
	   MP2.INT  := 0;   
         FOR I := 0 TO J-1 DO
            MP2 := MP2 + ( J - I ) *< MPTSIN[I] *< MPTF[J-I];
         MP1 := MP1 + ( K - J ) *< MPTF[K-J] *< MP2 / J;
      END;
	MP2 := COS(MPI);
      SIN := MP2 *< MPTF[K] - MP1 / K;
   END;
   MPFREE( MPI );
END;

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

GLOBAL FUNCTION COS( VAR MPTCOS : MPVECTOR; A , X_0 : REAL; K : INTEGER ) : MULPREC;

(*  K-TER TAYLORKOEFFIZIENT VON COS( A*X ) AN DER STELLE X_0 WIRD BERECHNET *)
(*  MPTCOS ENTHŽLT ERSTEN BIS (K-1)-TEN TAYLORKOEFFIZIENTEN VON COS( A*X )  *)
(*  GENAUE AUSWERTUNG MIT LANGZAHL-STANDARDFUNKTIONEN DES MODULS MPI_ARI    *)
(*  DIE GEWšNSCHTE GENAUIGKEIT WIRD IM HAUPTPROGRAMM FESTGELEGT             *)

VAR PREC                 : INTEGER;
    MPIA , MPIX , MPIRES : MPINTERVAL;
    RES                  : MULPREC;

BEGIN
   PREC := MPTCOS[K].PREC;
   IF K < 2 THEN BEGIN
      MPINIT( MPIA );
      MPINIT( MPIX );
      MPINIT( MPIRES );
      MPIA := A;
      MPIX := X_0;
      IF K = 0 THEN
         MPIRES := COS( MPIA * MPIX )
      ELSE
         MPIRES := - A * SIN( MPIA * MPIX );
      RES := CONVERT( MPIRES , PREC );
      MPFREE( MPIA );
      MPFREE( MPIX );
      MPFREE( MPIRES );
      END
   ELSE BEGIN
      RES := - ( MPTCOS[K-2] *> A ) *> A;
      IF RES.PREC < PREC THEN
         RES.PREC := PREC;
      RES := RES / ( K * ( K - 1 ) );
      RES := ROUND( RES , PREC );
   END;
   COS := RES;
END;

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

GLOBAL FUNCTION COS( VAR MPTCOS , MPTF : MPVECTOR; K : INTEGER) : MULPREC;

VAR MP1, MP2   : MULPREC;
    I, J       : INTEGER;
    MPI        : MPINTERVAL;

BEGIN
   MPINIT(MPI);
   IF K < 2 THEN
   BEGIN
      IF K = 0 THEN BEGIN
         MPI := MPTF[0];
         COS := COS(MPI)
         END
      ELSE BEGIN
         MPI := MPTF[0];
         MP1 := -SIN( MPI );
         COS := MP1 *< MPTF[1]
         END
      END
   ELSE BEGIN
      MP1.PREC := MPTCOS[0].PREC;
      MP1.STAG := 0;
      MP1.INT  := 0;
      FOR J := 1 TO K-1 DO BEGIN
          MP2.PREC := MPTCOS[0].PREC;
          MP2.STAG := 0;
          MP2.INT  := 0;
          FOR I := 0 TO J-1 DO
              MP2 := MP2 + (J-I) *< MPTCOS[I] *< MPTF[J-I];
          MP1 := MP1 + (K-J) *< MPTF[K-J] *< MP2/J;
      END;
      MP2 := -SIN(MPI);
      COS := MP2 *< MPTF[K] - MP1 / K;
   END;
   MPFREE(MPI);
END;

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

GLOBAL PROCEDURE SIN_COS( VAR MPTF, MPTSIN, MPTCOS : MPVECTOR; K : INTEGER);

(*  0-TE BIS K-TE TAYLORKOEFFIZIENTEN VON SIN UND COS WERDEN BERECHNET *)

VAR I              : INTEGER;
    MPSIN, MPCOS   : MULPREC;
    MPI            : MPINTERVAL;

BEGIN
   MPINIT( MPI );
   IF K = 0 THEN BEGIN
      MPI := MPTF[0];
      MPTSIN[K] := SIN(MPI);
      MPTCOS[K] := COS(MPI)
      END
   ELSE BEGIN
      MPSIN.PREC := MPTSIN[0].PREC;
      MPSIN.STAG := 0;
      MPSIN.INT  := 0;
      MPCOS.PREC := MPTCOS[0].PREC;
      MPCOS.STAG := 0;
      MPCOS.INT  := 0;
      FOR I := 0 TO K-1 DO BEGIN
          MPSIN := MPSIN + (K-I) *< MPTCOS[I] *< MPTF[K-I];
          MPCOS := MPCOS + (K-I) *< MPTSIN[I] *< MPTF[K-I];
      END;
      MPTSIN[K] := MPSIN / K;
      MPTCOS[K] := - MPCOS / K;
   END;
   MPFREE( MPI );
END;

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

GLOBAL FUNCTION SINH( VAR MPT , MPTSINH : MPVECTOR; K : INTEGER) : MULPREC;

VAR MP1, MP2   : MULPREC;
    I, J       : INTEGER;
    MPI        : MPINTERVAL;

BEGIN
   MPINIT( MPI );
   IF K < 2 THEN
      IF K = 0 THEN BEGIN
         MPI := MPT[0];
         SINH := SINH(MPI)
         END
      ELSE BEGIN
         MPI := MPT[0];
         MP1 := COSH(MPI);
         SINH := MP1 *< MPT[1]
      END
   ELSE BEGIN
      MP1.PREC := MPTSINH[0].PREC;
      MP1.STAG := 0;
      MP1.INT  := 0;
      FOR J := 1 TO K-1 DO BEGIN
          MP2.PREC := MPTSINH[0].PREC;
          MP2.STAG := 0;
          MP2.INT  := 0;
          FOR I := 0 TO J-1 DO
              MP2 := MP2 + (J-I) *< MPTSINH[I] *< MPT[J-I];
          MP1 := MP1 + (K-J) *< MPT[K-J] *< MP2/J;
      END;
      MP2 := COSH( MPI );
      SINH := MP2 *< MPT[K] + MP1 / K;
   END;
   MPFREE( MPI );
END;

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

GLOBAL FUNCTION COSH( VAR MPT , MPTCOSH : MPVECTOR; K : INTEGER) : MULPREC;

VAR MP1, MP2   : MULPREC;
    I, J       : INTEGER;
    MPI        : MPINTERVAL;

BEGIN
   MPINIT(MPI);
   IF K < 2 THEN
      IF K = 0 THEN BEGIN
         MPI := MPT[0];
         COSH := COSH(MPI)
         END
      ELSE BEGIN
         MPI := MPT[0];
         MP1 := SINH(MPI);
         COSH := MP1 *< MPT[1]
      END
   ELSE BEGIN
      MP1.PREC := MPTCOSH[0].PREC;
      MP1.STAG := 0;
      MP1.INT  := 0;
      FOR J := 1 TO K-1 DO BEGIN
          MP2.PREC := MPTCOSH[0].PREC;
          MP2.STAG := 0;
          MP2.INT  := 0;
          FOR I := 0 TO J-1 DO
              MP2 := MP2 + (J-I) *< MPTCOSH[I] *< MPT[J-I];
          MP1 := MP1 + (K-J) *< MPT[K-J] *< MP2/J;
      END;
      MP2 := SINH(MPI);
      COSH := MP2 *< MPT[K] + MP1/K;
   END;
   MPFREE(MPI);
END;

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

GLOBAL PROCEDURE SINH_COSH( VAR MPT , MPTSINH , MPTCOSH : MPVECTOR; K : INTEGER);

VAR I                : INTEGER;
    MPSINH, MPCOSH   : MULPREC;
    MPI              : MPINTERVAL;

BEGIN
   MPINIT(MPI);
   IF K = 0 THEN BEGIN
      MPI := MPT[0];
      MPTSINH[K] := SINH(MPI);
      MPTCOSH[K] := COSH(MPI)
      END
   ELSE BEGIN
      MPSINH.PREC := MPTSINH[0].PREC;
      MPSINH.STAG := 0;
      MPSINH.INT  := 0;
      MPCOSH.PREC := MPTCOSH[0].PREC;
      MPCOSH.STAG := 0;
      MPCOSH.INT  := 0;
      FOR I := 0 TO K-1 DO BEGIN
          MPSINH := MPSINH + (K-I) *< MPTCOSH[I] *< MPT[K-I];
          MPCOSH := MPCOSH + (K-I) *< MPTSINH[I] *< MPT[K-I];
      END;
      MPTSINH[K] := MPSINH / K;
      MPTCOSH[K] := MPCOSH / K;
   END;
   MPFREE( MPI );
END;

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

GLOBAL FUNCTION POWER( VAR MPT , MPTPOW : MPVECTOR; A : REAL; K : INTEGER)
                     : MULPREC;

VAR I   : INTEGER;
    MP  : MULPREC;
    MPI : MPINTERVAL;

BEGIN
   MPINIT( MPI );
   IF K = 0 THEN BEGIN
      MPI := MPT[0];
      POWER := EXP(A*LN(MPI))
      END
   ELSE BEGIN
      MP.PREC := MPTPOW[0].PREC;
      MP.STAG := 0;
      MP.INT  := 0;
      FOR I := 0 TO K-1 DO
          MP := MP + (A*(K-I)-I) *< MPTPOW[I] *< MPT[K-I];
      POWER := MP/(K *< MPT[0]);
   END;
   MPFREE( MPI );
END;

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

GLOBAL FUNCTION LN( VAR MPT , MPTLN : MPVECTOR; K : INTEGER) : MULPREC;

VAR I   : INTEGER;
    MP  : MULPREC;
    MPI : MPINTERVAL;

BEGIN
   MPINIT( MPI );
   IF K = 0 THEN BEGIN
      MPI := MPT[0];
      LN := LN( MPI )
      END
   ELSE BEGIN
      MP.PREC := MPTLN[0].PREC;
      MP.STAG := 0;
      MP.INT  := 0;
      FOR I := 1 TO K-1 DO
          MP := MP + I *< MPTLN[I] *< MPT[K-I];
      LN := ( MPT[K] - MP / K ) / MPT[0];
   END;
   MPFREE( MPI );
END;

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

END.

(*==========================================================================*)
(*                                                                          *)
(*  PASCAL-XSC - MODUL  MPS_TAYL                                            *)
(*                                                                          *)
(*  MODUL ZUR BERECHNUNG VON TAYLORKOEFFIZIENTEN                            *)
(*  MIT MEHRFACHER GENAUIGKEIT                                              *)
(*                                                                          *)
(*==========================================================================*)