MODULE ci_functions;
{ ************************************************************************** }
{ ************************************************************************** }
{ ***                                                                    *** }
{ ***        Weitere Standardfunktionen fuer komplexe Intervalle         *** }
{ ***                                                                    *** }
{ ************************************************************************** }
{ ************************************************************************** }
{ ************************************************************************** }
{ ***  Literatur :                                                       *** }
{ ***                                                                    *** }
{ ***   Braune, K.D. : 'Hochgenaue Standardfunktionen fuer reelle und    *** }
{ ***     und komplexe Punkte und Intervalle in beliebigen Gleitpunkt-   *** }
{ ***     rastern', Dissertation, Karlsruhe, 1987.                       *** }
{ ***   Buehler, G. : 'Standardfunktionen fuer komplexe Intervalle im    *** }
{ ***     64 Bit IEEE Datenformat', Diplomarbeit, Karlsruhe, 1993.       *** }
{ ***   Kraemer, W. : 'Inverse Standardfunktionen fuer reelle und        *** }
{ ***     komplexe Intervallargumente mit a priori Fehlerabschaetzungen  *** }
{ ***     fuer beliebige Datenformate', Dissertation, Karlsruhe, 1987.   *** }
{ ***                                                                    *** }
{ ************************************************************************** }
{ ***  Autor: Andreas Westphal                                           *** }
{ ***  Betreuer: Walter Kraemer                                          *** } 
{ ***     Institut fuer Wissenschaftliches Rechnen  und                  *** }  
{ ***     Mathematische Modellbildung, Universitaet Karlsruhe (TH)       *** }
{ ***                       Stand: Maerz 1999                            *** }
{ ************************************************************************** }

USE i_ari,ci_ari;

PROCEDURE Failed;
VAR Fehler : integer;
{ Was soll bei einem aufgetretenen Fehler (meist schlecht }
{ gewaehltes Eingabeintervall) geschehen?                 }
{ Hier: Abbruch mittels Division durch Null.              }
BEGIN
  Fehler:= 1 div 0;
END;


PROCEDURE z( vari : cinterval; VAR re_l,re_u,im_l,im_u : real );
BEGIN
  re_l:= inf(re(vari));
  re_u:= sup(re(vari));
  im_l:= inf(im(vari));
  im_u:= sup(im(vari));
END;

PROCEDURE z( x,y : interval; VAR x_l,x_u,y_l,y_u : real );
BEGIN
   x_l:= inf(x);
   x_u:= sup(x);
   y_l:= inf(y);
   y_u:= sup(y);
END;

FUNCTION max( x,y : real ) : real;
BEGIN
  IF x<y THEN max:=y
  ELSE max:= x;
END;   


{ ************************************************************************** }
{ ************************************************************************** }
{ ***                      Eindeutige Funktionen                         *** }
{ ************************************************************************** }
{ ************************************************************************** }


{ ************************************************************************** }
{ *** Der Operator 'Hoch'  pow  steht an spaeterer Stelle, da der natuer-*** }
{ *** liche Logarithmus verwendet wird, der weiter hinten bei den nicht  *** }
{ *** eindeutigen Funktionen steht.                                      *** }
{ ************************************************************************** }
 

{ ************************************************************************** }
{ *** Die hyperbolischen Funktionen exp, sin, cos, sinh, cosh besitzen   *** }
{ *** separable Formeln, d.h. ihr Real-bzw. Imaginaerteile lassen sich   *** }
{ *** als Produkt (reeller) Funktionen darstellen. Damit ist die Imple-  *** }
{ *** mentierung ohne Ueberschaetzung durch intervallmaessige Auswertung *** }
{ *** dieser Formeln moeglich.                                           *** }
{ ************************************************************************** }
{ ***   Formeln mit Re(z)=x und Im(z)=y :                                *** }
{ ***                                                                    *** }
{ ***        exp   :    Re(exp(z)) = exp(x) * cos(y)                     *** }
{ ***                   Im(exp(z)) = exp(x) * sin(y)                     *** }
{ ***                                                                    *** }
{ ***        sin   :    Re(sin(z)) = sin(x) * cosh(y)                    *** }
{ ***                   Im(sin(x)) = cos(x) * sinh(y)                    *** }
{ ***                                                                    *** }
{ ***        cos   :    Re(cos(z)) = cos(x) * cosh(y)                    *** }
{ ***                   Im(sin(x)) = -sin(x) * sinh(y)                   *** }
{ ***                                                                    *** }
{ ***       sinh   :   Re(sinh(z)) = sinh(x) * cos(y)                    *** }
{ ***                  Im(sinh(z)) = cosh(x) * sin(y)                    *** }
{ ***                                                                    *** }
{ ***       cosh   :   Re(cosh(z)) = cosh(x) * cos(y)                    *** }
{ ***                  Im(cosh(z)) = sinh(x) * sin(y)                    *** }
{ ***                                                                    *** }
{ ************************************************************************** }
 
GLOBAL FUNCTION exp( vari : cinterval ) : cinterval; 
VAR rewert,imwert: interval;
BEGIN
  rewert:= exp(re(vari))*cos(im(vari));
  imwert:= exp(re(vari))*sin(im(vari));
  exp:= compl(rewert,imwert);
END;

GLOBAL FUNCTION sin( vari : cinterval ) : cinterval;
VAR  rewert,imwert: interval;
BEGIN
  rewert:= sin(re(vari))*cosh(im(vari));
  imwert:= cos(re(vari))*sinh(im(vari));
  sin:= compl(rewert,imwert);
END;

GLOBAL FUNCTION cos( vari : cinterval ) : cinterval;
VAR rewert,imwert: interval;
BEGIN
  rewert:= cos(re(vari))*cosh(im(vari));
  imwert:= -sin(re(vari))*sinh(im(vari));
  cos:= compl(rewert,imwert); 
END;

GLOBAL FUNCTION sinh( vari : cinterval ) : cinterval;
VAR rewert,imwert: interval;
BEGIN
  rewert:= sinh(re(vari))*cos(im(vari));
  imwert:= cosh(re(vari))*sin(im(vari));
  sinh:= compl(rewert,imwert);
END;

GLOBAL FUNCTION cosh( vari : cinterval ) : cinterval;
VAR rewert,imwert: interval;
BEGIN
  rewert:= cosh(re(vari))*cos(im(vari));
  imwert:= sinh(re(vari))*sin(im(vari));
  cosh:= compl(rewert,imwert);
END;

{ ************************************************************************** }
{ *** Der Tangens ist NICHT separabel, so dass eine naiv ausgefuehrte    *** }
{ *** Auswertung im allgemeinen zu einer Ueberschaetzung fuehrt. Deshalb *** }
{ *** ist die folgende Implementierung aufwendiger und verwendet Fall-   *** }
{ *** unterscheidungen, die die Lage des Arguments betreffen. Um deren   *** }
{ *** Anzahl zu halbieren wird die Symmetrieeigenschaft                  *** }
{ ***                                                                    *** }
{ ***             tan(z) = transp( tan( transp(z) ) )                    *** }
{ ***                                                                    *** }
{ *** benutzt.                                                           *** }
{ ************************************************************************** }

PROCEDURE htan( x,y : interval;  VAR Re_tan,Im_tan : interval ; BOTH : boolean);
{ Ist BOTH=TRUE, so werden Einschliessungen von  Re(tan(x+i*y)) und  }
{ Im(tan(x+i*y)) berechnet und Re_tan bzw. Im_tan zugewiesen. Anson- }
{ sten wird nur eine Einschliessung von Im(tan(x+i*y)) berechnet und }
{ und Im_tan zugewiesen.                                             }
{ Es wird immer vorausgesetzt, das y positiv ist ( Nur obere Halbebene ). }
VAR argx,argy,hint  : interval;
    PiHalbe         : interval;
    xl,xu,yl,yu     : real;
    hint2           : interval; 
    Re_first,Im_first : boolean;
{}
PROCEDURE ReAdd( hx,hy : real );
BEGIN
  IF BOTH THEN
   IF Re_first THEN { Naive Formel fuer Realteil des Tangens }
   BEGIN
    Re_tan:= sin(intval(hx))*cos(intval(hx))/
              (sqr(cos(intval(hx)))+sqr(sinh(intval(hy))));
    Re_first:= FALSE;
   END ELSE { Interhuelle mit schon berechneten Werten }
    Re_tan:= Re_tan +* sin(intval(hx))*cos(intval(hx))/
                        (sqr(cos(intval(hx)))+sqr(sinh(intval(hy))));   
END;
{}
PROCEDURE ReAdd( hx : interval; hy : real );
BEGIN
  IF BOTH THEN
   IF Re_first THEN { Naive Formel fuer Realteil des Tangens }
   BEGIN
    Re_tan:= sin(hx)*cos(hx)/
              (sqr(cos(hx))+sqr(sinh(intval(hy))));
    Re_first:= FALSE;
   END ELSE { Interhuelle mit schon berechneten Werten }
    Re_tan:= Re_tan +* sin(hx)*cos(hx)/
                        (sqr(cos(hx))+sqr(sinh(intval(hy))));
END;
{}
PROCEDURE ImAdd( hx,hy : real );
BEGIN
   IF Im_first THEN
    BEGIN { Naive Formel fuer Imaginaerteil des Tangens }
      Im_tan:= sinh(intval(hy))*cosh(intval(hy))/
                 (sqr(cos(intval(hx)))+sqr(sinh(intval(hy))));
      Im_first:= FALSE;
    END
   ELSE { Intervallhuelle }
    Im_tan:= Im_tan +* sinh(intval(hy))*cosh(intval(hy))/
                         (sqr(cos(intval(hx)))+sqr(sinh(intval(hy))));
END;
{}
PROCEDURE ImAdd( hx : real; hy : interval );
BEGIN
   IF Im_first THEN
    BEGIN { Naive Formel fuer Imaginaerteil des Tangens }
      Im_tan:= sinh(hy)*cosh(hy)/
                 (sqr(cos(intval(hx)))+sqr(sinh(hy)));
      Im_first:= FALSE;
    END
   ELSE { Intervallhuelle }
    Im_tan:= Im_tan +* sinh(hy)*cosh(hy)/
                         (sqr(cos(intval(hx)))+sqr(sinh(hy)));
END;
{}
BEGIN
  Re_first:= TRUE;
  Im_first:= TRUE;
  PiHalbe:= arccos(intval(0.0));
  { Realteil des Arguments von links nach rechts }
  IF ( inf(x)< -inf(PiHalbe)/2 ) THEN { Realteil des Argument }
  BEGIN                               { schneidet [-pi/2,-pi/4] }
    IF (-inf(PiHalbe)/2<sup(x)) THEN hint:= intval(inf(x),-inf(PiHalbe)/2)
    ELSE hint:= x;              { hint = Schnitt von x mit [-pi/2,-pi/4] }
    ReAdd(inf(x),sup(y));     { evtl. Maximum }
    ReAdd(sup(hint),sup(y));  { evtl. Maximum }
    ImAdd(sup(hint),sup(y));  { evtl. Minimum }  
    ImAdd(sup(hint),inf(y));  { evtl. Minimum }
    IF BOTH THEN
    { Minimum des Realteils in hint }
     BEGIN
       IF (inf(y)>0) THEN 
        BEGIN
          argx:= arctan(coth(intval(inf(y)))); 
          IF ( argx >< hint ) THEN 
           BEGIN
             IF inf(x)< inf(argx) THEN ReAdd(inf(x),inf(y)) { Minimum }
             ELSE ReAdd(sup(hint),inf(y));                  { Minimum }
           END
          ELSE ReAdd(argx,inf(y));
        END
       ELSE ReAdd(inf(x),inf(y)); { Minimum }
     END;
    { Maximum des Imaginaerteils in hint }
    IF ( inf(x)< -sup(PiHalbe)/2 ) THEN
     BEGIN
       argy:= arcoth(-tan(intval(inf(x)))); { Maximum bei festem x }
                                            { und beliebigem y }
       IF ( argy >< y ) THEN { Maximum in einer der linken Eckpunkte }
        BEGIN
          IF inf(y)< sup(argy) THEN ImAdd(inf(x),inf(y))
          ELSE ImAdd(inf(x),sup(y));
        END
       ELSE ImAdd(inf(x),argy)
     END
    ELSE ImAdd(inf(x),y);
  END;
  hint:= intval(-inf(PiHalbe)/2 , 0);
  IF NOT ( hint >< x ) THEN { Realteil des Arguments schneidet [-pi/4,0] }
  BEGIN
    hint:= hint ** x;
    ReAdd(sup(hint),sup(y)); { Maximum } 
    ReAdd(inf(hint),inf(y)); { Minimum }
    ImAdd(inf(hint),sup(y)); { Maximum }
    ImAdd(sup(hint),inf(y)); { Minimum }
  END;
  hint:= intval( 0 ,inf(PiHalbe)/2 );
  IF NOT ( hint >< x ) THEN { Realteil des Argument schneidet [0,pi/4] }
  BEGIN
    hint:= hint ** x;
    ReAdd(sup(hint),inf(y)); { Maximum }
    ReAdd(inf(hint),sup(y)); { Minimum }
    ImAdd(sup(hint),sup(y)); { Maximum }
    ImAdd(inf(hint),inf(y)); { Minimum }
  END;
  IF ( sup(x)> inf(PiHalbe)/2 ) THEN { Realteil des Argument }
  BEGIN                              { schneidet [pi/4,pi/2] }
    IF (inf(x)<inf(PiHalbe)/2) THEN hint:= intval(inf(PiHalbe)/2,sup(x))
    ELSE hint:= x;
    IF BOTH THEN
     BEGIN
     { Maximum des Realteils in hint }
       IF ( inf(y)=0 ) THEN
        ReAdd(sup(x),inf(y))
       ELSE
        BEGIN
          argx:= arctan(coth(intval(inf(y)))); { Maximum fuer festes argy, }
                                               { beliebiges argx }
          IF ( argx >< x ) THEN { Maximum in einer der unteren Ecken }
           BEGIN
             IF sup(argx)< inf(hint) THEN ReAdd(inf(hint),inf(y))
             ELSE ReAdd(sup(x),inf(y));
           END;
        END;
     END;
    { Maximum des Imaginaerteils in hint }
    IF ( sup(x) < inf(PiHalbe)/2 ) THEN
     BEGIN
       argy:= arcoth(tan(intval(sup(x)))); { Maximum bei festem x }
                                           { und beliebigem y }
       IF ( argy >< y ) THEN { Maximum in einer der rechten Eckpunkte }
        IF sup(y)< inf(argy) THEN ImAdd(sup(x),sup(y))
        ELSE ImAdd(sup(x),inf(y));
     END
    ELSE ImAdd(sup(x),y);
    { In Frage kommende Minima : linke Eckpunkte }
    ImAdd(inf(hint),sup(y));
    ImAdd(inf(hint),inf(y));
    { In Frage kommende Minima: oberen Eckpunkte }
    ReAdd(inf(hint),sup(y));
    ReAdd(sup(x),sup(y));  
  END;
END;


PROCEDURE h_tan( x,y : interval; VAR re_tan,im_tan : interval;
                                          VAR Error : integer); 
VAR Pi,hx,hx2   : interval;
    Imtan,Retan : interval;
    BOTH        : boolean;
    Intervallaufteilung : boolean; { Liegt  Pi/2 mod Pi  in x ? }
BEGIN
  Error:= 0;
  IF ( x=intval(0) ) THEN
   BEGIN
     re_tan:= x;               { Re(tan) gleich NULL }
     im_tan:= tanh(y);
   END
  ELSE
   BEGIN
     IF NOT (intval(0)><y) AND NOT (intval(0)><cos(x)) THEN
      Error:=1
     ELSE
      BEGIN
        IF ( y=intval(0) ) THEN
         BEGIN
           re_tan:= tan(x);
           im_tan:= y;         { Im(tan) gleich NULL }
         END
        ELSE
         BEGIN 
           Pi:= arccos(intval(-1.0));
           { Wegen  vari equivalent zu (vari mod pi), suche }
           { entsprechendes Argument in [ -pi/2 , pi/2 ] .  }
           IF ( sign(inf(x)) = 1 ) THEN
            hx:= x - trunc( inf(x/Pi) + 0.5 ) * Pi
           ELSE
            hx:= x - trunc( inf(x/Pi) - 0.5 ) * Pi;
           IF ( inf(Pi)< 2*sup(hx) ) THEN
            {************}
            BEGIN
              Intervallaufteilung:= TRUE;
              hx2:= intval(-sup(Pi)/2,sup(hx)-inf(Pi)/2);
              hx:= intval(inf(hx),sup(Pi)/2);
            END; 
            {************}
           BOTH:= TRUE;
           { Nutze aus dass, tan(transp(z))=transp(tan(z)) fuer cinterval z. }
           { Damit sind nur Argumente in der oberen Halbebene noetig!!       }
           IF ( sup(y)<=0) THEN
            BEGIN
              IF Intervallaufteilung THEN
               BEGIN
                 htan(hx,-y,Re_tan,Im_tan,BOTH);
                 htan(hx2,-y,retan,imtan,BOTH);
                 Im_tan:= -( Im_tan +* imtan);
                 Re_tan:= Re_tan +* retan;
               END
              ELSE
               BEGIN
                 htan(hx,-y,Re_tan,Im_tan,BOTH);
                 Im_tan:= -Im_tan;
               END;
            END
           ELSE
            IF ( inf(y)>=0) THEN
             IF Intervallaufteilung THEN
              BEGIN
                htan(hx,y,Re_tan,Im_tan,BOTH);
                htan(hx2,y,retan,imtan,BOTH);
                Im_tan:= Im_tan +* imtan;
                Re_tan:= Re_tan +* retan;
              END
             ELSE  
              htan(hx,y,Re_tan,Im_tan,BOTH)
            ELSE
             IF -inf(y)<sup(y) THEN
              IF Intervallaufteilung THEN
               BEGIN
                 htan(hx,intval(0,sup(y)),Re_tan,Im_tan,BOTH);
                 htan(hx2,intval(0,sup(y)),retan,im_tan,BOTH);
                 Im_tan:= Im_tan +* imtan;
                 Re_tan:= Re_tan +* retan;
                 htan(hx,intval(0,-inf(y)),Re_tan,Imtan,NOT BOTH);
                 Im_tan:= Im_tan +* -Imtan; 
                 htan(hx2,intval(0,-inf(y)),Re_tan,Imtan,NOT BOTH);
                 Im_tan:= Im_tan +* -Imtan
               END
              ELSE
               BEGIN
                 htan(hx,intval(0,sup(y)),Re_tan,Im_tan,BOTH);
                 htan(hx,intval(0,-inf(y)),Re_tan,Imtan,NOT BOTH);
                 Im_tan:= Im_tan +* -Imtan;
               END
              ELSE
               IF Intervallaufteilung THEN
                BEGIN
                  htan(hx,intval(0,-inf(y)),Re_tan,Imtan,BOTH);
                  htan(hx,intval(0,sup(y)),Re_tan,Im_tan,NOT BOTH);
                  Im_tan:= Im_tan +* -Imtan;
                  htan(hx2,intval(0,-inf(y)),Retan,Imtan,BOTH);
                  Im_tan:= Im_tan +* imtan;
                  Re_tan:= Re_tan +* retan;
                  htan(hx2,intval(0,sup(y)),Re_tan,Im_tan,NOT BOTH);
                  Im_tan:= Im_tan +* -Imtan;
                END
               ELSE
                BEGIN
                  htan(hx,intval(0,-inf(y)),Re_tan,Imtan,BOTH);
                  htan(hx,intval(0,sup(y)),Re_tan,Im_tan,NOT BOTH);
                  Im_tan:= Im_tan +* -Imtan;
                END;
         END;
      END;
   END;
END;



{ ************************************************************************** }
{ ***  Implementierung von cot, tanh, coth basierend auf tan.            *** }
{ ************************************************************************** }
{ ************************************************************************** }
{ *** Formeln:    cot(z)  = tan( pi/2 - z )                              *** }
{ ***             tanh(z) = transp( i * tan( transp( i * z ) ) )         *** }
{ ***             coth(z) = i * cot( i * z )                             *** }
{ ***                     = i * tan( pi/2 - i * z )                      *** }
{ ************************************************************************** }


PROCEDURE h_cot( x,y : interval; VAR re_cot,im_cot : interval;
                                          VAR Error : integer);
VAR Pi,hx,hx2   : interval;
    Imcot,Recot : interval;
    BOTH        : boolean;
    Intervallaufteilung : boolean; { Liegt  0 mod Pi  in x ? }
            { Dann wird das Eingabeintervall in aufgeteilt }                    
BEGIN
  Error:= 0;
  IF NOT (intval(0)><y) AND NOT (intval(0)><sin(x)) THEN
   Error:=1
  ELSE
   IF ( y=intval(0) ) THEN
    BEGIN
      re_cot:= cot(x);
      im_cot:= y;               { Im(cot) gleich NULL } 
    END
   ELSE
    IF ( x=intval(0) ) THEN
     BEGIN
       re_cot:= x;              { Re(cot) gleich NULL }
       im_cot:= coth(y);
     END
    ELSE 
     BEGIN
       Pi:= arccos(intval(-1.0));
       { Wegen  vari equivalent zu (vari mod pi), suche entsprechendes }
       { Argument in [ 0 , Pi ] . Wegen der Verwendung der Beziehung   }
       {                    cot(z)=tan(Pi/2 - z)                       }
       { berechnet man gleich                                          }
       {       hz = Pi/2 - z    mod Pi    in  [ -Pi/2 , Pi/2 ]  ,      }
       { um die Procedure htan verwenden zu koennen.                   }
       IF ( sup(x)-inf(x)>inf(Pi) ) THEN
        hx:= intval(-sup(Pi),sup(Pi))/2
       ELSE
        BEGIN
          Pi:= arccos(intval(-1.0));
          { Wegen  vari equivalent zu (vari mod pi), suche }
          { entsprechendes Argument in [ -pi/2 , pi/2 ] .  }
          IF ( inf(x)<0 ) THEN
           
           hx:= - x + ( -1 - 2*trunc(-inf(x/Pi)) ) * Pi/2
          ELSE
           hx:= - x + ( 1 + 2*trunc(inf(x/Pi)) )* Pi/2;
          IF ( sup(hx)>inf(Pi/2) ) THEN 
           BEGIN
             Intervallaufteilung:= TRUE;
             hx2:= intval(-sup(Pi)/2,sup(hx)-inf(Pi)/2);
             hx:= intval(inf(hx),sup(Pi)/2);
           END
          ELSE Intervallaufteilung:= FALSE; 
        END;
       BOTH:= TRUE;
       { Nutze aus dass, tan(transp(z))=transp(tan(z)) fuer cinterval z. }
       { Damit sind nur Argumente in der oberen Halbebene noetig!! }
       IF ( sup(y)<=0) THEN
        BEGIN
          htan(hx,-y,Re_cot,Im_cot,BOTH);
          IF Intervallaufteilung THEN
           BEGIN
             htan(hx2,-y,Recot,Imcot,BOTH);
             Re_cot:= Re_cot +* Recot;
             Im_cot:= -(Im_cot +* Imcot);
           END
          ELSE 
           Im_cot:= -Im_cot;
        END
       ELSE
        IF ( inf(y)>=0) THEN
         IF Intervallaufteilung THEN
          BEGIN
            htan(hx,y,Re_cot,Im_cot,BOTH);
            htan(hx2,y,Recot,Imcot,BOTH);
            Re_cot:= Re_cot +* Recot;
            Im_cot:= Im_cot +* Imcot; 
          END
         ELSE
          htan(hx,y,Re_cot,Im_cot,BOTH)
        ELSE
         IF -inf(y)<sup(y) THEN
          BEGIN
            htan(hx,intval(0,sup(y)),Re_cot,Im_cot,BOTH);
            htan(hx,intval(0,-inf(y)),Recot,Imcot,NOT BOTH);
            Im_cot:= Im_cot +* -Imcot;
            IF Intervallaufteilung THEN
             BEGIN
               htan(hx2,intval(0,sup(y)),Recot,Imcot,BOTH);
               Re_cot:= Re_cot +* Recot;
               Im_cot:= Im_cot +* Imcot;
               htan(hx,intval(0,-inf(y)),Recot,Imcot,NOT BOTH);
               Im_cot:= Im_cot +* -Imcot;
             END;
          END
         ELSE
          BEGIN
            htan(hx,intval(0,-inf(y)),Re_cot,Imcot,BOTH);
            htan(hx,intval(0,sup(y)),Recot,Im_cot,NOT BOTH);
            Im_cot:= Im_cot +* -Imcot;
            IF Intervallaufteilung THEN
             BEGIN
               htan(hx2,intval(0,-inf(y)),Recot,Imcot,BOTH);
               Re_cot:= Re_cot +* Recot;
               Im_cot:= Im_cot +* -Imcot;  
               htan(hx2,intval(0,sup(y)),Recot,Imcot,NOT BOTH);
               Im_cot:= Im_cot +* Imcot;
             END;
          END;
     END;
END;


GLOBAL FUNCTION tan( vari : cinterval ) : cinterval;
VAR re_tan,im_tan : interval;
    Error    : integer;
BEGIN
  h_tan(re(vari),im(vari),re_tan,im_tan,Error);
  IF Error=0 THEN
   tan:= compl(re_tan,im_tan)
  ELSE
   BEGIN
     IF Error=1 THEN
      writeln('***Fehler : Eingabeintervall enthaelt Polstelle des Tangens***')
     ELSE 
      writeln('***Fehler : Realteil des Arguments zu weit von [-pi/2,pi/2] entfernt***');
     Failed;
   END;
END;

 
GLOBAL FUNCTION cot( vari : cinterval ) : cinterval;
VAR re_cot,im_cot : interval;
    Error    : integer;
BEGIN
  h_cot(re(vari),im(vari),re_cot,im_cot,Error);
  IF Error=0 THEN
   cot:= compl(re_cot,im_cot)
  ELSE
   BEGIN
     IF Error=1 THEN
      writeln('***Fehler : Eingabeintervall enthaelt Polstelle des Cotangens***')
     ELSE
      writeln('***Fehler : Realteil des Arguments zu weit von [0,pi] entfernt***');  
    Failed;    
  END;
END;


GLOBAL FUNCTION tanh( vari : cinterval ) : cinterval;
VAR re_tanh,im_tanh : interval;
    Error    : integer;
BEGIN
  h_tan(im(vari),re(vari),im_tanh,re_tanh,Error);
  IF Error=0 THEN
   tanh:= compl(re_tanh,im_tanh)
  ELSE
   BEGIN
     IF Error=1 THEN
      writeln('***Fehler : Polstelle des Tangens hyperbolicus im Eingabeintervall***')
     ELSE
      writeln('***Fehler : Imaginaerteil des Arguments zu weit von [-pi/2,pi/2] entfernt***');
    Failed;  
  END;
END;


GLOBAL FUNCTION coth( vari : cinterval ) : cinterval;
VAR m_re_coth,im_coth : interval;
    PiHalbe         : interval;
    Error           : integer;
BEGIN
  PiHalbe:= arccos(intval(0.0));
  { T:= tan(pi/2-i*z) = tan( pi/2 + im(z) - i * re(z) ) = cot( i*z ) } 
  { coth = i * T = - im(T) + i * re(T) }
  { h_tan(PiHalbe + im(vari),-re(vari),im_coth,m_re_coth,Error); }
  h_cot(-im(vari),re(vari),im_coth,m_re_coth,Error);
    IF Error=0 THEN
     coth:= compl(-m_re_coth,im_coth)
  ELSE
   BEGIN
     IF Error=1 THEN
      writeln('***Fehler : Polstelle des Kotangens hyperbolicus im Eingabeintervall***')
     ELSE
      writeln('***Fehler : Imaginaerteil des Arguments zu weit von [-pi/2,pi/2] entfernt***');
    Failed;
  END; 
END;

{ ************************************************************************** }


{ ************************************************************************** }
{ ************************************************************************** }
{ ***                      Mehrdeutige Funktionen                        *** }
{ ************************************************************************** }
{ ************************************************************************** }


GLOBAL FUNCTION arg( vari : cinterval ) : interval;
VAR PiHalbe : interval;
    wert    : interval; 
    hxl,hyl : interval;
    hxu,hyu : interval;
BEGIN
  IF (intval(0.0)=re(vari)) THEN
   IF (intval(0.0)=im(vari)) THEN
   BEGIN
     writeln('***Fehler: arg nicht fuer Null definiert***');
     Failed;
   END ELSE                         { vari rein imaginaer ungleich Null }
   BEGIN
     PiHalbe:= arccos(intval(0.0));
     IF ( 0<inf(im(vari)) ) THEN arg:= PiHalbe       { im(vari) positiv }
     ELSE IF ( 0>sup(im(vari)) ) THEN arg:= -PiHalbe { im(vari) negativ }
     ELSE arg:= intval(-sup(PiHalbe),sup(PiHalbe));
   END
  ELSE 
  BEGIN    
    PiHalbe:= arccos(intval(0.0));
    hxl:= intval(inf(re(vari))); hxu:= intval(sup(re(vari)));
    hyl:= intval(inf(im(vari))); hyu:= intval(sup(im(vari))); 
    IF ( intval(0.0) in re(vari) ) THEN
     IF ( intval(0.0) in im(vari) ) THEN           { alle Quadranten }
      arg:= 2 * intval(-sup(PiHalbe),sup(PiHalbe)) 
     ELSE 
     BEGIN
       IF ( inf(im(vari))>=0 ) THEN               { 1. und 2.Quadrant }
       arg:= intval(inf(arctan(hyl/hxu)),               { Minimum } 
                     sup(arctan(hyl/hxl)+2*PiHalbe))    { Maximum }
       ELSE                                       { 3. und 4.Quadrant }
       arg:= intval(inf(arctan(hyu/hxl)-2*PiHalbe),     { Minimum } 
                     sup(arctan(hyu/hxu)));             { Maximum }
     END
    ELSE
    BEGIN 
      IF ( intval(0.0) in im(vari) ) THEN
       IF ( inf(re(vari))>=0 ) THEN               { 1. und 4.Quadrant }
        IF ( inf(re(vari))=0 ) THEN
         arg:= intval(-sup(PiHalbe),                    { Minimum }
                       sup(PiHalbe))                    { Maximum }
        ELSE
         arg:= intval(inf(arctan(hyl/hxl)),             { Minimum }
                      sup(arctan(hyu/hxl)))             { Maximum }
       ELSE                              { 2. und 3.Quadrant(Schnitt) }
        IF ( sup(re(vari))=0 ) THEN
        arg:= intval(inf(PiHalbe),                      { Minimum }
                     sup(3*PiHalbe))                    { Maximum }
        ELSE
        arg:= intval(inf(arctan(hyu/hxu)),              { Minimum }
                     sup(arctan(hyl/hxu))) + 2*PiHalbe  { Maximum }
      ELSE                 { Argument liegt in genau einem Quadranten } 
      IF ( inf(im(vari))>=0 ) THEN                  { obere Halbebene }   
       IF ( inf(re(vari))>=0 ) THEN                    { 1.Quadrant }
        IF ( inf(re(vari))=0 ) THEN
        arg:= intval(inf(arctan(hyl/hxu)),              { Minimum }
                      sup(PiHalbe))                     { Maximum }
        ELSE 
        arg:= intval(inf(arctan(hyl/hxu)),              { Minimum }
                      sup(arctan(hyu/hxl)))             { Maximum }
       ELSE                                            { 2.Quadrant }
        IF ( sup(re(vari))=0 ) THEN 
        arg:= intval(inf(PiHalbe),                      { Minimum }
                     sup(arctan(hyl/hxl)+2*PiHalbe))    { Maximum }
        ELSE
        arg:= intval(inf(arctan(hyu/hxu)),              { Minimum }
                      sup(arctan(hyl/hxl))) + 2*PiHalbe { Maximum }
      ELSE                                          { untere Halbebene }
       IF ( inf(re(vari))>=0 ) THEN                    { 4.Quadrant }
        IF ( inf(re(vari))=0 ) THEN
        arg:= intval(-sup(PiHalbe),                     { Minimum }
                      sup(arctan(hyu/hxu)))             { Maximum }
        ELSE
        arg:= intval(inf(arctan(hyl/hxl)),              { Minimum }
                      sup(arctan(hyu/hxu)))             { Maximum }
       ELSE                                            { 3.Quadrant }
        IF ( sup(re(vari))=0 ) THEN
        arg:= intval(-sup(PiHalbe),                     { Minimum }
                      sup(arctan(hxu/hyl)-2*PiHalbe))   { Maximum }
        ELSE
        arg:= intval(inf(arctan(hyu/hxl)),              { Minimum }
                     sup(arctan(hyl/hxu))) - 2*PiHalbe; { Maximum }
    END;
  END;
END;   
    
GLOBAL FUNCTION ln( vari : cinterval ) : cinterval;
VAR re_min,im_min : real;
    re_max,im_max : real;
    re_ln   : interval;
FUNCTION Betrag( x,y : interval ) : interval;
BEGIN
{  Betrag:= sqrt( sqr(x) + sqr(y) ); }
   Betrag:= abs(compl(x,y));
END; 

BEGIN
  IF sup(re(vari))<0 THEN re_min:= -sup(re(vari))
  ELSE IF inf(re(vari))<0 THEN re_min:= 0
  ELSE re_min:=inf(re(vari));
  IF sup(im(vari))<0 THEN im_min:= -sup(im(vari))
  ELSE IF inf(im(vari))<0 THEN im_min:= 0
  ELSE im_min:=inf(im(vari));
  IF ( re_min=0 ) AND ( im_min=0 ) THEN
  BEGIN 
    writeln('*** Fehler: Singularitaet des Logarithmus im Eingabeintervall ***');
    Failed;
  END ELSE
  BEGIN 
    re_max:= max(sup(re(vari)),-inf(re(vari)));
    im_max:= max(sup(im(vari)),-inf(im(vari)));
    re_ln:= intval(inf(ln(Betrag(intval(re_min),intval(im_min)))),
                   sup(ln(Betrag(intval(re_max),intval(im_max)))));
    ln:= compl(re_ln,arg(vari));
  END;
END;

GLOBAL FUNCTION log2( vari : cinterval ) : cinterval;
BEGIN
  log2:= log2(exp(intval(1.0)))*ln(vari);
END;

GLOBAL FUNCTION log10( vari : cinterval ) : cinterval;
BEGIN
  log10:= log10(exp(intval(1.0)))*ln(vari);
END;

FUNCTION sign( x : interval ) : interval;
BEGIN
  sign:= intval(sign(inf(x)),sign(sup(x))); 
END;

FUNCTION re_sqrt( x,y : interval ) : interval;
{ Formeln fuer einzelne Quadranten }
BEGIN
  IF sup(x)<0 THEN
   IF (y=intval(0.0)) THEN re_sqrt:= intval(0.0)
   ELSE re_sqrt:= 1/sqrt(intval(2.0))*abs(y)/sqrt(abs(compl(x,y))-x)
  ELSE 
   IF (y=intval(0.0)) THEN re_sqrt:= sqrt(x)
   ELSE re_sqrt:= 1/sqrt(intval(2.0))*sqrt(abs(compl(x,y))+x);
END;

FUNCTION im_sqrt( x,y : interval ) : interval;
{ Formeln fuer einzelne Quadranten }
BEGIN
  IF sup(x)<0 THEN
   IF (y=intval(0.0)) THEN im_sqrt:= 1/sqrt(-x)
   ELSE im_sqrt:= 1/sqrt(intval(2.0))*sign(y)*sqrt(abs(compl(x,y))-x)
  ELSE
   IF (y=intval(0.0)) THEN im_sqrt:= intval(0.0) 
   ELSE im_sqrt:= 1/sqrt(intval(2.0))*y/sqrt(abs(compl(x,y))+x);
END;

GLOBAL FUNCTION sqrt( vari : cinterval ) : cinterval;
VAR rwert,iwert : interval;
    hxl,hyl : interval;
    hxu,hyu : interval;
BEGIN
  hxl:= intval(inf(re(vari))); hxu:= intval(sup(re(vari)));
  hyl:= intval(inf(im(vari))); hyu:= intval(sup(im(vari)));
  IF ( sup(im(vari))<0 ) THEN { untere Halbebene (ohne Schnitt) }
  BEGIN
    rwert:= re_sqrt(hxl,hyu); { Minimum }
    rwert:= intval(inf(rwert),sup(re_sqrt(hxu,hyl)));{ mit Maximum }
    iwert:= im_sqrt(hxl,hyl); { Minimum }
    iwert:= intval(inf(iwert),sup(im_sqrt(hxu,hyu)));{ mit Maximum }
  END ELSE
  IF ( inf(im(vari))>=0 ) THEN { obere Halbebene }
  BEGIN
    rwert:= re_sqrt(hxl,hyl); { Minimum }
    rwert:= intval(inf(rwert),sup(re_sqrt(hxu,hyu)));{ mit Maximum }
    iwert:= im_sqrt(hxu,hyl); { Minimum }
    iwert:= intval(inf(iwert),sup(im_sqrt(hxl,hyu)));{ mit Maximum }
  END ELSE  { Null liegt im Imaginaerteil des Eingabeintervalls ! }
  IF ( 0<=inf(re(vari)) ) THEN { rechte Halbebene (d.h kein Schnitt) }
  BEGIN
    rwert:= sqrt(hxl); { Minimum } 
    rwert:= intval(inf(rwert)  ,  { mit moeglichen Maxima : }
            max(sup(re_sqrt(hxu,hyu)),sup(re_sqrt(hxu,hyl))));
    iwert:= im_sqrt(hxl,hyl); { Minimum }
    iwert:= intval(inf(iwert),sup(im_sqrt(hxl,hyu)));{ mit Maximum }
  END ELSE 
  IF ( 0>sup(re(vari)) ) THEN { linke Halbebene (mit Schnitt) }
  BEGIN
    rwert:= -re_sqrt(hxu,hyl); { Minima (anderes Vorzeichen, da Schnitt) }
    rwert:= intval(inf(rwert),sup(re_sqrt(hxu,hyu))); { mit Maximum }
    iwert:= sqrt(-hxu); { Minimum, beachte Null in im(vari) }
    iwert:= intval(inf(iwert),            { mit moeglichen Maxima : } 
            max(sup(im_sqrt(hxl,hyu)),   
                -inf(im_sqrt(hxl,hyl)))); { anderes Vorzeichen, wegen Schnitt } 
  END ELSE { alle Quadranten }
  BEGIN                                       { Minimum ( anderes Vor-  }
     rwert:= -1/sqrt(intval(2.0))*sqrt(-hyl); {    zeichen, da Schnitt) }
     { mit moeglichen Maxima }
     rwert:= intval(inf(rwert),
             max(sup(re_sqrt(hxu,hyu)),sup(re_sqrt(hxu,hyl)))); 
     iwert:= -1/sqrt(intval(2.0))*sqrt(-hyl); { Minimum }
     iwert:= intval(inf(iwert),            { mit moeglichen Maxima : }
             max(sup(im_sqrt(hxl,hyu)),   
                 -inf(im_sqrt(hxl,hyl)))); { anderes Vorzeichen, wegen Schnitt }
  END;
  sqrt:= compl(rwert,iwert);
END;

{ ***************************************************************** }
{ *** arcsin uebernohmen aus Diplomarbeit von Gabriele Buehler. *** }
{ *** Werden im Original Langzahlen verwendet, so stehen hier   *** }
{ *** entsprechend Zahlen einfacher Laenge.                     *** } 
{ ***************************************************************** }

FUNCTION pi : interval;
BEGIN
  pi:= 2.0 * arccos(intval(0.0));
END;

function g( s : interval ):interval;
begin
  g:= sqrt( 1.0 + s ) - 1.0;
end;

function s_re( ix,iy : interval ):interval;
begin
  s_re:= 2.0 * iy / ( ix * ix + iy * iy  - 1.0 );
end;


function re_arcsin (x,y:real) :interval;
var ix,iy,hilf1,hilf2,hilf3,nenner,zaehler,sqrs_re:interval;
    null,eins:real;
begin
  ix:=intval(x);
  iy:=intval(y);
  null:= 0.0;
  eins:=1.0;
  { Realteil }
  if x = null then
    re_arcsin:= null
  else
    if y = null then
      if x >= eins then
        re_arcsin:= pi / 2.0
      else
        re_arcsin:= arcsin(ix)
    else
      if (x*x - y*y) < 0.5 then
        begin
          hilf1:= sqrt(sqr(ix+1.0)+sqr(iy));
          hilf2:= sqrt(sqr(ix-1.0)+sqr(iy));
          nenner:= hilf1 + hilf2;
          hilf3:= (2.0 * ix)/nenner;
          re_arcsin:= arcsin(hilf3);
        end
      else                    { x*x - y*y > 0.5 }
        begin
          hilf1:=ix*ix + iy*iy;
          hilf2:=hilf1 - eins;
          nenner:=hilf1 + eins + sqrt( hilf2 * hilf2 + 4.0 * iy * iy);
          sqrs_re:=sqr( s_re(ix,iy) );
          if hilf1.sup < eins then
            begin
              zaehler:=2.0 - 2.0 * ix * ix - hilf2 *g(sqrs_re);
              hilf3:=sqrt(zaehler / nenner);
              re_arcsin:= pi /2.0 - arcsin(hilf3);
            end
          else                { x*x - y*y > 1.0 }
            begin
              zaehler:=2.0 * iy * iy + hilf2 * g(sqrs_re);
              hilf3:=sqrt(zaehler / nenner);
              re_arcsin:= pi /2.0 - arcsin(hilf3);
            end;
        end;
end;

function im_arcsin (x, y:real) :interval;
{ Berechnung des Imaginaerteilwertes des arcsin(z) intervallmaessig. d.h, }
{ Eingabeargumente sind Punktintervalle mit den Formeln nach Kraemer      }
{ Seite 126 }
var ix,iy,hilf1,hilf2,hilf3,hilf4,hilf5,hilf6,nenner,zaehler:interval;
    t,r:interval;
    null,eins:real;
begin
  ix:=intval(x);
  iy:=intval(y);
  null:=0.0;
  eins:=1.0;
  x:=abs(x);
  ix:=   abs(ix);
  { Imaginaerteil }
  if y = null then           { y = 0.0 }
    if x<= eins then
      im_arcsin:= null
    else                           { x> 1.0 }
      if x > 1.1 then       { x> 1.1 }
        begin
          t:= 0.5 * ( (ix + eins) + (ix - eins) );
          im_arcsin := ln(t + sqrt(sqr(t) - eins));
        end
      else                          { 1< x < 1.1 }
        begin
          t:= ix;
          r:= t - eins;
          im_arcsin := ln(eins + (r + sqrt(sqr(r) + 2.0 * r)));
        end
  else                       {y <> 0.0 }
    begin
      hilf1:=ix*ix + iy*iy;
      hilf2:=hilf1 + eins;
      if x = null then         { x= 0.0 }
        im_arcsin := ln(sqrt(eins + iy * iy) + iy)
      else                   { x <> 0.0 }
        begin
          t := 0.5 * ( sqrt( hilf2 + 2.0 * ix) + sqrt( hilf2 - 2.0 * ix ));
          if t.sup <= 1.1 then   { t <= 1.1 }
            begin
              if x = eins then   { x = 1.0 }
                r:=  g(iy * iy /4.0) + 0.5 * iy
              else             { x <>1.0 }
                begin
                  hilf3:=iy/(ix + eins);
                  hilf4:= (ix + eins) * g(hilf3 * hilf3);
                  if x < eins then { x < 1.0 }
                    begin
                      hilf5 := iy/( eins -ix);
                      hilf6 := (eins - ix) * g(hilf5*hilf5);
                      r:= 0.5 * ( hilf4 + hilf6);
                    end
                  else              { x > 1.0 }
                    r := 0.5 * (hilf4 + (ix - eins) + sqrt(hilf2 -2.0 * ix));
                end;
              im_arcsin := ln(eins + (r + sqrt(sqr(r) + 2.0 * r)));
            end
          else                    { t > 1.1 }
            im_arcsin := ln(t + sqrt(sqr(t) - eins));
        end;
    end;
end;

function fortsetz_asin(x,y:real):interval;
var hilf,null:interval;
begin
  null:=0.0;
  hilf:=re_arcsin(x,y);
  if hilf = null then
    fortsetz_asin := pi
  else
    fortsetz_asin := pi - hilf;
end;

global function real_asin(c:cinterval):interval;
var xl,xu,yl,yu,maxy,max:real;
    re_spiegel:boolean;
    null,eins:real;
    ergxl,ergxu,ergx:interval;
    c1:cinterval;
begin
  c1:=c;
  z(c1,xl,xu,yl,yu);
  { Die Funktion z weisst den reellen Zahlen xl,xu,yl,yu die entsprechenden }
  { Zahlen des komplexen Intervalls c1 zu }

  null:=0.0;
  eins:=1.0;
  { Falls das Eingaberechteck in der linken Halbebene liegt wird es }
  { nach rechts gespiegelt }
  if xu <= null then
    begin
      c1.re := -c1.re;
      re_spiegel:=true;
    end
  else
    re_spiegel:=false;
  { Falls das Eingaberechteck in der unteren Halbebene liegt wird es }
  { nach oben gespiegelt}
  if yu <=null then
    c1.im:=-c1.im;
  z(c1,xl,xu,yl,yu);
  if -yl >= yu then  { Berechnet das Maximum von yl und yu (Betrag) }
    maxy:=-yl
  else
    maxy:=yu;

  { Hier beginnt die Fallunterscheidung }
  if yl >= null then      { obere Halbebene }
    if xl >= null then    { 1. Quadrant }
      begin
        ergxl:=re_arcsin(xl,yu);
        ergxu:=re_arcsin(xu,yl);
      end
    else                  { 0 in Innern das Realteil von c }
      begin
        ergxl:=-re_arcsin(-xl,yl);
        ergxu:=re_arcsin(xu,yl);
      end
  else                   { 0 im Innern des Imaginaerteils von c }
    if xl >= eins then   { Schnitt ohne Windungspunkt im Tntervall }
      begin
        ergxl:= re_arcsin(xl,yl);
        ergxu:= fortsetz_asin(xl,yu);
      end
    else                  { xl < 1.0 }
      begin
        if xl < -eins then { xl kleiner als unterer Windungspunkt }
          if xu < eins then  { nur Windungspunkt (-1,0) im Intervall }
            begin
              ergxl:= -pi / 2.0;
              ergxu:= re_arcsin(xu,null);
            end
          else                { beide Windungspunkte im Intervall }
            begin
              ergxl:= -pi / 2.0;
              ergxu:= pi / 2.0;
            end
        else                   { -1 <= xl < 1 }
          if xl <= null then   { xl <= 0 }
            if xu <= eins then { enthaelt keinen Schnitt aber Ursprung }
              begin
                ergxl:= -re_arcsin(-xl,null);
                ergxu:= re_arcsin(xu,null);
              end
            else               { enthält windungspunkt (1.0) und Ursprung }
              begin
                ergxl:= -re_arcsin(-xl,null);
                ergxu:= pi / 2.0;
              end
          else               { 0 <= xl < 1 }
            begin
              if xu < eins then { enthaelt keinen WP und Ursprung }
                begin
                  ergxl:= re_arcsin(xl,maxy);
                  ergxu:= re_arcsin(xu,null);
                end
              else            { enthaelt WP (1,0) }
                begin
                  ergxl:= re_arcsin(xl,max);
                  ergxu:= pi / 2.0;
                end
            end;
      end;
  ergx:=intval(inf(ergxl),sup(ergxu));
  if re_spiegel then
    real_asin:=-ergx
  else
    real_asin:=ergx;
end;

global function imag_asin(c:cinterval):interval;
var xl,xu,yl,yu,maxx,maxy,max:real;
    im_spiegel:boolean;
    null,eins:real;
    ergyl,ergyu,ergy:interval;
    c1:cinterval;
begin
  c1:=c;
  z(c1,xl,xu,yl,yu);
  null:=0.0;
  eins:=1.0;
  { Falls das Eingaberechteck in der linken Halbebene liegt wird es }
  { nach rechts gespiegelt }
  if xu <= null then
    c1.re:= -c1.re;
    { Falls das Eingaberechteck in der unteren Halbebene liegt wird es }
    { nach oben gespiegelt }
   if yu <= null then
     begin
       c1.im:=-c1.im;
       im_spiegel:=true;
     end
   else
     im_spiegel:=false;
   z(c1,xl,xu,yl,yu);
   if -xl>= yu then  { Berechnet das Maximum von xl und xu (Betrag) }
     maxx:=-xl
   else
     maxx :=xu;
   if -yl >= yu then { Berechnet das Maximum von yl und yu (Betrag) }
     maxy:=-yl
   else
     maxy:=yu;

   { Hier beginnt die Fallunterscheidung }
   if yl >= null then            { obere Halbebene }
     if xl >null then            { 1. Quadrant }
       begin
         ergyl:=im_arcsin(xl,yl);
         ergyu:=im_arcsin(xu,yu);
       end
     else                   { 0 im Innern des Realteil von c }
       begin
         ergyl:=im_arcsin(null,yl);
         ergyu:=im_arcsin(maxx,yu);
       end
   else                     { 0 in Innern des Imaginaerteils von c }
     if xl >= eins then     { Schnitt ohne Windungspunkt im Intervall }
       begin
         ergyl:=-im_arcsin(xu,maxy);
         ergyu:= im_arcsin(xl,null);
       end
     else                   { xl < 1.0 }
       begin
         ergyl:=-im_arcsin(maxx,-yl);
         ergyu:= im_arcsin(maxx,yu);
       end;

  ergy:=intval(inf(ergyl) ,sup(ergyu));

  if im_spiegel then
    imag_asin:=-ergy
  else
    imag_asin:=ergy;
end;


{ ************************************************************************** }
{ *** Implementierung von arccos, arsinh und arcosh basierend auf arcsin *** }
{ ************************************************************************** }
{ ***     Aus den obigen Formeln zur Berechnung von sin, cos, sinh       *** }
{ ***     und cosh kann man folgende Beziehungen ermitteln               *** }
{ ***                                                                    *** }
{ ***  (1)   sin(z)  = - sin(-z)        ,                                *** }
{ ***  (2)   cos(z)  =   cos(-z)                                         *** }
{ ***  (3)           =  sin( pi/2 -z )   ,                               *** }
{ ***  (4)   sinh( i * z ) = i * sin(z)    und                           *** }
{ ***  (5)   cosh(z) =  cos( i * z )     .                               *** }
{ ***                                                                    *** }
{ ***     Hieraus lassen sich Beziehungen der zugehoerigen ( nicht       *** }
{ ***     eindeutigen ) Umkehrfunktionen arcsin, arccos, arsinh und      *** }
{ ***     arcosh herleiten, so das mit einer dieser Funktionen alle      *** }
{ ***     weiteren einfach zu berechnen sind:                            *** }
{ ***                                                                    *** }
{ ***     Es gilt                                                        *** }
{ ***                          (2)                                       *** }
{ ***        z = cos(Arccos(z)) = cos(-Arccos(z))                        *** }
{ ***         (3)                                                        *** }
{ ***          = sin(pi/2 + Arccos(z))                                   *** }
{ ***                                                                    *** }
{ ***          = sin(Arcsin(z))       ,                                  *** }
{ ***                                                                    *** }
{ ***     also  pi/2 +/- Arccos(z) = Arcsin(z) ( mod 2*pi )              *** }
{ ***     oder                                                           *** }
{ ***                                                                    *** }
{ ***  (6)    Arccos(z) = -/+ ( Arcsin(z) - pi/2 )   ,                   *** }
{ ***         ------------------------------------                       *** }
{ ***                                                                    *** }
{ ***                                                                    *** }
{ ***        z = sinh(Arsinh(z))                                         *** }
{ ***                                                                    *** }
{ ***          = i * -i * z = i * sin(Arcsin( -i * z ))                  *** }
{ ***         (4)                                                        *** }
{ ***          = sinh( i * Arcsin( -i * z ) )    ,                       *** }
{ ***                                                                    *** }
{ ***     also                                                           *** }
{ ***                                                                    *** }
{ ***  (7)   Arsinh(z) = i * Arcsin( -i * z ) ( mod i*2*pi )             *** }
{ ***        --------------------------------                            *** }
{ ***                                                                    *** }
{ ***     und                                                            *** }
{ ***                           (6)                                      *** }
{ ***        z = cosh(Arcosh(z)) = cos(i*Arcosh(z))                      *** }
{ ***                           (2)                                      *** }
{ ***          = cos(Arccos(z))  = cos(-Arccos(z))                       *** }
{ ***                                                                    *** }
{ ***                                                                    *** }
{ ***     also  i * Arcosh(z) = +/- Arccos(z) ( mod 2*i )                *** }
{ ***     oder                                                           *** }
{ ***                                                                    *** }
{ ***  (8)    Arcosh(z) = -/+ i * Arccos(z)     .                        *** }
{ ***         -----------------------------                              *** }
{ ***                                                                    *** }
{ ************************************************************************** }
{ ***    Die Funktionen Arccos und Arcosh sind bis auf den Faktor -1     *** }
{ ***    und   mod 2*pi   bzw.   mod i*2*pi   eindeutig bestimmt.        *** }
{ ***    Die Funktionen Arcsin und Arsinh sind bis auf    mod 2*pi       *** }
{ ***    bzw.   mod i*2*pi   eindeutig bestimmt. Deshalb werden nur      *** }
{ ***    sogenannte Hauptwerte berechnet. Diese sind eindeutig bestimmte *** }
{ ***    Vertreter der moeglichen Werte. Liegen die Hauptwerte eines     *** }
{ ***    komplexen Intervalls auf disjunkten Mengen, so versucht man     *** }
{ ***    ein moeglichst kleines Intervall zu bestimmen, das u.U. andere  *** }
{ ***    Vertreter als die Hauptwerte enthaelt.                          *** }
{ ************************************************************************** }
{ *** verwendete Formeln:   arccos(z) = +/- (pi/2 - arcsin(z))           *** }
{ ***                       arcosh(z) = -/+ i * arccos(z)                *** }
{ ***                       arsinh(z) =  -i * arcsin(i*z)                *** }
{ ************************************************************************** }

GLOBAL FUNCTION arcsin( c : cinterval): cinterval;
BEGIN
   arcsin:= compl(real_asin(c),imag_asin(c));
END;

GLOBAL FUNCTION arccos( c : cinterval): cinterval;
VAR PiHalbe : interval;
BEGIN
   PiHalbe:= arccos(intval(0.0));
   arccos:= compl(real_asin(c)-PiHalbe,imag_asin(c));
END;

GLOBAL FUNCTION arsinh( c : cinterval): cinterval;
VAR hc : cinterval;
BEGIN
   hc:= compl(-im(c),re(c));                     { hc:= i*c }
   arsinh:= compl(imag_asin(hc),-real_asin(hc)); { arsinh(c):= -i*arcsin(i*c) }
END;

GLOBAL FUNCTION arcosh( c : cinterval): cinterval;
VAR hc : cinterval;
BEGIN
   hc:= arccos(c);
   arcosh:= compl(-im(hc),re(hc));               { arcosh:= i * arccos(c) }
END;

{ ************************************************************************** }


{ ************************************************************************** }
{ ***  Hilfsprozeduren zur Berechnung einer Umkehrfunktion des Tangens   *** }
{ ************************************************************************** }

FUNCTION qbetrag( lx,ux,ly,uy : real ) : interval;
VAR infqbetr,supqbetr : real;
BEGIN
  IF ( lx>=0 ) THEN
   BEGIN
     IF ( ly>0 ) THEN
      BEGIN
        infqbetr:=#<( lx*lx + ly*ly );
        supqbetr:=#>( ux*ux + uy*uy );
      END
     ELSE
      IF ( uy<0 )THEN
       BEGIN
         infqbetr:=#<( lx*lx + uy*uy );
         supqbetr:=#>( ux*ux + ly*ly );
       END
      ELSE
       BEGIN
         infqbetr:=#<( lx*lx );
         IF (-ly<uy) THEN
          supqbetr:=#>( ux*ux +  uy*uy )
         ELSE
          supqbetr:=#>( ux*ux + ly*ly );
       END;
     qbetrag:= intval(infqbetr,supqbetr);
   END
  ELSE
   IF ( ux<=0 ) THEN 
    qbetrag:= qbetrag( -ux,-lx,ly,uy )
   ELSE
    IF ( -lx<ux ) THEN
     qbetrag:= qbetrag( 0.0,ux,ly,uy )
    ELSE
     qbetrag:= qbetrag( 0.0,-lx,ly,uy );
END; 

PROCEDURE harctan( x,y : interval; VAR re_arct,im_arct : interval );
VAR Pi        : interval; 
    h_re,h_im : interval;
{}
FUNCTION re_fun( hx,hy : real ) : interval;
VAR PiHalbe : interval;
    qbetr   : interval;
BEGIN
  qbetr:= qbetrag(hx,hx,hy,hy);
  IF 1>sup(qbetr) THEN
   re_fun:= arctan( intval(2*hx)/(1-qbetr) ) / 2
  ELSE
   IF 1<inf(qbetr) THEN
    re_fun:= (arctan( intval(2*hx)/(1-qbetr) ) - Pi )/2 
   ELSE
    re_fun:= -Pi/4;
END;
{}
FUNCTION re_fun( hx : interval; hy : real ) : interval;
VAR PiHalbe : interval;
    qbetr   : interval;
BEGIN
  qbetr:= qbetrag(inf(hx),sup(hx),hy,hy);
  IF 1>sup(qbetr) THEN
   re_fun:= arctan( 2*hx/(1-qbetr) ) / 2
  ELSE
   IF 1<inf(qbetr) THEN 
    re_fun:= (arctan( 2*hx/(1-qbetr) ) - Pi )/2
   ELSE 
    re_fun:= -Pi/4 ;
END;
{}
FUNCTION im_fun( hx,hy : real ) : interval;
VAR nenner : interval;
    ly,uy  : real;
BEGIN
  ly:= #< ( 1 - hy );
  uy:= #> ( 1 - hy );
  nenner:= qbetrag(hx,hx,ly,uy);
  IF ( 0 < inf(nenner) ) THEN
   im_fun:= ln( abs(1+4*intval(hy)/nenner) ) / 4
  ELSE
   BEGIN
     writeln('*** Fehler : Eingabeintervall zu nah an Singularitaet ***');
     Failed;  
   END;
END;
{}
FUNCTION im_fun( hx : real; hy : interval ) : interval;
VAR nenner : interval;
    ly,uy  : real;
BEGIN
  ly:= #< ( 1 - sup(hy) );
  uy:= #> ( 1 - inf(hy) );
  nenner:= qbetrag(hx,hx,ly,uy);
  IF ( 0 < inf(nenner) ) THEN  
   im_fun:= ln(abs(1+4*hy/(sqr(intval(hx))+sqr(1-hy)))) / 4
   ELSE
   BEGIN
     writeln('*** Fehler : Eingabeintervall zu nah an Singularitaet ***');
     Failed;
   END;
END;
{}
PROCEDURE left_side( hx,hy : interval; VAR re_a,im_a : interval);
BEGIN
  IF sup(sqr(intval(sup(hy))))<= inf( 1.0 + sqr(intval(sup(hx)))) THEN
  { Eingabeintervall unter Extremalkurve y=sqrt(1+sqr(x)) }
   BEGIN
     re_a:= intval( inf(re_fun(inf(hx),sup(hy))),   { Minimum }
                    sup(re_fun(sup(hx),inf(hy))) ); { Maximum }
     im_a:= intval( inf(im_fun(inf(hx),inf(hy))),   { Minimum } 
                    sup(im_fun(sup(hx),sup(hy))) ); { Maximum }
   END
  ELSE
   IF inf(sqr(intval(inf(hy))))>= sup( 1.0 + sqr(intval(inf(hx)))) THEN
   { Eingabeintervall ueber Extremalkurve y=sqrt(1+sqr(x)) }
    BEGIN
      re_a:= intval( inf(re_fun(sup(hx),sup(hy))),   { Minimum }
                     sup(re_fun(inf(hx),inf(hy))) ); { Maximum }
      im_a:= intval( inf(im_fun(inf(hx),sup(hy))),   { Minimum } 
                     sup(im_fun(sup(hx),inf(hy))) ); { Maximum }
    END
   ELSE { Eingabeintervall schneidet Extremalkurve y=sqrt(1+sqr(x)) }
    BEGIN
      re_a:= re_fun(sup(hx),sup(hy));          { evtl. Minimum }
      re_a:= re_a +* re_fun(inf(hx),sup(hy));  { evtl. Minimum }
      re_a:= re_a +* re_fun(hx,inf(hy));       { Maximum }
      im_a:= im_fun(inf(hx),inf(hy));          { evtl. Minimum }
      im_a:= im_a +* im_fun(inf(hx),sup(hy));  { evtl. Minimum }
      im_a:= im_fun(sup(hx),hy);               { Maximum }
    END;
END;
{}
PROCEDURE right_side( hx,hy : interval; VAR re_a,im_a : interval );
{ Verwende transp(arctan(z))=arctan(transp(z)) }
BEGIN
  left_side(-x,y,re_a,im_a);                       
  re_a:=-re_a;
END;
{}
BEGIN
  Pi:= arccos(intval(-1));
  IF inf(x)>=0 THEN { 1.Quadrant } 
   right_side(x,y,re_arct,im_arct)
  ELSE
   IF sup(x)<=0 THEN { 2.Quadrant }
    left_side(x,y,re_arct,im_arct)
   ELSE
    IF sup(y)<1.0 THEN
     BEGIN
       re_arct:= intval( inf(re_fun(inf(x),sup(y))),
                         sup(re_fun(sup(x),sup(y))) ) ;
       im_arct:= intval( inf(im_fun(max(sup(x),-inf(x)),inf(y))),
                         sup(im_fun(0.0,sup(y))) ) ;
     END
    ELSE { Schnitt im Eingabeintervall }
     BEGIN
       right_side(intval(0,sup(x)),y,re_arct,im_arct);
       left_side(intval(inf(x),0),y,h_re,h_im);
       re_arct:= re_arct +* ( h_re + Pi );
       im_arct:= im_arct +* h_im; 
     END;
END;


PROCEDURE h_arctan( x,y : interval; VAR re_arct,im_arct : interval;
                                             VAR singulaer : boolean );
VAR z      : cinterval;
    imarct,rearct : interval;
BEGIN
  IF ( (inf(x)<=0) AND (sup(x)>=0) ) AND 
   ( ((inf(y)<=-1) AND (sup(y)>=-1)) OR ((inf(y)<=1) AND (sup(y)>=1)) ) THEN
   singulaer:= TRUE
  ELSE
   BEGIN
     singulaer:= FALSE;
     IF (y=intval(0) ) THEN
      BEGIN
        re_arct:= arctan(x);
        im_arct:= y;  { im(arctan) gleich NULL }
      END
     ELSE
      IF ( x=intval(0) ) AND ( inf(y)>-1 ) AND ( sup(y)<1 ) THEN
       BEGIN
         re_arct:= x; { re(arctan) gleich NULL }
         im_arct:= artanh(y); { Nur fuer  -1 < y < 1  definiert }
       END
      ELSE
       BEGIN

         { Verwende Arctan(-z) = - Arctan(z) } 
         IF inf(y)>=0 THEN
          harctan(x,y,re_arct,im_arct)
         ELSE  
          IF sup(y)<=0 THEN
           BEGIN
             harctan(-x,-y,re_arct,im_arct); 
             re_arct:= - re_arct;
             im_arct:= - im_arct;
           END
          ELSE
           BEGIN
             harctan(-x,intval(0,-inf(y)),rearct,imarct); 
             harctan(x,intval(0,sup(y)),re_arct,im_arct);
             re_arct:= re_arct +* -rearct;
             im_arct:= im_arct +* -imarct;
           END; 
       END;
   END;
END;

{ ************************************************************************** }
{ *** Implementierung von arccot, artanh und arcoth basierend auf arctan *** }
{ ************************************************************************** }
{ ***    Wie im vorigen Fall gelten gewisse Beziehungen zwischen den     *** }
{ ***    Funktionen tan, cot, tanh, coth :                               *** }
{ ***                                                                    *** }
{ ***  (1)      tan(z)  = - tan(-z)                                      *** }
{ ***  (2)      cot(z)  =  tan( pi/2 - z )                               *** }
{ ***  (3)              = - cot(-z)                                      *** }
{ ***  (4)      tanh(z) = -i * tan( i * z )                              *** }
{ ***  (5)      coth(z) = i * cot( i * z )                               *** }
{ ***                                                                    *** }
{ ***    Wieder lassen sich entsprechende Beziehungen fuer die Umkehr-   *** }
{ ***    funktionen arctan, arccot, artanh und arcoth herleiten:         *** }
{ ***                                                                    *** }
{ ***         z = Arccot(cot(z))                                         *** }
{ ***                                                                    *** }
{ ***           = z - pi/2 + pi/2 = -Arctan(tan( pi/2 -z )) + pi/2       *** }
{ ***          (2)                                                       *** }
{ ***           = -Arctan(cot(z)) + pi/2     ,                           *** }
{ ***                                                                    *** }
{ ***    also ( ersetze cot(z) durch z )                                 *** }
{ ***                                                                    *** }
{ ***  (6)      Arccot(z) = pi/2 - Arctan(z)   ,                         *** }
{ ***           ----------------------------                             *** }
{ ***                                                                    *** }
{ ***                                                                    *** }
{ ***         z = Artanh(tanh(z))                                        *** }
{ ***                                                                    *** }
{ ***           = -i * i * z = -i * Arctan(tan( i * z ))                 *** }
{ ***          (4)                                                       *** }
{ ***           = -i * Arctan( i * tanh(z) )    ,                        *** }
{ ***                                                                    *** }
{ ***    also ( ersetze tanh(z) durch z )                                *** }
{ ***                                                                    *** }
{ ***  (7)      Artanh(z) = -i * Arctan( i * z )    ,                    *** }
{ ***           --------------------------------                         *** }
{ ***                                                                    *** }
{ ***                                                                    *** }
{ ***         z = coth(Arcoth(z))                                        *** }
{ ***                                                                    *** }
{ ***           = i * -i * z = i * cot(Arccot( -i * z ))                 *** }
{ ***          (5)                                                       *** }
{ ***           = coth( -i * Arccot( -i * z ) )                          *** }
{ ***                                                                    *** }
{ ***    also  ( mit (3) gilt Arccot(z) = - Arccot(-z) )                 *** }
{ ***                                                                    *** }
{ ***  (8)      Arcoth(z) = i * Arccot( i * z )     .                    *** }
{ ***           -------------------------------                          *** }
{ ***                                                                    *** }
{ ************************************************************************** }
{ ***   Der Bildraum des Tangens ueber den komplexen Zahlen nimmt die    *** }
{ ***   Werte +/- i nicht an. Entsprechend sind dies singulaere Werte    *** }
{ ***   der Umkehrfunktionen arctan, sowie arccot ( siehe (2) ).         *** }
{ ***   Nach (4) und (5) sind +/- 1 singulaere Werte von den Umkehr-     *** }
{ ***   funktionen arctanh und arccot.                                   *** }
{ ************************************************************************** }
{ *** Formeln:     arccot(z) = arctan(-z) +/-  pi/2  , z ungleich +/- i  *** }
{ ***              artanh(z) = -i * arctan( i * z )  , z ungleich +/- 1  *** }
{ ***              arcoth(z) = i * arccot( i * z )   ,  ----- ''-----    *** }
{ ************************************************************************** }

GLOBAL FUNCTION arctan( vari : cinterval ) : cinterval;
VAR im_arct,re_arct : interval;
    singulaer       : boolean;
BEGIN
  h_arctan(re(vari),im(vari),re_arct,im_arct,singulaer);
  IF singulaer THEN
   BEGIN
     writeln('***Fehler : Singularitaet des arctan im Eingabeintervall ***');
     Failed;
   END
  ELSE
   arctan:= compl(re_arct,im_arct);
END;

PROCEDURE h_arccot( x,y : interval; VAR re_arcc,im_arcc : interval;
                       VAR singulaer : boolean );
VAR PiHalbe         : interval;
BEGIN
  PiHalbe:= arccos(intval(0.0));
  h_arctan(x,y,re_arcc,im_arcc,singulaer);
  IF NOT singulaer THEN
   IF inf(re_arcc) >= 0.0 THEN
    re_arcc:= -re_arcc - PiHalbe
   ELSE
    re_arcc:= PiHalbe - re_arcc;
END;


GLOBAL FUNCTION arccot( vari : cinterval ) : cinterval;
VAR im_arcc,re_arcc : interval;
    singulaer       : boolean;
BEGIN
  h_arccot(re(vari),im(vari),re_arcc,im_arcc,singulaer);
  IF singulaer THEN
   BEGIN
     writeln('***Fehler : Singularitaet des arccot im Eingabeintervall ***');
     Failed;
   END
  ELSE
   arccot:= compl(re_arcc,im_arcc);
END;


GLOBAL FUNCTION artanh( vari : cinterval ) : cinterval;
VAR im_art,m_re_art : interval;
    singulaer       : boolean;
BEGIN
  h_arctan(im(vari),-re(vari),im_art,m_re_art,singulaer);
  IF singulaer THEN
   BEGIN
     writeln('***Fehler : Singularitaet des artanh im Eingabeintervall ***');
     Failed;
   END
  ELSE
   artanh:= compl(-m_re_art,im_art);
END;

GLOBAL FUNCTION arcoth( vari : cinterval ) : cinterval;
VAR im_arco,m_re_arco : interval;
    singulaer         : boolean;
BEGIN
  h_arccot(-im(vari),re(vari),im_arco,m_re_arco,singulaer); 
  IF singulaer THEN
   BEGIN
     writeln('***Fehler : Singularitaet des arcoth im Eingabeintervall ***');
     Failed;
   END
  ELSE
   arcoth:= compl(-m_re_arco,im_arco); 
END;

{ ************************************************************************** }
{ *** Die eigentlich eindeutige Operation 'Hoch'  pow  wird erst hier de-*** }
{ *** finiert, da die Implementierung auf  ln  zurueckgreift.            *** }
{ ************************************************************************** }


GLOBAL FUNCTION power( bas : cinterval; expo : integer ) : cinterval;
BEGIN
  IF ( inf(re(bas))>0 ) OR ( inf(im(bas))>0 )
   OR ( sup(re(bas))<0 ) OR ( sup(im(bas))<0 ) THEN
   IF (expo=0) THEN
    power:= compl(intval(1),intval(0))
   ELSE
    BEGIN
      IF odd(expo) THEN
       IF (expo<0) THEN
        power:= 1/sqr(power(bas,-expo div 2)) / bas
       ELSE
        power:= sqr(power(bas,expo div 2)) * bas
      ELSE
       IF (expo<0) THEN
        power:= 1/sqr(power(bas,-expo div 2))
       ELSE
        power:= sqr(power(bas,expo div 2));
    END
  ELSE
   BEGIN
     writeln('*** Fehler : Null liegt in Basis ***');
     Failed;
   END;
END;

GLOBAL FUNCTION power( bas : cinterval; expo : interval ) : cinterval;
BEGIN
  IF ( inf(re(bas))>0 ) OR ( inf(im(bas))>0 )
     OR ( sup(re(bas))<0 ) OR ( sup(im(bas))<0 ) THEN
   power:= exp(expo*ln(bas))
  ELSE
   BEGIN
     writeln('*** Fehler : Null liegt in Basis ***');
     Failed;
   END;
END;

GLOBAL FUNCTION power( bas : cinterval; expo : cinterval ) : cinterval;
BEGIN
  IF ( inf(re(bas))>0 ) OR ( inf(im(bas))>0 )
     OR ( sup(re(bas))<0 ) OR ( sup(im(bas))<0 ) THEN
   power:= exp(expo*ln(bas))
  ELSE
   BEGIN
     writeln('*** Fehler : Null liegt in Basis ***');
     Failed;
   END;
END;

PRIORITY pow = *;

GLOBAL OPERATOR pow(bas : cinterval; expo : integer) icomplPow : cinterval;
BEGIN
  icomplPow:= power( bas,expo);
END;

GLOBAL OPERATOR pow(bas : cinterval; expo : interval) Value : cinterval;
BEGIN
  Value:= power( bas,expo);
END;

GLOBAL OPERATOR pow(bas : cinterval; expo : cinterval) Value : cinterval;
BEGIN
  Value:= power( bas,expo);
END;

END.