(* CAPO - Computational Analysis Platform for Oberon - by Alan Freed and Felix Friedrich. *)
(* Version 1, Update 2 *)

MODULE MathCplx;   (** AUTHOR "adf"; PURPOSE "Complex math functions"; *)

(*	Refs: M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical
				Tables, National Bureau of Standards, Applied Mathematics Series, Vol. 55, 1964.
			R. V. Churchill, J. W. Brown and R. F. Verhey, Complex Variables and Applications, 3rd ed., McGraw-Hill, 1974.
			P. Midy and Y. Yakovlev, Mathematics and Computers in Simulation, Vol. 33, 1991, 33-49. *)

IMPORT NbrInt, NbrRe, NbrCplx, DataErrors, MathRe;

VAR
	ln2, ln2Inv, ln10, ln10Inv: NbrRe.Real;

	(** Returns a pseudo-random number  z = r exp(if)  uniformly distributed over the unit circle,  r N (0,1),  f N (-p,p)). *)
	PROCEDURE Random*( ): NbrCplx.Complex;
	VAR abs, arg: NbrRe.Real;  z: NbrCplx.Complex;
	BEGIN
		abs := MathRe.Sqrt( MathRe.Random() );  arg := NbrRe.Pi * (2 * MathRe.Random() - 1);
		NbrCplx.Set( abs, arg, z );  RETURN z
	END Random;

(** Power Functions *)
	PROCEDURE Sqrt*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR abs, arg, im, re: NbrRe.Real;  sqrt: NbrCplx.Complex;
	BEGIN
		im := NbrCplx.Im( z );
		IF im = 0 THEN
			re := NbrCplx.Re( z );
			(* Handle the branch cut along the negative real axis. *)
			IF re < 0 THEN NbrCplx.Set( 0, MathRe.Sqrt( -re ), sqrt )
			ELSIF re = 0 THEN sqrt := 0
			ELSE NbrCplx.Set( MathRe.Sqrt( re ), 0, sqrt )
			END
		ELSE abs := NbrCplx.Abs( z );  arg := NbrCplx.Arg( z );  NbrCplx.SetPolar( MathRe.Sqrt( abs ), arg / 2, sqrt )
		END;
		RETURN sqrt
	END Sqrt;

	PROCEDURE IntPower*( z: NbrCplx.Complex;  n: NbrInt.Integer ): NbrCplx.Complex;
	VAR abs, arg: NbrRe.Real;  power: NbrCplx.Complex;
	BEGIN
		IF n = 0 THEN
			IF z # 0 THEN power := 1 ELSE DataErrors.Error( "Both argument and exponent cannot be zero." );  power := 0 END
		ELSIF z = 0 THEN
			IF n > 0 THEN power := 0
			ELSE DataErrors.IntError( n, "Exponent cannot be negative when argument is zero." );  power := 0
			END
		ELSE
			abs := MathRe.IntPower( NbrCplx.Abs( z ), n );  arg := n * NbrCplx.Arg( z );
			NbrCplx.SetPolar( abs, arg, power )
		END;
		RETURN power
	END IntPower;

	PROCEDURE LnMod( z: NbrCplx.Complex ): NbrRe.Real;
	(* Algorithm of Midy and Yakovlev. *)
	VAR r, r2, im, im2, lower, upper, result: NbrRe.Real;
	BEGIN
		lower := 0.82;  upper := 1.22;  r := NbrRe.Abs( NbrCplx.Re( z ) );  r2 := r * r;
		im := NbrRe.Abs( NbrCplx.Im( z ) );  im2 := im * im;
		IF r > im THEN
			IF (2 * im2 < r2) & (lower < r) & (r < upper) THEN
				result := MathRe.Ln( r ) + MathRe.ArcTanh( im2 / (2 * r2 + im2) )
			ELSE result := MathRe.Ln( r2 + im2 ) / 2
			END
		ELSE
			IF (2 * r2 < im2) & (lower < im) & (im < upper) THEN
				result := MathRe.Ln( im ) + MathRe.ArcTanh( r2 / (2 * im2 + r2) )
			ELSE result := MathRe.Ln( r2 + im2 ) / 2
			END
		END;
		RETURN result
	END LnMod;

	PROCEDURE RealPower*( z: NbrCplx.Complex;  x: NbrRe.Real ): NbrCplx.Complex;
	VAR abs, arg: NbrRe.Real;  power: NbrCplx.Complex;
	BEGIN
		abs := NbrCplx.Abs( z );
		IF x = 0 THEN
			IF abs > 0 THEN power := 1 ELSE DataErrors.Error( "Both argument and exponent cannot be zero." );  power := 0 END
		ELSIF abs = 0 THEN
			IF x > 0 THEN power := 0
			ELSE DataErrors.ReError( x, "Exponent cannot be negative when argument is zero." );  power := 0
			END
		ELSE abs := MathRe.Exp( x * LnMod( z ) );  arg := x * NbrCplx.Arg( z );  NbrCplx.SetPolar( abs, arg, power )
		END;
		RETURN power
	END RealPower;

	PROCEDURE Power*( zc, ze: NbrCplx.Complex ): NbrCplx.Complex;
	VAR abs, absC, absE, arg, im, lnRho, phi, re: NbrRe.Real;  power: NbrCplx.Complex;
		(* Formula listed in Midy and Yakovlev. *)
	BEGIN
		absC := NbrCplx.Abs( zc );  absE := NbrCplx.Abs( ze );
		IF absE = 0 THEN
			IF absC > 0 THEN power := 1
			ELSE DataErrors.Error( "Both argument and exponent cannot be zero." );  power := 0
			END
		ELSIF absC = 0 THEN power := 0
		ELSE
			re := NbrCplx.Re( ze );  im := NbrCplx.Im( ze );  lnRho := LnMod( zc );  phi := NbrCplx.Arg( zc );
			abs := MathRe.Exp( lnRho * re - phi * im );  arg := lnRho * im + phi * re;  NbrCplx.SetPolar( abs, arg, power )
		END;
		RETURN power
	END Power;

	(** Logarithmic Functions *)
(** Exp(z) = Exp(z 1 2pi), i.e., it is periodic with imaginary period  2pi. *)
	PROCEDURE Exp*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR coef, re, real, im, imag: NbrRe.Real;  exp: NbrCplx.Complex;
	BEGIN
		(* Euler's formula *)
		re := NbrCplx.Re( z );  im := NbrCplx.Im( z );  coef := MathRe.Exp( re );  real := coef * MathRe.Cos( im );
		imag := coef * MathRe.Sin( im );  NbrCplx.Set( real, imag, exp );  RETURN exp
	END Exp;

	PROCEDURE Exp2*( z: NbrCplx.Complex ): NbrCplx.Complex;
	BEGIN
		RETURN Exp( ln2 * z )
	END Exp2;

	PROCEDURE Exp10*( z: NbrCplx.Complex ): NbrCplx.Complex;
	BEGIN
		RETURN Exp( ln10 * z )
	END Exp10;

(** Ln(z) returns its principal value where Im(Ln(z)) N [-p, p]. *)
	PROCEDURE Ln*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR re, im: NbrRe.Real;  ln: NbrCplx.Complex;
	BEGIN
		re := LnMod( z );  im := NbrCplx.Arg( z );  NbrCplx.Set( re, im, ln );  RETURN ln
	END Ln;

	PROCEDURE Log2*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR log: NbrCplx.Complex;
	BEGIN
		log := ln2Inv * Ln( z );  RETURN log
	END Log2;

	PROCEDURE Log*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR log: NbrCplx.Complex;
	BEGIN
		log := ln10Inv * Ln( z );  RETURN log
	END Log;

	(* Using the Exp and Ln functions from above to define all functions below
		implies that they will return principal value(s) in their periodic component(s). *)

(** Trigonometric Functions *)
	PROCEDURE Sin*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR sin: NbrCplx.Complex;
	BEGIN
		sin := (Exp( NbrCplx.I * z ) - Exp( -NbrCplx.I * z )) / (2 * NbrCplx.I);  RETURN sin
	END Sin;

	PROCEDURE Cos*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR cos: NbrCplx.Complex;
	BEGIN
		cos := (Exp( NbrCplx.I * z ) + Exp( -NbrCplx.I * z )) / 2;  RETURN cos
	END Cos;

	PROCEDURE Tan*( z: NbrCplx.Complex ): NbrCplx.Complex;
	(* Algorithm of Midy and Yakovlev. *)
	VAR d, w, w2, x, x2, y, y2, re, real, im, imag, four, limit: NbrRe.Real;  tan: NbrCplx.Complex;
	BEGIN
		limit := 0.65;  re := NbrCplx.Re( z );  im := NbrCplx.Im( z );  x := MathRe.Tan( re );  x2 := x * x;
		IF NbrRe.Abs( im ) < limit THEN
			y := MathRe.Tanh( im );  y2 := y * y;  d := 1 + x2 * y2;  real := x * (1 - y2) / d;  imag := y * (1 + x2) / d
		ELSE
			four := 4;  y := MathRe.Exp( MathRe.Ln( four ) - 2 * NbrRe.Abs( im ) );  w := y / four;  w2 := w * w;
			d := 1 + 2 * w + w2 + (1 - 2 * w + w2) * x2;  real := x * y / d;
			imag := NbrRe.Sign( im ) * (1 - w2) * (1 + x2) / d
		END;
		NbrCplx.Set( real, imag, tan );  RETURN tan
	END Tan;

	PROCEDURE ArcSin*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR asin: NbrCplx.Complex;
	BEGIN
		IF NbrCplx.Im( z ) = 0 THEN asin := MathRe.ArcSin( NbrCplx.Re( z ) )
		ELSE asin := -NbrCplx.I * Ln( NbrCplx.I * z + Sqrt( 1 - z * z ) )
		END;
		RETURN asin
	END ArcSin;

	PROCEDURE ArcCos*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR acos: NbrCplx.Complex;
	BEGIN
		IF NbrCplx.Im( z ) = 0 THEN acos := MathRe.ArcCos( NbrCplx.Re( z ) )
		ELSE acos := -NbrCplx.I * Ln( z + NbrCplx.I * Sqrt( 1 - z * z ) )
		END;
		RETURN acos
	END ArcCos;

	PROCEDURE ArcTan*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR im: NbrRe.Real;  atan: NbrCplx.Complex;
	BEGIN
		IF NbrCplx.Re( z ) = 0 THEN
			(* Address the branch cuts along the imaginary axis. *)
			im := NbrCplx.Im( z );
			IF NbrRe.Abs( im ) < 1 THEN atan := NbrCplx.I * (MathRe.Ln( (1 + im) / (1 - im) ) / 2)
			ELSE DataErrors.CplxError( z, "Argument lies outside the admissible range for this function." );  atan := 0
			END
		ELSE atan := NbrCplx.I * Ln( (NbrCplx.I + z) / (NbrCplx.I - z) ) / 2
		END;
		RETURN atan
	END ArcTan;

(** Hyperbolic functions *)
	PROCEDURE Sinh*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR sinh: NbrCplx.Complex;
	BEGIN
		sinh := (Exp( z ) - Exp( -z )) / 2;  RETURN sinh
	END Sinh;

	PROCEDURE Cosh*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR cosh: NbrCplx.Complex;
	BEGIN
		cosh := (Exp( z ) + Exp( -z )) / 2;  RETURN cosh
	END Cosh;

	PROCEDURE Tanh*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR tanh: NbrCplx.Complex;
	BEGIN
		tanh := -NbrCplx.I * Tan( NbrCplx.I * z );  RETURN tanh
	END Tanh;

	PROCEDURE ArcSinh*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR abs, im: NbrRe.Real;  asinh: NbrCplx.Complex;
	BEGIN
		IF NbrCplx.Re( z ) = 0 THEN
			(* Address the branch cut along the imaginary axis. *)
			im := NbrCplx.Im( z );  abs := NbrRe.Abs( im );
			IF abs < 1 THEN asinh := Ln( z + MathRe.Sqrt( 1 - im * im ) )
			ELSIF abs = 1 THEN asinh := 0
			ELSE DataErrors.CplxError( z, "Argument lies outside the admissible range for this function." );  asinh := 0
			END
		ELSE asinh := Ln( z + Sqrt( z * z + 1 ) )
		END;
		RETURN asinh
	END ArcSinh;

	PROCEDURE ArcCosh*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR acosh: NbrCplx.Complex;
	BEGIN
		IF NbrCplx.Im( z ) = 0 THEN acosh := MathRe.ArcCosh( NbrCplx.Re( z ) ) ELSE acosh := Ln( z + Sqrt( z * z - 1 ) ) END;
		RETURN acosh
	END ArcCosh;

	PROCEDURE ArcTanh*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR atanh: NbrCplx.Complex;
	BEGIN
		IF NbrCplx.Im( z ) = 0 THEN atanh := MathRe.ArcTanh( NbrCplx.Re( z ) ) ELSE atanh := Ln( (1 + z) / (1 - z) ) / 2 END;
		RETURN atanh
	END ArcTanh;

BEGIN
	ln2 := MathRe.Ln( 2 );  ln2Inv := 1/ln2;  ln10 := MathRe.Ln( 10 );  ln10Inv := 1/ln10
END MathCplx.