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

MODULE MathErf;   (** AUTHOR "adf"; PURPOSE "The error function:  erf(x) = (2/Vp) x0x exp(-t2) dt"; *)

(* To change to 64-bit reals, address the code fragments written in light red. *)

(* Ref:  J.F. Hart, E.W. Cheney, C.L. Lawson, H.J. Maehly, C.K. Mesztenyi, J.R. Rice, H.G. Thacher, Jr., and C. Witzgall,
		"Computer Approximations," in: The SIAM Series in Applied MathLematics, Wiley, New York, 1968. *)

IMPORT NbrInt, NbrRe, DataErrors, MathRe, MathReSeries;

VAR
	maxIterations: NbrInt.Integer;  twoBySqrtPi: NbrRe.Real;
	(* Whenever NbrRe.Real is a 32-bit real, define the following arrays. *)
	erfcP1: ARRAY 3 OF NbrRe.Real;
	erfcP, erfcQ1: ARRAY 4 OF NbrRe.Real;
	erfcQ: ARRAY 5 OF NbrRe.Real;
	(* Or whenever NbrRe.Real is a 64-bit real, define the following arrays. *)
	(*  erfcP1: ARRAY 6 OF NbrRe.Real;
	erfcP, erfcQ1: ARRAY 7 OF NbrRe.Real;
	erfcQ: ARRAY 8 OF NbrRe.Real;  *)

TYPE
	ErfP = OBJECT (MathReSeries.Coefficient)

		PROCEDURE Evaluate*;
		VAR m: NbrInt.Integer;
		BEGIN
			IF n = 0 THEN coef := 0
			ELSIF n = 1 THEN coef := 1
			ELSE
				coef := 2 * (n - 1) * x;  m := n + 1;
				IF NbrInt.Odd( m ) THEN coef := -coef END
			END;
			IF n > maxIterations THEN eos := TRUE;  DataErrors.ReWarning( x, "Did not converge -  timed out." ) END
		END Evaluate;

	END ErfP;

	ErfQ = OBJECT (MathReSeries.Coefficient)

		PROCEDURE Evaluate*;
		BEGIN
			IF n = 0 THEN coef := 0
			ELSIF n = 1 THEN coef := 1
			ELSE coef := 1 + 2 * (n - 1)
			END
		END Evaluate;

	END ErfQ;

	PROCEDURE Erfc( x: NbrRe.Real ): NbrRe.Real;
	VAR max: NbrInt.Integer;  erfc: NbrRe.Real;
	BEGIN
		max := 26;   (* Asymptote reached by this value. *)
		IF x > max THEN erfc := 0
		ELSIF x < -max THEN erfc := 2
		ELSE
			IF NbrRe.Abs( x ) > 9 THEN
				erfc := MathReSeries.TruncatedRationalFunction( erfcP1, erfcQ1, NbrRe.Abs( x ) )
			ELSE erfc := MathReSeries.TruncatedRationalFunction( erfcP, erfcQ, NbrRe.Abs( x ) )
			END;
			erfc := erfc * NbrRe.Exp( -x * x );
			IF x < 0 THEN erfc := 2 - erfc END
		END;
		RETURN erfc
	END Erfc;

	PROCEDURE Fn*( x: NbrRe.Real ): NbrRe.Real;
	VAR erf: NbrRe.Real;  p: ErfP;  q: ErfQ;
	BEGIN
		IF NbrRe.Abs( x ) > 0.1 THEN erf := 1 - Erfc( x )
		ELSE
			NEW( p );  NEW( q );
			erf := twoBySqrtPi * MathRe.Exp( -x * x ) * MathReSeries.ContinuedFraction( p, q, x )
		END;
		RETURN erf
	END Fn;

(*
	PROCEDURE CplxFn*( z: NbrCplx.Complex ): NbrCplx.Complex;
	VAR
	BEGIN
	END CplxFn;
*)


BEGIN
	maxIterations := 1000;  twoBySqrtPi := 2 / MathRe.Sqrt( NbrRe.Pi );
	(* Whenever NbrRe.Real is a 32-bit real, use the following eonstants. *)
	(* Constants from Table ERFC 5663 from "Computer Approximations".   *)
	erfcP[0] := 1.000464117E1;  erfcP[1] := 8.426552865;  erfcP[2] := 3.460259332;  erfcP[3] := 5.623536121E-1;
	erfcQ[0] := 1.000464117E1;  erfcQ[1] := 1.971558074E1;  erfcQ[2] := 1.570228809E1;  erfcQ[3] := 6.090748787;
	erfcQ[4] := 1.0;
	(* Constants from Table ERFC 5722 from "Computer Approximations".  *)
	erfcP1[0] := 6.141050179E-1;  erfcP1[1] := 3.295899049E-1;  erfcP1[2] := 5.641895902E-1;
	erfcQ1[0] := 2.953483222E-1;  erfcQ1[1] := 1.588352760;  erfcQ1[2] := 5.841849230E-1;  erfcQ1[3] := 1.0
	(* Or, whenever NbrRe.Real is a 64-bit real, use the following eonstants. *)
	(*  (* Constants from Table ERFC 5666 from "Computer Approximations". *)
	erfcP[0] := 4.4041373582475223325269D2;  erfcP[1] := 6.25686535769683006135193D2;
	erfcP[2] := 4.483171659914834429345063D2;  erfcP[3] := 1.9184014058796685752390301D2;
	erfcP[4] := 5.0991697628253203795669523D1;  erfcP[5] := 7.9239298287415448527194039D0;
	erfcP[6] := 0.56419994559825748305038673D0;  erfcQ[0] := 4.4041373582475221430599D2;
	erfcQ[1] := 1.122640220177046578853681D3;  erfcQ[2] := 1.2746672667576622866580268D3;
	erfcQ[3] := 8.3881036547840644197941778D2;  erfcQ[4] := 3.47122928811784331429078057D2;
	erfcQ[5] := 9.08727112035370362595507371D1;  erfcQ[6] := 1.404533730546113829847322587D1;
	erfcQ[7] := 1.0D0;
	(* Constants from Table ERFC 5725 from "Computer Approximations". *)
	erfcP1[0] := 2.9788656263939928862D0;  erfcP1[1] := 7.409740605964741794425D0;
	erfcP1[2] := 6.1602098531096305440906D0;  erfcP1[3] := 5.019049726784267463450058D0;
	erfcP1[4] := 1.275366644729965952479585264D0;  erfcP1[5] := 0.5641895835477550741253201704D0;
	erfcQ1[0] := 3.3690752069827527677D0;  erfcQ1[1] := 9.608965327192787870698D0;
	erfcQ1[2] := 1.708144074746600431571095D1;  erfcQ1[3] := 1.20489519278551290360340491D1;
	erfcQ1[4] := 9.396034016235054150430579648D0;  erfcQ1[5] := 2.260528520767326969591866945D0;
	erfcQ1[6] := 1.0D0  *)
END MathErf.