```(* CAPO - Computational Analysis Platform for Oberon - by Alan Freed and Felix Friedrich. *) (* Version 1, Update 2 *) MODULE MathReSeries; (** AUTHOR "adf"; PURPOSE "Procedures for computing real series"; *) IMPORT NbrInt, NbrRe, 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. *) x-: NbrRe.Real; (** The End Of Series - to be updated by Evaluate. *) eos*: BOOLEAN; (** The derived coefficinet - to be updated by Evaluate. *) coef*: NbrRe.Real; (** 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(x) y = b0(x) + ---------- b1(x) + a2(x) ---------- b2(x) + a3(x) -------- b3(x) + ... *) PROCEDURE ContinuedFraction*( a, b: Coefficient; x: NbrRe.Real ): NbrRe.Real; (* Based on an algorithm from: Lentz, W.J., Applied Optics, 15, 1976, 668-671. *) VAR convergedLast, convergedThis: BOOLEAN; c, d, delta, f: NbrRe.Real; BEGIN (* Get the leading coefficient, b0. *) a.n := 0; a.x := x; a.eos := FALSE; b.n := 0; b.x := x; 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 * x / c; d := 1 / (b.coef + a.coef * x * d); delta := c * d; f := delta * f; delta := 1 - delta; convergedLast := convergedThis; convergedThis := NbrRe.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; x: NbrRe.Real ): NbrRe.Real; VAR i, aLen, bLen: NbrInt.Integer; c, d, f: NbrRe.Real; 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] * x / c; d := 1 / (b[i] + a[i] * x * d); f := c * d * f; NbrInt.Inc( i ) UNTIL i = aLen; RETURN f END TruncatedContinuedFraction; (** y = a0 + a1x + a2x2 + a3x3 + ... *) PROCEDURE PowerSeries*( a: Coefficient; x: NbrRe.Real ): NbrRe.Real; VAR convergedLast, convergedThis: BOOLEAN; sum, update, xx: NbrRe.Real; BEGIN a.x := x; a.n := 0; a.eos := FALSE; a.Evaluate; sum := a.coef; xx := 1; convergedThis := FALSE; REPEAT NbrInt.Inc( a.n ); a.Evaluate; xx := x * xx; update := a.coef * xx; sum := sum + update; convergedLast := convergedThis; convergedThis := NbrRe.Abs( update ) < (epsilon * NbrRe.Abs( sum )) UNTIL (convergedLast & convergedThis) OR a.eos; RETURN sum END PowerSeries; PROCEDURE TruncatedPowerSeries*( a: ARRAY OF NbrRe.Real; x: NbrRe.Real ): NbrRe.Real; VAR i, len: NbrInt.Integer; prod: NbrRe.Real; BEGIN len := LEN( a ); prod := a[len - 1] * x; FOR i := len - 2 TO 1 BY -1 DO prod := (a[i] + prod) * x END; prod := a[0] + prod; RETURN prod END TruncatedPowerSeries; (** a0 + a1x + a2x2 + a3x3 + ... y = -------------------- b0 + b1x + b2x2 + b3x3 + ... *) PROCEDURE RationalFunction*( a, b: Coefficient; x: NbrRe.Real ): NbrRe.Real; VAR denom, num, ratio: NbrRe.Real; BEGIN num := PowerSeries( a, x ); denom := PowerSeries( b, x ); ratio := num / denom; RETURN ratio END RationalFunction; PROCEDURE TruncatedRationalFunction*( a, b: ARRAY OF NbrRe.Real; x: NbrRe.Real ): NbrRe.Real; VAR denom, num, ratio: NbrRe.Real; BEGIN num := TruncatedPowerSeries( a, x ); denom := TruncatedPowerSeries( b, x ); ratio := num / denom; RETURN ratio END TruncatedRationalFunction; BEGIN epsilon := 10 * NbrRe.Epsilon END MathReSeries.```