(* ***********************************************
    Модуль работы с комплексными числами.
    Вадим Исаев, 2020
    Module for complex numbers.
    Vadim Isaev, 2020
*************************************************** *)

MODULE CMath;

IMPORT Math, Out;

TYPE
  complex* = POINTER TO RECORD
    re*: REAL;
    im*: REAL
  END;

VAR
  result: complex;

  i* : complex;
  _0*: complex;

(* Инициализация комплексного числа.
   Init complex number. *)
PROCEDURE CInit* (re : REAL; im: REAL): complex;
VAR
  temp: complex;
BEGIN
  NEW(temp);
  temp.re:=re;
  temp.im:=im;

  RETURN temp
END CInit;


(* Четыре основных арифметических операций.
   Four base operations  +, -, * , / *)

(* Сложение
   addition : z := z1 + z2 *)
PROCEDURE CAdd* (z1, z2: complex): complex;
BEGIN
  result.re := z1.re + z2.re;
  result.im := z1.im + z2.im;

  RETURN result
END CAdd;

(* Сложение с REAL.
   addition : z := z1 + r1 *)
PROCEDURE CAdd_r* (z1: complex; r1: REAL): complex;
BEGIN
  result.re := z1.re + r1;
  result.im := z1.im;

  RETURN result
END CAdd_r;

(* Сложение с INTEGER.
   addition : z := z1 + i1 *)
PROCEDURE CAdd_i* (z1: complex; i1: INTEGER): complex;
BEGIN
  result.re := z1.re + FLT(i1);
  result.im := z1.im;

  RETURN result
END CAdd_i;

(* Смена знака.
   substraction : z := - z1 *)
PROCEDURE CNeg (z1 : complex): complex;
BEGIN
  result.re := -z1.re;
  result.im := -z1.im;

  RETURN result
END CNeg;

(* Вычитание.
   substraction : z := z1 - z2 *)
PROCEDURE CSub* (z1, z2 : complex): complex;
BEGIN
  result.re := z1.re - z2.re;
  result.im := z1.im - z2.im;

  RETURN result
END CSub;

(* Вычитание REAL.
   substraction : z := z1 - r1 *)
PROCEDURE CSub_r1* (z1 : complex; r1 : REAL): complex;
BEGIN
  result.re := z1.re - r1;
  result.im := z1.im;

  RETURN result
END CSub_r1;

(* Вычитание из REAL.
   substraction : z := r1 - z1 *)
PROCEDURE CSub_r2* (r1 : REAL; z1 : complex): complex;
BEGIN
  result.re := r1 - z1.re;
  result.im := - z1.im;

  RETURN result
END CSub_r2;

(* Вычитание INTEGER.
   substraction : z := z1 - i1 *)
PROCEDURE CSub_i* (z1 : complex; i1 : INTEGER): complex;
BEGIN
  result.re := z1.re - FLT(i1);
  result.im := z1.im;

  RETURN result
END CSub_i;

(* Умножение.
   multiplication : z := z1 * z2 *)
PROCEDURE CMul (z1, z2 : complex): complex;
BEGIN
  result.re := (z1.re * z2.re) - (z1.im * z2.im);
  result.im := (z1.re * z2.im) + (z1.im * z2.re);

  RETURN result
END CMul;

(* Умножение с REAL.
   multiplication : z := z1 * r1 *)
PROCEDURE CMul_r (z1 : complex; r1 : REAL): complex;
BEGIN
  result.re := z1.re * r1;
  result.im := z1.im * r1;

  RETURN result
END CMul_r;

(* Умножение с INTEGER.
   multiplication : z := z1 * i1 *)
PROCEDURE CMul_i (z1 : complex; i1 : INTEGER): complex;
BEGIN
  result.re := z1.re * FLT(i1);
  result.im := z1.im * FLT(i1);

  RETURN result
END CMul_i;

(* Деление.
   division : z := znum / zden *)
PROCEDURE CDiv (z1, z2 : complex): complex;
    (* The following algorithm is used to properly handle
      denominator overflow:

                 |  a + b(d/c)   c - a(d/c)
                 |  ---------- + ---------- I     if |d| < |c|
      a + b I    |  c + d(d/c)   a + d(d/c)
      -------  = |
      c + d I    |  b + a(c/d)   -a+ b(c/d)
                 |  ---------- + ---------- I     if |d| >= |c|
                 |  d + c(c/d)   d + c(c/d)
    *)
VAR
  tmp, denom : REAL;
BEGIN
   IF ( ABS(z2.re) > ABS(z2.im) ) THEN
     tmp := z2.im / z2.re;
     denom := z2.re + z2.im * tmp;
     result.re := (z1.re + z1.im * tmp) / denom;
     result.im := (z1.im - z1.re * tmp) / denom;
   ELSE
     tmp := z2.re / z2.im;
     denom := z2.im + z2.re * tmp;
     result.re := (z1.im + z1.re * tmp) / denom;
     result.im := (-z1.re + z1.im * tmp) / denom;
   END;

   RETURN result
END CDiv;

(* Деление на REAL.
   division : z := znum / r1 *)
PROCEDURE CDiv_r* (z1 : complex; r1 : REAL): complex;
BEGIN
  result.re := z1.re / r1;
  result.im := z1.im / r1;

  RETURN result
END CDiv_r;

(* Деление на INTEGER.
   division : z := znum / i1 *)
PROCEDURE CDiv_i* (z1 : complex; i1 : INTEGER): complex;
BEGIN
  result.re := z1.re / FLT(i1);
  result.im := z1.im / FLT(i1);

  RETURN result
END CDiv_i;

(* fonctions elementaires *)

(* Вывод на экран.
   out complex number *)
PROCEDURE CPrint* (z: complex; width: INTEGER);
BEGIN
  Out.Real(z.re, width);
  IF z.im>=0.0 THEN
    Out.String("+");
  END;
  Out.Real(z.im, width);
  Out.String("i");
END CPrint;

PROCEDURE CPrintLn* (z: complex; width: INTEGER);
BEGIN
  CPrint(z, width);
  Out.Ln;
END CPrintLn;

(* Вывод на экран с фиксированным кол-вом знаков
   после запятой (p) *)
PROCEDURE CPrintFix* (z: complex; width, p: INTEGER);
BEGIN
  Out.FixReal(z.re, width, p);
  IF z.im>=0.0 THEN
    Out.String("+");
  END;
  Out.FixReal(z.im, width, p);
  Out.String("i");
END CPrintFix;

PROCEDURE CPrintFixLn* (z: complex; width, p: INTEGER);
BEGIN
  CPrintFix(z, width, p);
  Out.Ln;
END CPrintFixLn;

(* Модуль числа.
   module : r = |z| *)
PROCEDURE CMod* (z1 : complex): REAL;
BEGIN
  RETURN Math.sqrt((z1.re * z1.re) + (z1.im * z1.im))
END CMod;

(* Квадрат числа.
   square : r := z*z *)
PROCEDURE CSqr* (z1: complex): complex;
BEGIN
  result.re := z1.re * z1.re - z1.im * z1.im;
  result.im := 2.0 * z1.re * z1.im;

  RETURN result
END CSqr;

(* Квадратный корень числа.
   square root : r := sqrt(z) *)
PROCEDURE CSqrt* (z1: complex): complex;
VAR
  root, q: REAL;
BEGIN
  IF (z1.re#0.0) OR (z1.im#0.0) THEN
    root := Math.sqrt(0.5 * (ABS(z1.re) + CMod(z1)));
    q := z1.im / (2.0 * root);
    IF z1.re >= 0.0 THEN
      result.re := root;
      result.im := q;
    ELSE
      IF z1.im < 0.0 THEN
        result.re := - q;
        result.im := - root
      ELSE
        result.re :=  q;
        result.im :=  root
      END
    END
  ELSE
    result := z1;
  END;

  RETURN result
END CSqrt;

(* Экспонента.
   exponantial : r := exp(z) *)
(* exp(x + iy) = exp(x).exp(iy) = exp(x).[cos(y) + i sin(y)] *)
PROCEDURE CExp* (z: complex): complex;
VAR
  expz : REAL;
BEGIN
  expz := Math.exp(z.re);
  result.re := expz * Math.cos(z.im);
  result.im := expz * Math.sin(z.im);

  RETURN result
END CExp;

(* Натуральный логарифм.
   natural logarithm : r := ln(z) *)
(* ln( p exp(i0)) = ln(p) + i0 + 2kpi *)
PROCEDURE CLn* (z: complex): complex;
BEGIN
  result.re := Math.ln(CMod(z));
  result.im := Math.arctan2(z.im, z.re);

  RETURN result
END CLn;

(* Число в степени.
   exp : z := z1^z2 *)
PROCEDURE CPower* (z1, z2 : complex): complex;
VAR
  a: complex;
BEGIN
   a:=CLn(z1);
   a:=CMul(z2, a);
   result:=CExp(a);

   RETURN result
END CPower;

(* Число в степени REAL.
   multiplication : z := z1^r *)
PROCEDURE CPower_r* (z1: complex; r: REAL): complex;
VAR
  a: complex;
BEGIN
  a:=CLn(z1);
  a:=CMul_r(a, r);
  result:=CExp(a);

  RETURN result
END CPower_r;

(* Обратное число.
   inverse : r := 1 / z *)
PROCEDURE CInv* (z: complex): complex;
VAR
  denom : REAL;
BEGIN
  denom := (z.re * z.re) + (z.im * z.im);
  (* generates a fpu exception if denom=0 as for reals *)
  result.re:=z.re/denom;
  result.im:=-z.im/denom;

  RETURN result
END CInv;

(* direct trigonometric functions *)

(* Косинус.
   complex cosinus *)
(* cos(x+iy) = cos(x).cos(iy) - sin(x).sin(iy) *)
(* cos(ix) = cosh(x) et sin(ix) = i.sinh(x) *)
PROCEDURE CCos* (z: complex): complex;
BEGIN
  result.re := Math.cos(z.re) * Math.cosh(z.im);
  result.im := - Math.sin(z.re) * Math.sinh(z.im);

  RETURN result
END CCos;

(* Синус.
   sinus complex *)
(* sin(x+iy) = sin(x).cos(iy) + cos(x).sin(iy) *)
(* cos(ix) = cosh(x) et sin(ix) = i.sinh(x) *)
PROCEDURE CSin (z: complex): complex;
BEGIN
  result.re := Math.sin(z.re) * Math.cosh(z.im);
  result.im := Math.cos(z.re) * Math.sinh(z.im);

  RETURN result
END CSin;

(* Тангенс.
   tangente *)
PROCEDURE CTg* (z: complex): complex;
VAR
  temp1, temp2: complex;
BEGIN
  temp1:=CSin(z);
  temp2:=CCos(z);
  result:=CDiv(temp1, temp2);

  RETURN result
END CTg;

(* inverse complex hyperbolic functions *)

(* Гиперболический арккосинус.
   hyberbolic arg cosinus *)
(*                          _________  *)
(* argch(z) = -/+ ln(z + i.V 1 - z.z)  *)
PROCEDURE CArcCosh* (z : complex): complex;
BEGIN
  result:=CNeg(CLn(CAdd(z, CMul(i, CSqrt(CSub_r2(1.0, CMul(z, z)))))));

  RETURN result
END CArcCosh;

(* Гиперболический арксинус.
   hyperbolic arc sinus       *)
(*                    ________  *)
(* argsh(z) = ln(z + V 1 + z.z) *)
PROCEDURE CArcSinh* (z : complex): complex;
BEGIN
  result:=CLn(CAdd(z, CSqrt(CAdd_r(CMul(z, z), 1.0))));

  RETURN result
END CArcSinh;

(* Гиперболический арктангенс.
   hyperbolic arc tangent *)
(* argth(z) = 1/2 ln((z + 1) / (1 - z)) *)
PROCEDURE CArcTgh (z : complex): complex;
BEGIN
  result:=CDiv_r(CLn(CDiv(CAdd_r(z, 1.0), CSub_r2(1.0, z))), 2.0);

  RETURN result
END CArcTgh;

(* trigonometriques inverses *)

(* Арккосинус.
   arc cosinus complex *)
(* arccos(z) = -i.argch(z) *)
PROCEDURE CArcCos* (z: complex): complex;
BEGIN
  result := CNeg(CMul(i, CArcCosh(z)));

  RETURN result
END CArcCos;

(* Арксинус.
   arc sinus complex *)
(* arcsin(z) = -i.argsh(i.z) *)
PROCEDURE CArcSin* (z : complex): complex;
BEGIN
  result := CNeg(CMul(i, CArcSinh(z)));

  RETURN result
END CArcSin;

(* Арктангенс.
   arc tangente complex *)
(* arctg(z) = -i.argth(i.z) *)
PROCEDURE CArcTg* (z : complex): complex;
BEGIN
  result := CNeg(CMul(i, CArcTgh(CMul(i, z))));

  RETURN result
END CArcTg;

BEGIN

  result:=CInit(0.0, 0.0);
  i :=CInit(0.0, 1.0);
  _0:=CInit(0.0, 0.0);

END CMath.