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