MODULE CalcGrunwald;
IMPORT NbrInt, NbrRe, MathRe, CalcFn;
CONST
OKay* = 0; MaxSubDivReached* = 1; Oscillation* = 2;
VAR
MaxIntervals-: NbrInt.Integer;
PROCEDURE Grunwald( f: CalcFn.ReArg; a, b, order: NbrRe.Real; VAR error, result: NbrRe.Real; VAR res: NbrInt.Integer );
VAR i, im1, ip1, n: NbrInt.Integer; absError, absX, dx, errorLast, fn, fnm1, last, tol: NbrRe.Real;
fns, save: POINTER TO ARRAY OF NbrRe.Real;
BEGIN
result := 0; n := 4; tol := NbrRe.Max( error, 100 * NbrRe.Epsilon ); NEW( fns, 3 ); dx := (b - a) / 2;
FOR i := 0 TO 2 DO fns[i] := f( b - i * dx ) END;
error := NbrRe.MaxNbr;
LOOP
save := fns; fns := NIL; NEW( fns, n + 1 ); dx := (b - a) / n;
FOR i := 0 TO n DO
IF NbrInt.Odd( i ) THEN fns[i] := f( b - i * dx ) ELSE fns[i] := save[i DIV 2] END
END;
last := result; result := 0;
FOR i := n - 1 TO 1 BY -1 DO
im1 := i - 1; ip1 := i + 1;
fn := order * (fns[ip1] - 2 * fns[i] + fns[im1]) / 2; fn := order * ((fns[ip1] - fns[im1]) - fn) / 4;
fn := fns[i] - fn; result := (result + fn) * (im1 - order) / i
END;
fnm1 := f( b + dx ); fn := order * (fns[1] - 2 * fns[0] + fnm1) / 2; fn := order * ((fns[1] - fnm1) - fn) / 4;
result := (result + fns[0] - fn) * MathRe.Power( n / b, order ); n := 2 * n; errorLast := error;
absError := NbrRe.Abs( result - last ); absX := NbrRe.Abs( result );
IF absX > 1 THEN error := absError / absX ELSE error := NbrRe.Max( absError, 50 * NbrRe.Epsilon * absX ) END;
IF error < tol THEN res := OKay; EXIT END;
IF n = MaxIntervals THEN res := MaxSubDivReached; EXIT END;
IF error > errorLast THEN result := last; res := Oscillation; EXIT END
END
END Grunwald;
PROCEDURE SolveD*( f: CalcFn.ReArg; x, order: NbrRe.Real; VAR error: NbrRe.Real; VAR res: NbrInt.Integer ): NbrRe.Real;
VAR result, zero: NbrRe.Real;
BEGIN
error := 0; result := 0; res := OKay; zero := 0;
IF f # NIL THEN
IF zero # x THEN
IF order < zero THEN order := -order END;
IF zero < x THEN Grunwald( f, zero, x, order, error, result, res )
ELSE Grunwald( f, x, zero, order, error, result, res ); result := -result
END
END
ELSE result := 0
END;
RETURN result
END SolveD;
PROCEDURE SolveI*( f: CalcFn.ReArg; x, order: NbrRe.Real; VAR error: NbrRe.Real; VAR res: NbrInt.Integer ): NbrRe.Real;
VAR result, zero: NbrRe.Real;
BEGIN
error := 0; result := 0; res := OKay; zero := 0;
IF f # NIL THEN
IF zero # x THEN
IF order > zero THEN order := -order END;
IF zero < x THEN Grunwald( f, zero, x, order, error, result, res )
ELSE Grunwald( f, x, zero, order, error, result, res ); result := -result
END
END
END;
RETURN result
END SolveI;
BEGIN
MaxIntervals := NbrInt.MaxNbr
END CalcGrunwald.