MODULE numericalIntegrationModula2;

(* ISO Modula-2 libraries. *)
IMPORT LongMath, SLongIO, STextIO;

TYPE functionRealToReal = PROCEDURE (LONGREAL) : LONGREAL;

PROCEDURE leftRule (f : functionRealToReal;
                    a : LONGREAL;
                    b : LONGREAL;
                    n : INTEGER) : LONGREAL;
  VAR sum : LONGREAL;
      h   : LONGREAL;
      i   : INTEGER;
BEGIN
  sum := 0.0;
  h := (b - a) / LFLOAT (n);
  FOR i := 1 TO n DO
    sum := sum + f (a + (h * LFLOAT (i - 1)))
  END;
  RETURN (sum * h)
END leftRule;

PROCEDURE rightRule (f : functionRealToReal;
                     a : LONGREAL;
                     b : LONGREAL;
                     n : INTEGER) : LONGREAL;
  VAR sum : LONGREAL;
      h   : LONGREAL;
      i   : INTEGER;
BEGIN
  sum := 0.0;
  h := (b - a) / LFLOAT (n);
  FOR i := 1 TO n DO
    sum := sum + f (a + (h * LFLOAT (i)))
  END;
  RETURN (sum * h)
END rightRule;

PROCEDURE midpointRule (f : functionRealToReal;
                        a : LONGREAL;
                        b : LONGREAL;
                        n : INTEGER) : LONGREAL;
  VAR sum    : LONGREAL;
      h      : LONGREAL;
      half_h : LONGREAL;
      i      : INTEGER;
BEGIN
  sum := 0.0;
  h := (b - a) / LFLOAT (n);
  half_h := 0.5 * h;
  FOR i := 1 TO n DO
    sum := sum + f (a + (h * LFLOAT (i)) - half_h)
  END;
  RETURN (sum * h)
END midpointRule;

PROCEDURE trapeziumRule (f : functionRealToReal;
                         a : LONGREAL;
                         b : LONGREAL;
                         n : INTEGER) : LONGREAL;
  VAR sum : LONGREAL;
      y0  : LONGREAL;
      y1  : LONGREAL;
      h   : LONGREAL;
      i   : INTEGER;
BEGIN
  sum := 0.0;
  h := (b - a) / LFLOAT (n);
  y0 := f (a);
  FOR i := 1 TO n DO
    y1 := f (a + (h * LFLOAT (i)));
    sum := sum + 0.5 * (y0 + y1);
    y0 := y1
  END;
  RETURN (sum * h)
END trapeziumRule;


PROCEDURE simpsonRule (f : functionRealToReal;
                       a : LONGREAL;
                       b : LONGREAL;
                       n : INTEGER) : LONGREAL;
  VAR sum1   : LONGREAL;
      sum2   : LONGREAL;
      h      : LONGREAL;
      half_h : LONGREAL;
      x      : LONGREAL;
      i      : INTEGER;
BEGIN
  h := (b - a) / LFLOAT (n);
  half_h := 0.5 * h;
  sum1 := f (a + half_h);
  sum2 := 0.0;
  FOR i := 2 TO n DO
    x := a + (h * LFLOAT (i - 1));
    sum1 := sum1 + f (x + half_h);
    sum2 := sum2 + f (x);
  END;
  RETURN (h / 6.0) * (f (a) + f (b) + (4.0 * sum1) + (2.0 * sum2));
END simpsonRule;

PROCEDURE cube (x : LONGREAL) : LONGREAL;
BEGIN
  RETURN x * x * x;
END cube;

PROCEDURE reciprocal (x : LONGREAL) : LONGREAL;
BEGIN
  RETURN 1.0 / x;
END reciprocal;

PROCEDURE identity (x : LONGREAL) : LONGREAL;
BEGIN
  RETURN x;
END identity;

PROCEDURE printResults (f       : functionRealToReal;
                        a       : LONGREAL;
                        b       : LONGREAL;
                        n       : INTEGER;
                        nominal : LONGREAL);
  PROCEDURE printOneResult (y : LONGREAL);
  BEGIN
    SLongIO.WriteFloat (y, 16, 20);
    STextIO.WriteString ('  (nominal + ');
    SLongIO.WriteFloat (y - nominal, 6, 0);
    STextIO.WriteString (')');
    STextIO.WriteLn;
  END printOneResult;
BEGIN
  STextIO.WriteString ('  left rule       ');
  printOneResult (leftRule (f, a, b, n));

  STextIO.WriteString ('  right rule      ');
  printOneResult (rightRule (f, a, b, n));

  STextIO.WriteString ('  midpoint rule   ');
  printOneResult (midpointRule (f, a, b, n));

  STextIO.WriteString ('  trapezium rule  ');
  printOneResult (trapeziumRule (f, a, b, n));

  STextIO.WriteString ('  Simpson rule    ');
  printOneResult (simpsonRule (f, a, b, n));
END printResults;

BEGIN
  STextIO.WriteLn;

  STextIO.WriteString ('x³ in [0,1] with n = 100');
  STextIO.WriteLn;
  printResults (cube, 0.0, 1.0, 100, 0.25);

  STextIO.WriteLn;

  STextIO.WriteString ('1/x in [1,100] with n = 1000');
  STextIO.WriteLn;
  printResults (reciprocal, 1.0, 100.0, 1000, LongMath.ln (100.0));

  STextIO.WriteLn;

  STextIO.WriteString ('x in [0,5000] with n = 5000000');
  STextIO.WriteLn;
  printResults (identity, 0.0, 5000.0, 5000000, 12500000.0);

  STextIO.WriteLn;

  STextIO.WriteString ('x in [0,6000] with n = 6000000');
  STextIO.WriteLn;
  printResults (identity, 0.0, 6000.0, 6000000, 18000000.0);

  STextIO.WriteLn
END numericalIntegrationModula2.