integrals: procedure options (main);   /* 1 September 2019 */

f: procedure (x, function) returns (float(18));
   declare x float(18), function fixed binary;
   select (function);
      when (1) return (x**3);
      when (2) return (1/x);
      when (3) return (x);
      when (4) return (x);
   end;
end f;

   declare (a, b) fixed decimal (10);
   declare (rect_area, trap_area, Simpson) float(18);
   declare (d, dx) float(18);
   declare (S1, S2) float(18);
   declare N fixed decimal (15), function fixed binary;
   declare k fixed decimal (7,2);

   put ('     Rectangle-left           Rectangle-mid            Rectangle-right' ||
        '        Trapezoid                 Simpson');
   do function = 1 to 4;
      select(function);
         when (1) do; N = 100;     a = 0; b = 1;    end;
         when (2) do; N = 1000;    a = 1; b = 100;  end;
         when (3) do; N = 5000000; a = 0; b = 5000; end;
         when (4) do; N = 6000000; a = 0; b = 6000; end;
      end;

      dx = (b-a)/float(N);

      /* Rectangle method, left-side */
      rect_area = 0;
      do d = 0 to N-1;
         rect_area = rect_area + dx*f(a + d*dx, function);
      end;
      put skip edit (rect_area) (E(25, 15));

      /* Rectangle method, mid-point */
      rect_area = 0;
      do d = 0 to N-1;
         rect_area = rect_area + dx*f(a + d*dx + dx/2, function);
      end;
      put edit (rect_area) (E(25, 15));

      /* Rectangle method, right-side */
      rect_area = 0;
      do d = 1 to N;
         rect_area = rect_area + dx*f(a + d*dx, function);
      end;
      put edit (rect_area) (E(25, 15));

      /* Trapezoid method */
      trap_area = 0;
      do d = 0 to N-1;
         trap_area = trap_area + dx*(f(a+d*dx, function) + f(a+(d+1)*dx, function))/2;
      end;
      put edit (trap_area) (X(1), E(25, 15));

      /* Simpson's Rule */
      S1 = f(a+dx/2, function);
      S2 = 0;
      do d = 1 to N-1;
         S1 = S1 + f(a+d*dx+dx/2, function);
         S2 = S2 + f(a+d*dx, function);
      end;
      Simpson = dx * (f(a, function) + f(b, function) + 4*S1 + 2*S2) / 6;
      put edit (Simpson) (X(1), E(25, 15));
   end;

end integrals;