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

Взято из http://www.ibrtses.com/delphi/complexrec.html
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.

Copyright (99,2000) Ing.BЁ№ro R.Tschaggelar
Взято с Vingrad.ru http://forum.vingrad.ru

В процедурах C_Div и C_DivN имеются грубейшие ошибки! Неверно вычисляется результат деления комплексных чисел! Эти процедуры разрабатывал неуч, не знающий элементарной математики! По-видимому, аналогичные ошибки есть и в последующих процедурах...

В процедуре C_Div должно быть: t := (a.a*b.a+a.b*b.b)/(b.a*b.a+b.b*b.b); a.b := (a.b*b.a-a.a*b.b)/(b.a*b.a+b.b*b.b); a.a := t;

Отправить комментарий

Проверка
Антиспам проверка
Image CAPTCHA
...