Как работать с комплексными числами?
{ 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. |
