GENERIC MODULE ChebyPolynomialFast(R, RT);

Arithmetic for Modula-3, see doc for details

Abstract: Implementation of Modula-3 version of NR92, ch 5.

IMPORT Arithmetic AS Arith;
FROM RT IMPORT Cos, Pi;

CONST Module = "ChebyPolynomialFast.";

PROCEDURE Expand (func: Ftn; n: CARDINAL; ): T =
  VAR
    nr := FLOAT(n, R.T);
    f  := NEW(T, n + 1);         (* we skip f[0] *)
    x, factor, sum, jr, kr: R.T;
    z                           := NEW(T, n);
  BEGIN
    (*---load up the function values---*)
    FOR k := 1 TO n DO
      kr := FLOAT(k, R.T);
      x := Cos(Pi * (kr - RT.Half) / nr);
      f[k] := func(x);
    END;

    (*---compute coeffs---*)
    factor := R.Two / nr;
    FOR j := 0 TO n - 1 DO
      jr := FLOAT(j, R.T);
      sum := R.Zero;
      FOR k := 1 TO n DO
        kr := FLOAT(k, R.T);
        sum := sum + f[k] * Cos(Pi * jr * (kr - RT.Half) / nr);
      END;
      z[j] := factor * sum;
    END;
    RETURN z;
  END Expand;

PROCEDURE FindUpperExp (x: T; prec: R.T; ): CARDINAL =
  <* UNUSED *>
  CONST
    ftn = Module & "FindUpperExp";
  VAR error: R.T := R.Zero;
  BEGIN
    FOR i := LAST(x^) TO 0 BY -1 DO
      error := error + ABS(x[i]);
      IF error > prec THEN RETURN i + 1; END;
    END;
    RETURN 0;
  END FindUpperExp;

PROCEDURE Abort (x: T;           (* abort the expansion *)
                 m: CARDINAL;    (* before the m-th term *)
  ): T =
  VAR z := NEW(T, m);
  BEGIN
    z^ := SUBARRAY(x^, 0, m);
    RETURN z;
  END Abort;

PROCEDURE Eval (x: T; xi: R.T; ): R.T RAISES {Arith.Error} =
  <* UNUSED *>
  CONST
    ftn = Module & "Eval";
  VAR
    dj   := R.Zero;
    djp1 := R.Zero;
    djp2 := R.Zero;
  BEGIN
    IF xi < -R.One OR xi > R.One THEN
      (* need -1<x<+1 *)
      RAISE Arith.Error(NEW(Arith.ErrorOutOfRange).init());
    END;
    FOR j := LAST(x^) TO 1 BY -1 DO
      dj := R.Two * xi * djp1 - djp2 + x[j];
      djp2 := djp1;
      djp1 := dj;
    END;
    RETURN xi * djp1 - djp2 + RT.Half * x[0];
  END Eval;

PROCEDURE Derive (x: T; ): T =
  VAR
    n := NUMBER(x^);
    z := NEW(T, n - 1);
  BEGIN
    IF n >= 2 THEN
      z[n - 2] := FLOAT(2 * (n - 1), R.T) * x[n - 1];
      IF n >= 3 THEN
        z[n - 3] := FLOAT(2 * (n - 2), R.T) * x[n - 2];
        FOR j := n - 4 TO 0 BY -1 DO
          z[j] := z[j + 2] + FLOAT(2 * (j + 1), R.T) * x[j + 1];
        END;
      END;
    END;
    RETURN z;
  END Derive;

PROCEDURE Integrate (x: T; ): T =
  VAR
    n := NUMBER(x^);
    z := NEW(T, n + 1);
  BEGIN
    IF 0 <= n THEN
      z[0] := R.Zero;
      FOR j := 1 TO n - 2 DO
        z[j] := (x[j - 1] - x[j + 1]) / FLOAT(2 * j, R.T);
      END;
      IF 2 <= n THEN z[n - 1] := x[n - 2] / FLOAT(2 * (n - 1), R.T); END;
      IF 1 <= n THEN z[n] := x[n - 1] / FLOAT(2 * (n), R.T); END;
    END;
    RETURN z;
  END Integrate;

BEGIN
END ChebyPolynomialFast.