// -----------------------------------------------------------------
// Complex functions for Formulae Compilator 2.02 and greater
// Fractal Explorer
// <c> 2000-2004, Sirotinsky A, Fedorenko O.
//
// =================================================================
// Please, do not modify this code without authors permission !
// =================================================================
Unit Complex;

Interface

  Uses Math;

  Type
    TComplex = record
      real: Extended;
      imag: Extended;
    end;

  Function  MakeComplex(real, imag: Extended): TComplex;
  Procedure SetResult(var x,y: Extended; complex: TComplex);
  Function  CModule(Cmp1: TComplex): Extended;

  Procedure CAddV(var Cmp1: TComplex; Cmp2: TComplex);     //  Z:=Z+A
  Function  CAdd (Cmp1, Cmp2: TComplex): TComplex;         //  V:=Z+A
  Procedure CSubV(var Cmp1: TComplex; Cmp2: TComplex);     //  Z:=Z-A
  Function  CSub (Cmp1, Cmp2: TComplex): TComplex;         //  V:=Z-A
  Procedure CMulV(var Cmp1: TComplex; Cmp2: TComplex);     //  Z:=Z*A
  Function  CMul (Cmp1, Cmp2: TComplex): TComplex;         //  V:=Z*A
  Procedure CDivV(var Cmp1: TComplex; Cmp2: TComplex);     //  Z:=Z/A
  Function  CDiv (Cmp1, Cmp2: TComplex): TComplex;         //  V:=Z/A
  // ++19/08/2000
     Function  CAddR(Cmp1: TComplex; t: Extended): TComplex;  //  V:=Z+var
     Function  CSubR(Cmp1: TComplex; t: Extended): TComplex;  //  V:=Z-var
     Function  CMulR(Cmp1: TComplex; t: Extended): TComplex;  //  V:=Z*var
     Function  CDivR(Cmp1: TComplex; t: Extended): TComplex;  //  V:=Z/var


  Procedure CSqrV(var Cmp1: TComplex);                           //  Z:=Z*Z
  Function  CSqr (Cmp1: TComplex): TComplex;               //  V:=Z*Z
  Procedure CTripleV(var Cmp1: TComplex);                        //  Z:=Z*Z*Z
  Function  CTriple (Cmp1: TComplex): TComplex;            //  V:=Z*Z*Z
  Procedure CFourV(var Cmp1: TComplex);                          //  Z:=Z*Z*Z*Z
  Function  CFour (Cmp1: TComplex): TComplex;              //  V:=Z*Z*Z*Z

  Procedure CFlipV(var Cmp1: TComplex);
  Function  CFlip (Cmp1: TComplex): TComplex;
  Procedure CRevV (var Cmp1: TComplex);                          //  Z:=1/Z
  Function  CRev  (Cmp1: TComplex): TComplex;              //  V:=1/Z
  Procedure CRev2V(var Cmp1: TComplex; Cmp2: TComplex);    //  Z:=1/(Z-A)
  Function  CRev2 (Cmp1,Cmp2: TComplex): TComplex;         //  V:=1/(Z-A)
  Function  CConj (Cmp1: TComplex): TComplex;
  Procedure CConjV(var Cmp1: TComplex);

  Function  CAbs  (Cmp1: TComplex): TComplex;              //  ++19/12/2001
  Procedure CAbsV (var Cmp1: TComplex);
  Function  CAbs2 (Cmp1: TComplex): TComplex;              //  ++19/12/2001
  Procedure CAbs2V(var Cmp1: TComplex);


  Function  CReal (Cmp1: TComplex): TComplex;              //  CReal ++10/12/2001
  Procedure CRealV(var Cmp1: TComplex);                          //  suggested by Pete
  Function  CImag (Cmp1: TComplex): TComplex;              //  CImag ++10/12/2001
  Procedure CImagV(var Cmp1: TComplex);                          //  suggested by Pete
  Procedure CZeroV(var Cmp1: TComplex);


  Procedure CSqrtV(var Cmp1: TComplex);                          //  Z:=Sqrt(Z)
  Function  CSqrt (Cmp1: TComplex): TComplex;                    //  V:=Sqrt(Z)
  Procedure CExpV (var Cmp1: TComplex);                          //  Z:=Exp(Z)
  Function  CExp  (Cmp1: TComplex): TComplex;                    //  V:=Exp(Z)
  Procedure CLnV  (var Cmp1: TComplex);                          //  Z:=Ln(Z)
  Function  CLn   (Cmp1: TComplex): TComplex;                    //  V:=Ln(Z)
  Procedure CPowerV(var Cmp1: TComplex; Cmp2: TComplex);         //  Z:=Z^A
  Function  CPower (Cmp1,Cmp2: TComplex): TComplex;              //  V:=Z^A
  // ++19/08/2000
     Procedure CPowerRV(var Cmp1: TComplex; t: Extended);        //  Z:=Z^var
     Function  CPowerR (Cmp1:TComplex;t: Extended): TComplex;    //  V:=Z^var

  // ++12.01.2005 fixed versions, suggested by RaDaR3D
     Procedure CPowerV2(var Cmp1: TComplex; Cmp2: TComplex);     //  Z:=Z^A
     Function  CPower2 (Cmp1,Cmp2: TComplex): TComplex;          //  V:=Z^A
     Procedure CPowerRV2(var Cmp1: TComplex; t: Extended);        //  Z:=Z^var
     Function  CPowerR2 (Cmp1:TComplex;t: Extended): TComplex;    //  V:=Z^var


  Procedure CSinV  (var Cmp1: TComplex);                         //  Z:=Sin(Z)
  Function  CSin   (Cmp1: TComplex): TComplex;             //  V:=Sin(Z)
  Procedure CCosV  (var Cmp1: TComplex);                         //  Z:=Cos(Z)
  Function  CCos   (Cmp1: TComplex): TComplex;             //  V:=Cos(Z)
  Procedure CTanV  (var Cmp1: TComplex);                         //  Z:=Tan(Z)
  Function  CTan   (Cmp1: TComplex): TComplex;             //  V:=Tan(Z)
  Procedure CSinhV (var Cmp1: TComplex);                         //  Z:=Sinh(Z)
  Function  CSinh  (Cmp1: TComplex): TComplex;             //  V:=Sinh(Z)
  Procedure CCoshV (var Cmp1: TComplex);                         //  Z:=Cosh(Z)
  Function  CCosh  (Cmp1: TComplex): TComplex;             //  V:=Cosh(Z)
  Procedure CCosh2V(var Cmp1: TComplex);                         //  Z:=Cosh2(Z)
  Function  CCosh2 (Cmp1: TComplex): TComplex;             //  V:=Cosh2(Z)


  // ++20/08/2000
     Procedure CCotanV(var Cmp1: TComplex);                         //  Z:=Cotan(Z)
     Function  CCotan (Cmp1: TComplex): TComplex;             //  V:=Cotan(Z)
     Procedure CTanhV (var Cmp1: TComplex);                         //  Z:=Tanh(Z)
     Function  CTanh  (Cmp1: TComplex): TComplex;             //  V:=Tanh(Z)
     Procedure CCotanhV(var Cmp1: TComplex);                        //  Z:=Cotanh(Z)
     Function  CCotanh(Cmp1: TComplex): TComplex;             //  V:=Cotanh(Z)

  // ++19/08/2000
     Procedure CASinV (var Cmp1: TComplex);                         //  Z:=ArcSin(Z)
     Function  CASin  (Cmp1: TComplex): TComplex;             //  V:=ArcSin(Z)
     Procedure CACosV (var Cmp1: TComplex);                         //  Z:=ArcCos(Z)
     Function  CACos  (Cmp1: TComplex): TComplex;             //  V:=ArcCos(Z)
     Procedure CATanV (var Cmp1: TComplex);                         //  Z:=ArcTan(Z)
     Function  CATan  (Cmp1: TComplex): TComplex;             //  V:=ArcTan(Z)

  // ++19/08/2000
     Procedure CASinhV(var Cmp1: TComplex);                         //  Z:=ArcSinh(Z)
     Function  CASinh (Cmp1: TComplex): TComplex;             //  V:=ArcSinh(Z)
     Procedure CACoshV(var Cmp1: TComplex);                         //  Z:=ArcCosh(Z)
     Function  CACosh (Cmp1: TComplex): TComplex;             //  V:=ArcCosh(Z)
     Procedure CATanhV(var Cmp1: TComplex);                         //  Z:=ArcTanh(Z)
     Function  CATanh (Cmp1: TComplex): TComplex;             //  V:=ArcTanh(Z)

  // ++22/08/2000
     Procedure FuncDisp(Fn: Integer; var Cmp1: TComplex);     //  functions dispatcher
  // ++26.09.2004
     Function  FuncDispF(Fn: Integer; Cmp1: TComplex): TComplex;  //  another functions dispatcher

  // ++16/06/2001
     Procedure CSqrRecipV(var Cmp1: TComplex);       // #33
     Procedure CTripRecipV(var Cmp1: TComplex);      // #34
     Procedure CFourRecipV(var Cmp1: TComplex);      // #35
     Procedure CSin_ZZV(var Cmp1: TComplex);         // #36
     Procedure CCos_ZZV(var Cmp1: TComplex);         // #37


Implementation

Const
        SmallTol: Extended = 1E-25;
        BigTol  : Single   = 1E9;
        ExpTol  : Single   = 64;

Function  MakeComplex(real, imag: Extended): TComplex;
Begin
  Result.real:=real;
  Result.imag:=imag;
End;

Procedure SetResult(var x,y: Extended; complex: TComplex);
Begin
  x:=complex.real;
  y:=complex.imag;
End;


Function  CModule(Cmp1: TComplex): Extended;
Begin
  Result:=Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag;
End;


{ ------------------------------------------------------------------- }
Procedure CAddV(var Cmp1: TComplex; Cmp2: TComplex);
Begin
  Cmp1.real:=Cmp1.real+Cmp2.real;
  Cmp1.imag:=Cmp1.imag+Cmp2.imag;
End;

Function  CAdd (Cmp1, Cmp2: TComplex): TComplex;
Begin
  Result.real:=Cmp1.real+Cmp2.real;
  Result.imag:=Cmp1.imag+Cmp2.imag;
End;

Function  CAddR(Cmp1: TComplex; t: Extended): TComplex;  //  V:=Z+var
Begin
  Result.real:=Cmp1.real + t;
  Result.imag:=Cmp1.imag;
End;




Procedure CSubV(var Cmp1: TComplex; Cmp2: TComplex);
Begin
  Cmp1.real:=Cmp1.real-Cmp2.real;
  Cmp1.imag:=Cmp1.imag-Cmp2.imag;
End;

Function  CSub (Cmp1, Cmp2: TComplex): TComplex;
Begin
  Result.real:=Cmp1.real-Cmp2.real;
  Result.imag:=Cmp1.imag-Cmp2.imag;
End;

Function  CSubR(Cmp1: TComplex; t: Extended): TComplex;  //  V:=Z-var
Begin
  Result.real:=Cmp1.real-t;
  Result.imag:=Cmp1.imag;
End;




Procedure CMulV(var Cmp1: TComplex; Cmp2: TComplex);
  var tmp: Extended;
Begin
  tmp      :=Cmp1.real*Cmp2.real - Cmp1.imag*Cmp2.imag;
  Cmp1.imag:=Cmp1.real*Cmp2.imag + Cmp1.imag*Cmp2.real;
  Cmp1.real:=tmp;
End;

Function  CMul (Cmp1, Cmp2: TComplex): TComplex;
Begin
  Result.real:=Cmp1.real*Cmp2.real - Cmp1.imag*Cmp2.imag;
  Result.imag:=Cmp1.real*Cmp2.imag + Cmp1.imag*Cmp2.real;
End;

Function  CMulR(Cmp1: TComplex; t: Extended): TComplex;  //  V:=Z*var
Begin
  Result.real:=Cmp1.real*t;
  Result.imag:=Cmp1.imag*t;
End;




Procedure CDivV(var Cmp1: TComplex; Cmp2: TComplex);
  var tmp1,tmp2: Extended;
Begin
  tmp1       := Cmp2.real*Cmp2.real + Cmp2.imag*Cmp2.imag + SmallTol;
  tmp2       :=(Cmp1.real*Cmp2.real + Cmp1.imag*Cmp2.imag)/tmp1;
  Cmp1.imag  :=(Cmp1.imag*Cmp2.real - Cmp1.real*Cmp2.imag)/tmp1;
  Cmp1.real  := tmp2;
End;

Function  CDiv (Cmp1, Cmp2: TComplex): TComplex;
  var tmp: Extended;
Begin
  tmp        := Cmp2.real*Cmp2.real + Cmp2.imag*Cmp2.imag + SmallTol;
  Result.real:=(Cmp1.real*Cmp2.real + Cmp1.imag*Cmp2.imag)/tmp;
  Result.imag:=(Cmp1.imag*Cmp2.real - Cmp1.real*Cmp2.imag)/tmp;
End;

Function  CDivR(Cmp1: TComplex; t: Extended): TComplex;  //  V:=Z/var
Begin
  Result.real:=Cmp1.real/t;
  Result.imag:=Cmp1.imag/t;
End;



{ ------------------------------------------------------------------- }
Procedure CSqrV(var Cmp1: TComplex);
  var tmp: Extended;
Begin
  tmp      :=(Cmp1.real+Cmp1.imag)*(Cmp1.real-Cmp1.imag);
  Cmp1.imag:= Cmp1.real*Cmp1.imag * 2;
  Cmp1.real:=tmp;
End;

Function  CSqr (Cmp1: TComplex): TComplex;
Begin
  Result.real:=(Cmp1.real+Cmp1.imag)*(Cmp1.real-Cmp1.imag);
  Result.imag:= Cmp1.real*Cmp1.imag * 2;
End;




Procedure CTripleV(var Cmp1: TComplex);
  var tmp: Extended;
Begin
  tmp      :=Cmp1.real*(Cmp1.real*Cmp1.real - 3*Cmp1.imag*Cmp1.imag);
  Cmp1.imag:=Cmp1.imag*(3*Cmp1.real*Cmp1.real - Cmp1.imag*Cmp1.imag);
  Cmp1.real:=tmp;
End;

Function  CTriple (Cmp1: TComplex): TComplex;
Begin
  Result.real:=Cmp1.real*(Cmp1.real*Cmp1.real - 3*Cmp1.imag*Cmp1.imag);
  Result.imag:=Cmp1.imag*(3*Cmp1.real*Cmp1.real - Cmp1.imag*Cmp1.imag);
End;



Procedure CFourV(var Cmp1: TComplex);
  var tmpR, tmpI: Extended;
Begin
  tmpR:=(Cmp1.real-Cmp1.imag)*(Cmp1.real+Cmp1.imag);
  tmpI:= Cmp1.real*Cmp1.imag * 2;

  Cmp1.real:=(tmpR-tmpI)*(tmpR+tmpI);
  Cmp1.imag:= tmpR*tmpI * 2;
End;

Function  CFour (Cmp1: TComplex): TComplex;
  var tmpR, tmpI: Extended;
Begin
  tmpR:=(Cmp1.real-Cmp1.imag)*(Cmp1.real+Cmp1.imag);
  tmpI:= Cmp1.real*Cmp1.imag * 2;

  Result.real:=(tmpR-tmpI)*(tmpR+tmpI);
  Result.imag:= tmpR*tmpI * 2;
End;



{ ------------------------------------------------------------------- }
Procedure CFlipV(var Cmp1: TComplex);
  var tmp: Extended;
Begin
  tmp:=Cmp1.real;
  Cmp1.real:=Cmp1.imag;
  Cmp1.imag:=tmp;
End;

Function  CFlip (Cmp1: TComplex): TComplex;
Begin
  Result.real:=Cmp1.imag;
  Result.imag:=Cmp1.real;
End;




Procedure CRevV(var Cmp1: TComplex);
  var tmp: Extended;
Begin
  tmp       := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Cmp1.real := Cmp1.real/tmp;
  Cmp1.imag :=-Cmp1.imag/tmp;
End;

Function  CRev (Cmp1: TComplex): TComplex;
  var tmp: Extended;
Begin
  tmp        := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Result.real:= Cmp1.real/tmp;
  Result.imag:=-Cmp1.imag/tmp;
End;




Procedure CRev2V(var Cmp1: TComplex; Cmp2: TComplex);
  var tmp: Extended;
Begin
  tmp       := (Cmp1.real-Cmp2.real)*(Cmp1.real-Cmp2.real) +
               (Cmp1.imag-Cmp2.imag)*(Cmp1.imag-Cmp2.imag) + SmallTol;
  Cmp1.real := (Cmp1.real-Cmp2.real)/tmp;
  Cmp1.imag := (Cmp2.imag-Cmp1.imag)/tmp;
End;

Function  CRev2 (Cmp1,Cmp2: TComplex): TComplex;
  var tmp: Extended;
Begin
  tmp        := (Cmp1.real-Cmp2.real)*(Cmp1.real-Cmp2.real) +
                (Cmp1.imag-Cmp2.imag)*(Cmp1.imag-Cmp2.imag) + SmallTol;
  Result.real:= (Cmp1.real-Cmp2.real)/tmp;
  Result.imag:= (Cmp2.imag-Cmp1.imag)/tmp;
End;


Function  CReal (Cmp1: TComplex): TComplex;              //  CReal ++10/12/2001
Begin                                                          //  suggested by Pete
  Result.real:=Cmp1.real;
  Result.imag:=0;
End;

Procedure CRealV(var Cmp1: TComplex);
Begin
  Cmp1.imag:=0;
End;

Function  CImag (Cmp1: TComplex): TComplex;              //  CImag ++10/12/2001
Begin                                                          //  suggested by Pete
  Result.real:=Cmp1.imag;
  Result.imag:=0;
End;

Procedure CImagV(var Cmp1: TComplex);
Begin
  Cmp1.real:=Cmp1.imag;
  Cmp1.imag:=0;
End;


Function  CAbs  (Cmp1: TComplex): TComplex;              //  ++19/12/2001
Begin
  Result.real:=Abs(Cmp1.real);
  Result.imag:=Cmp1.imag;
End;

Procedure CAbsV (var Cmp1: TComplex);
Begin
  Cmp1.real:=Abs(Cmp1.real);
End;

Function  CAbs2 (Cmp1: TComplex): TComplex;              //  ++19/12/2001
Begin
  Result.real:=Abs(Cmp1.real);
  Result.imag:=Abs(Cmp1.imag);
End;

Procedure CAbs2V(var Cmp1: TComplex);
Begin
  Cmp1.real:=Abs(Cmp1.real);
  Cmp1.imag:=Abs(Cmp1.imag);
End;

{ ------------------------------------------------------------------- }
Procedure CSqrtV(var Cmp1: TComplex);                          //  Z:=Sqrt(Z)
  var a,b  : Extended;
      ck,sk: Extended;
Begin
  With Cmp1 do begin
    a:=Sqrt(Sqrt(real*real+imag*imag));
    b:=ArcTan2(imag,real)/2;
    asm
      FLD   b
      FSINCOS
      FNSTSW  AX
      SAHF
      JPO     @@a1
      FSTP    ST(0)
      FLD1    // FLDZ
      FLD1    // FLDZ
    @@a1:
      FSTP  ck
      FSTP  sk
      FWAIT
    end;
    real:=a*ck;
    imag:=a*sk;
  end;
End;

Function  CSqrt (Cmp1: TComplex): TComplex;              //  V:=Sqrt(Z)
  var a,b: Extended;
      ck,sk: Extended;
Begin
  With Cmp1 do begin
    a:=Sqrt(Sqrt(real*real+imag*imag));
    b:=ArcTan2(imag,real)/2;
  end;
  asm
    FLD   b
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  ck
    FSTP  sk
    FWAIT
  end;
  Result.real:=a*ck;
  Result.imag:=a*sk;
End;



Procedure CExpV (var Cmp1: TComplex);                          //  Z:=Exp(Z)
  var tmp: Extended;
      ck,sk: Extended;
Begin
  tmp:=Cmp1.imag;
  asm
    FLD     tmp
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  ck
    FSTP  sk
    FWAIT
  end;
  If  Cmp1.real> ExpTol Then Cmp1.real:= ExpTol;
  If  Cmp1.real<-ExpTol Then Cmp1.real:=-ExpTol;
  tmp :=Exp(Cmp1.real);
  Cmp1.real:=tmp*ck;
  Cmp1.imag:=tmp*sk;
End;

Function  CExp  (Cmp1: TComplex): TComplex;              //  V:=Exp(Z)
  var tmp: Extended;
      ck,sk: Extended;
Begin
  tmp:=Cmp1.imag;
  asm
    FLD     tmp
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  ck
    FSTP  sk
    FWAIT
  end;
  tmp:=Cmp1.real;
  If  tmp> ExpTol Then tmp:= ExpTol;
  If  tmp<-ExpTol Then tmp:=-ExpTol;
  tmp :=Exp(tmp);
  Result.real:=tmp*ck;
  Result.imag:=tmp*sk;
End;




Procedure CLnV  (var Cmp1: TComplex);                          //  Z:=Ln(Z)
  var tmp: Extended;
Begin
  tmp :=Log2(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol)/2.7182818285;
  Cmp1.imag:=ArcTan2(Cmp1.imag, Cmp1.real);
  Cmp1.real:=tmp;
End;

Function  CLn   (Cmp1: TComplex): TComplex;              //  V:=Ln(Z)
Begin
  Result.real:=Log2(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol)/2.7182818285;
  Result.imag:=ArcTan2(Cmp1.imag, Cmp1.real);
End;




Procedure CPowerV(var Cmp1: TComplex; Cmp2: TComplex);         //  Z:=Z^A
  var h1x,h1y: Extended;
      h2x,h2y: Extended;
      f      : Extended;
      sXk,cXk: Extended;
Begin
  h1x:=Log2(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol)/2.7182818285;
  h1y:=ArcTan2(Cmp1.imag, Cmp1.real);

  h2x:=h1x*Cmp2.real-h1y*Cmp2.imag;
  h2y:=h1y*Cmp2.real+h1x*Cmp2.imag;
  asm
    FLD     h2y
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;
  If h2x> ExpTol Then h2x:= ExpTol;
  If h2x<-ExpTol Then h2x:=-ExpTol;
  f  :=Exp(h2x);
  Cmp1.real :=f*cXk;
  Cmp1.imag :=f*sXk;
End;

Function  CPower (Cmp1,Cmp2: TComplex): TComplex;        //  V:=Z^A
  var h1x,h1y: Extended;
      h2x,h2y: Extended;
      f      : Extended;
      sXk,cXk: Extended;
Begin
  h1x:=Log2(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol)/2.7182818285;
  h1y:=ArcTan2(Cmp1.imag, Cmp1.real);

  h2x:=h1x*Cmp2.real-h1y*Cmp2.imag;
  h2y:=h1y*Cmp2.real+h1x*Cmp2.imag;
  asm
    FLD     h2y
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;

  If h2x> ExpTol Then h2x:= ExpTol;
  If h2x<-ExpTol Then h2x:=-ExpTol;
  f  :=Exp(h2x);
  Result.real :=f*cXk;
  Result.imag :=f*sXk;
End;


Procedure CPowerV2(var Cmp1: TComplex; Cmp2: TComplex);         //  Z:=Z^A
  var h1x,h1y: Extended;
      h2x,h2y: Extended;
      f      : Extended;
      sXk,cXk: Extended;
Begin
  h1x:=0.5*Ln (Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol);
  h1y:=ArcTan2(Cmp1.imag, Cmp1.real);

  h2x:=h1x*Cmp2.real-h1y*Cmp2.imag;
  h2y:=h1y*Cmp2.real+h1x*Cmp2.imag;
  asm
    FLD     h2y
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;
  If h2x> ExpTol Then h2x:= ExpTol;
  If h2x<-ExpTol Then h2x:=-ExpTol;
  f  :=Exp(h2x);
  Cmp1.real :=f*cXk;
  Cmp1.imag :=f*sXk;
End;

Function  CPower2(Cmp1,Cmp2: TComplex): TComplex;        //  V:=Z^A
  var h1x,h1y: Extended;
      h2x,h2y: Extended;
      f      : Extended;
      sXk,cXk: Extended;
Begin
  h1x:=0.5*Ln (Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol);
  h1y:=ArcTan2(Cmp1.imag, Cmp1.real);

  h2x:=h1x*Cmp2.real-h1y*Cmp2.imag;
  h2y:=h1y*Cmp2.real+h1x*Cmp2.imag;
  asm
    FLD     h2y
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;

  If h2x> ExpTol Then h2x:= ExpTol;
  If h2x<-ExpTol Then h2x:=-ExpTol;
  f  :=Exp(h2x);
  Result.real :=f*cXk;
  Result.imag :=f*sXk;
End;




Procedure CPowerRV(var Cmp1: TComplex; t: Extended);           //  Z:=Z^var
  var tr,ti,f: Extended;
      sXk,cXk: Extended;
Begin
  tr:=t*Log2(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol)/2.7182818285;
  ti:=t*ArcTan2(Cmp1.imag,Cmp1.real);
  asm
    FLD     ti
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;
  If tr> ExpTol Then tr:= ExpTol;
  If tr<-ExpTol Then tr:=-ExpTol;
  f  :=Exp(tr);
  Cmp1.real:=f*cXk;
  Cmp1.imag:=f*sXk;
End;

Function  CPowerR (Cmp1:TComplex;t: Extended): TComplex; //  Result:=Z^var
  var tr,ti,f: Extended;
      sXk,cXk: Extended;
Begin
  tr:=t*Log2(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol)/2.7182818285;
  ti:=t*ArcTan2(Cmp1.imag,Cmp1.real);
  asm
    FLD     ti
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;
  If tr> ExpTol Then tr:= ExpTol;
  If tr<-ExpTol Then tr:=-ExpTol;
  f  :=Exp(tr);
  Result.real:=f*cXk;
  Result.imag:=f*sXk;
End;


Procedure CPowerRV2(var Cmp1: TComplex; t: Extended);           //  Z:=Z^var
  var tr,ti,f: Extended;
      sXk,cXk: Extended;
Begin
  tr:=0.5*t*Ln(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol);
  ti:=t*ArcTan2(Cmp1.imag,Cmp1.real);
  asm
    FLD     ti
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;
  If tr> ExpTol Then tr:= ExpTol;
  If tr<-ExpTol Then tr:=-ExpTol;
  f  :=Exp(tr);
  Cmp1.real:=f*cXk;
  Cmp1.imag:=f*sXk;
End;

Function  CPowerR2(Cmp1:TComplex;t: Extended): TComplex; //  Result:=Z^var
  var tr,ti,f: Extended;
      sXk,cXk: Extended;
Begin
  tr:=0.5*t*Ln(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol);
  ti:=t*ArcTan2(Cmp1.imag,Cmp1.real);
  asm
    FLD     ti
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;
  If tr> ExpTol Then tr:= ExpTol;
  If tr<-ExpTol Then tr:=-ExpTol;
  f  :=Exp(tr);
  Result.real:=f*cXk;
  Result.imag:=f*sXk;
End;




{ ------------------------------------------------------------------- }
Procedure CSinV (var Cmp1: TComplex);                          //  Z:=Sin(Z)
  var a1,a2  : Extended;
      sXk,cXk: Extended;
Begin
  a1  := Cmp1.real;
  asm
    FLD     a1
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;
  a1:=Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;
  Cmp1.real:=sXk*(a1+a2);
  Cmp1.imag:=cXk*(a1-a2);
End;

Function  CSin  (Cmp1: TComplex): TComplex;              //  Result:=Sin(Z)
  var a1,a2  : Extended;
      sXk,cXk: Extended;
Begin
  a1  := Cmp1.real;
  asm
    FLD     a1
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP    cXk
    FSTP    sXk
    FWAIT
  end;
  a1:=Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;
  Result.real:=sXk*(a1+a2);
  Result.imag:=cXk*(a1-a2); 
End;


Procedure CCosV (var Cmp1: TComplex);                          //  Z:=Cos(Z)
  var a1,a2  : Extended;
      sXk,cXk: Extended;
Begin
//  cos(x)cosh(y) - i sin(x)sinh(y)
  a1  := Cmp1.real;
  asm
    FLD   a1
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;
  Cmp1.real:= cXk*(a1+a2);
  Cmp1.imag:=-sXk*(a1-a2);
End;

Function  CCos  (Cmp1: TComplex): TComplex;              //  Result:=Cos(Z)
  var a1,a2  : Extended;
      sXk,cXk: Extended;
Begin
//  cos(x)cosh(y) - i sin(x)sinh(y)
  a1  := Cmp1.real;
  asm
    FLD   a1
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag;
  If a1> ExpTol Then a1:= ExpTol;
  If a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;
  Result.real:= cXk*(a1+a2);
  Result.imag:=-sXk*(a1-a2);
End;

Procedure CCosxxV(var Cmp1: TComplex);
  var a1 ,a2 : Extended;
      sXk,cXk: Extended;
Begin
//  cos(x)cosh(y) + i sin(x)sinh(y)
  a1:=Cmp1.real;
  asm
    FLD     a1
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@m1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@m1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag;
  If a1> ExpTol Then a1:= ExpTol;
  If a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2  :=0.25/a1;

  Cmp1.real:=cXk*(a1+a2);
  Cmp1.imag:=sXk*(a1-a2);
End;

Function  CCosxx(Cmp1: TComplex): TComplex;
  var a1 ,a2 : Extended;
      sXk,cXk: Extended;
Begin
//  cos(x)cosh(y) + i sin(x)sinh(y)
  a1:=Cmp1.real;
  asm
    FLD     a1
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@m1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@m1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag;
  If a1> ExpTol Then a1:= ExpTol;
  If a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2  :=0.25/a1;

  Result.real:=cXk*(a1+a2);
  Result.imag:=sXk*(a1-a2);
End;






Procedure CTanV (var Cmp1: TComplex);                          //  Z:=Tan(Z)
  var b,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
//        sin(2x)               sinh(2y)
//  ------------------  + i------------------
//  cos(2x) + cosh(2y)     cos(2x) + cosh(2y)
  a1:=Cmp1.real;
  asm
    FLD     a1
    FADD    ST,ST(0)
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag+Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;

  b:=cXk+(a1+a2);
  If b=0 Then b:=SmallTol;
  Cmp1.real:=sXk/b;
  Cmp1.imag:=(a1-a2)/b;
End;

Function  CTan  (Cmp1: TComplex): TComplex;              //  Result:=Tan(Z)
  var b,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
  a1:=Cmp1.real;
  asm
    FLD     a1
    FADD    ST,ST(0)
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag+Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;

  b:=cXk+(a1+a2);
  If b=0 Then b:=SmallTol;
  Result.real:=sXk/b;
  Result.imag:=(a1-a2)/b;
End;




Procedure CCotanV(var Cmp1: TComplex);                         //  Z:=Cotan(Z)
  var b,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
//  sin(2x) - i*sinh(2y)
//  --------------------
//   cosh(2y) - cos(2x)
  a1:=Cmp1.real;
  asm
    FLD     a1
    FADD    ST,ST(0)
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag+Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;

  b:=(a1+a2)-cXk;
  If b=0 Then b:=SmallTol;
  Cmp1.real:= sXk/b;
  Cmp1.imag:=(a2-a1)/b;
End;

Function  CCotan (Cmp1: TComplex): TComplex;             //  V:=Cotan(Z)
  var b,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
  a1:=Cmp1.real;
  asm
    FLD     a1
    FADD    ST,ST(0)
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag+Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;

  b:=(a1+a2)-cXk;
  If b=0 Then b:=SmallTol;
  Result.real:= sXk/b;
  Result.imag:=(a2-a1)/b;
End;




Procedure CSinhV(var Cmp1: TComplex);                          //  Z:=Sinh(Z)
  var a,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
//  Sinh(real)*Cos(imag) + j*Cosh(real)*Sin(imag)
  a  :=Cmp1.imag;
  asm
    FLD     a
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.real;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1  :=Exp(a1)/2;
  a2  :=0.25/a1;

  Cmp1.real:=(a1-a2)*cXk;
  Cmp1.imag:=(a1+a2)*sXk;
End;

Function  CSinh (Cmp1: TComplex): TComplex;              //  Result:=Sinh(Z)
  var a,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
//  Sinh(real)*Cos(imag) + j*Cosh(real)*Sin(imag)
  a  :=Cmp1.imag;
  asm
    FLD     a
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk  { Cos(imag) }
    FSTP  sXk  { Sin(imag) }
    FWAIT
  end;
  a1:=Cmp1.real;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1  :=Exp(a1)/2;
  a2  :=0.25/a1;

  Result.real:=(a1-a2)*cXk;
  Result.imag:=(a1+a2)*sXk;
End;


Procedure CCoshV(var Cmp1: TComplex);                          //  Z:=Cosh(Z)
  var a,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
//  Cosh(real)*Cos(imag) + j*Sinh(real)*Sin(imag);
  a  :=Cmp1.imag;
  asm
    FLD     a
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLD1    // FLDZ
    FLD1    // FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  If Cmp1.real> ExpTol Then Cmp1.real:= ExpTol;
  If Cmp1.real<-ExpTol Then Cmp1.real:=-ExpTol;
  a1  :=Exp(Cmp1.real)/2;
  If a1=0 Then a1:=SmallTol;
  a2  :=0.25/a1;

  Cmp1.real:=(a1+a2)*cXk;
  Cmp1.imag:=(a1-a2)*sXk;
End;


Function  CCosh (Cmp1: TComplex): TComplex;              //  Result:=Cosh(Z)
  var a,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
//Cosh(real)*Cos(imag) + j*Sinh(real)*Sin(imag);
  a  :=Cmp1.imag;
  asm
    FLD     a
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLD1    // FLDZ
    FLD1    // FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a:=Cmp1.real;
  If a> ExpTol Then a:= ExpTol;
  If a<-ExpTol Then a:=-ExpTol;
  a1 :=Exp(a)/2;
  a2 :=0.25/a1;

  Result.real:=(a1+a2)*cXk;
  Result.imag:=(a1-a2)*sXk;
End;


// ++17.06.2002
Procedure CCosh2V(var Cmp1: TComplex);
  var a,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
//  cosh(real)*cos(imag) + i cosh(real)sin(imag)
  a  :=Cmp1.imag;
  asm
    FLD     a
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a_1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a_1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a:=Cmp1.real;
  If a> ExpTol Then a:= ExpTol;
  If a<-ExpTol Then a:=-ExpTol;
  a1  :=Exp(a)/2;
  a2  :=0.25/a1;

  Cmp1.real:=(a1+a2)*cXk;
  Cmp1.imag:=(a1+a2)*sXk;
End;

Function CCosh2(Cmp1: TComplex): TComplex;
  var a,a1,a2: Extended;
      sXk,cXk: Extended;
Begin
//  cosh(real)*cos(imag) + i cosh(real)sin(imag)
  a  :=Cmp1.imag;
  asm
    FLD     a
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a_1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a_1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a:=Cmp1.real;
  If a> ExpTol Then a:= ExpTol;
  If a<-ExpTol Then a:=-ExpTol;
  a1  :=Exp(a)/2;
  a2  :=0.25/a1;

  Result.real:=(a1+a2)*cXk;
  Result.imag:=(a1+a2)*sXk;
End;



Procedure CTanhV (var Cmp1: TComplex);                         //  Z:=Tanh(Z)
  var b: Extended;
Begin
  b:=Cosh(2*Cmp1.real)+Cos(2*Cmp1.imag);
  If b=0 Then b:=SmallTol;
  Cmp1.real:=Sinh(2*Cmp1.real)/b;
  Cmp1.imag:=Sin( 2*Cmp1.imag)/b;
End;

Function  CTanh  (Cmp1: TComplex): TComplex;             //  Result:=Tanh(Z)
  var b: Extended;
Begin
  b:=Cosh(2*Cmp1.real)+Cos(2*Cmp1.imag);
  If b=0 Then b:=SmallTol;
  Result.real:=Sinh(2*Cmp1.real)/b;
  Result.imag:=Sin( 2*Cmp1.imag)/b;
End;




Procedure CCotanhV(var Cmp1: TComplex);                        //  Z:=Cotanh(Z)
  var b: Extended;
Begin
  b:=Cosh(2*Cmp1.real)-Cos(2*Cmp1.imag); If b=0 Then b:=SmallTol;
  Cmp1.real:= Sinh(2*Cmp1.real)/b;
  Cmp1.imag:=-Sin( 2*Cmp1.imag)/b;
End;

Function  CCotanh(Cmp1: TComplex): TComplex;             //  V:=Cotanh(Z)
  var b: Extended;
Begin
  b:=Cosh(2*Cmp1.real)-Cos(2*Cmp1.imag); If b=0 Then b:=SmallTol;
  Result.real:= Sinh(2*Cmp1.real)/b;
  Result.imag:=-Sin( 2*Cmp1.imag)/b;
End;

{ ------------------------------------------------------------------- }
Procedure CASinV (var Cmp1: TComplex);                         //  Z:=Sin(Z)
  var a,b  : Extended;
      xr,xi: Extended;    //  [-i * log(i*z+sqrt(1-z*z))]
Begin
  With Cmp1 do begin
    a :=(real - imag)*(real + imag);
    xi:= -2*real*imag;
    xr:=1-a;

    a :=Sqrt(Sqrt(xr*xr+xi*xi));
    b :=ArcTan2(xi,xr)/2;
    xr:=a*cos(b)-imag;
    xi:=a*sin(b)+real;

    imag:=-Ln(xr*xr + xi*xi + SmallTol)/2;
    real:= ArcTan2(xi, xr);
 end;
End;

Function  CASin  (Cmp1: TComplex): TComplex;             //  Result:=ArcSin(Z)
  var a,b  : Extended;
      xr,xi: Extended;    //  [-i * log(i*z+sqrt(1-z*z))]
Begin
  With Cmp1 do begin
    a :=(real - imag)*(real + imag);
    xi:= -2*real*imag;
    xr:=1-a;

    a :=Sqrt(Sqrt(xr*xr+xi*xi));
    b :=ArcTan2(xi,xr)/2;
    xr:=a*cos(b)-imag;
    xi:=a*sin(b)+real;

    Result.imag:=-Ln(xr*xr + xi*xi + SmallTol)/2;
    Result.real:=ArcTan2(xi, xr);
 end;
End;




Procedure CACosV (var Cmp1: TComplex);                         //  Z:=Cos(Z)
  var a,b  : Extended;
      xr,xi: Extended;    //  [-i * log(z+sqrt(z*z-1))]
Begin
  With Cmp1 do begin
    a :=(real - imag)*(real + imag);
    xi:=2*real*imag;
    xr:=a-1;

    a :=Sqrt(Sqrt(xr*xr+xi*xi));
    b :=ArcTan2(xi,xr)/2;
    xr:=a*cos(b)+real;
    xi:=a*sin(b)+imag;

    imag:=-Ln(xr*xr + xi*xi + SmallTol)/2;
    real:=ArcTan2(xi, xr);
  end;
End;

Function  CACos  (Cmp1: TComplex): TComplex;             //  Result:=ArcCos(Z)
  var a,b  : Extended;
      xr,xi: Extended;    //  [-i * log(z+sqrt(z*z-1))]
Begin
  With Cmp1 do begin
    a :=(real - imag)*(real + imag);
    xi:=2*real*imag;
    xr:=a-1;

    a :=Sqrt(Sqrt(xr*xr+xi*xi));
    b :=ArcTan2(xi,xr)/2;
    xr:=a*cos(b)+real;
    xi:=a*sin(b)+imag;

    Result.imag:=-Ln(xr*xr + xi*xi + SmallTol)/2;
    Result.real:=ArcTan2(xi, xr);
  end;
End;




Procedure CATanV (var Cmp1: TComplex);                         //  Z:=ArcTan(Z)
  var a,b  : Extended;
      xr,xi: Extended;
      tr,ti: Extended;
      f    : Extended;
Begin
//  i/2 * log((1-i*z)/(1+i*z))
  a :=1+Cmp1.imag; b := -Cmp1.real;
  xr:=1-Cmp1.imag; xi:=  Cmp1.real;

  f := xr*xr + xi*xi + SmallTol;
  tr:=(a*xr + b*xi)/f;
  ti:=(b*xr - a*xi)/f;

  xr:=Ln(tr*tr + ti*ti + SmallTol)/4;
  xi:=ArcTan2(ti, tr)/2;
  Cmp1.real:=-xi;
  Cmp1.imag:= xr;
End;

Function  CATan  (Cmp1: TComplex): TComplex;             //  V:=ArcTan(Z)
  var a,b  : Extended;
      xr,xi: Extended;
      tr,ti: Extended;
      f    : Extended;
Begin
//  i/2 * log((1-i*z)/(1+i*z))
  a :=1+Cmp1.imag; b := -Cmp1.real;
  xr:=1-Cmp1.imag; xi:=  Cmp1.real;

  f := xr*xr + xi*xi + SmallTol;
  tr:=(a*xr + b*xi)/f;
  ti:=(b*xr - a*xi)/f;

  xr:=Ln(tr*tr + ti*ti + SmallTol)/4;
  xi:=ArcTan2(ti, tr)/2;
  Result.real:=-xi;
  Result.imag:= xr;
End;



{ ------------------------------------------------------------------- }
Procedure CASinhV(var Cmp1: TComplex);                         //  Z:=ArcSinh(Z)
  var a,b  : Extended;
      xr,xi: Extended;
Begin
//  log(z+sqrt(z*z+1))
  a :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  xi:=2*Cmp1.real*Cmp1.imag;
  xr:=a+1;

  a :=Sqrt(Sqrt(xr*xr+xi*xi));
  b :=ArcTan2(xi,xr)/2;
  xr:=a*cos(b)+Cmp1.real;
  xi:=a*sin(b)+Cmp1.imag;

  Cmp1.real:=Ln(xr*xr + xi*xi + SmallTol)/2;
  Cmp1.imag:=ArcTan2(xi, xr);
End;

Function  CASinh (Cmp1: TComplex): TComplex;             //  V:=ArcSinh(Z)
  var a,b  : Extended;      //  log(z+sqrt(z*z+1))
      xr,xi: Extended;
Begin
  a :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  xi:=2*Cmp1.real*Cmp1.imag;
  xr:=a+1;

  a :=Sqrt(Sqrt(xr*xr+xi*xi));
  b :=ArcTan2(xi,xr)/2;
  xr:=a*cos(b)+Cmp1.real;
  xi:=a*sin(b)+Cmp1.imag;

  Result.real:=Ln(xr*xr + xi*xi + SmallTol)/2;
  Result.imag:=ArcTan2(xi, xr);
End;




Procedure CACoshV(var Cmp1: TComplex);                         //  Z:=ArcCosh(Z)
  var a,b  : Extended;    //  log(z+sqrt(z*z-1))
      xr,xi: Extended;
Begin
  a :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  xi:=2*Cmp1.real*Cmp1.imag;
  xr:=a-1;

  a :=Sqrt(Sqrt(xr*xr+xi*xi));
  b :=ArcTan2(xi,xr)/2;
  xr:=a*cos(b)+Cmp1.real;
  xi:=a*sin(b)+Cmp1.imag;

  Cmp1.real:=Ln(xr*xr + xi*xi + SmallTol)/2;
  Cmp1.imag:=ArcTan2(xi, xr);
End;

Function  CACosh (Cmp1: TComplex): TComplex;             //  V:=ArcCosh(Z)
  var a,b  : Extended;    //  log(z+sqrt(z*z-1))
      xr,xi: Extended;
Begin
  a :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  xi:=2*Cmp1.real*Cmp1.imag;
  xr:=a-1;

  a :=Sqrt(Sqrt(xr*xr+xi*xi));
  b :=ArcTan2(xi,xr)/2;
  xr:=a*cos(b)+Cmp1.real;
  xi:=a*sin(b)+Cmp1.imag;

  Result.real:=Ln(xr*xr + xi*xi + SmallTol)/2;
  Result.imag:=ArcTan2(xi, xr);
End;




Procedure CATanhV(var Cmp1: TComplex);                         //  Z:=ArcTanh(Z)
  var a,b  : Extended;    //  log((1+z)/(1-z))/2
      xr,xi: Extended;
      tr,ti: Extended;
      f    : Extended;
Begin
  a :=1+Cmp1.real; b := Cmp1.imag;
  xr:=1-Cmp1.real; xi:=-Cmp1.imag;

  f:= xr*xr + xi*xi + SmallTol;
  tr:=(a*xr + b*xi)/f;
  ti:=(b*xr - a*xi)/f;

  Cmp1.real:=Ln(tr*tr + ti*ti + SmallTol)/4;
  Cmp1.imag:=ArcTan2(ti, tr)/2;
End;

Function  CATanh (Cmp1: TComplex): TComplex;             //  V:=ArcTanh(Z)
  var a,b  : Extended;    //  log((1+z)/(1-z))/2
      xr,xi: Extended;
      tr,ti: Extended;
      f    : Extended;
Begin
  a :=1+Cmp1.real; b := Cmp1.imag;
  xr:=1-Cmp1.real; xi:=-Cmp1.imag;

  f:= xr*xr + xi*xi + SmallTol;
  tr:=(a*xr + b*xi)/f;
  ti:=(b*xr - a*xi)/f;

  Result.real:=Ln(tr*tr + ti*ti + SmallTol)/4;
  Result.imag:=ArcTan2(ti, tr)/2;
End;




{ =================================================================== }
Procedure CCabsV(var Cmp1: TComplex);
Begin
  Cmp1.real:=Sqrt(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag);
  Cmp1.imag:=0;
End;

Function  CCabs(Cmp1: TComplex): TComplex;
Begin
  Result.real:=Sqrt(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag);
  Result.imag:=0;
End;


Procedure CFloorV(var Cmp1: TComplex);
Begin
  Cmp1.real:=Floor(Cmp1.real);
  Cmp1.imag:=Floor(Cmp1.imag);
End;

Function  CFloor(Cmp1: TComplex): TComplex;
Begin
  Result.real:=Floor(Cmp1.real);
  Result.imag:=Floor(Cmp1.imag);
End;


Procedure CCeilV(var Cmp1: TComplex);
Begin
  Cmp1.real:=Trunc(Cmp1.real);
  Cmp1.imag:=Trunc(Cmp1.imag);
End;

Function  CCeil(Cmp1: TComplex): TComplex;
Begin
  Result.real:=Trunc(Cmp1.real);
  Result.imag:=Trunc(Cmp1.imag);
End;


Procedure CTruncV(var Cmp1: TComplex);
Begin
  Cmp1.real:=Trunc(Cmp1.real);
  Cmp1.imag:=Trunc(Cmp1.imag);
End;

Function  CTrunc(Cmp1: TComplex): TComplex;
Begin
  Result.real:=Trunc(Cmp1.real);
  Result.imag:=Trunc(Cmp1.imag);
End;


Procedure CRoundV(var Cmp1: TComplex);
Begin
  Cmp1.real:=Round(Cmp1.real);
  Cmp1.imag:=Round(Cmp1.imag);
End;

Function  CRound(Cmp1: TComplex): TComplex;
Begin
  Result.real:=Round(Cmp1.real);
  Result.imag:=Round(Cmp1.imag);
End;



Procedure COneV(var Cmp1: TComplex);
Begin
  Cmp1.real:=1;
  Cmp1.imag:=0;
End;

Function  COne(Cmp1: TComplex): TComplex;
Begin
  Result.real:=1;
  Result.imag:=0;
End;




Procedure CRecipV(var Cmp1: TComplex);
  var a: Extended;
Begin
//  (x-iy) / (x^2+y^2)
  With Cmp1 do begin
    a   := real*real + imag*imag + SmallTol;
    real:= real/a;
    imag:=-imag/a
  end;
End;

Function  CRecip(Cmp1: TComplex): TComplex;
  var a: Extended;
Begin
  a   := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Result.real:= Cmp1.real/a;
  Result.imag:=-Cmp1.imag/a
End;



Function  CConj(Cmp1: TComplex): TComplex;
Begin
  Result.real:= Cmp1.real;
  Result.imag:=-Cmp1.imag;
End;

Procedure CConjV(var Cmp1: TComplex);
Begin
  Cmp1.imag:=-Cmp1.imag;
End;




Procedure CZeroV(var Cmp1: TComplex);
Begin
  Cmp1.real:=0;
  Cmp1.imag:=0;
End;

Function  CZero(Cmp1: TComplex): TComplex;
Begin
  Result.real:=0;
  Result.imag:=0;
End;



Procedure CLogV(var Cmp1: TComplex);
  var a: Extended;
Begin
//  (1/2)ln(x^2 + y^2) + i*arctan2(y/x)
  a   :=Ln(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol)/2;
  Cmp1.imag:=ArcTan2(Cmp1.imag, Cmp1.real);
  Cmp1.real:=a;
End;

Function  CLog(Cmp1: TComplex): TComplex;
  var a: Extended;
Begin
//  (1/2)ln(x^2 + y^2) + i*arctan2(y/x)
  a   :=Ln(Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol)/2;
  Result.imag:=ArcTan2(Cmp1.imag, Cmp1.real);
  Result.real:=a;
End;


{ ------------------------------------------------------------------- }
Procedure CSqrRecipV(var Cmp1: TComplex);       // #33
  var a: Extended;
Begin
  a        :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  Cmp1.imag:= 2*Cmp1.real*Cmp1.imag;
  Cmp1.real:= a;
  a        := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Cmp1.real:= Cmp1.real/a;
  Cmp1.imag:=-Cmp1.imag/a
End;

Function  CSqrRecip(Cmp1: TComplex): TComplex;
  var a: Extended;
Begin
  a        :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  Cmp1.imag:= 2*Cmp1.real*Cmp1.imag;
  Cmp1.real:= a;
  a        := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Result.real:= Cmp1.real/a;
  Result.imag:=-Cmp1.imag/a
End;



Procedure CTripRecipV(var Cmp1: TComplex);      // #34
  var a: Extended;
Begin
  a        := Cmp1.real*(Cmp1.real*Cmp1.real   - Cmp1.imag*Cmp1.imag*3);
  Cmp1.imag:= Cmp1.imag*(Cmp1.real*Cmp1.real*3 - Cmp1.imag*Cmp1.imag  );
  Cmp1.real:= a;
  a        := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Cmp1.real:= Cmp1.real/a;
  Cmp1.imag:=-Cmp1.imag/a
End;

Function  CTripRecip(Cmp1: TComplex): TComplex;      // #34
  var a: Extended;
Begin
  a        := Cmp1.real*(Cmp1.real*Cmp1.real   - Cmp1.imag*Cmp1.imag*3);
  Cmp1.imag:= Cmp1.imag*(Cmp1.real*Cmp1.real*3 - Cmp1.imag*Cmp1.imag  );
  Cmp1.real:= a;
  a        := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Result.real:= Cmp1.real/a;
  Result.imag:=-Cmp1.imag/a
End;



Procedure CFourRecipV(var Cmp1: TComplex);      // #35
  var a: Extended;
Begin
  a        :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  Cmp1.imag:= 2*Cmp1.real*Cmp1.imag;
  Cmp1.real:= a;
  Cmp1.real:=(a - Cmp1.imag)*(a + Cmp1.imag);
  Cmp1.imag:= 2*a*Cmp1.imag;
  a        := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Cmp1.real:= Cmp1.real/a;
  Cmp1.imag:=-Cmp1.imag/a
End;

Function  CFourRecip(Cmp1: TComplex): TComplex;    
  var a: Extended;
Begin
  a        :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  Cmp1.imag:= 2*Cmp1.real*Cmp1.imag;
  Cmp1.real:= a;
  Cmp1.real:=(a - Cmp1.imag)*(a + Cmp1.imag);
  Cmp1.imag:= 2*a*Cmp1.imag;
  a        := Cmp1.real*Cmp1.real + Cmp1.imag*Cmp1.imag + SmallTol;
  Result.real:= Cmp1.real/a;
  Result.imag:=-Cmp1.imag/a
End;



Procedure CSin_ZZV(var Cmp1: TComplex);         // #36
  var a1 ,a2 : Extended;
      sXk,cXk: Extended;
Begin
  a1  :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  Cmp1.imag:=2*Cmp1.real*Cmp1.imag;
  Cmp1.real:=a1;

  a1  := Cmp1.real;
  asm
    FLD   a1
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;
  Cmp1.real:=sXk*(a1+a2);
  Cmp1.imag:=cXk*(a1-a2);
End;

Function  CSin_ZZ(Cmp1: TComplex): TComplex;         // #36
Begin
  Result:=Cmp1;
  CSin_ZZV(Result);
End;

Procedure CCos_ZZV(var Cmp1: TComplex);         // #37
  var a1 ,a2 : Extended;
      sXk,cXk: Extended;
Begin
  a1  :=(Cmp1.real - Cmp1.imag)*(Cmp1.real + Cmp1.imag);
  Cmp1.imag:=2*Cmp1.real*Cmp1.imag;
  Cmp1.real:=a1;

  //  cos(x)cosh(y) - i sin(x)sinh(y)
  a1  := Cmp1.real;
  asm
    FLD   a1
    FSINCOS
    FNSTSW  AX
    SAHF
    JPO     @@a1
    FSTP    ST(0)
    FLD1    // FLDZ
    FLD1    // FLDZ
  @@a1:
    FSTP  cXk
    FSTP  sXk
    FWAIT
  end;
  a1:=Cmp1.imag;
  If  a1> ExpTol Then a1:= ExpTol;
  If  a1<-ExpTol Then a1:=-ExpTol;
  a1:=Exp(a1)/2;
  a2:=0.25/a1;
  Cmp1.real:= cXk*(a1+a2);
  Cmp1.imag:=-sXk*(a1-a2);
End;

Function  CCos_ZZ(Cmp1: TComplex): TComplex;         // #37
Begin
  Result:=Cmp1;
  CCos_ZZV(Result);
End;


{ ------------------------------------------------------------------- }
Procedure FuncDisp(Fn: Integer; var Cmp1: TComplex);     //  dispatcher
Begin
  If Cmp1.real> BigTol Then Cmp1.real:= BigTol;
  If Cmp1.imag> BigTol Then Cmp1.imag:= BigTol;
  If Cmp1.real<-BigTol Then Cmp1.real:=-BigTol;
  If Cmp1.imag<-BigTol Then Cmp1.imag:=-BigTol;
  Case Fn of
    0: ;                     // CIdent(Cmp1);
    1: CCosV(Cmp1);
    2: CTanV(Cmp1);
    3: CTanhV(Cmp1);
    4: CCotanV(Cmp1);
    5: CCotanhV(Cmp1);
    6: CFlipV(Cmp1);
    7: CConjV(Cmp1);
    8: CZeroV(Cmp1);
    9: CASinV(Cmp1);
   10: CASinhV(Cmp1);
   11: CACosV(Cmp1);
   12: CACoshV(Cmp1);
   13: CATanV(Cmp1);
   14: CATanhV(Cmp1);
   15: CCabsV(Cmp1);
   16: CAbsV(Cmp1);
   17: CSqrtV(Cmp1);
   18: CFloorV(Cmp1);
   19: CCeilV(Cmp1);
   20: CTruncV(Cmp1);
   21: CRoundV(Cmp1);
   22: COneV(Cmp1);
   23: CSinV(Cmp1);
   24: CCosxxV(Cmp1);
   25: CSinhV(Cmp1);
   26: CCoshV(Cmp1);
   27: CExpV(Cmp1);
   28: CLogV(Cmp1);
   29: CSqrV(Cmp1);
   30: CRecipV(Cmp1);

   31: CTripleV(Cmp1);
   32: CFourV(Cmp1);
   33: CSqrRecipV(Cmp1);
   34: CTripRecipV(Cmp1);
   35: CFourRecipV(Cmp1);
   36: CSin_ZZV(Cmp1);
   37: CCos_ZZV(Cmp1);
   38: CRealV(Cmp1);
   39: CImagV(Cmp1);
   40: CAbs2V(Cmp1);
   41: CCosh2V(Cmp1);
  end;
End;


Function FuncDispF(Fn: Integer; Cmp1: TComplex): TComplex;     //  new dispatcher  ++26.09.2004
Begin
  If Cmp1.real> BigTol Then Cmp1.real:= BigTol;
  If Cmp1.imag> BigTol Then Cmp1.imag:= BigTol;
  If Cmp1.real<-BigTol Then Cmp1.real:=-BigTol;
  If Cmp1.imag<-BigTol Then Cmp1.imag:=-BigTol;
  Case Fn of
    0: Result:=Cmp1;                    // CIdent(Cmp1);
    1: Result:=CCos(Cmp1);
    2: Result:=CTan(Cmp1);
    3: Result:=CTanh(Cmp1);
    4: Result:=CCotan(Cmp1);
    5: Result:=CCotanh(Cmp1);
    6: Result:=CFlip(Cmp1);
    7: Result:=CConj(Cmp1);
    8: Result:=CZero(Cmp1);
    9: Result:=CASin(Cmp1);
   10: Result:=CASinh(Cmp1);
   11: Result:=CACos(Cmp1);
   12: Result:=CACosh(Cmp1);
   13: Result:=CATan(Cmp1);
   14: Result:=CATanh(Cmp1);
   15: Result:=CCabs(Cmp1);
   16: Result:=CAbs(Cmp1);
   17: Result:=CSqrt(Cmp1);
   18: Result:=CFloor(Cmp1);
   19: Result:=CCeil(Cmp1);
   20: Result:=CTrunc(Cmp1);
   21: Result:=CRound(Cmp1);
   22: Result:=COne(Cmp1);
   23: Result:=CSin(Cmp1);
   24: Result:=CCosxx(Cmp1);
   25: Result:=CSinh(Cmp1);
   26: Result:=CCosh(Cmp1);
   27: Result:=CExp(Cmp1);
   28: Result:=CLog(Cmp1);
   29: Result:=CSqr(Cmp1);
   30: Result:=CRecip(Cmp1);

   31: Result:=CTriple(Cmp1);
   32: Result:=CFour(Cmp1);
   33: Result:=CSqrRecip(Cmp1);
   34: Result:=CTripRecip(Cmp1);
   35: Result:=CFourRecip(Cmp1);
   36: Result:=CSin_ZZ(Cmp1);
   37: Result:=CCos_ZZ(Cmp1);
   38: Result:=CReal(Cmp1);
   39: Result:=CImag(Cmp1);
   40: Result:=CAbs2(Cmp1);
   41: Result:=CCosh2(Cmp1);
  end;
End;

END.
