(*==========================================================================*)
(*                                                                          *)
(*  PASCAL-XSC - MODUL  ITAYLOR                                     000221  *)
(*                                                                          *)
(*  MODUL ZUR BERECHNUNG VON TAYLORKOEFFIZIENTEN MIT EINFACHER GENAUIGKEIT  *)
(*                                                                          *)
(*==========================================================================*)

MODULE itaylor;

USE i_ari;

GLOBAL TYPE ITAYLOR = GLOBAL DYNAMIC ARRAY [*] OF INTERVAL;

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

GLOBAL FUNCTION PLUS( IT1, IT2 : ITAYLOR; K : INTEGER ) : INTERVAL;

BEGIN
   PLUS := IT1[K] + IT2[K];
END;

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

GLOBAL FUNCTION MINUS( IT1 , IT2 : ITAYLOR; K : INTEGER ) : INTERVAL;

BEGIN
   MINUS := IT1[K] - IT2[K];
END;

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

GLOBAL FUNCTION TIMES( IT1 , IT2 : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR I : INTEGER;

BEGIN
   TIMES := ##(FOR I := 0 TO K SUM(IT1[I]*IT2[K-I]));
END;

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

GLOBAL FUNCTION DIVIDE( IT1 , IT2 , IT1_2 : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR I : INTEGER;

BEGIN
   DIVIDE := ( IT1[K] -
               ##( FOR I := 1 TO K SUM( IT1_2[K-I] * IT2[I] ) ) ) / IT2[0];
END;

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

GLOBAL FUNCTION SQR( IT : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR I : INTEGER;

BEGIN
   IF K = 0 THEN
	SQR := IT[0] * IT[0]
   ELSE
   BEGIN
	IF ODD( K ) THEN
	   SQR := 2 * ##( FOR I := 0 TO ( ( K - 1 ) DIV 2 )
                              SUM ( IT[I] * IT[K-I] ) )
	ELSE
	   SQR := 2 * ##( FOR I := 0 TO ( ( K - 2 ) DIV 2 )
                                SUM ( IT[I] * IT[K-I] ) )
                                    + IT[K DIV 2] * IT[K DIV 2];
   END;	
END;

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

GLOBAL FUNCTION SQRT( IT , ITSQ : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR I : INTEGER;

BEGIN
   IF K = 0 THEN
      SQRT := SQRT( IT[0] )
   ELSE
   BEGIN
	IF ODD( K ) THEN
	   SQRT := ( IT[K] - 2 * ##( FOR I := 1 TO ( ( K-1 ) DIV 2 )
	             SUM ( ITSQ[I] * ITSQ[K-I] ) ) ) / ( 2 * ITSQ[0] )
	ELSE
	   SQRT := ( IT[K] - 2 * ##( FOR I := 1 TO (( K-2 ) DIV 2 )
	             SUM ( ITSQ[I] * ITSQ[K-I] ) )
		           - ITSQ[K DIV 2] * ITSQ[K DIV 2] ) / ( 2 * ITSQ[0] );
   END;
END;

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

GLOBAL FUNCTION EXP( IT , ITEXP : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR I   : INTEGER;
    INT : INTERVAL;

BEGIN
   IF K = 0 THEN
      EXP := EXP( IT[0] )
   ELSE
   BEGIN
      INT := 0;
      FOR I := 0 TO K-1 DO
         INT := INT + ( K - I ) * ITEXP[I] * IT[K-I];
      EXP := INT / K;
   END;
END;

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

GLOBAL FUNCTION SIN( IT , ITSIN : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR INT1, INT2 : INTERVAL;
    I , J      : INTEGER;

BEGIN
   IF K < 2 THEN
   BEGIN
      IF K = 0 THEN
         SIN := SIN( IT[0] )
      ELSE
         SIN := COS( IT[0] ) * IT[1]
      END
   ELSE BEGIN
      INT1 := 0;
      FOR J := 1 TO K - 1 DO
      BEGIN
         INT2 := 0;
         FOR I := 0 TO J - 1 DO
            INT2 := INT2 + ( J - I ) * ITSIN[I] * IT[J-I];
         INT1 := INT1 + ( K - J ) * IT[K-J] * INT2 / J;
      END;
      SIN := COS( IT[0] ) * IT[K] - INT1 / K;
   END;
END;

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

GLOBAL FUNCTION COS( IT , ITCOS : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR INT1 , INT2 : INTERVAL;
    I , J       : INTEGER;

BEGIN
   IF K < 2 THEN
   BEGIN
      IF K = 0 THEN
         COS := COS( IT[0] )
      ELSE
         COS := - SIN( IT[0] ) * IT[1]
      END
   ELSE
   BEGIN
      INT1 := 0;
      FOR J := 1 TO K-1 DO
      BEGIN
         INT2 := 0;
         FOR I := 0 TO J - 1 DO
            INT2 := INT2 + ( J - I ) * ITCOS[I] * IT[J-I];
         INT1 := INT1 + ( K - J ) * IT[K-J] * INT2 / J;
      END;
      COS := - SIN( IT[0] ) * IT[K] - INT1 / K;
   END;
END;

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

GLOBAL PROCEDURE SIN_COS( VAR IT , ITSIN , ITCOS : ITAYLOR; K : INTEGER);

VAR I               : INTEGER;
    INTSIN , INTCOS : INTERVAL;

BEGIN
   IF K = 0 THEN
   BEGIN
      ITSIN[K] := SIN( IT[0] );
      ITCOS[K] := COS( IT[0] )
      END
   ELSE BEGIN
      INTSIN := 0;
      INTCOS := 0;
      FOR I := 0 TO K-1 DO
      BEGIN
         INTSIN := INTSIN + ( K - I ) * ITCOS[I] * IT[K-I];
         INTCOS := INTCOS + ( K - I ) * ITSIN[I] * IT[K-I];
      END;
      ITSIN[K] := INTSIN / K;
      ITCOS[K] := - INTCOS / K;
   END;
END;

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

GLOBAL FUNCTION SINH( IT , ITSINH : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR INT1 , INT2 : INTERVAL;
    I , J       : INTEGER;

BEGIN
   IF K < 2 THEN
   BEGIN
      IF K = 0 THEN
         SINH := SINH( IT[0] )
      ELSE
         SINH := COSH( IT[0] ) * IT[1]
      END
   ELSE
   BEGIN
      INT1 := 0;
      FOR J := 1 TO K - 1 DO
      BEGIN
         INT2 := 0;
         FOR I := 0 TO J - 1 DO
            INT2 := INT2 + ( J - I ) * ITSINH[I] * IT[J-I];
         INT1 := INT1 + ( K - J ) * IT[K-J] * INT2 / J;
      END;
      SINH := COSH( IT[0] ) * IT[K] + INT1 / K;
   END;
END;

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

GLOBAL FUNCTION COSH( IT , ITCOSH : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR INT1 , INT2 : INTERVAL;
    I , J       : INTEGER;

BEGIN
   IF K < 2 THEN
   BEGIN
      IF K = 0 THEN
         COSH := COSH( IT[0] )
      ELSE
         COSH := SINH( IT[0] ) * IT[1]
      END
   ELSE
   BEGIN
      INT1 := 0;
      FOR J := 1 TO K - 1 DO
      BEGIN
         INT2 := 0;
         FOR I := 0 TO J - 1 DO
            INT2 := INT2 + ( J - I ) * ITCOSH[I] * IT[J-I];
         INT1 := INT1 + ( K - J ) * IT[K-J] * INT2 / J;
      END;
      COSH := SINH( IT[0] ) * IT[K] + INT1 / K;
   END;
END;

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

GLOBAL PROCEDURE SINH_COSH( VAR IT , ITSINH , ITCOSH : ITAYLOR; K : INTEGER );

VAR I                 : INTEGER;
    INTSINH , INTCOSH : INTERVAL;

BEGIN
   IF K = 0 THEN
   BEGIN
      ITSINH[K] := SINH( IT[0] );
      ITCOSH[K] := COSH( IT[0] )
   END
   ELSE
   BEGIN
      INTSINH := 0;
      INTCOSH := 0;
      FOR I := 0 TO K - 1 DO
      BEGIN
         INTSINH := INTSINH + ( K - I ) * ITCOSH[I] * IT[K-I];
         INTCOSH := INTCOSH + ( K - I ) * ITSINH[I] * IT[K-I];
      END;
      ITSINH[K] := INTSINH / K;
      ITCOSH[K] := INTCOSH / K;
   END;
END;

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

GLOBAL FUNCTION POWER( IT, ITPOW: ITAYLOR; A: REAL; K: INTEGER ) : INTERVAL;

VAR I   : INTEGER;
    INT : INTERVAL;

BEGIN
   IF K = 0 THEN
      POWER := EXP( A * LN( IT[0] ) )
   ELSE
   BEGIN
      INT := 0;
      FOR I := 0 TO K - 1 DO
         INT := INT + ( A * ( K - I ) - I ) * ITPOW[I] * IT[K-I];
      POWER := INT / ( K * IT[0] );
   END;
END;

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

GLOBAL FUNCTION LN ( IT , ITLN : ITAYLOR; K : INTEGER ) : INTERVAL;

VAR I   : INTEGER;
    INT : INTERVAL;

BEGIN
   IF K = 0 THEN
      LN := LN( IT[0] )
   ELSE
   BEGIN
      INT := 0;
      FOR I := 1 TO K - 1 DO
         INT := INT + I * ITLN[I] * IT[K-I];
      LN := ( IT[K] - INT / K) / IT[0];
   END;
END;

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

END.

(*==========================================================================*)
(*                                                                          *)
(*  PASCAL-XSC - MODUL  ITAYLOR                                             *)
(*                                                                          *)
(*  MODUL ZUR BERECHNUNG VON TAYLORKOEFFIZIENTEN MIT EINFACHER GENAUIGKEIT  *)
(*                                                                          *)
(*==========================================================================*)