(* ***********************************************
Модуль работы с комплексными числами.
Вадим Исаев, 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.