mirror of
https://gitlab.com/freepascal.org/fpc/source.git
synced 2025-10-31 22:11:12 +01:00
578 lines
15 KiB
ObjectPascal
578 lines
15 KiB
ObjectPascal
{
|
|
$Id$
|
|
This file is part of the Free Pascal run time library.
|
|
Copyright (c) 1993,97 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 COMPLEX;
|
|
|
|
(* Bibliotheque mathematique pour type complexe *)
|
|
(* Version a fonctions et pointeurs *)
|
|
(* JD GAYRARD mai 95 *)
|
|
|
|
(* modified for FPK by Pierre Muller *)
|
|
(* FPK supports operator overloading *)
|
|
(* This library is based on functions instead of procedures.
|
|
To allow a function to return complex type, the trick is
|
|
is to use a pointer on the result of the function. All
|
|
functions are of Pcomplex type (^complexe).
|
|
In the main program the function computation is accessed
|
|
by z := function_name(param1, param2)^ *)
|
|
|
|
{$G+}
|
|
{$N+}
|
|
{$E-}
|
|
|
|
interface
|
|
|
|
uses MATHLIB, HYPERBOL;
|
|
|
|
const author = 'GAYRARD J-D';
|
|
version = 'ver 0.0 - 05/95';
|
|
|
|
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);
|
|
|
|
|
|
(* quatre operations : +, -, * , / and = *)
|
|
operator + (z1, z2 : complex) z : complex; (* addition *)
|
|
{ assignment overloading is not implement }
|
|
{ if it was the next two would not be unnecessary }
|
|
operator + (z1 : complex; r : real) z : complex; (* addition *)
|
|
operator + (r : real; z1 : complex) z : complex; (* addition *)
|
|
|
|
operator - (z1, z2 : complex) z : complex; (* soustraction *)
|
|
operator - (z1 : complex;r : real) z : complex; (* soustraction *)
|
|
operator - (r : real; z1 : complex) z : complex; (* soustraction *)
|
|
|
|
operator * (z1, z2 : complex) z : complex; (* multiplication *)
|
|
operator * (z1 : complex; r : real) z : complex; (* multiplication *)
|
|
operator * (r : real; z1 : complex) z : complex; (* multiplication *)
|
|
|
|
operator / (znum, zden : complex) z : complex; (* division znum / zden *)
|
|
operator / (znum : complex; r : real) z : complex; (* division znum / r *)
|
|
operator / (r : real; zden : complex) z : complex; (* division r / zden *)
|
|
|
|
{ ^ is not built in FPK, but can be overloaded also }
|
|
operator ^ (r1, r2 : real) r : real;
|
|
operator ^ (z1, z2 : complex) z : complex;
|
|
operator ^ (z1 : complex; r2 : real) z : complex;
|
|
operator ^ (r : real; z1 : complex) z : complex;
|
|
|
|
operator = (z1, z2 : complex) b : boolean;
|
|
operator = (z1 : complex;r : real) b : boolean;
|
|
operator = (r : real; z1 : complex) b : boolean;
|
|
|
|
(* fonctions complexes particulieres *)
|
|
function cneg (z : complex) : complex; (* negatif *)
|
|
function ccong (z : complex) : complex; (* conjuge *)
|
|
function crcp (z : complex) : complex; (* inverse *)
|
|
function ciz (z : complex) : complex; (* multiplication par i *)
|
|
function c_iz (z : complex) : complex; (* multiplication par -i *)
|
|
function czero : complex; (* return zero *)
|
|
|
|
(* fonctions complexes a retour non complex *)
|
|
function cmod (z : complex) : real; (* module *)
|
|
function cequal (z1, z2 : complex) : boolean; (* compare deux complexes *)
|
|
function carg (z : complex) : real; (* argument : a / z = p.e^ia *)
|
|
|
|
(* fonctions elementaires *)
|
|
function cexp (z : complex) : complex; (* exponantielle *)
|
|
function cln (z : complex) : complex; (* logarithme naturel *)
|
|
function csqrt (z : complex) : complex; (* racine carre *)
|
|
|
|
(* fonctions trigonometrique directe *)
|
|
function ccos (z : complex) : complex; (* cosinus *)
|
|
function csin (z : complex) : complex; (* sinus *)
|
|
function ctg (z : complex) : complex; (* tangente *)
|
|
|
|
(* fonctions trigonometriques inverses *)
|
|
function carc_cos (z : complex) : complex; (* arc cosinus *)
|
|
function carc_sin (z : complex) : complex; (* arc sinus *)
|
|
function carc_tg (z : complex) : complex; (* arc tangente *)
|
|
|
|
(* fonctions trigonometrique hyperbolique *)
|
|
function cch (z : complex) : complex; (* cosinus hyperbolique *)
|
|
function csh (z : complex) : complex; (* sinus hyperbolique *)
|
|
function cth (z : complex) : complex; (* tangente hyperbolique *)
|
|
|
|
(* fonctions trigonometrique hyperbolique inverse *)
|
|
function carg_ch (z : complex) : complex; (* arc cosinus hyperbolique *)
|
|
function carg_sh (z : complex) : complex; (* arc sinus hyperbolique *)
|
|
function carg_th (z : complex) : complex; (* arc tangente hyperbolique *)
|
|
|
|
|
|
|
|
implementation
|
|
|
|
(* quatre operations de base +, -, * , / *)
|
|
|
|
operator + (z1, z2 : complex) z : complex;
|
|
(* 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 *)
|
|
begin
|
|
z.re := z1.re + r;
|
|
z.im := z1.im;
|
|
end;
|
|
|
|
operator + (r : real; z1 : complex) z : complex;
|
|
(* addition : z := r + z1 *)
|
|
begin
|
|
z.re := z1.re + r;
|
|
z.im := z1.im;
|
|
end;
|
|
|
|
operator - (z1, z2 : complex) z : complex;
|
|
(* 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;
|
|
(* substraction : z := z1 - r *)
|
|
begin
|
|
z.re := z1.re - r;
|
|
z.im := z1.im;
|
|
end;
|
|
|
|
operator - (r : real; z1 : complex) z : complex;
|
|
(* 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 *)
|
|
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;
|
|
(* multiplication : z := z1 * r *)
|
|
begin
|
|
z.re := z1.re * r;
|
|
z.im := z1.im * r;
|
|
end;
|
|
|
|
operator * (r : real; z1 : complex) z : complex;
|
|
(* multiplication : z := r + z1 *)
|
|
begin
|
|
z.re := z1.re * r;
|
|
z.im := z1.im * r;
|
|
end;
|
|
|
|
operator / (znum, zden : complex) z : complex;
|
|
(* division : z := znum / zden *)
|
|
var denom : real;
|
|
begin
|
|
with zden do denom := (re * re) + (im * im);
|
|
if denom = 0.0
|
|
then begin
|
|
writeln('******** function Cdiv ********');
|
|
writeln('******* DIVISION BY ZERO ******');
|
|
runerror(200);
|
|
end
|
|
else begin
|
|
z.re := ((znum.re * zden.re) + (znum.im * zden.im)) / denom;
|
|
z.im := ((znum.im * zden.re) - (znum.re * zden.im)) / denom;
|
|
end;
|
|
end;
|
|
|
|
operator / (znum : complex; r : real) z : complex;
|
|
(* division : z := znum / r *)
|
|
begin
|
|
if r = 0.0
|
|
then begin
|
|
writeln('******** function Cdiv ********');
|
|
writeln('******* DIVISION BY ZERO ******');
|
|
runerror(200);
|
|
end
|
|
else begin
|
|
z.re := znum.re / r;
|
|
z.im := znum.im / r;
|
|
end;
|
|
end;
|
|
|
|
operator / (r : real; zden : complex) z : complex;
|
|
(* division : z := r / zden *)
|
|
var denom : real;
|
|
begin
|
|
with zden do denom := (re * re) + (im * im);
|
|
if denom = 0.0
|
|
then begin
|
|
writeln('******** function Cdiv ********');
|
|
writeln('******* DIVISION BY ZERO ******');
|
|
runerror(200);
|
|
end
|
|
else begin
|
|
z.re := (r * zden.re) / denom;
|
|
z.im := - (r * zden.im) / denom;
|
|
end;
|
|
end;
|
|
|
|
(* fonctions complexes particulieres *)
|
|
|
|
function cneg (z : complex) : complex;
|
|
(* negatif : cneg = - z *)
|
|
begin
|
|
cneg.re := - z.re;
|
|
cneg.im := - z.im;
|
|
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 ccong (z : complex) : complex;
|
|
(* conjuge : z := x + i.y alors r = x - i.y *)
|
|
begin
|
|
ccong.re := z.re;
|
|
ccong.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);
|
|
if denom = 0.0
|
|
then begin
|
|
writeln('******** function Cinv ********');
|
|
writeln('******* DIVISION BY ZERO ******');
|
|
runerror(200);
|
|
end
|
|
else begin
|
|
cinv.re := z.re / denom;
|
|
cinv.im := - z.im / denom
|
|
end;
|
|
end;
|
|
|
|
function ciz (z : complex) : complex;
|
|
(* multiplication par i *)
|
|
(* z = x + i.y , r = i.z = - y + i.x *)
|
|
begin
|
|
ciz.re := - z.im;
|
|
ciz.im := z.re;
|
|
end;
|
|
|
|
function c_iz (z : complex) : complex;
|
|
(* multiplication par -i *)
|
|
(* z = x + i.y , r = -i.z = y - i.x *)
|
|
begin
|
|
c_iz.re := z.im;
|
|
c_iz.im := - z.re;
|
|
end;
|
|
|
|
function czero : complex;
|
|
(* return a zero complex *)
|
|
begin
|
|
czero.re := 0.0;
|
|
czero.im := 0.0;
|
|
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;
|
|
(* exponantielle : 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;
|
|
(* logarithme naturel : r := ln(z) *)
|
|
(* ln( p exp(i0)) = ln(p) + i0 + 2kpi *)
|
|
var modz : real;
|
|
begin
|
|
with z do modz := (re * re) + (im * im);
|
|
if modz = 0.0
|
|
then begin
|
|
writeln('********* function Cln *********');
|
|
writeln('****** LOGARITHME OF ZERO ******');
|
|
runerror(200)
|
|
end
|
|
else begin
|
|
cln.re := ln(modz);
|
|
cln.im := arctan2(z.re, z.im);
|
|
end
|
|
end;
|
|
|
|
function csqrt (z : complex) : complex;
|
|
(* racine carre : 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 result := z;
|
|
end;
|
|
|
|
operator ^ (z1, z2 : complex) z : complex;
|
|
(* multiplication : z := z1 * z2 *)
|
|
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;
|
|
(* multiplication : z := z1 * r *)
|
|
begin
|
|
z.re := z1.re * r;
|
|
z.im := z1.im * r;
|
|
end;
|
|
|
|
operator * (r : real; z1 : complex) z : complex;
|
|
(* multiplication : z := r + z1 *)
|
|
begin
|
|
z.re := z1.re * r;
|
|
z.im := z1.im * r;
|
|
end;
|
|
|
|
(* fonctions trigonometriques directes *)
|
|
|
|
function ccos (z : complex) : complex;
|
|
(* cosinus complex *)
|
|
(* cos(x+iy) = cos(x).cos(iy) - sin(x).sin(iy) *)
|
|
(* cos(ix) = ch(x) et sin(ix) = i.sh(x) *)
|
|
begin
|
|
ccos.re := cos(z.re) * ch(z.im);
|
|
ccos.im := - sin(z.re) * sh(z.im);
|
|
end;
|
|
|
|
function csin (z : complex) : complex;
|
|
(* sinus complex *)
|
|
(* sin(x+iy) = sin(x).cos(iy) + cos(x).sin(iy) *)
|
|
(* cos(ix) = ch(x) et sin(ix) = i.sh(x) *)
|
|
begin
|
|
csin.re := sin(z.re) * ch(z.im);
|
|
csin.im := cos(z.re) * sh(z.im);
|
|
end;
|
|
|
|
function ctg (z : complex) : complex;
|
|
(* tangente *)
|
|
var ccosz, temp : complex;
|
|
begin
|
|
ccosz := ccos(z);
|
|
if (ccosz.re = 0.0) and (ccosz.im = 0.0)
|
|
then begin
|
|
writeln('********* function Ctg *********');
|
|
writeln('******* DIVISION BY ZERO ******');
|
|
runerror(200)
|
|
end
|
|
else begin
|
|
temp := csin(z);
|
|
ctg := temp / ccosz;
|
|
end
|
|
end;
|
|
|
|
(* fonctions trigonometriques inverses *)
|
|
|
|
function carc_cos (z : complex) : complex;
|
|
(* arc cosinus complex *)
|
|
(* arccos(z) = -i.argch(z) *)
|
|
begin
|
|
z := carg_ch(z);
|
|
carc_cos := c_iz(z);
|
|
end;
|
|
|
|
function carc_sin (z : complex) : complex;
|
|
(* arc sinus complex *)
|
|
(* arcsin(z) = -i.argsh(i.z) *)
|
|
begin
|
|
z := ciz(z);
|
|
z := carg_sh(z);
|
|
carc_sin := c_iz(z);
|
|
end;
|
|
|
|
function carc_tg (z : complex) : complex;
|
|
(* arc tangente complex *)
|
|
(* arctg(z) = -i.argth(i.z) *)
|
|
begin
|
|
z := ciz(z);
|
|
z := carg_th(z);
|
|
carc_tg := c_iz(z);
|
|
end;
|
|
|
|
(* fonctions trigonometriques hyperboliques *)
|
|
|
|
function cch (z : complex) : complex;
|
|
(* cosinus hyperbolique *)
|
|
(* ch(x+iy) = ch(x).ch(iy) + sh(x).sh(iy) *)
|
|
(* ch(iy) = cos(y) et sh(iy) = i.sin(y) *)
|
|
begin
|
|
cch.re := ch(z.re) * cos(z.im);
|
|
cch.im := sh(z.re) * sin(z.im);
|
|
end;
|
|
|
|
function csh (z : complex) : complex;
|
|
(* sinus hyperbolique *)
|
|
(* sh(x+iy) = sh(x).ch(iy) + ch(x).sh(iy) *)
|
|
(* ch(iy) = cos(y) et sh(iy) = i.sin(y) *)
|
|
begin
|
|
csh.re := sh(z.re) * cos(z.im);
|
|
csh.im := ch(z.re) * sin(z.im);
|
|
end;
|
|
|
|
function cth (z : complex) : complex;
|
|
(* tangente hyperbolique complex *)
|
|
(* th(x) = sh(x) / ch(x) *)
|
|
(* ch(x) > 1 qq x *)
|
|
var temp : complex;
|
|
begin
|
|
temp := cch(z);
|
|
z := csh(z);
|
|
cth := z / temp;
|
|
end;
|
|
|
|
(* fonctions trigonometriques hyperboliques inverses *)
|
|
|
|
function carg_ch (z : complex) : complex;
|
|
(* arg cosinus hyperbolique *)
|
|
(* _________ *)
|
|
(* argch(z) = -/+ ln(z + i.V 1 - z.z) *)
|
|
var temp : complex;
|
|
begin
|
|
with temp do begin
|
|
re := 1 - z.re * z.re + z.im * z.im;
|
|
im := - 2 * z.re * z.im
|
|
end;
|
|
temp := csqrt(temp);
|
|
temp := ciz(temp);
|
|
temp := temp + z;
|
|
temp := cln(temp);
|
|
carg_ch := cneg(temp);
|
|
end;
|
|
|
|
function carg_sh (z : complex) : complex;
|
|
(* arc sinus hyperbolique *)
|
|
(* ________ *)
|
|
(* argsh(z) = ln(z + V 1 + z.z) *)
|
|
var temp : complex;
|
|
begin
|
|
with temp do begin
|
|
re := 1 + z.re * z.re - z.im * z.im;
|
|
im := 2 * z.re * z.im
|
|
end;
|
|
temp := csqrt(temp);
|
|
temp := cadd(temp, z);
|
|
carg_sh := cln(temp);
|
|
end;
|
|
|
|
function carg_th (z : complex) : complex;
|
|
(* arc tangente hyperbolique *)
|
|
(* argth(z) = 1/2 ln((z + 1) / (1 - z)) *)
|
|
var temp : complex;
|
|
begin
|
|
with temp do begin
|
|
re := 1 + z.re;
|
|
im := z.im
|
|
end;
|
|
with result do begin
|
|
re := 1 - re;
|
|
im := - im
|
|
end;
|
|
result := temp / result;
|
|
carg_th.re := 0.5 * re;
|
|
carg_th.im := 0.5 * im
|
|
end;
|
|
|
|
end.
|
|
{
|
|
$Log$
|
|
Revision 1.1 1998-03-25 11:18:43 root
|
|
Initial revision
|
|
|
|
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
|
|
=============================================================================
|
|
}
|