Unit PCOMPLEX;

(* Bibliotheque mathematique pour type complexe *)
(* Version a fonctions et pointeurs *)
(* JD GAYRARD mai 95 *)

(* 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 complexe = record
                reel : float;
                imag : float
                end;

pcomplexe = ^complexe;

const _i : complexe = (reel : 0.0; imag : 1.0);
      _0 : complexe = (reel : 0.0; imag : 0.0);

var result : complexe; { all functions result points on this varaible }

(* quatre operations : +, -, * , / *)
function cadd (z1, z2 : complexe) : pcomplexe;      (* addition *)
function csub (z1, z2 : complexe) : pcomplexe;      (* soustraction *)
function cmul (z1, z2 : complexe) : pcomplexe;      (* multiplication *)
function cdiv (znum, zden : complexe) : pcomplexe;  (* division znum / zden *)

(* fonctions complexes particulieres *)
function cneg (z : complexe) : pcomplexe;       (* negatif *)
function ccong (z : complexe) : pcomplexe;      (* conjuge *)
function crcp (z : complexe) : pcomplexe;       (* inverse *)
function ciz (z : complexe) : pcomplexe;        (* multiplication par i *)
function c_iz (z : complexe) : pcomplexe;       (* multiplication par -i *)
function czero : pcomplexe;                     (* return zero *)

(* fonctions complexes a retour non complexe *)
function cmod (z : complexe) : float;           (* module *)
function cequal (z1, z2 : complexe) : boolean;  (* compare deux complexes *)
function carg (z : complexe) : float;           (* argument : a / z = p.e^ia *)

(* fonctions elementaires *)
function cexp (z : complexe) : pcomplexe;       (* exponantielle *)
function cln (z : complexe) : pcomplexe;        (* logarithme naturel *)
function csqrt (z : complexe) : pcomplexe;      (* racine carre *)

(* fonctions trigonometrique directe *)
function ccos (z : complexe) : pcomplexe;       (* cosinus *)
function csin (z : complexe) : pcomplexe;       (* sinus *)
function ctg  (z : complexe) : pcomplexe;       (* tangente *)

(* fonctions trigonometriques inverses *)
function carc_cos (z : complexe) : pcomplexe;   (* arc cosinus *)
function carc_sin (z : complexe) : pcomplexe;   (* arc sinus *)
function carc_tg  (z : complexe) : pcomplexe;   (* arc tangente *)

(* fonctions trigonometrique hyperbolique *)
function cch (z : complexe) : pcomplexe;        (* cosinus hyperbolique *)
function csh (z : complexe) : pcomplexe;        (* sinus hyperbolique *)
function cth (z : complexe) : pcomplexe;        (* tangente hyperbolique *)

(* fonctions trigonometrique hyperbolique inverse *)
function carg_ch (z : complexe) : pcomplexe;    (* arc cosinus hyperbolique *)
function carg_sh (z : complexe) : pcomplexe;    (* arc sinus hyperbolique *)
function carg_th (z : complexe) : pcomplexe;    (* arc tangente hyperbolique *)



implementation

(* quatre operations de base +, -, * , / *)

function cadd (z1, z2 : complexe) : pcomplexe;
(* addition : r := z1 + z2 *)
begin
result.reel := z1.reel + z2.reel;
result.imag := z1.imag + z2.imag;
cadd := @result
end;

function csub (z1, z2 : complexe) : pcomplexe;
(* soustraction : r :=  z1 - z2 *)
begin
result.reel := z1.reel - z2.reel;
result.imag := z1.imag - z2.imag;
csub := @result
end;

function cmul (z1, z2 : complexe) : pcomplexe;
(* multiplication : r := z1 * z2 *)
begin
result.reel := (z1.reel * z2.reel) - (z1.imag * z2.imag);
result.imag := (z1.reel * z2.imag) + (z1.imag * z2.reel);
cmul := @result
end;

function cdiv (znum, zden : complexe) : pcomplexe;
(* division : r := znum / zden *)
var denom : float;
begin
with zden do denom := (reel * reel) + (imag * imag);
if denom = 0.0
   then begin
        writeln('******** function Cdiv ********');
        writeln('******* DIVISION PAR ZERO ******');
        halt
        end
   else begin
        result.reel := ((znum.reel * zden.reel) + (znum.imag * zden.imag)) / denom;
        result.imag := ((znum.imag * zden.reel) - (znum.reel * zden.imag)) / denom
        end;
cdiv := @result
end;

(* fonctions complexes particulieres *)

function cneg (z : complexe) : pcomplexe;
(* negatif : r = - z *)
begin
result.reel := - z.reel;
result.imag := - z.imag;
cneg := @result
end;

function cmod (z : complexe): float;
(* module : r = |z| *)
begin
with z do cmod := sqrt((reel * reel) + (imag * imag))
end;

function carg (z : complexe): float;
(* argument : 0 / z = p ei0 *)
begin
carg := arctan2(z.reel, z.imag)
end;

function ccong (z : complexe) : pcomplexe;
(* conjuge : z := x + i.y alors r = x - i.y *)
begin
result.reel := z.reel;
result.imag := - z.imag;
ccong := @result
end;

function crcp (z : complexe) : pcomplexe;
(* inverse : r := 1 / z *)
var denom : float;
begin
with z do denom := (reel * reel) + (imag * imag);
if denom = 0.0
   then begin
        writeln('******** function Crcp ********');
        writeln('******* DIVISION PAR ZERO ******');
        halt
        end
   else begin
        result.reel := z.reel / denom;
        result.imag := - z.imag / denom
        end;
crcp := @result
end;

function ciz (z : complexe) : pcomplexe;
(* multiplication par i *)
(* z = x + i.y , r = i.z = - y + i.x *)
begin
result.reel := - z.imag;
result.imag := z.reel;
ciz := @result
end;

function c_iz (z : complexe) : pcomplexe;
(* multiplication par -i *)
(* z = x + i.y , r = i.z = y - i.x *)
begin
result.reel := z.imag;
result.imag := - z.reel;
c_iz := @result
end;

function czero : pcomplexe;
(* return a zero complexe *)
begin
result.reel := 0.0;
result.imag := 0.0;
czero := @result
end;

function cequal (z1, z2 : complexe) : boolean;
(* retourne TRUE si z1 = z2 *)
begin
cequal := (z1.reel = z2.reel) and (z1.imag = z2.imag)
end;

(* fonctions elementaires *)

function cexp (z : complexe) : pcomplexe;
(* exponantielle : r := exp(z) *)
(* exp(x + iy) = exp(x).exp(iy) = exp(x).[cos(y) + i sin(y)] *)
var expz : float;
begin
expz := exp(z.reel);
result.reel := expz * cos(z.imag);
result.imag := expz * sin(z.imag);
cexp := @result
end;

function cln (z : complexe) : pcomplexe;
(* logarithme naturel : r := ln(z) *)
(* ln( p exp(i0)) = ln(p) + i0 + 2kpi *)
var modz : float;
begin
with z do modz := (reel * reel) + (imag * imag);
if modz = 0.0
   then begin
        writeln('********* function Cln *********');
        writeln('****** LOGARITHME DE ZERO ******');
        halt
        end
   else begin
   result.reel := ln(modz);
   result.imag := arctan2(z.reel, z.imag);
   cln := @result
        end
end;

function csqrt (z : complexe) : pcomplexe;
(* racine carre : r := sqrt(z) *)
var root, q : float;
begin
if (z.reel <> 0.0) or (z.imag <> 0.0)
   then begin
        root := sqrt(0.5 * (abs(z.reel) + cmod(z)));
        q := z.imag / (2.0 * root);
        if z.reel >= 0.0
           then with result do
                begin
                reel := root;
                imag := q
                end
           else if z.imag < 0.0
                   then with result do
                        begin
                        reel := - q;
                        imag := - root
                        end
                   else with result do
                        begin
                        reel :=  q;
                        imag :=  root
                        end
        end
   else result := z;
csqrt := @result
end;

(* fonctions trigonometriques directes *)

function ccos (z : complexe) : pcomplexe;
(* cosinus complexe *)
(* cos(x+iy) = cos(x).cos(iy) - sin(x).sin(iy) *)
(* cos(ix) = ch(x) et sin(ix) = i.sh(x) *)
begin
result.reel := cos(z.reel) * ch(z.imag);
result.imag := - sin(z.reel) * sh(z.imag);
ccos := @result
end;

function csin (z : complexe) : pcomplexe;
(* sinus complexe *)
(* sin(x+iy) = sin(x).cos(iy) + cos(x).sin(iy) *)
(* cos(ix) = ch(x) et sin(ix) = i.sh(x) *)
begin
result.reel := sin(z.reel) * ch(z.imag);
result.imag := cos(z.reel) * sh(z.imag);
csin := @result
end;

function ctg (z : complexe) : pcomplexe;
(* tangente *)
var ccosz, temp : complexe;
begin
ccosz := ccos(z)^;
if (ccosz.reel = 0.0) and (ccosz.imag = 0.0)
   then begin
        writeln('********* function Ctg *********');
        writeln('******* DIVISION PAR ZERO ******');
        halt
        end
   else begin
        temp := csin(z)^;
        result := cdiv(temp, ccosz)^;
        ctg := @result
        end
end;

(* fonctions trigonometriques inverses *)

function carc_cos (z : complexe) : pcomplexe;
(* arc cosinus complexe *)
(* arccos(z) = -i.argch(z) *)
begin
z := carg_ch(z)^;
result := c_iz(z)^;
carc_cos := @result
end;

function carc_sin (z : complexe) : pcomplexe;
(* arc sinus complexe *)
(* arcsin(z) = -i.argsh(i.z) *)
begin
z := ciz(z)^;
z := carg_sh(z)^;
result := c_iz(z)^;
carc_sin := @result
end;

function carc_tg (z : complexe) : pcomplexe;
(* arc tangente complexe *)
(* arctg(z) = -i.argth(i.z) *)
begin
z := ciz(z)^;
z := carg_th(z)^;
result := c_iz(z)^;
carc_tg := @result
end;

(* fonctions trigonometriques hyperboliques *)

function cch (z : complexe) : pcomplexe;
(* cosinus hyperbolique *)
(* ch(x+iy) = ch(x).ch(iy) + sh(x).sh(iy) *)
(* ch(iy) = cos(y) et sh(iy) = i.sin(y) *)
begin
result.reel := ch(z.reel) * cos(z.imag);
result.imag := sh(z.reel) * sin(z.imag);
cch := @result
end;

function csh (z : complexe) : pcomplexe;
(* sinus hyperbolique *)
(* sh(x+iy) = sh(x).ch(iy) + ch(x).sh(iy) *)
(* ch(iy) = cos(y) et sh(iy) = i.sin(y) *)
begin
result.reel := sh(z.reel) * cos(z.imag);
result.imag := ch(z.reel) * sin(z.imag);
csh := @result
end;

function cth (z : complexe) : pcomplexe;
(* tangente hyperbolique complexe *)
(* th(x) = sh(x) / ch(x) *)
(* ch(x) > 1 qq x *)
var temp : complexe;
begin
temp := cch(z)^;
z := csh(z)^;
result := cdiv(z, temp)^;
cth := @result
end;

(* fonctions trigonometriques hyperboliques inverses *)

function carg_ch (z : complexe) : pcomplexe;
(*   arg cosinus hyperbolique    *)
(*                          _________  *)
(* argch(z) = -/+ ln(z + i.V 1 - z.z)  *)
var temp : complexe;
begin
with temp do begin
             reel := 1 - z.reel * z.reel + z.imag * z.imag;
             imag := - 2 * z.reel * z.imag
             end;
temp := csqrt(temp)^;
temp := ciz(temp)^;
temp := cadd(temp, z)^;
temp := cln(temp)^;
result := cneg(temp)^;
carg_ch := @result
end;

function carg_sh (z : complexe) : pcomplexe;
(*   arc sinus hyperbolique    *)
(*                    ________  *)
(* argsh(z) = ln(z + V 1 + z.z) *)
var temp : complexe;
begin
with temp do begin
             reel := 1 + z.reel * z.reel - z.imag * z.imag;
             imag := 2 * z.reel * z.imag
             end;
temp := csqrt(temp)^;
temp := cadd(temp, z)^;
result := cln(temp)^;
carg_sh := @result
end;

function carg_th (z : complexe) : pcomplexe;
(* arc tangente hyperbolique *)
(* argth(z) = 1/2 ln((z + 1) / (1 - z)) *)
var temp : complexe;
begin
with temp do begin
             reel := 1 + z.reel;
             imag := z.imag
             end;
with result do begin
          reel := 1 - reel;
          imag := - imag
          end;
result := cdiv(temp, result)^;
with result do begin
          reel := 0.5 * reel;
          imag := 0.5 * imag
          end;
carg_th := @result
end;

end.