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

MODULE MathCplxSeries;   (** AUTHOR "adf"; PURPOSE "Procedures for computing complex series"; *)

IMPORT NbrInt, NbrRe, NbrCplx, DataErrors;

TYPE
	Coefficient* = OBJECT
	VAR
		(** The index of a coefficient to be determined - supplied by the calling function, e.g., PowerSeries. *)
		n-: NbrInt.Integer;
		(** The argument of the function being evaluated - supplied by the calling function, e.g., PowerSeries. *)
		z-: NbrCplx.Complex;
		(** The End Of Series - to be updated by Evaluate. *)
		eos*: BOOLEAN;
		(** The derived coefficinet - to be updated by Evaluate. *)
		coef*: NbrCplx.Complex;

		(** Abstract - must be extended. Used to evaluate coefficients in user-defined series. *)
		PROCEDURE Evaluate*;
		BEGIN
			DataErrors.Error( "This function is abstract. It must be extended." )
		END Evaluate;

	END Coefficient;

VAR
	epsilon: NbrRe.Real;

	(**                   a1(z)
		y = b0(z) + ----------
		                  b1(z) + a2(z)
		                              ----------
		                              b2(z) + a3(z)
		                                          --------
		                                          b3(z) + ...
	 *)
	PROCEDURE ContinuedFraction*( a, b: Coefficient;  z: NbrCplx.Complex ): NbrCplx.Complex;
	(* Based on an algorithm from: Lentz, W.J., Applied Optics, 15, 1976, 668-671. *)
	VAR convergedLast, convergedThis: BOOLEAN;  c, d, delta, f: NbrCplx.Complex;
	BEGIN
		(* Get the leading coefficient, b0. *)
		a.n := 0;  a.z := z;  a.eos := FALSE;  b.n := 0;  b.z := z;  b.eos := FALSE;  b.Evaluate;
		(* Initialize Lentz's recursive algorithm. *)
		IF b.coef = 0 THEN b.coef := 1/NbrRe.MaxNbr END;
		f := b.coef;  c := f;  d := 0;
		(* The recursive algorithm of Lentz. *)
		convergedThis := FALSE;
		REPEAT
			NbrInt.Inc( a.n );  a.Evaluate;  NbrInt.Inc( b.n );  b.Evaluate;  c := b.coef + a.coef * z / c;
			d := 1 / (b.coef + a.coef * z * d);  delta := c * d;  f := delta * f;  convergedLast := convergedThis;
			convergedThis := NbrRe.Abs( 1 - NbrCplx.Abs( delta ) ) < epsilon
		UNTIL (convergedLast & convergedThis) OR (a.eos OR b.eos);
		RETURN f
	END ContinuedFraction;

(** The lengths of the supplied arrays  a & b  must be the same. *)
	PROCEDURE TruncatedContinuedFraction*( a, b: ARRAY OF NbrRe.Real;  z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR i, aLen, bLen: NbrInt.Integer;  c, d, f: NbrCplx.Complex;
	BEGIN
		aLen := LEN( a );  bLen := LEN( b );
		IF aLen # bLen THEN
			DataErrors.Error( "Lengths of supplied arrays must be equal." );  f := 0;  RETURN f END;
		IF b[0] = 0 THEN b[0] := 1/NbrRe.MaxNbr END;
		f := b[0];  c := f;  d := 0;  i := 1;
		REPEAT c := b[i] + a[i] * z / c;  d := 1 / (b[i] + a[i] * z * d);  f := c * d * f;  NbrInt.Inc( i ) UNTIL i = aLen;
		RETURN f
	END TruncatedContinuedFraction;

(** y = a0 + a1z + a2z2 + a3z3 + ... *)
	PROCEDURE PowerSeries*( a: Coefficient;  z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR convergedLast, convergedThis: BOOLEAN;  sum, update, zz: NbrCplx.Complex;
	BEGIN
		a.z := z;  a.n := 0;  a.eos := FALSE;  a.Evaluate;  sum := a.coef;  zz := 1;  convergedThis := FALSE;
		REPEAT
			NbrInt.Inc( a.n );  a.Evaluate;  zz := z * zz;  update := a.coef * zz;  sum := sum + update;
			convergedLast := convergedThis;  convergedThis := NbrCplx.Abs( update ) < (epsilon * NbrCplx.Abs( sum ))
		UNTIL (convergedLast & convergedThis) OR a.eos;
		RETURN sum
	END PowerSeries;

	PROCEDURE TruncatedPowerSeries*( a: ARRAY OF NbrRe.Real;  z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR i, len: NbrInt.Integer;  prod: NbrCplx.Complex;
	BEGIN
		len := LEN( a );  prod := a[len - 1] * z;
		FOR i := len - 2 TO 1 BY -1 DO prod := (a[i] + prod) * z END;
		prod := a[0] + prod;  RETURN prod
	END TruncatedPowerSeries;

(**         a0 + a1z + a2z2 + a3z3 + ...
	 	y = --------------------
		    b0 + b1z + b2z2 + b3z3 + ...
	 *)
	PROCEDURE RationalFunction*( a, b: Coefficient;  z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR denom, num, ratio: NbrCplx.Complex;
	BEGIN
		num := PowerSeries( a, z );  denom := PowerSeries( b, z );  ratio := num / denom;  RETURN ratio
	END RationalFunction;

	PROCEDURE TruncatedRationalFunction*( a, b: ARRAY OF NbrRe.Real;  z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR denom, num, ratio: NbrCplx.Complex;
	BEGIN
		num := TruncatedPowerSeries( a, z );  denom := TruncatedPowerSeries( b, z );  ratio := num / denom;
		RETURN ratio
	END TruncatedRationalFunction;

BEGIN
	epsilon := 10 * NbrRe.Epsilon
END MathCplxSeries.