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

Previous  Top  Next

    
 

 

Взято из 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.

Code:

{ 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

©Drkb::04064

Взято с Vingrad.ru http://forum.vingrad.ru