﻿ Delphi World - Как работать с комплексными числами
 Как работать с комплексными числами Автор: R.Tschaggelar WEB-сайт: http://forum.vingrad.ru Complex numbers Complex numbers have two representations : rectanglar : Z = a + i * b, a being the real part, and b being the imaginary part polar : Z = r * exp(i * phi), r being the absolute value, and phi being the argument(angle) a reason to demotivate compiler writers to have it as native type. Here is a unit that approaches the complex as record. the used record is of dual use, either rectangular or polar, one just has to keep in mind what in is at the moment. ```{ unit for complex numbers based on C_reords ----------------------------------------- they are efficient on arrays } unit ComplexRec; interface type float = extended; ComplexPtr = ^Complex; Complex = record // C_record without rectangular/polar discrimination a, b: float; // (re,im) or (abs,arg) end; function C_Copy(a: ComplexPtr): ComplexPtr; // result:=a function C_One: ComplexPtr; // result:=1 BOTH function C_I: ComplexPtr; // result:=i RECTANGULAR function C_IP: ComplexPtr; // result:=i POLAR procedure C_P2R(a: ComplexPtr); // polar to rectangular procedure C_R2P(a: ComplexPtr); // rectangular to polar function C_abs(a: ComplexPtr): float; // RECTANGULAR function C_arg(a: ComplexPtr): float; // RECTANGULAR function C_re(a: ComplexPtr): float; // POLAR function C_im(a: ComplexPtr): float; // POLAR procedure C_Inv(a: ComplexPtr); // a:=-a RECTANGULAR procedure C_InvP(a: ComplexPtr); // a:=-a POLAR procedure C_Conj(a: ComplexPtr); // a:=konjug(a) BOTH function C_ConjN(a: ComplexPtr): ComplexPtr; //result:=konjug(a) BOTH procedure C_Scale(a: ComplexPtr; u: float); // a:=a*u; procedure C_ScaleP(a: ComplexPtr; u: float); // a:=a*u; procedure C_Add(a, b: ComplexPtr); //a:=a+b RECTANGULAR function C_AddN(a, b: ComplexPtr): ComplexPtr; //result:=a+b RECTANGULAR procedure C_Sub(a, b: ComplexPtr); //a:=a-b RECTANGULAR function C_SubN(a, b: ComplexPtr): ComplexPtr; //result:=a-b RECTANGULAR procedure C_Mul(a, b: ComplexPtr); //a:=a*b RECTANGULAR function C_MulN(a, b: ComplexPtr): ComplexPtr; //result:=a*b RECTANGULAR procedure C_MulP(a, b: ComplexPtr); //a:=a*b POLAR function C_MulNP(a, b: ComplexPtr): ComplexPtr; //result:=a*b POLAR procedure C_DivP(a, b: ComplexPtr); //a:=a/b POLAR function C_DivNP(a, b: ComplexPtr): ComplexPtr; //result:=a/b POLAR procedure C_Div(a, b: ComplexPtr); //a:=a/b POLAR function C_DivN(a, b: ComplexPtr): ComplexPtr; //result:=a/b POLAR function C_ExpN(a: ComplexPtr): ComplexPtr; // RECTANGLE function C_LogN(a: ComplexPtr): ComplexPtr; // POLAR function C_SinN(a: ComplexPtr): ComplexPtr; function C_CosN(a: ComplexPtr): ComplexPtr; function C_TanN(a: ComplexPtr): ComplexPtr; function C_SinhN(a: ComplexPtr): ComplexPtr; function C_CoshN(a: ComplexPtr): ComplexPtr; function C_TanhN(a: ComplexPtr): ComplexPtr; function C_IntPowerN(a: ComplexPtr; n: integer): ComplexPtr; // RECTANGLE function C_IntPowerNP(a: ComplexPtr; n: integer): ComplexPtr; // POLAR function C_ParallelN(a, b: ComplexPtr): ComplexPtr; // result:=a//b =(a*b)/(a+b) RECTANGULAR // electronic parallel circuit implementation uses math; const AlmostZero = 1E-30; // test for zero function C_Copy(a: ComplexPtr): ComplexPtr; // result:=a begin result := new(ComplexPtr); result.a := a.a; result.b := a.b; end; function C_One: ComplexPtr; // result:=1 begin result := new(ComplexPtr); result.a := 1; result.b := 0; end; function C_I: ComplexPtr; // result:=i RECTANGULAR begin result := new(ComplexPtr); result.a := 0; result.b := 1; end; function C_IP: ComplexPtr; // result:=i POLAR begin result := new(ComplexPtr); result.a := 1; result.b := pi / 2; end; procedure C_P2R(a: ComplexPtr); var t, u, v: float; begin t := a.a; sincos(a.b, u, v); a.a := t * v; a.b := t * u; end; procedure C_R2P(a: ComplexPtr); var t: float; begin t := a.a; a.a := sqrt(sqr(a.a) + sqr(a.b)); if (abs(t)0 then a.b := pi / 2 else a.b := -pi / 2; end else begin a.b := arctan(a.b / t); if (t < 0) then a.b := a.b + pi; end; end; function C_abs(a: ComplexPtr): float; begin result := sqrt(sqr(a.a) + sqr(a.b)); end; function C_arg(a: ComplexPtr): float; begin if (abs(a.a)0 then result := pi / 2 else result := -pi / 2; end else begin result := arctan(a.b / a.a); if (a.a < 0) then result := result + pi; end; end; function C_re(a: ComplexPtr): float; // POLAR begin result := a.a * cos(a.b); end; function C_im(a: ComplexPtr): float; // POLAR begin result := a.a * sin(a.b); end; procedure C_Inv(a: ComplexPtr); // a:=-a RECTANGULAR begin a.a := -a.a; a.b := -a.b; end; procedure C_InvP(a: ComplexPtr); // a:=-a POLAR begin a.b := a.b + pi; end; procedure C_Conj(a: ComplexPtr); // a:=konjug(a) BOTH begin a.b := -a.b; end; function C_ConjN(a: ComplexPtr): ComplexPtr; //result:=konjug(a) BOTH begin result := new(ComplexPtr); result.a := a.a; result.b := -a.b; end; procedure C_Scale(a: ComplexPtr; u: float); // a:=a*u; begin a.a := a.a * u; a.b := a.b * u; end; procedure C_ScaleP(a: ComplexPtr; u: float); // a:=a*u; begin a.a := a.a * u; end; procedure C_Add(a, b: ComplexPtr); //a:=a+b RECTANGULAR begin a.a := a.a + b.a; a.b := a.b + b.b; end; function C_AddN(a, b: ComplexPtr): ComplexPtr; //result:=a+b RECTANGULAR begin result := new(ComplexPtr); result.a := a.a + b.a; result.b := a.b + b.b; end; procedure C_Sub(a, b: ComplexPtr); //a:=a-b RECTANGULAR begin a.a := a.a - b.a; a.b := a.b - b.b; end; function C_SubN(a, b: ComplexPtr): ComplexPtr; //result:=a-b RECTANGULAR begin result := new(ComplexPtr); result.a := a.a - b.a; result.b := a.b - b.b; end; procedure C_Mul(a, b: ComplexPtr); //a:=a*b RECTANGULAR var u, v: float; begin u := a.a * b.a - a.b * b.b; v := a.a * b.b + a.b * b.a; a.a := u; a.b := v; end; function C_MulN(a, b: ComplexPtr): ComplexPtr; //result:=a*b RECTANGULAR begin result := new(ComplexPtr); result.a := a.a * b.a - a.b * b.b; result.b := a.a * b.b + a.b * b.a; end; procedure C_MulP(a, b: ComplexPtr); //a:=a*b POLAR begin a.a := a.a * b.a; a.b := a.b + b.b; end; function C_MulNP(a, b: ComplexPtr): ComplexPtr; //result:=a*b POLAR begin result := new(ComplexPtr); result.a := a.a * b.a; result.b := a.b + b.b; end; procedure C_Div(a, b: ComplexPtr); //a:=a/b RECTANGULAR var t: float; begin t := a.a / b.a + a.b / b.b; a.b := -a.a / b.b + a.b / b.a; a.a := t; end; function C_DivN(a, b: ComplexPtr): ComplexPtr; //result:=a/b RECTANGULAR begin result := new(ComplexPtr); result.a := a.a / b.a + a.b / b.b; result.b := -a.a / b.b + a.b / b.a; end; procedure C_DivP(a, b: ComplexPtr); //a:=a/b POLAR begin a.a := a.a / b.a; a.b := a.b - b.b; end; function C_DivNP(a, b: ComplexPtr): ComplexPtr; //result:=a/b POLAR begin result := new(ComplexPtr); result.a := a.a / b.a; result.b := a.b - b.b; end; function C_ExpN(a: ComplexPtr): ComplexPtr; // RECTANGLE begin result := new(ComplexPtr); result.a := exp(a.a); result.b := a.b; C_P2R(result); end; function C_LogN(a: ComplexPtr): ComplexPtr; // POLAR begin result := new(ComplexPtr); result.a := ln(a.a); result.b := a.b; C_R2P(result); end; function C_SinN(a: ComplexPtr): ComplexPtr; var z, n, v, t: ComplexPtr; begin t := C_I; v := C_MulN(a, t); // i*a z := C_expN(a); // exp(i*a) t := C_Copy(v); C_Inv(t); // -i*a t := C_ExpN(v); // exp(-i*a) C_Sub(z, t); n := C_I; C_Scale(n, 2); result := C_DivN(z, n); dispose(z); dispose(n); dispose(v); dispose(t); end; function C_CosN(a: ComplexPtr): ComplexPtr; var z, n, v, t: ComplexPtr; begin t := C_I; v := C_MulN(a, t); // i*a z := C_expN(a); // exp(i*a) t := C_Copy(v); C_Inv(t); // -i*a t := C_ExpN(v); // exp(-i*a) C_Add(z, t); n := C_One; C_Scale(n, 2); result := C_DivN(z, n); dispose(z); dispose(n); dispose(v); dispose(t); end; function C_TanN(a: ComplexPtr): ComplexPtr; begin end; function C_SinhN(a: ComplexPtr): ComplexPtr; var u, v, t: ComplexPtr; begin u := C_ExpN(a); t := C_Copy(a); C_inv(t); v := C_ExpN(t); result := C_SubN(u, v); C_Scale(result, 1 / 2); dispose(u); dispose(v); dispose(t); end; function C_CoshN(a: ComplexPtr): ComplexPtr; var u, v, t: ComplexPtr; begin u := C_ExpN(a); t := C_Copy(a); C_inv(t); v := C_ExpN(t); result := C_AddN(u, v); C_Scale(result, 1 / 2); dispose(u); dispose(v); dispose(t); end; function C_TanhN(a: ComplexPtr): ComplexPtr; begin end; function C_IntPowerN(a: ComplexPtr; n: integer): ComplexPtr; var j: integer; u, v: float; begin if n = 0 then result := C_One else begin result := C_Copy(a); if n > 1 then begin C_R2P(result); u := result.a; v := result.b; for j := 2 to n do begin u := u * result.a; v := v + result.b; end; result.a := u; result.b := v; C_P2R(result); end; if n < 0 then begin end; end; end; function C_IntPowerNP(a: ComplexPtr; n: integer): ComplexPtr; var j: integer; u, v: float; begin result := C_Copy(a); u := result.a; v := result.b; for j := 2 to n do begin u := u * result.a; v := v + result.b; end; result.a := u; result.b := v; end; function C_ParallelN(a, b: ComplexPtr): ComplexPtr; // result:=a//b = (a*b)/(a+b) var z, n: ComplexPtr; begin z := C_MulN(a, b); n := C_AddN(a, b); C_R2P(n); C_R2P(z); result := C_DivNP(z, n); C_P2R(result); dispose(n); dispose(z); end; end.```
Проект Delphi World © Выпуск 2002 - 2017
Автор проекта: Эксклюзивные курсы программирования