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

MODULE NbrCplx;   (** AUTHOR "adf"; PURPOSE "Defines a base Complex type for scientific computing."; *)

(** Cartesian componets:  z = x + i y  where  x, y N !.
	Polar components:  z = r exp(if)  where  r, f N !; r 3 0.
	Imaginary unit:  i = V-1. *)

IMPORT Streams, NbrInt, NbrRat, NbrRe;

TYPE
	Complex* = RECORD
		re, im: NbrRe.Real
	END;

VAR
	zero, one: NbrRe.Real;
	(**  The imaginary unit:  i = V-1. *)
	I-: Complex;

	(* Local procedure *)
	PROCEDURE Reciprocal( x: Complex ): Complex;
	VAR cplx: Complex;  denom, ratio: NbrRe.Real;
	BEGIN
		IF NbrRe.Abs( x.re ) > NbrRe.Abs( x.im ) THEN
			ratio := x.im / x.re;  denom := x.re + x.im * ratio;  cplx.re := 1 / denom;  cplx.im := -ratio / denom
		ELSE ratio := x.re / x.im;  denom := x.re * ratio + x.im;  cplx.re := ratio / denom;  cplx.im := -1 / denom
		END;
		RETURN cplx
	END Reciprocal;

	(** Monadic Operators *)
(** Negative, i.e., -z = -x - i y *)
	PROCEDURE "-"*( x: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := -x.re;  cplx.im := -x.im;  RETURN cplx
	END "-";

(** Complex conjugate, i.e., ~z = x - i y  *)
	PROCEDURE "~"*( x: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := x.re;  cplx.im := -x.im;  RETURN cplx
	END "~";

	(** Dyadic Operators *)
(** Type Conversions *)
	PROCEDURE ":="*( VAR l: Complex;  r: NbrInt.Integer );
	BEGIN
		l.re := r;  l.im := 0
	END ":=";

	PROCEDURE ":="*( VAR l: Complex;  r: NbrRat.Rational );
	BEGIN
		l.re := r;  l.im := 0
	END ":=";

	PROCEDURE ":="*( VAR l: Complex;  r: NbrRe.Real );
	BEGIN
		l.re := r;  l.im := 0
	END ":=";

(** Comparison Operators *)
	PROCEDURE "="*( l, r: Complex ): BOOLEAN;
	BEGIN
		RETURN (l.re = r.re) & (l.im = r.im)
	END "=";

	PROCEDURE "="*( l: Complex;  r: NbrRe.Real ): BOOLEAN;
	BEGIN
		RETURN (l.re = r) & (l.im = 0)
	END "=";

	PROCEDURE "="*( l: Complex;  r: NbrRat.Rational ): BOOLEAN;
	BEGIN
		RETURN (l.re = r) & (l.im = 0)
	END "=";

	PROCEDURE "="*( l: Complex;  r: NbrInt.Integer ): BOOLEAN;
	BEGIN
		RETURN (l.re = r) & (l.im = 0)
	END "=";

	PROCEDURE "="*( l: NbrRe.Real;  r: Complex ): BOOLEAN;
	BEGIN
		RETURN (l = r.re) & (r.im = 0)
	END "=";

	PROCEDURE "="*( l: NbrRat.Rational;  r: Complex ): BOOLEAN;
	BEGIN
		RETURN (l = r.re) & (r.im = 0)
	END "=";

	PROCEDURE "="*( l: NbrInt.Integer;  r: Complex ): BOOLEAN;
	BEGIN
		RETURN (l = r.re) & (r.im = 0)
	END "=";

	PROCEDURE "#"*( l, r: Complex ): BOOLEAN;
	BEGIN
		RETURN (l.re # r.re) OR (l.im # r.im)
	END "#";

	PROCEDURE "#"*( l: Complex;  r: NbrRe.Real ): BOOLEAN;
	BEGIN
		RETURN (l.re # r) OR (l.im # 0)
	END "#";

	PROCEDURE "#"*( l: Complex;  r: NbrRat.Rational ): BOOLEAN;
	BEGIN
		RETURN (l.re # r) OR (l.im # 0)
	END "#";

	PROCEDURE "#"*( l: Complex;  r: NbrInt.Integer ): BOOLEAN;
	BEGIN
		RETURN (l.re # r) OR (l.im # 0)
	END "#";

	PROCEDURE "#"*( l: NbrRe.Real;  r: Complex ): BOOLEAN;
	BEGIN
		RETURN (l # r.re) OR (r.im # 0)
	END "#";

	PROCEDURE "#"*( l: NbrRat.Rational;  r: Complex ): BOOLEAN;
	BEGIN
		RETURN (l # r.re) OR (r.im # 0)
	END "#";

	PROCEDURE "#"*( l: NbrInt.Integer;  r: Complex ): BOOLEAN;
	BEGIN
		RETURN (l # r.re) OR (r.im # 0)
	END "#";

(** Arithmetic *)
	PROCEDURE "+"*( l, r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re + r.re;  cplx.im := l.im + r.im;  RETURN cplx
	END "+";

	PROCEDURE "+"*( l: Complex;  r: NbrRe.Real ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re + r;  cplx.im := l.im;  RETURN cplx
	END "+";

	PROCEDURE "+"*( l: Complex;  r: NbrRat.Rational ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re + r;  cplx.im := l.im;  RETURN cplx
	END "+";

	PROCEDURE "+"*( l: Complex;  r: NbrInt.Integer ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re + r;  cplx.im := l.im;  RETURN cplx
	END "+";

	PROCEDURE "+"*( l: NbrRe.Real;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l + r.re;  cplx.im := r.im;  RETURN cplx
	END "+";

	PROCEDURE "+"*( l: NbrRat.Rational;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l + r.re;  cplx.im := r.im;  RETURN cplx
	END "+";

	PROCEDURE "+"*( l: NbrInt.Integer;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l + r.re;  cplx.im := r.im;  RETURN cplx
	END "+";

	PROCEDURE "-"*( l, r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re - r.re;  cplx.im := l.im - r.im;  RETURN cplx
	END "-";

	PROCEDURE "-"*( l: Complex;  r: NbrRe.Real ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re - r;  cplx.im := l.im;  RETURN cplx
	END "-";

	PROCEDURE "-"*( l: Complex;  r: NbrRat.Rational ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re - r;  cplx.im := l.im;  RETURN cplx
	END "-";

	PROCEDURE "-"*( l: Complex;  r: NbrInt.Integer ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re - r;  cplx.im := l.im;  RETURN cplx
	END "-";

	PROCEDURE "-"*( l: NbrRe.Real;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l - r.re;  cplx.im := -r.im;  RETURN cplx
	END "-";

	PROCEDURE "-"*( l: NbrRat.Rational;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l - r.re;  cplx.im := -r.im;  RETURN cplx
	END "-";

	PROCEDURE "-"*( l: NbrInt.Integer;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l - r.re;  cplx.im := -r.im;  RETURN cplx
	END "-";

	PROCEDURE "*"*( l, r: Complex ): Complex;
	VAR cplx: Complex;  left, right: NbrRe.Real;
	BEGIN
		left := l.re * r.re;  right := l.im * r.im;  cplx.re := left - right;  left := l.re * r.im;  right := l.im * r.re;
		cplx.im := left + right;  RETURN cplx
	END "*";

	PROCEDURE "*"*( l: Complex;  r: NbrRe.Real ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re * r;  cplx.im := l.im * r;  RETURN cplx
	END "*";

	PROCEDURE "*"*( l: Complex;  r: NbrRat.Rational ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re * r;  cplx.im := l.im * r;  RETURN cplx
	END "*";

	PROCEDURE "*"*( l: Complex;  r: NbrInt.Integer ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l.re * r;  cplx.im := l.im * r;  RETURN cplx
	END "*";

	PROCEDURE "*"*( l: NbrRe.Real;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l * r.re;  cplx.im := l * r.im;  RETURN cplx
	END "*";

	PROCEDURE "*"*( l: NbrRat.Rational;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l * r.re;  cplx.im := l * r.im;  RETURN cplx
	END "*";

	PROCEDURE "*"*( l: NbrInt.Integer;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx.re := l * r.re;  cplx.im := l * r.im;  RETURN cplx
	END "*";

	PROCEDURE "/"*( l, r: Complex ): Complex;
	(* Algorithm can be found in Press et al. *)
	VAR cplx: Complex;  denom, ratio: NbrRe.Real;
	BEGIN
		IF NbrRe.Abs( r.re ) > NbrRe.Abs( r.im ) THEN
			ratio := r.im / r.re;  denom := r.re + r.im * ratio;  cplx.re := (l.re + l.im * ratio) / denom;
			cplx.im := (l.im - l.re * ratio) / denom
		ELSE
			ratio := r.re / r.im;  denom := r.re * ratio + r.im;  cplx.re := (l.re * ratio + l.im) / denom;
			cplx.im := (l.im * ratio - l.re) / denom
		END;
		RETURN cplx
	END "/";

	PROCEDURE "/"*( l: Complex;  r: NbrRe.Real ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx := l * (1 / r);  RETURN cplx
	END "/";

	PROCEDURE "/"*( l: Complex;  r: NbrRat.Rational ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx := l * (1 / r);  RETURN cplx
	END "/";

	PROCEDURE "/"*( l: Complex;  r: NbrInt.Integer ): Complex;
	VAR re: NbrRe.Real;  cplx: Complex;
	BEGIN
		re := r;  cplx := l * (1 / re);  RETURN cplx
	END "/";

	PROCEDURE "/"*( l: NbrRe.Real;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx := l * Reciprocal( r );  RETURN cplx
	END "/";

	PROCEDURE "/"*( l: NbrRat.Rational;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx := l * Reciprocal( r );  RETURN cplx
	END "/";

	PROCEDURE "/"*( l: NbrInt.Integer;  r: Complex ): Complex;
	VAR cplx: Complex;
	BEGIN
		cplx := l * Reciprocal( r );  RETURN cplx
	END "/";

	(** Basic functions *)
(** x  in  z = x + i y *)
	PROCEDURE Re*( x: Complex ): NbrRe.Real;
	BEGIN
		RETURN x.re
	END Re;

(** y  in  z = x + i y *)
	PROCEDURE Im*( x: Complex ): NbrRe.Real;
	BEGIN
		RETURN x.im
	END Im;

(** r  in  z = r exp(if) *)
	PROCEDURE Abs*( x: Complex ): NbrRe.Real;
	VAR abs, absRe, absIm, ratio: NbrRe.Real;
	BEGIN
		absRe := NbrRe.Abs( x.re );  absIm := NbrRe.Abs( x.im );
		IF absRe > absIm THEN ratio := absIm / absRe;  abs := absRe * NbrRe.Sqrt( 1 + ratio * ratio )
		ELSIF absIm = 0 THEN abs := 0
		ELSE ratio := absRe / absIm;  abs := absIm * NbrRe.Sqrt( 1 + ratio * ratio )
		END;
		RETURN abs
	END Abs;

(** f  in  z = r exp(if) *)
	PROCEDURE Arg*( x: Complex ): NbrRe.Real;
	VAR arg: NbrRe.Real;

		PROCEDURE ArcTan( xn, xd: NbrRe.Real ): NbrRe.Real;
		VAR sn, sd: NbrInt.Integer;  atan, ratio: NbrRe.Real;
		BEGIN
			IF xn < 0 THEN sn := -1
			ELSIF xn = 0 THEN sn := 0
			ELSE sn := 1
			END;
			IF xd < 0 THEN sd := -1
			ELSIF xd = 0 THEN sd := 0
			ELSE sd := 1
			END;
			IF xd = 0 THEN atan := sn * NbrRe.Pi / 2
			ELSIF xn = 0 THEN atan := (1 - sd) * NbrRe.Pi / 2
			ELSE ratio := xn / xd;  atan := NbrRe.ArcTan( ratio ) + sn * (1 - sd) * NbrRe.Pi / 2
			END;
			RETURN atan
		END ArcTan;

	BEGIN
		IF x = 0 THEN arg := 0 ELSE arg := ArcTan( x.im, x.re ) END;
		RETURN arg
	END Arg;

(** re = x  and  im = y  in  z = x + i y. *)
	PROCEDURE Get*( x: Complex;  VAR re, im: NbrRe.Real );
	BEGIN
		re := x.re;  im := x.im
	END Get;

	PROCEDURE Set*( re, im: NbrRe.Real;  VAR x: Complex );
	BEGIN
		x.re := re;  x.im := im
	END Set;

(** abs = r  and  arg = f  in  z = r exp(if). *)
	PROCEDURE GetPolar*( x: Complex;  VAR abs, arg: NbrRe.Real );
	BEGIN
		abs := Abs( x );  arg := Arg( x )
	END GetPolar;

	PROCEDURE SetPolar*( abs, arg: NbrRe.Real;  VAR x: Complex );
	VAR cos, sin: NbrRe.Real;
	BEGIN
		cos := NbrRe.Cos( arg );  sin := NbrRe.Sin( arg );  x.re := abs * cos;  x.im := abs * sin
	END SetPolar;

	(** String conversions. *)
(** Admissible characters: {" ", "+", "-", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "E", "i", ",", "."}. *)
	PROCEDURE StringToCplx*( string: ARRAY OF CHAR;  VAR x: Complex );
	VAR negative: BOOLEAN;  i, j: NbrInt.Integer;
		reString, imString: ARRAY 32 OF CHAR;
	BEGIN
		i := 0;
		WHILE (string[i] # 0X) & (string[i] # "i") DO reString[i] := string[i];  NbrInt.Inc( i ) END;
		j := i;
		LOOP
			NbrInt.Dec( j );
			IF j < 1 THEN x.re := 0;  EXIT END;
			(* Pass over any tailing white space and a plus or minus sign. *)
			IF (reString[j] # CHR( 20H )) & ((reString[j] # CHR( 2BH )) OR (reString[j] = CHR( 2DH ))) THEN
				reString[j + 1] := 0X;  NbrRe.StringToRe( reString, x.re );  EXIT
			END
		END;
		IF string[i] = "i" THEN
			j := i;
			LOOP
				NbrInt.Dec( j );
				(* Determine the sign of the imaginary part. *)
				IF j < 1 THEN negative := FALSE;  EXIT END;
				IF string[j] = CHR( 2BH ) THEN negative := FALSE;  EXIT END;
				IF string[j] = CHR( 2DH ) THEN negative := TRUE;  EXIT END
			END;
			NbrInt.Inc( i );  j := 0;
			REPEAT imString[j] := string[i];  NbrInt.Inc( i );  NbrInt.Inc( j ) UNTIL string[i] = 0X;
			NbrRe.StringToRe( imString, x.im );
			IF negative THEN x.im := -x.im END
		ELSE x.im := 0
		END
	END StringToCplx;

(** LEN(string) >= 2*significantFigures+18 *)
	PROCEDURE CplxToString*( x: Complex;  significantFigures: NbrInt.Integer;  VAR string: ARRAY OF CHAR );
	VAR i, j: NbrInt.Integer;
		reString, imString: ARRAY 32 OF CHAR;
	BEGIN
		NbrRe.ReToString( x.re, significantFigures, reString );  i := 0;
		REPEAT string[i] := reString[i];  NbrInt.Inc( i ) UNTIL reString[i] = 0X;
		IF x.im > 0 THEN
			string[i] := CHR( 20H );  NbrInt.Inc( i );  string[i] := "+";  NbrInt.Inc( i );  string[i] := CHR( 20H );  NbrInt.Inc( i );  string[i] := "i";
			NbrInt.Inc( i );  string[i] := CHR( 20H );  NbrInt.Inc( i );  NbrRe.ReToString( x.im, significantFigures, imString );  j := 0;
			REPEAT string[i] := imString[j];  NbrInt.Inc( i );  NbrInt.Inc( j ) UNTIL imString[j] = 0X
		ELSIF x.im = 0 THEN  (* write nothing *)
		ELSE
			string[i] := CHR( 20H );  NbrInt.Inc( i );  string[i] := "-";  NbrInt.Inc( i );  string[i] := CHR( 20H );  NbrInt.Inc( i );  string[i] := "i";
			NbrInt.Inc( i );  string[i] := CHR( 20H );  NbrInt.Inc( i );  NbrRe.ReToString( -x.im, significantFigures, imString );  j := 0;
			REPEAT string[i] := imString[j];  NbrInt.Inc( i );  NbrInt.Inc( j ) UNTIL imString[j] = 0X
		END;
		string[i] := 0X
	END CplxToString;

(** Admissible characters: {" ", "+", "-", "0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "E", "exp", "(", ")", "i", ",", "."}. *)
	PROCEDURE PolarStringToCplx*( string: ARRAY OF CHAR;  VAR x: Complex );
	VAR negative: BOOLEAN;  i, j: NbrInt.Integer;  abs, arg: NbrRe.Real;
		absString, argString: ARRAY 32 OF CHAR;
	BEGIN
		i := 0;
		WHILE (string[i] # 0X) & (string[i] # "e") DO absString[i] := string[i];  NbrInt.Inc( i ) END;
		IF i = 0 THEN x.re := 0;  x.im := 0;  RETURN END;
		absString[i] := 0X;  NbrRe.StringToRe( absString, abs );
		IF string[i] = "e" THEN
			REPEAT NbrInt.Inc( i ) UNTIL string[i] = "(";
			NbrInt.Inc( i );
			(* Pass over leading white space. *)
			WHILE string[i] = CHR( 20H ) DO NbrInt.Inc( i ) END;
			(* Determine the sign. *)
			IF string[i] = CHR( 2DH ) THEN negative := TRUE;  NbrInt.Inc( i ) ELSE negative := FALSE END;
			(* Move to the first digit. *)
			WHILE (string[i] < CHR( 30H )) OR (string[i] > CHR( 39H )) DO NbrInt.Inc( i ) END;
			j := 0;
			WHILE string[i] # ")" DO argString[j] := string[i];  NbrInt.Inc( i );  NbrInt.Inc( j ) END;
			argString[j] := 0X;  NbrRe.StringToRe( argString, arg );
			IF negative THEN arg := -arg END
		ELSE arg := 0
		END;
		SetPolar( abs, arg, x )
	END PolarStringToCplx;

(** LEN(string) >= 2*significantFigures+22 *)
	PROCEDURE CplxToPolarString*( x: Complex;  significantFigures: NbrInt.Integer;  VAR string: ARRAY OF CHAR );
	VAR i, j: NbrInt.Integer;  abs, arg: NbrRe.Real;
		absString, argString: ARRAY 32 OF CHAR;
	BEGIN
		GetPolar( x, abs, arg );  NbrRe.ReToString( abs, significantFigures, absString );  i := 0;
		REPEAT string[i] := absString[i];  NbrInt.Inc( i ) UNTIL absString[i] = 0X;
		IF arg > 0 THEN
			string[i] := CHR( 20H );  NbrInt.Inc( i );  string[i] := "e";  NbrInt.Inc( i );  string[i] := "x";  NbrInt.Inc( i );  string[i] := "p";
			NbrInt.Inc( i );  string[i] := "(";  NbrInt.Inc( i );  string[i] := "i";  NbrInt.Inc( i );  string[i] := CHR( 20H );  NbrInt.Inc( i );
			NbrRe.ReToString( arg, significantFigures, argString );  j := 0;
			REPEAT string[i] := argString[j];  NbrInt.Inc( i );  NbrInt.Inc( j ) UNTIL argString[j] = 0X;
			string[i] := ")";  NbrInt.Inc( i )
		ELSIF arg = 0 THEN
			(* write nothing *)
		ELSE
			string[i] := CHR( 20H );  NbrInt.Inc( i );  string[i] := "e";  NbrInt.Inc( i );  string[i] := "x";  NbrInt.Inc( i );  string[i] := "p";
			NbrInt.Inc( i );  string[i] := "(";  NbrInt.Inc( i );  string[i] := "-";  NbrInt.Inc( i );  string[i] := "i";  NbrInt.Inc( i );
			string[i] := CHR( 20H );  NbrInt.Inc( i );  NbrRe.ReToString( -arg, significantFigures, argString );  j := 0;
			REPEAT string[i] := argString[j];  NbrInt.Inc( i );  NbrInt.Inc( j ) UNTIL argString[j] = 0X;
			string[i] := ")";  NbrInt.Inc( i )
		END;
		string[i] := 0X
	END CplxToPolarString;

(** Persistence: file IO *)
	PROCEDURE Load*( R: Streams.Reader;  VAR x: Complex );
	BEGIN
		NbrRe.Load( R, x.re );  NbrRe.Load( R, x.im )
	END Load;

	PROCEDURE Store*( W: Streams.Writer;  x: Complex );
	BEGIN
		NbrRe.Store( W, x.re );  NbrRe.Store( W, x.im )
	END Store;

BEGIN
	zero := 0;  one := 1;  Set( zero, one, I )
END NbrCplx.