complextype.pas

Delphi math processing compononets and sources. Release.

{
@abstract(EBK&NVS Pascal-Delphi Math Library: Complex Math Types and constants)
@author(Nikolai V. Shokhirev <[e-mail]> <[e-mail]>)
@author(Eugene B. Krissinel <[e-mail]> <[e-mail]>)
@created(02.02.1994)
@lastmod(04.04.2002)
This is a temporary publication (reduced variant), will be updated later
check: http://www.shokhirev.com/nikolai/programs/samplecode.html
㎞ikolai V. Shokhirev, 2002
}
unit ComplexType;

interface

uses
  MathTypes{, SpecFunc};
  
type

  Complex = record
              Re,Im: RealType;
            end;
{ The dimension in declarations can be of any value }
  CVector   =  array [1..999] of Complex;
  CVcPtr    =  ^CVector;
  CMatrix   =  array [1..999] of CVcPtr;
  CMtPtr    =  ^CMatrix;

procedure  GetCVectorMemory (var Vc: CVcPtr; N: integer );

procedure  FreeCVectorMemory(var Vc: CVcPtr; N: integer );

procedure  GetCMatrixMemory (var Mt: CMtPtr; M, N: integer );

procedure  FreeCMatrixMemory(var Mt: CMtPtr; M, N: integer );

{ Conversion functions for complex }

function CmplxToStr(z: complex): string;
function StrToCmplx(s: string): complex;

{ Basic functions for complex Calculations }

function cmplx(x: RealType; y: RealType = 0.0): complex; { result := x + i*y }
function Cabs (z: complex): RealType; { sqrt(sqr(Re)+sqr(Im))    }
function Cabs0(z: complex): RealType; { max(abs(z.re),abs(z.im)) }
function Cabs1(z: complex): RealType; { abs(z.Re) + abs(z.Im)    }
function Cabs2(z: complex): RealType; { sqr(z.Re)+ sqr(z.Im)     }
function Cneg(z: complex): complex; { result := -z }
function Cinv(z: complex): complex; { result := 1/z }
function conjug(z: complex): complex;  { result := z.Re -i*z.Im }
function Re(z: complex): RealType;  { result := z.Re }
function Im(z: complex): RealType;  { result := z.Im }
function Cunit(z: complex; var r: RealType): complex; { result := z/r r := Cabs (z) }
function Cadd(z: complex;  w: complex): complex; overload; { result := z + w }
function Cadd(z: complex;  x: RealType): complex; overload;
function Cadd(x: RealType; w: complex): complex; overload;
function Csub(z: complex;  w: complex): complex; overload; { result := z - w }
function Csub(z: complex;  x: RealType): complex; overload;
function Csub(x: RealType; w: complex): complex; overload;
function Cmul(z: complex;  w: complex): complex; overload; { result := z * w }
function Cmul(x: RealType; w: complex): complex; overload;
function Cmul(z: complex;  x: RealType): complex; overload;
function Cdiv(z: complex;  w: complex): complex; overload; { result := z / w }
function Cdiv(x: RealType; w: complex): complex; overload;
function Cdiv(z: complex;  x: RealType): complex; overload;
function Cexp(z: complex): complex;         { z := exp(v)  }
function Csqrt(z: complex): complex;        { result := sqrt(z), result.Re > 0 }
function Croot(x: RealType): complex;        { result := sqrt(x) }
procedure Csum(var s: complex; z: complex); overload; { s := s + z }
procedure Csum(var s: complex; x: RealType); overload; { s := s + x }

{ Special functions for Complex Calculations }

function cer(z: complex): complex; { result := er(z) }

const
  Cmplx1: complex = (Re: 1.0; Im: 0.0);  // cmplx(1.0, 0.0)
  Cmplx0: complex = (Re: 0.0; Im: 0.0);  // cmplx(0.0, 0.0)

{ ===================================================================== }

implementation

uses
  SpecFunc, sysutils, Dialogs;

var
  S: pointer;

procedure  GetCVectorMemory(  var  Vc: CVcPtr;  N : integer );
begin
  GetMem( S, N*SizeOf(complex) );
  Vc := S;
end;

procedure  FreeCVectorMemory(  var  Vc: CVcPtr;  N : integer );
begin
  S := Vc;
  FreeMem( S, N*SizeOf(complex) );
  Vc := nil;
end;

procedure  GetCMatrixMemory(  var  Mt: CMtPtr;  M,N : integer );
var  i : integer;
begin
  GetMem( S, SizeOf(CVcPtr)*M );
  Mt := S;
  for i:=1 to M  do
    begin
      GetMem( S, SizeOf(complex)*N );
      Mt^[i] := S;
    end;
end;

procedure  FreeCMatrixMemory(  var  Mt: CMtPtr;  M,N : integer );
var  i : integer;
begin
  for i:=M downto 1  do
    begin
      S := Mt^[i];
      FreeMem( S, SizeOf(complex)*N );
    end;
  S := Mt;
  FreeMem( S, SizeOf(CVcPtr)*M );
  Mt := nil;
end;

{ Conversion functions for complex }

function CmplxToStr(z: complex): string;
var
  s: string;
begin
  if z.Im = 0.0 then
    result := FloatToStr(z.Re)
  else
  begin
    if z.Im > 0 then s := ' + i*' else s := ' - i*';
    result := FloatToStr(z.Re)+s+FloatToStr(abs(z.Im));
  end;
end;

function StrToCmplx(s: string): complex;
//  2.2 + i*1.1
// +2.2 - I*1.1
// -2.2 + j 1.1
//  2.2 + J*1.1
var
  r, i, ss: string;
  n: IntType;
  sig: RealType;
begin
  ss := StringReplace(s, ' ', '', [rfReplaceAll, rfIgnoreCase]);
//  L := Length(ss);
  n := Pos('i',ss);
  if n = 0 then n := Pos('I',ss);
  if n = 0 then n := Pos('j',ss);
  if n = 0 then n := Pos('J',ss);
  if n = 0 then // real
  begin
    r := ss;
    i := '';
  end else
  begin
    case n of
    1: begin // imaginary
         sig := 1.0;
         r := '';
         if ss[n+1] = '*' then i := copy(ss, n+2,512)
         else i := copy(ss, n+1,512);
       end;
    2: begin // imaginary
         r := '';
         if ss[1] = '-' then sig := -1.0 else sig := 1.0;
         if ss[n+1] = '*' then i := copy(ss, n+2,512)
         else i := copy(ss, n+1,512);
       end;
    else
      sig := 1.0;
      if ss[n-1] = '-' then
      begin
        sig := -1.0;
        r := copy(ss, 1, n-2);
      end else
        if ss[n-1] = '+' then
          r := copy(ss, 1, n-2)
        else // no sign = +
          r := copy(ss, 1, n-1);
      if n = Length(ss) then i := '1.0' else
        if ss[n+1] = '*' then i := copy(ss, n+2,512)
        else i := copy(ss, n+1,512);
    end;
  end;

  if r ='' then
    result.Re := 0.0
  else
    result.Re := StrToFloat(r);
    
  if i ='' then
    result.Im := 0.0 
  else
    result.Im := sig*StrToFloat(i);
end;

{ Basic functions for complex Calculations }

function cmplx(x: RealType; y: RealType = 0.0): complex;
{ result := x + i*y }
begin
  result.Re := x;
  result.Im := y;
end;

function Cabs (z: complex ): RealType;
{ sqrt(sqr(Re)+sqr(Im)) }
begin
  result := sqroot2(z.Re, z.Im);
end;

function Cabs0(z: complex): RealType;
{ max(abs(z.re),abs(z.im)) }
begin
  if (abs(z.re)>abs(z.im))
  then result := abs(z.re) else result := abs(z.im);
end;

function Cabs1(z: complex ): RealType;
{ abs(z.Re) + abs(z.Im)    }
begin
  result := abs(z.re) + abs(z.im);
end;

function Cabs2(z: complex ): RealType;
{ sqr(z.Re)+ sqr(z.Im)     }
begin
  result := sqr(z.re) + sqr(z.im);
end;

function Re(z: complex):RealType;
begin
  result := z.re;
end;

function Im(z: complex):RealType;
begin
  result := z.im;
end;

function Cneg(z: complex): complex;
{ result := -z }
begin
  result.Re := -z.Re;
  result.Im := -z.Im;
end;

function Cinv(z: complex): complex;
{ result := 1/z }
var
  r: RealType;
begin
  if Cabs0(z) > MinReal then
  begin
    r := Cabs2(z);
    result.Re :=  z.Re/r;
    result.Im := -z.Im/r;
  end else
  begin
    result.Re :=  sign(z.Re)*MaxReal;
    result.Im := -sign(z.Im)*MaxReal;
  end;
end;

function conjug(z: complex): complex;
{ result := z.Re -i*z.Im }
begin
  result.Re :=  z.Re;
  result.Im := -z.Im;
end;

function Cunit(z: complex; var r: RealType): complex;
{ result := z/r r := Cabs (z) }
begin
  r := Cabs(z);
  if r > MinReal then
  begin
    result.Re := z.Re/r;
    result.Im := z.Im/r;
  end else
  begin
    r := 0.0;
    result.Re := 0.0;
    result.Im := 0.0;
  end;
end;

function Cadd(z: complex;  w: complex): complex;
{ result := z + w }
begin
  result.Re :=  z.Re + w.Re;
  result.Im :=  z.Im + w.Im;
end;

function Cadd(z: complex;  x: RealType): complex;
begin
  result.Re :=  z.Re + x;
  result.Im :=  z.Im;
end;

function Cadd(x: RealType; w: complex): complex;
begin
  result.Re :=  w.Re + x;
  result.Im :=  w.Im;
end;

function Csub(z: complex;  w: complex): complex;
{ result := z - w }
begin
  result.Re :=  z.Re - w.Re;
  result.Im :=  z.Im - w.Im;
end;

function Csub(z: complex;  x: RealType): complex;
begin
  result.Re :=  z.Re - x;
  result.Im :=  z.Im;
end;

function Csub(x: RealType; w: complex): complex;
begin
  result.Re :=  x - w.Re;
  result.Im :=    - w.Im;
end;

function Cmul(z: complex;  w: complex): complex;
{ result := z * w }
begin
  result.Re :=  z.Re*w.Re - z.Im*w.Im;
  result.Im :=  z.Im*w.Re + z.Re*w.Im;
end;

function Cmul(x: RealType; w: complex): complex;
begin
  result.Re :=  x*w.Re;
  result.Im :=  x*w.Im;
end;

function Cmul(z: complex;  x: RealType): complex;
begin
  result.Re :=  z.Re*x;
  result.Im :=  z.Im*x;
end;

function Cdiv(z: complex;  w: complex): complex;
{ result := z / w }
var
  rz, rw: RealType;
  zu, wu: complex;
begin
  wu := Cunit(w, rw);
  if rw > MinReal then
  begin
    zu := Cunit(z, rz);
    try
      result := Cmul(Cmul(zu, conjug(wu)),rz/rw);
    except
      on E:Exception do
        ShowMessage(Format('Complex division Error: %s',[E.Message]));
    end;
  end;
end;

function Cdiv(x: RealType; w: complex): complex;
var
  rw: RealType;
  wu: complex;
begin
  wu := Cunit(w, rw);
  if rw > MinReal then
  begin
    try
      result := Cmul(x/rw, conjug(wu));
    except
      on E:Exception do
        ShowMessage(Format('Complex division Error: %s',[E.Message]));
    end;
  end;
end;

function Cdiv(z: complex;  x: RealType): complex;
begin
  try
    result.Re :=  z.Re/x;
    result.Im :=  z.Im/x;
  except
    on E:Exception do
      ShowMessage(Format('Complex division Error: %s',[E.Message]));
  end;
end;

function Cexp(z: complex): complex;
{ z := exp(v)  }
begin
  result.Re := exp1(z.re)*cos(z.im);
  result.Im := exp1(z.re)*sin(z.im);
end;

function Csqrt(z: complex): complex;
{ result := sqrt(z), result.Re > 0 }
begin
  result.Re :=  sqrt(abs(0.5*(Cabs(z) + z.re)));
  result.Im :=  sqrt(abs(0.5*(Cabs(z) - z.re)));
end;

function Croot(x: RealType): complex;
{ result := sqrt(x) }
begin
  if x < 0.0 then result := cmplx(0.0, sqrt(-x))
             else result := cmplx(sqrt(x));
end;

procedure Csum(var s: complex; z: complex);
{ s := s + z }
begin
  s := Cadd(s, z);
end;
procedure Csum(var s: complex; x: RealType);
{ s := s + x }
begin
  s := Cadd(s, x);
end;

{ Special functions for Complex Calculations }

function cer(z: complex): complex;
{ result := er(z) }

  function cerser(z: complex): complex;
  { series expansion of er(z) }
  const
    epsilon = 1.0E-14;
         sp = 1.7724538509055160;{ sqrt(pi) }
  var
    s, t, m, zz, zz2: complex;
  begin
    m := cmplx1;
    t := cmplx1;
    s := cmplx0;
    zz := Cmul(z, z);
    zz2 := Cmul(zz, 2.0);
    repeat
      Csum(s, t);
      Csum(m, 2.0);
      t := Cmul(t, zz2);
      t := Cdiv(t, m);
    until (Cabs0(t) < epsilon);
    s := Cmul(s, z);
    s := Cmul(s, 2.0/sp);
    t := Cexp(zz);
    result := Csub(t, s);
  end;

  function cerabr(z: complex): complex;
  { approximation of er(z): asimptotics from Abramowitz }
  const
    a1 = 0.4613135;    b1 = 0.1901635;
    a2 = 0.09999216;   b2 = 1.7844927;
    a3 = 0.002883894;  b3 = 5.5253437;
  var
    u, v, zz                : complex;
  begin
    zz := Cmul(z, z);
    u := Cadd(zz, b1);    v := Cdiv(a1, u);
    u := Cadd(zz, b2);    u := Cdiv(a2, u);
    Csum(v, u);
    u := Cadd(zz, b3);    u := Cdiv(a3, u);
    Csum(v, u);
    result := Cmul(v, z);
  end;

begin
  if (abs(z.re) > 3.0) or (abs(z.im) > 3.9)
     then result := cerabr(z)
     else result := cerser(z);
end;{ of cer }

end.