mirror of
https://gitlab.com/freepascal.org/fpc/source.git
synced 2025-10-31 20:31:51 +01:00
622 lines
15 KiB
ObjectPascal
622 lines
15 KiB
ObjectPascal
{
|
|
$Id$
|
|
This file is part of the Free Pascal run time library.
|
|
Copyright (c) 1999-2000 by Pierre Muller,
|
|
member of the Free Pascal development team.
|
|
|
|
See the file COPYING.FPC, included in this distribution,
|
|
for details about the copyright.
|
|
|
|
This program is distributed in the hope that it will be useful,
|
|
but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
|
|
|
|
**********************************************************************}
|
|
Unit UComplex;
|
|
|
|
{ created for FPC by Pierre Muller }
|
|
{ inpired from the complex unit from JD GAYRARD mai 95 }
|
|
{ FPC supports operator overloading }
|
|
|
|
|
|
interface
|
|
|
|
uses math;
|
|
|
|
type complex = record
|
|
re : real;
|
|
im : real;
|
|
end;
|
|
|
|
pcomplex = ^complex;
|
|
|
|
const i : complex = (re : 0.0; im : 1.0);
|
|
_0 : complex = (re : 0.0; im : 0.0);
|
|
|
|
|
|
{ assignment overloading is also used in type conversions
|
|
(beware also in implicit type conversions)
|
|
after this operator any real can be passed to a function
|
|
as a complex arg !! }
|
|
|
|
operator := (r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
{ operator := (i : longint) z : complex;
|
|
not needed because longint can be converted to real }
|
|
|
|
|
|
{ four operator : +, -, * , / and comparison = }
|
|
operator + (z1, z2 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
{ these ones are created because the code
|
|
is simpler and thus faster }
|
|
operator + (z1 : complex; r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator + (r : real; z1 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
|
|
operator - (z1, z2 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator - (z1 : complex;r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator - (r : real; z1 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
|
|
operator * (z1, z2 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator * (z1 : complex; r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator * (r : real; z1 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
|
|
operator / (znum, zden : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator / (znum : complex; r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator / (r : real; zden : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
{ ** is the exponentiation operator }
|
|
operator ** (z1, z2 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator ** (z1 : complex; r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator ** (r : real; z1 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
|
|
operator = (z1, z2 : complex) b : boolean;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator = (z1 : complex;r : real) b : boolean;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator = (r : real; z1 : complex) b : boolean;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
operator - (z1 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
|
|
{ complex functions }
|
|
function cong (z : complex) : complex; { conjuge }
|
|
|
|
{ inverse function 1/z }
|
|
function cinv (z : complex) : complex;
|
|
|
|
{ complex functions with real return values }
|
|
function cmod (z : complex) : real; { module }
|
|
function carg (z : complex) : real; { argument : a / z = p.e^ia }
|
|
|
|
{ fonctions elementaires }
|
|
function cexp (z : complex) : complex; { exponential }
|
|
function cln (z : complex) : complex; { natural logarithm }
|
|
function csqrt (z : complex) : complex; { square root }
|
|
|
|
{ complex trigonometric functions }
|
|
function ccos (z : complex) : complex; { cosinus }
|
|
function csin (z : complex) : complex; { sinus }
|
|
function ctg (z : complex) : complex; { tangent }
|
|
|
|
{ inverse complex trigonometric functions }
|
|
function carc_cos (z : complex) : complex; { arc cosinus }
|
|
function carc_sin (z : complex) : complex; { arc sinus }
|
|
function carc_tg (z : complex) : complex; { arc tangent }
|
|
|
|
{ hyperbolic complex functions }
|
|
function cch (z : complex) : complex; { hyperbolic cosinus }
|
|
function csh (z : complex) : complex; { hyperbolic sinus }
|
|
function cth (z : complex) : complex; { hyperbolic tangent }
|
|
|
|
{ inverse hyperbolic complex functions }
|
|
function carg_ch (z : complex) : complex; { hyperbolic arc cosinus }
|
|
function carg_sh (z : complex) : complex; { hyperbolic arc sinus }
|
|
function carg_th (z : complex) : complex; { hyperbolic arc tangente }
|
|
|
|
{ functions to write out a complex value }
|
|
function cstr(z : complex) : string;
|
|
function cstr(z:complex;len : integer) : string;
|
|
function cstr(z:complex;len,dec : integer) : string;
|
|
|
|
implementation
|
|
|
|
operator := (r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
begin
|
|
z.re:=r;
|
|
z.im:=0.0;
|
|
end;
|
|
|
|
{ four base operations +, -, * , / }
|
|
|
|
operator + (z1, z2 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
{ addition : z := z1 + z2 }
|
|
begin
|
|
z.re := z1.re + z2.re;
|
|
z.im := z1.im + z2.im;
|
|
end;
|
|
|
|
operator + (z1 : complex; r : real) z : complex;
|
|
{ addition : z := z1 + r }
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
begin
|
|
z.re := z1.re + r;
|
|
z.im := z1.im;
|
|
end;
|
|
|
|
operator + (r : real; z1 : complex) z : complex;
|
|
{ addition : z := r + z1 }
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
|
|
begin
|
|
z.re := z1.re + r;
|
|
z.im := z1.im;
|
|
end;
|
|
|
|
operator - (z1, z2 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
{ substraction : z := z1 - z2 }
|
|
begin
|
|
z.re := z1.re - z2.re;
|
|
z.im := z1.im - z2.im;
|
|
end;
|
|
|
|
operator - (z1 : complex; r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
{ substraction : z := z1 - r }
|
|
begin
|
|
z.re := z1.re - r;
|
|
z.im := z1.im;
|
|
end;
|
|
|
|
operator - (z1 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
{ substraction : z := - z1 }
|
|
begin
|
|
z.re := -z1.re;
|
|
z.im := -z1.im;
|
|
end;
|
|
|
|
operator - (r : real; z1 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
{ substraction : z := r - z1 }
|
|
begin
|
|
z.re := r - z1.re;
|
|
z.im := - z1.im;
|
|
end;
|
|
|
|
operator * (z1, z2 : complex) z : complex;
|
|
{ multiplication : z := z1 * z2 }
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
begin
|
|
z.re := (z1.re * z2.re) - (z1.im * z2.im);
|
|
z.im := (z1.re * z2.im) + (z1.im * z2.re);
|
|
end;
|
|
|
|
operator * (z1 : complex; r : real) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
{ multiplication : z := z1 * r }
|
|
begin
|
|
z.re := z1.re * r;
|
|
z.im := z1.im * r;
|
|
end;
|
|
|
|
operator * (r : real; z1 : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
{ multiplication : z := r * z1 }
|
|
begin
|
|
z.re := z1.re * r;
|
|
z.im := z1.im * r;
|
|
end;
|
|
|
|
operator / (znum, zden : complex) z : complex;
|
|
{$ifdef TEST_INLINE}
|
|
inline;
|
|
{$endif TEST_INLINE}
|
|
{ division : z := znum / zden }
|
|
var
|
|
denom : real;
|
|
begin
|
|
with zden do denom := (re * re) + (im * im);
|
|
{ generates a fpu exception if denom=0 as for reals }
|
|
z.re := ((znum.re * zden.re) + (znum.im * zden.im)) / denom;
|
|
z.im := ((znum.im * zden.re) - (znum.re * zden.im)) / denom;
|
|
end;
|
|
|
|
operator / (znum : complex; r : real) z : complex;
|
|
{ division : z := znum / r }
|
|
begin
|
|
z.re := znum.re / r;
|
|
z.im := znum.im / r;
|
|
end;
|
|
|
|
operator / (r : real; zden : complex) z : complex;
|
|
{ division : z := r / zden }
|
|
var denom : real;
|
|
begin
|
|
with zden do denom := (re * re) + (im * im);
|
|
{ generates a fpu exception if denom=0 as for reals }
|
|
z.re := (r * zden.re) / denom;
|
|
z.im := - (r * zden.im) / denom;
|
|
end;
|
|
|
|
function cmod (z : complex): real;
|
|
{ module : r = |z| }
|
|
begin
|
|
with z do
|
|
cmod := sqrt((re * re) + (im * im));
|
|
end;
|
|
|
|
function carg (z : complex): real;
|
|
{ argument : 0 / z = p ei0 }
|
|
begin
|
|
carg := arctan2(z.re, z.im);
|
|
end;
|
|
|
|
function cong (z : complex) : complex;
|
|
{ complex conjugee :
|
|
if z := x + i.y
|
|
then cong is x - i.y }
|
|
begin
|
|
cong.re := z.re;
|
|
cong.im := - z.im;
|
|
end;
|
|
|
|
function cinv (z : complex) : complex;
|
|
{ inverse : r := 1 / z }
|
|
var
|
|
denom : real;
|
|
begin
|
|
with z do denom := (re * re) + (im * im);
|
|
{ generates a fpu exception if denom=0 as for reals }
|
|
cinv.re:=z.re/denom;
|
|
cinv.im:=-z.im/denom;
|
|
end;
|
|
|
|
operator = (z1, z2 : complex) b : boolean;
|
|
{ returns TRUE if z1 = z2 }
|
|
begin
|
|
b := (z1.re = z2.re) and (z1.im = z2.im);
|
|
end;
|
|
|
|
operator = (z1 : complex; r :real) b : boolean;
|
|
{ returns TRUE if z1 = r }
|
|
begin
|
|
b := (z1.re = r) and (z1.im = 0.0)
|
|
end;
|
|
|
|
operator = (r : real; z1 : complex) b : boolean;
|
|
{ returns TRUE if z1 = r }
|
|
begin
|
|
b := (z1.re = r) and (z1.im = 0.0)
|
|
end;
|
|
|
|
|
|
{ fonctions elementaires }
|
|
|
|
function cexp (z : complex) : complex;
|
|
{ exponantial : r := exp(z) }
|
|
{ exp(x + iy) = exp(x).exp(iy) = exp(x).[cos(y) + i sin(y)] }
|
|
var expz : real;
|
|
begin
|
|
expz := exp(z.re);
|
|
cexp.re := expz * cos(z.im);
|
|
cexp.im := expz * sin(z.im);
|
|
end;
|
|
|
|
function cln (z : complex) : complex;
|
|
{ natural logarithm : r := ln(z) }
|
|
{ ln( p exp(i0)) = ln(p) + i0 + 2kpi }
|
|
var modz : real;
|
|
begin
|
|
with z do
|
|
modz := (re * re) + (im * im);
|
|
cln.re := ln(modz);
|
|
cln.im := arctan2(z.re, z.im);
|
|
end;
|
|
|
|
function csqrt (z : complex) : complex;
|
|
{ square root : r := sqrt(z) }
|
|
var
|
|
root, q : real;
|
|
begin
|
|
if (z.re<>0.0) or (z.im<>0.0) then
|
|
begin
|
|
root := sqrt(0.5 * (abs(z.re) + cmod(z)));
|
|
q := z.im / (2.0 * root);
|
|
if z.re >= 0.0 then
|
|
begin
|
|
csqrt.re := root;
|
|
csqrt.im := q;
|
|
end
|
|
else if z.im < 0.0 then
|
|
begin
|
|
csqrt.re := - q;
|
|
csqrt.im := - root
|
|
end
|
|
else
|
|
begin
|
|
csqrt.re := q;
|
|
csqrt.im := root
|
|
end
|
|
end
|
|
else csqrt := z;
|
|
end;
|
|
|
|
|
|
operator ** (z1, z2 : complex) z : complex;
|
|
{ exp : z := z1 ** z2 }
|
|
begin
|
|
z := cexp(z2*cln(z1));
|
|
end;
|
|
|
|
operator ** (z1 : complex; r : real) z : complex;
|
|
{ multiplication : z := z1 * r }
|
|
begin
|
|
z := cexp( r *cln(z1));
|
|
end;
|
|
|
|
operator ** (r : real; z1 : complex) z : complex;
|
|
{ multiplication : z := r + z1 }
|
|
begin
|
|
z := cexp(z1*ln(r));
|
|
end;
|
|
|
|
{ direct trigonometric functions }
|
|
|
|
function ccos (z : complex) : complex;
|
|
{ complex cosinus }
|
|
{ cos(x+iy) = cos(x).cos(iy) - sin(x).sin(iy) }
|
|
{ cos(ix) = cosh(x) et sin(ix) = i.sinh(x) }
|
|
begin
|
|
ccos.re := cos(z.re) * cosh(z.im);
|
|
ccos.im := - sin(z.re) * sinh(z.im);
|
|
end;
|
|
|
|
function csin (z : complex) : complex;
|
|
{ sinus complex }
|
|
{ sin(x+iy) = sin(x).cos(iy) + cos(x).sin(iy) }
|
|
{ cos(ix) = cosh(x) et sin(ix) = i.sinh(x) }
|
|
begin
|
|
csin.re := sin(z.re) * cosh(z.im);
|
|
csin.im := cos(z.re) * sinh(z.im);
|
|
end;
|
|
|
|
function ctg (z : complex) : complex;
|
|
{ tangente }
|
|
var ccosz, temp : complex;
|
|
begin
|
|
ccosz := ccos(z);
|
|
temp := csin(z);
|
|
ctg := temp / ccosz;
|
|
end;
|
|
|
|
{ fonctions trigonometriques inverses }
|
|
|
|
function carc_cos (z : complex) : complex;
|
|
{ arc cosinus complex }
|
|
{ arccos(z) = -i.argch(z) }
|
|
begin
|
|
carc_cos := -i*carg_ch(z);
|
|
end;
|
|
|
|
function carc_sin (z : complex) : complex;
|
|
{ arc sinus complex }
|
|
{ arcsin(z) = -i.argsh(i.z) }
|
|
begin
|
|
carc_sin := -i*carg_sh(i*z);
|
|
end;
|
|
|
|
function carc_tg (z : complex) : complex;
|
|
{ arc tangente complex }
|
|
{ arctg(z) = -i.argth(i.z) }
|
|
begin
|
|
carc_tg := -i*carg_th(i*z);
|
|
end;
|
|
|
|
{ hyberbolic complex functions }
|
|
|
|
function cch (z : complex) : complex;
|
|
{ hyberbolic cosinus }
|
|
{ cosh(x+iy) = cosh(x).cosh(iy) + sinh(x).sinh(iy) }
|
|
{ cosh(iy) = cos(y) et sinh(iy) = i.sin(y) }
|
|
begin
|
|
cch.re := cosh(z.re) * cos(z.im);
|
|
cch.im := sinh(z.re) * sin(z.im);
|
|
end;
|
|
|
|
function csh (z : complex) : complex;
|
|
{ hyberbolic sinus }
|
|
{ sinh(x+iy) = sinh(x).cosh(iy) + cosh(x).sinh(iy) }
|
|
{ cosh(iy) = cos(y) et sinh(iy) = i.sin(y) }
|
|
begin
|
|
csh.re := sinh(z.re) * cos(z.im);
|
|
csh.im := cosh(z.re) * sin(z.im);
|
|
end;
|
|
|
|
function cth (z : complex) : complex;
|
|
{ hyberbolic complex tangent }
|
|
{ th(x) = sinh(x) / cosh(x) }
|
|
{ cosh(x) > 1 qq x }
|
|
var temp : complex;
|
|
begin
|
|
temp := cch(z);
|
|
z := csh(z);
|
|
cth := z / temp;
|
|
end;
|
|
|
|
{ inverse complex hyperbolic functions }
|
|
|
|
function carg_ch (z : complex) : complex;
|
|
{ hyberbolic arg cosinus }
|
|
{ _________ }
|
|
{ argch(z) = -/+ ln(z + i.V 1 - z.z) }
|
|
begin
|
|
carg_ch:=-cln(z+i*csqrt(z*z-1.0));
|
|
end;
|
|
|
|
function carg_sh (z : complex) : complex;
|
|
{ hyperbolic arc sinus }
|
|
{ ________ }
|
|
{ argsh(z) = ln(z + V 1 + z.z) }
|
|
begin
|
|
carg_sh:=cln(z+csqrt(z*z+1.0));
|
|
end;
|
|
|
|
function carg_th (z : complex) : complex;
|
|
{ hyperbolic arc tangent }
|
|
{ argth(z) = 1/2 ln((z + 1) / (1 - z)) }
|
|
begin
|
|
carg_th:=cln((z+1.0)/(z-1.0))/2.0;
|
|
end;
|
|
|
|
{ functions to write out a complex value }
|
|
function cstr(z : complex) : string;
|
|
var
|
|
istr : string;
|
|
begin
|
|
str(z.im,istr);
|
|
str(z.re,cstr);
|
|
while istr[1]=' ' do
|
|
delete(istr,1,1);
|
|
if z.im<0 then
|
|
cstr:=cstr+istr+'i'
|
|
else if z.im>0 then
|
|
cstr:=cstr+'+'+istr+'i';
|
|
end;
|
|
|
|
function cstr(z:complex;len : integer) : string;
|
|
var
|
|
istr : string;
|
|
begin
|
|
str(z.im:len,istr);
|
|
while istr[1]=' ' do
|
|
delete(istr,1,1);
|
|
str(z.re:len,cstr);
|
|
if z.im<0 then
|
|
cstr:=cstr+istr+'i'
|
|
else if z.im>0 then
|
|
cstr:=cstr+'+'+istr+'i';
|
|
end;
|
|
|
|
function cstr(z:complex;len,dec : integer) : string;
|
|
var
|
|
istr : string;
|
|
begin
|
|
str(z.im:len:dec,istr);
|
|
while istr[1]=' ' do
|
|
delete(istr,1,1);
|
|
str(z.re:len:dec,cstr);
|
|
if z.im<0 then
|
|
cstr:=cstr+istr+'i'
|
|
else if z.im>0 then
|
|
cstr:=cstr+'+'+istr+'i';
|
|
end;
|
|
|
|
|
|
end.
|
|
{
|
|
$Log$
|
|
Revision 1.2 2000-07-13 11:33:46 michael
|
|
+ removed logs
|
|
|
|
}
|