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

Delphi , Синтаксис , Математика

Автор: 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 и FreePascal.

Комментарии и вопросы

:: Главная :: Математика ::