Как работать с комплексными числами?

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_records
-----------------------------------------
they are efficient on arrays
}
unit ComplexRec;
interface
type
 float = extended;
 PComplex = ^TComplex;
 TComplex = record // C_record without rectangular/polar discrimination
 a, b: float; // (re,im) or (abs,arg)
 end;
function C_Copy(const a: TComplex): TComplex; // Result := a
function C_One: TComplex; // Result := 1 BOTH
function C_I: TComplex; // Result := i RECTANGULAR
function C_IP: TComplex; // Result := i POLAR
procedure C_P2R(var a: TComplex); // polar to rectangular
procedure C_R2P(var a: TComplex); // rectangular to polar
function C_abs(const a: TComplex): float; // RECTANGULAR
function C_arg(const a: TComplex): float; // RECTANGULAR
function C_re(const a: TComplex): float; // POLAR
function C_im(const a: TComplex): float; // POLAR
procedure C_Inv(var a: TComplex); // a := -a RECTANGULAR
procedure C_InvP(var a: TComplex); // a := -a POLAR
procedure C_Conj(var a: TComplex); // a := konjug(a) BOTH
function C_ConjN(const a: TComplex): TComplex; // Result := konjug(a) BOTH
procedure C_Scale(var a: TComplex; u: float); // a := a * u;
procedure C_ScaleP(var a: TComplex; u: float); // a := a * u;
procedure C_Add(var a, b: TComplex); // a := a + b RECTANGULAR
function C_AddN(const a, b: TComplex): TComplex; // Result := a + b RECTANGULAR
procedure C_Sub(var a, b: TComplex); // a := a - b RECTANGULAR
function C_SubN(const a, b: TComplex): TComplex; // Result := a - b RECTANGULAR
procedure C_Mul(var a, b: TComplex); //a:=a*b RECTANGULAR
function C_MulN(const a, b: TComplex): TComplex; //result:=a*b RECTANGULAR
procedure C_MulP(var a, b: TComplex); //a:=a*b POLAR
function C_MulNP(const a, b: TComplex): TComplex; //result:=a*b POLAR
procedure C_DivP(var a, b: TComplex); //a:=a/b POLAR
function C_DivNP(const a, b: TComplex): TComplex; //result:=a/b POLAR
procedure C_Div(var a, b: TComplex); //a:=a/b POLAR
function C_DivN(const a, b: TComplex): TComplex; //result := a / b RECTANGULAR
function C_ExpN(const a: TComplex): TComplex; // RECTANGLE
function C_LogN(const a: TComplex): TComplex; // POLAR
function C_SinN(const a: TComplex): TComplex;
function C_CosN(const a: TComplex): TComplex;
function C_TanN(a: TComplex): TComplex;
function C_SinhN(const a: TComplex): TComplex;
function C_CoshN(const a: TComplex): TComplex;
function C_TanhN(a: TComplex): TComplex;
function C_IntPowerN(const a: TComplex; n: Integer): TComplex; // RECTANGLE
function C_IntPowerNP(const a: TComplex; n: Integer): TComplex; // POLAR
function C_ParallelN(const a, b: TComplex): TComplex;
// result:=a//b =(a*b)/(a+b) RECTANGULAR
// electronic parallel circuit
implementation
uses Math;
const
 AlmostZero = 1E-30; // test for zero
function C_Copy(const a: TComplex): TComplex; // result:=a
begin
 Result.a := a.a;
 Result.b := a.b;
end;
function C_One: TComplex; // result:=1
begin
 Result.a := 1;
 Result.b := 0;
end;
function C_I: TComplex; // result:=i RECTANGULAR
begin
 Result.a := 0;
 Result.b := 1;
end;
function C_IP: TComplex; // result:=i POLAR
begin
 Result.a := 1;
 Result.b := pi / 2;
end;
procedure C_P2R(var a: TComplex);
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(var a: TComplex);
var
 t: float;
begin
 t := a.a;
 a.a := sqrt(sqr(a.a) + sqr(a.b));
 a.b := arctan(a.b / t);
 if (t < 0) then a.b := a.b + pi;
end;
function C_abs(const a: TComplex): float;
begin
 Result := sqrt(sqr(a.a) + sqr(a.b));
end;
function C_arg(const a: TComplex): float;
begin
 Result := arctan2(a.b, a.a);
end;
function C_re(const a: TComplex): float; // POLAR
begin
 Result := a.a * cos(a.b);
end;
function C_im(const a: TComplex): float; // POLAR
begin
 Result := a.a * sin(a.b);
end;
procedure C_Inv(var a: TComplex); // a := -a RECTANGULAR
begin
 a.a := -a.a;
 a.b := -a.b;
end;
procedure C_InvP(var a: TComplex); // a := -a POLAR
begin
 a.b := a.b + pi;
end;
procedure C_Conj(var a: TComplex); // a := konjug(a) BOTH
begin
 a.b := -a.b;
end;
function C_ConjN(const a: TComplex): TComplex; // result := konjug(a) BOTH
begin
 Result.a := a.a;
 Result.b := -a.b;
end;
procedure C_Scale(var a: TComplex; u: float); // a := a * u;
begin
 a.a := a.a * u;
 a.b := a.b * u;
end;
procedure C_ScaleP(var a: TComplex; u: float); // a := a * u;
begin
 a.a := a.a * u;
end;
procedure C_Add(var a, b: TComplex); // a := a + b RECTANGULAR
begin
 a.a := a.a + b.a;
 a.b := a.b + b.b;
end;
function C_AddN(const a, b: TComplex): TComplex; // Result := a + b RECTANGULAR
begin
 Result.a := a.a + b.a;
 Result.b := a.b + b.b;
end;
procedure C_Sub(var a, b: TComplex); // a := a - b RECTANGULAR
begin
 a.a := a.a - b.a;
 a.b := a.b - b.b;
end;
function C_SubN(const a, b: TComplex): TComplex; // Result := a - b RECTANGULAR
begin
 Result.a := a.a - b.a;
 Result.b := a.b - b.b;
end;
procedure C_Mul(var a, b: TComplex); // 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(const a, b: TComplex): TComplex; // Result := a * b RECTANGULAR
begin
 Result.a := a.a * b.a - a.b * b.b;
 Result.b := a.a * b.b + a.b * b.a;
end;
procedure C_MulP(var a, b: TComplex); // a := a * b POLAR
begin
 a.a := a.a * b.a;
 a.b := a.b + b.b;
end;
function C_MulNP(const a, b: TComplex): TComplex; // Result := a * b POLAR
begin
 Result.a := a.a * b.a;
 Result.b := a.b + b.b;
end;
procedure C_Div(var a, b: TComplex); // 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(const a, b: TComplex): TComplex; // Result := a / b RECTANGULAR
begin
 Result.a := a.a / b.a + a.b / b.b;
 Result.b := -a.a / b.b + a.b / b.a;
end;
procedure C_DivP(var a, b: TComplex); // a := a / b POLAR
begin
 a.a := a.a / b.a;
 a.b := a.b - b.b;
end;
function C_DivNP(const a, b: TComplex): TComplex; // Result := a / b POLAR
begin
 Result.a := a.a / b.a;
 Result.b := a.b - b.b;
end;
function C_ExpN(const a: TComplex): TComplex; // RECTANGLE
begin
 Result.a := exp(a.a);
 Result.b := a.b;
 C_P2R(Result);
end;
function C_LogN(const a: TComplex): TComplex; // POLAR
begin
 Result.a := ln(a.a);
 Result.b := a.b;
 C_R2P(Result);
end;
function C_SinN(const a: TComplex): TComplex;
var
 z, n, v, t: TComplex;
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);
end;
function C_CosN(const a: TComplex): TComplex;
var
 z, n, v, t: TComplex;
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);
end;
function C_TanN(a: TComplex): TComplex;
begin
end;
function C_SinhN(const a: TComplex): TComplex;
var
 u, v, t: TComplex;
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);
end;
function C_CoshN(const a: TComplex): TComplex;
var
 u, v, t: TComplex;
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);
end;
function C_TanhN(a: TComplex): TComplex;
begin
end;
function C_IntPowerN(const a: TComplex; n: Integer): TComplex;
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(const a: TComplex; n: Integer): TComplex;
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(const a, b: TComplex): TComplex; // result:=a//b = (a*b)/(a+b)
var
 z, n: TComplex;
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);
end;
end.
Copyright (99,2000) Ing.BЁ№ro R.Tschaggelar
Source: Взято с Vingrad.ru http://forum.vingrad.ru
ID: 04064
Category:Алгоритмы
Category:Математика