fpc/rtl/inc/ucomplex.pp
2000-01-07 16:41:28 +00:00

655 lines
16 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.4 2000-01-07 16:41:37 daniel
* copyright 2000
Revision 1.3 2000/01/07 16:32:25 daniel
* copyright 2000 added
Revision 1.2 1999/12/20 22:24:48 pierre
+ cinv in interface
Revision 1.1 1998/06/15 15:45:42 pierre
+ complex.pp replaced by ucomplex.pp
complex operations working
Revision 1.1.1.1 1998/03/25 11:18:43 root
* Restored version
Revision 1.3 1998/01/26 11:59:25 michael
+ Added log at the end
Working file: rtl/inc/complex.pp
description:
----------------------------
revision 1.2
date: 1997/12/01 15:33:30; author: michael; state: Exp; lines: +14 -0
+ added copyright reference in header.
----------------------------
revision 1.1
date: 1997/11/27 08:33:46; author: michael; state: Exp;
Initial revision
----------------------------
revision 1.1.1.1
date: 1997/11/27 08:33:46; author: michael; state: Exp; lines: +0 -0
FPC RTL CVS start
=============================================================================
}