(****************************************************************)
(*                                                              *)
(*       TMT Pascal 3 Runtime Library                           *)
(*       Complex Mathematics Unit                               *)
(*                                                              *)
(*       Copyright (c) 1995,2000 TMT Development Corporation    *)
(*       Portions Copyright 1990, by J. W. Rider (CMATH unit)   *)
(*       Authors: Anton Moscal, Vadim Bodrov                    *)
(*                                                              *)
(****************************************************************)

unit Comp;

interface

type
  CReal = Extended;

type Complex = record
  re, im: CReal
end;
TComplex = Complex;

function Compl_RR(const re, im: CReal): Complex;
overload Complex = Compl_RR;
function Compl_R(const re: CReal): Complex;
overload Complex = Compl_R;
function Conj(const z: Complex): Complex; // complex conjugate
function Cjx(const z: Complex): Complex;  // result := j*z, j = sqrt(-1)
function Cnjx(const z: Complex): Complex; // result :=-j*z, j = sqrt(-1)

function add_cc(const a: Complex; const b: Complex): Complex;
overload + = add_cc;
function add_cr(const a: Complex; const b: CReal): Complex;
overload + = add_cr;
function add_rc(const a: CReal; const b: Complex): Complex;
overload + = add_rc;
function sub_cc(const a: Complex; const b: Complex): Complex;
overload - = sub_cc;
function sub_cr(const a: Complex; const b: CReal): Complex;
overload - = sub_cr;
function sub_rc(const a: CReal; const b: Complex): Complex;
overload - = sub_rc;
function mul_cc(const a: Complex; const b: Complex): Complex;
overload * = mul_cc;
function mul_cr(const a: Complex; const b: CReal): Complex;
overload * = mul_cr;
function mul_rc(const a: CReal; const b: Complex): Complex;
overload * = mul_rc;
function div_cc(const a: Complex; const b: Complex): Complex;
overload / = div_cc;
function div_cr(const a: Complex; const b: CReal): Complex;
overload / = div_cr;
function div_rc(const a: CReal; const b: Complex): Complex;
overload / = div_rc;

procedure addab_cc(var a: Complex; const b: Complex);
overload +:= = addab_cc;
procedure addab_cr(var a: Complex; const b: CReal);
overload +:= = addab_cr;
procedure subab_cc(var a: Complex; const b: Complex);
overload -:= = subab_cc;
procedure subab_cr(var a: Complex; const b: CReal);
overload -:= = subab_cr;
procedure mulab_cc(var a: Complex; const b: Complex);
overload *:= = mulab_cc;
procedure mulab_cr(var a: Complex; const b: CReal);
overload *:= = mulab_cr;
procedure divab_cc(var a: Complex; const b: Complex);
overload /:= = divab_cc;
procedure divab_cr(var a: Complex; const b: CReal);
overload /:= = divab_cr;

function eq_cc(const a: Complex; const b: Complex): Boolean;
overload = = eq_cc;
function eq_cr(const a: Complex; const b: CReal): Boolean;
overload = = eq_cr;
function eq_rc(const a: CReal; const b: Complex): Boolean;
overload = = eq_rc;
function ne_cc(const a: Complex; const b: Complex): Boolean;
overload <> = ne_cc;
function ne_cr(const a: Complex; const b: CReal): Boolean;
overload <> = ne_cr;
function ne_rc(const a: CReal; const b: Complex): Boolean;
overload <> = ne_rc;

function sqr_cc(const a: Complex): Complex;
overload Sqr = sqr_cc;
function sqrt_cc(const a: Complex): Complex;
overload Sqrt = sqrt_cc;
function ln_cc(const a: Complex): Complex;
overload Ln = ln_cc;
function exp_cc(const a: Complex): Complex;
overload Exp = exp_cc;
function power_cc(const a: Complex; const b: Complex): Complex;
overload Power = power_cc;
function sin_cc(const a: Complex): Complex;
overload Sin = sin_cc;
function sinh_cc(const a: Complex): Complex;
overload SinH = sinh_cc;
function cos_cc(const a: Complex): Complex;
overload Cos = cos_cc;
function cosh_cc(const a: Complex): Complex;
overload CosH = cosh_cc;
function tan_cc(const a: Complex): Complex;
overload Tan = tan_cc;
function tanh_cc(const a: Complex): Complex;
overload TanH = tanh_cc;
function cotan_cc(const a: Complex): Complex;
overload Cotan = cotan_cc;
function cotanh_cc(const a: Complex): Complex;
overload CotanH = cotanh_cc;
function arcsin_cc(const a: Complex): Complex;
overload ArcSin = arcsin_cc;
function arcsinh_cc(const a: Complex): Complex;
overload ArcSinH = arcsinh_cc;
function arccos_cc(const a: Complex): Complex;
overload ArcCos = arccos_cc;
function arccosh_cc(const a: Complex): Complex;
overload ArcCosH = arccosh_cc;
function arctan_cc(const a: Complex): Complex;
overload ArcTan = arctan_cc;
function arctanh_cc(const a: Complex): Complex;
overload ArcTanH = arctanh_cc;

implementation

uses Math;

function Compl_RR;
begin
  Result.re := re;
  Result.im := im
end;

function Compl_R;
begin
  Result.re := re;
  Result.im := 0
end;

function Conj;
with Result do
begin
  re := z.re;
  im :=-z.im
end;

function Cjx;
var
  tmp: CReal;
begin
  tmp := -z.im;
  with Result do
  begin
    im := z.re;
    re := tmp;
  end;
end;

function Cnjx;
var
  tmp: CReal;
begin
  tmp := z.im;
  with Result do
  begin
    im := -z.re;
    re := tmp;
  end;
end;

function add_cc;
with Result do
begin
  re := b.re + a.re;
  im := a.im + b.im
end;

function add_cr;
with Result do
begin
  re := b + a.re;
  im := a.im;
end;

function add_rc;
with Result do
begin
  re := a + b.re;
  im := b.im;
end;

function sub_cc;
with Result do
begin
  re := a.re - b.re;
  im := a.im - b.im
end;

function sub_cr;
with Result do
begin
  re := a.re - b;
  im := a.im;
end;

function sub_rc;
with Result do
begin
  re := a - b.re;
  im :=   - b.im;
end;

function mul_cc;
with Result do
begin
  re := a.re * b.re - a.im * b.im;
  im := a.re * b.im + a.im * b.re
end;

function mul_cr;
with Result do
begin
  re := b * a.re;
  im := b * a.im;
end;

function mul_rc;
with Result do
begin
  re := a * b.re;
  im := a * b.im;
end;

function div_cc;
var
  x: CReal;
begin
  x := sqr(b.re) + sqr(b.im);
  with Result do begin
    re := (a.re * b.re + a.im * b.im) / x;
    im := (a.im * b.re - a.re * b.im) / x;
  end;
end;

function div_cr;
with Result do
begin
  re := b / a.re;
  im := b / a.im;
end;

function div_rc;
with Result do
begin
  re := a / b.re;
  im := a / b.im;
end;

procedure addab_cc;
with a do begin
  re +:= b.re;
  im +:= b.im
end;

procedure addab_cr;
begin
  a.re +:= b;
end;

procedure subab_cc;
with a do begin
  re -:= b.re;
  im -:= b.im
end;

procedure subab_cr;
begin
  a.re -:= b;
end;

procedure mulab_cc;
var
  c: Complex;
with c do
begin
  re := a.re * b.re - a.im * b.im;
  im := a.re * b.im + a.im * b.re;
  a := c
end;

procedure mulab_cr;
with a do
begin
  re *:= b;
  im *:= b;
end;

procedure divab_cc;
var
  c: Complex;
  x: CReal;
begin
  x := sqr(b.re) + sqr(b.im);
  with c do
  begin
    re := (a.re * b.re + a.im * b.im) / x;
    im := (a.im * b.re - a.re * b.im) / x;
  end;
  a := c
end;

procedure divab_cr;
with a do
begin
  re /:= b;
  im /:= b;
end;

function eq_cc;
begin
  Result := (a.re  = b.re) and (a.im =  b.im);
end;

function eq_cr;
begin
  Result := (a.re = b) and (a.im = 0);
end;

function eq_rc;
begin
  Result := (a = b.re) and (b.im = 0);
end;

function ne_cc;
begin
  Result := (a.re <> b.re) or (a.im <> b.im);
end;

function ne_cr;
begin
  Result := (a.re <> b) OR (a.im <> 0);
end;

function ne_rc;
begin
  Result := (a <> b.re) or (b.im <> 0);
end;

function sqr_cc;
with Result do
begin
  im := 2 * a.re * a.im;
  re := sqr(a.re) - sqr(a.im);
end;

function sqrt_cc;
var
  absx: CReal;
begin
   absx := Hypot(a.re, a.im);
   with Result do begin
     im := sgn(a.im) * sqrt((absx - a.re) / 2);
     re := sqrt((absx + a.re) / 2);
   end;
end;

function ln_cc;
with Result do
begin
  im := ArcTan2(a.re, a.im);
  re := ln(Hypot(a.re, a.im));
end;

function exp_cc;
var
  tmp: CReal;
begin
  with Result do
  begin
    tmp := Exp(a.re);
    re := tmp * Cos(a.im);
    im := tmp * Sin(a.im);
  end;
end;

function power_cc;
begin
  Result := exp_cc(mul_cc(ln_cc(a), b));
end;

function sin_cc;
with Result do
begin
  im := Cos(a.re)* Sinh(a.im);
  re := Sin(a.re)* Cosh(a.im);
end;

function sinh_cc;
with Result do
begin
  im := Cosh(a.re) * Sin(a.im);
  re := Sinh(a.re) * Cos(a.im);
end;

function cos_cc;
with Result do
begin
  im := -Sin(a.re) * Sinh(a.im);
  re := Cos(a.re) * Cosh(a.im);
end;

function cosh_cc;
with Result do
begin
  im := Sinh(a.re) * Sin(a.im);
  re := Cosh(a.re) * Cos(a.im);
end;

function tan_cc;
var
  tmp: CReal;
begin
  with Result do
  begin
    tmp := Cos(2 * a.re) + Cosh(2 * a.im);
    re := Sin(2 * a.re) / tmp;
    im := Sinh(2 * a.re) / tmp;
  end;
end;

function tanh_cc;
var
  tmp: CReal;
begin
  with Result do
  begin
    tmp := Cosh(2 * a.re) + Cos(2 * a.im);
    re := Sinh(2 * a.re) / tmp;
    im := Sin(2 * a.im) / tmp;
  end;
end;

function cotan_cc;
var
  tmp: CReal;
begin
  with Result do
  begin
    tmp := Cos(2 * a.re) + Cosh(2 * a.im);
    re := Sin(2 * a.re) / tmp;
    im := Sinh(2 * a.re) / tmp;
  end;
end;

function cotanh_cc;
var
  tmp: CReal;
begin
  with Result do
  begin
    tmp := Cosh(2 * a.re) + Cos(2 * a.im);
    re := Sinh(2 * a.re) / tmp;
    im := Sin(2 * a.im) / tmp;
  end;
end;

function alphbeta(a: Complex): Complex;
var
  tmp: CReal;
  res: Complex;
begin
  tmp := Sqrt(Sqr(a.re + 1) + Sqr(a.im));
  res.re := Sqrt(Sqr(a.re - 1) + Sqr(a.im));
  res.im := (tmp + a.re) / 2;
  with Result do begin
    re := res.im - res.re;
    im := ln(res.im + sqrt(sqr(res.im) - 1))
  end;
end;

function arccos_cc;
var
  tmp: Complex;
begin
  tmp := alphbeta(a);
  with Result do
  begin
     re := ArcCos(tmp.re);
    im := tmp.im;
  end;
end;

function arccosh_cc;
var
  tmp: Complex;
begin
  tmp := arccos_cc(a);
  Result := Cjx(tmp);
end;

function arcsin_cc;
var
  tmp: Complex;
begin
  tmp := alphbeta(a);
  with Result do
  begin
     re := ArcSin(tmp.re);
    im := tmp.im;
  end;
end;

function arcsinh_cc;
var
  tmp: Complex;
begin
  tmp := arcsinh_cc(Cjx(a));
  Result := Cnjx(tmp);
end;

function arctan_cc;
var
  re2, im2: CReal;
begin
  re2 := Sqr(a.re);
  im2 := Sqr(a.im);
  with Result do
  begin
    im := Ln((re2 + im2 + 2 * a.im + 1) / (re2 + im2 - 2 * a.im + 1)) / 4;
    re := ArcTan(2 * a.re / (1 - re2 - re2)) / 2;
  end;
end;

function arctanh_cc;
var
  tmp: Complex;
begin
  tmp := arctanh_cc(Cjx(a));
  Result := Cnjx(tmp);
end;

end.