Numerical integration/Adaptive Simpson's method

From Rosetta Code
Numerical integration/Adaptive Simpson's method is a draft programming task. It is not yet considered ready to be promoted as a complete task, for reasons that should be found in its talk page.

Lyness's (1969) Modified Adaptive Simpson's method (doi:10.1145/321526.321537) is a numerical quadrature method that recursively bisects the interval until the precision is high enough.

Write an implementation of quadrature by adaptive Simpson’s method. Use the implementation to estimate the definite integral of sin(x) from 0 to 1. Show your output.

You can use the following pseudocode, which includes Lyness's modifications 1, 2, and 3.

Pseudocode: Simpson's method, adaptive
; Lyness's Adaptive Simpson's Rule, Modifications 1, 2, 3
procedure simpson_rule(f, a, fa, b, fb)
    m := (a + b) / 2
    fm := f(m)
    h := b - a
    return multiple [m, fm, (h / 6) * (fa + 4*fm + fb)]

procedure recursive_simpson(f, a, fa, b, fb, tol, whole, m, fm, depth)
    lm, flm, left  := simpson_rule(f, a, fa, m, fm)
    rm, frm, right := simpson_rule(f, m, fm, b, fb)
    delta := left + right - whole
    tol' := tol / 2
    if depth <= 0 or tol' == tol or abs(delta) <= 15 * tol:
        return left + right + (delta / 15)
    else:
        return recursive_simpson(f, a, fa, m, fm, tol', left , lm, flm, depth - 1) +
               recursive_simpson(f, m, fm, b, fb, tol', right, rm, frm, depth - 1)

procedure quad_asr(f, a, b, tol, depth)
   fa := f(a)
   fb := f(b)
   m, fm, whole := simpson_rule(f, a, fa, b, fb)
   return recursive_simpson(f, a, fa, b, fb, tol, whole, m, fm, depth)

11l

Translation of: Python
T Triple
   Float m, fm, simp
   F (m, fm, simp)
      .m = m
      .fm = fm
      .simp = simp

F _quad_simpsons_mem(f, Float a, fa, b, fb)
   ‘
     Evaluates Simpson's Rule, also returning m and f(m) to reuse.
   ’
   V m = a + (b - a) / 2
   V fm = f(m)
   V simp = abs(b - a) / 6 * (fa + 4 * fm + fb)
   R Triple(m, fm, simp)

F _quad_asr(f, Float a, fa, b, fb, eps, whole, m, fm) -> Float
   ‘
     Efficient recursive implementation of adaptive Simpson's rule.
     Function values at the start, middle, end of the intervals are retained.
   ’
   V lt = _quad_simpsons_mem(f, a, fa, m, fm)
   V rt = _quad_simpsons_mem(f, m, fm, b, fb)
   V delta = lt.simp + rt.simp - whole
   R (I (abs(delta) <= eps * 15) {lt.simp + rt.simp + delta / 15}
      E _quad_asr(f, a, fa, m, fm, eps / 2, lt.simp, lt.m, lt.fm)
      + _quad_asr(f, m, fm, b, fb, eps / 2, rt.simp, rt.m, rt.fm))

F quad_asr(f, Float a, b, eps) -> Float
   ‘
     Integrate f from a to b using ASR with max error of eps.
   ’
   V fa = f(a)
   V fb = f(b)
   V t = _quad_simpsons_mem(f, a, fa, b, fb)
   R _quad_asr(f, a, fa, b, fb, eps, t.simp, t.m, t.fm)

V (a, b) = (0.0, 1.0)
V sinx = quad_asr(x -> sin(x), a, b, 1e-09)
print(‘Simpson's integration of sine from #. to #. = #.’.format(a, b, sinx))
Output:
Simpson's integration of sine from 0 to 1 = 0.459697694

Ada

with ada.text_io; use ada.text_io;
with ada.numerics; use ada.numerics;

with ada.numerics.elementary_functions;
use ada.numerics.elementary_functions;

procedure adaptive_simpson_task is

  subtype flt is float;
  type function_flt_to_flt is
    access function (x : in flt) return flt;

  procedure simpson_rule (f : in function_flt_to_flt;
                          a, fa, b, fb : in flt;
                          m, fm, quadval : out flt) is
  begin
    m := 0.5 * (a + b);
    fm := f(m);
    quadval := ((b - a) / 6.0) * (fa + (4.0 * fm) + fb);
  end;

  function recursive_simpson (f : in function_flt_to_flt;
                              a, fa, b, fb : in flt;
                              tol, whole, m, fm : in flt;
                              depth : in integer) return flt is
    lm, flm, left : flt;
    rm, frm, right : flt;
    diff, tol2, quadval : flt;
  begin
    simpson_rule (f, a, fa, m, fm, lm, flm, left);
    simpson_rule (f, m, fm, b, fb, rm, frm, right);
    diff := left + right - whole;
    tol2 := 0.5 * tol;
    if depth <= 0 or tol2 = tol or abs (diff) <= 15.0 * tol then
      quadval := left + right + (diff / 15.0);
    else
      quadval := recursive_simpson (f, a, fa, m, fm, tol2,
                                    left, lm, flm, depth - 1)
                    + recursive_simpson (f, m, fm, b, fb, tol2,
                                         right, rm, frm, depth - 1);
    end if;
    return quadval;
  end;

  function quad_asr (f : in function_flt_to_flt;
                     a, b, tol : in flt;
                     depth : in integer) return flt is
    fa, fb, m, fm, whole : flt;
  begin
    fa := f(a);
    fb := f(b);
    simpson_rule (f, a, fa, b, fb, m, fm, whole);
    return recursive_simpson (f, a, fa, b, fb, tol,
                              whole, m, fm, depth);
  end;

  function sine (x : in flt) return flt is
  begin
    return sin (x);
  end;

  quadval : flt;

begin
  quadval := quad_asr (sine'access, 0.0, 1.0, 1.0e-5, 100);
  put ("estimate of ∫ sin x dx from 0 to 1:");
  put_line (quadval'image);
end;

-- local variables:
-- mode: indented-text
-- tab-width: 2
-- end:
Output:
$ gnatmake -q adaptive_simpson_task.adb && ./adaptive_simpson_task
estimate of ∫ sin x dx from 0 to 1: 4.59698E-01

ALGOL 60

Works with: GNU MARST version 2.7
procedure simpsonrule (f, a, fa, b, fb, m, fm, quadval);
  value a, fa, b, fb;
  real procedure f;
  real a, fa, b, fb, m, fm, quadval;
begin
  m := (a + b) / 2;
  fm := f(m);
  quadval := ((b - a) / 6) * (fa + (4 * fm) + fb)
end simpsonrule;

real procedure recursion (f, a, fa, b, fb, tol, whole, m, fm, depth);
  value a, fa, b, fb, tol, whole, m, fm, depth;
  real procedure f;
  real a, fa, b, fb, tol, whole, m, fm;
  integer depth;
begin
  real lm, flm, left;
  real rm, frm, right;
  real delta, tol2;
  simpsonrule (f, a, fa, m, fm, lm, flm, left);
  simpsonrule (f, m, fm, b, fb, rm, frm, right);
  delta := left + right - whole;
  tol2 := tol / 2;
  if depth <= 0 | tol2 = tol | abs (delta) <= 15 * tol then
    recursion := left + right + (delta / 15)
  else
    recursion :=
      recursion (f, a, fa, m, fm, tol2, left , lm, flm, depth - 1)
        + recursion (f, m, fm, b, fb, tol2, right, rm, frm, depth - 1)
end recursion;

real procedure quadASR (f, a, b, tol, depth);
  value a, b, tol, depth;
  real procedure f;
  real a, b, tol;
  integer depth;
begin
  real fa, fb, m, fm, whole;
  fa := f(a);
  fb := f(b);
  simpsonrule (f, a, fa, b, fb, m, fm, whole);
  quadASR := recursion (f, a, fa, b, fb, tol, whole, m, fm, depth)
end quadASR;

begin
  real q;
  q := quadASR (sin, 0, 1, 0.000000001, 1000);
  outstring (1, "estimated definite integral of sin(x) ");
  outstring (1, "for x from 0 to 1: ");
  outreal (1, q);
  outstring (1, "\n")
end
Output:
$ marst adaptive_simpson_task.algol60 -o adaptive_simpson_task_algol60.c && gcc -g adaptive_simpson_task_algol60.c -lalgol -lm && ./a.out
estimated definite integral of sin(x) for x from 0 to 1: 0.459697694132

ALGOL 68

Translation of: C
...with the helper routines nested inside the main quadrature routine.
BEGIN # Numerical integration using an Adaptive Simpson's method #

    PROC quad asr = ( PROC(REAL)REAL f, REAL a, b, eps )REAL:
         BEGIN
            MODE TRIPLE = STRUCT( REAL m, fm, simp );

            # "structured" adaptive version, translated from Racket #
            PROC quad simpsons mem = ( PROC(REAL)REAL f, REAL a, fa, b, fb )TRIPLE:
                 BEGIN # Evaluates Simpson's Rule, also returning m and f(m) to reuse. #
                    REAL m    = ( a + b ) / 2;
                    REAL fm   = f( m );
                    REAL simp = ABS ( b - a ) / 6 * ( fa + 4*fm + fb );
                    ( m, fm, simp )
                 END # quad simpsons mem # ;

            PROC quad asr rec = ( PROC(REAL)REAL f, REAL a, fa, b, fb, eps, whole, m, fm )REAL:
                 IF # Efficient recursive implementation of adaptive Simpson's rule. #
                    # Function values at the start, middle, end of the intervals are retained. #
                    TRIPLE lt    = quad simpsons mem( f, a, fa, m, fm );
                    TRIPLE rt    = quad simpsons mem( f, m, fm, b, fb );
                    REAL   delta = simp OF lt + simp OF rt - whole;
                    ABS delta <= eps * 15
                 THEN simp OF lt + simp OF rt + delta / 15
                 ELSE quad asr rec( f, a, fa, m, fm, eps / 2, simp OF lt, m OF lt, fm OF lt )
                    + quad asr rec( f, m, fm, b, fb, eps / 2, simp OF rt, m OF rt, fm OF rt )
                 FI # quad asr # ;

             # Integrate f from a to b using ASR with max error of eps. #
             REAL   fa = f( a );
             REAL   fb = f( b );
             TRIPLE t  = quad simpsons mem( f, a, fa, b, fb );
             quad asr rec( f, a, fa, b, fb, eps, simp OF t, m OF t, fm OF t )
         END # quad asr # ;

    REAL a = 0.0, b = 1.0;
    REAL sinx = quad asr( sin, a, b, 1e-09 );
    print( ( "Simpson's integration of sine from ", fixed( a,    -6, 3 )
           , " to ",                                fixed( b,    -6, 3 )
           , " = ",                                 fixed( sinx, -9, 6 )
           , newline
           )
         )
END
Output:
Simpson's integration of sine from  0.000 to  1.000 =  0.459698

ALGOL W

Translation of: ALGOL 68
which is
Translation of: C
begin % Numerical integration using an Adaptive Simpson's method %

    real procedure quad_asr ( real procedure f; real value a, b, eps ) ;
    begin
        record triple ( real mOf, fmOf, simpOf );

        % "structured" adaptive version, translated from Racket %
        reference(triple) procedure quad_simpsons_mem ( real procedure f
                                                      ; real value a, fa, b, fb
                                                      ) ;
        begin % Evaluates Simpson's Rule, also returning m and f(m) to reuse. %
            real m, fm, simp;
            m    := ( a + b ) / 2;
            fm   := f( m );
            simp := abs ( b - a ) / 6 * ( fa + 4*fm + fb );
            triple( m, fm, simp )
        end quad_simpsons_mem ;

        real procedure quad_asr_rec ( real procedure f
                                    ; real value a, fa, b, fb, eps, whole, m, fm
                                    ) ;
        begin
            % Efficient recursive implementation of adaptive Simpson's rule. %
            % Function values at the start, middle, end of the intervals are retained. %
            reference(triple) lt, rt;
            real              delta;
            lt    := quad_simpsons_mem( f, a, fa, m, fm );
            rt    := quad_simpsons_mem( f, m, fm, b, fb );
            delta := simpOf(lt) + simpOf(rt) - whole;
            if abs delta <= eps * 15
            then simpOf(lt) + simpOf(rt) + delta / 15
            else quad_asr_rec( f, a, fa, m, fm, eps / 2, simpOf(lt), mOf(lt), fmOf(lt) )
               + quad_asr_rec( f, m, fm, b, fb, eps / 2, simpOf(rt), mOf(rt), fmOf(rt) )
        end quad_asr ;

        real              fa, fb;
        reference(triple) t;

        % Integrate f from a to b using ASR with max error of eps. %
        fa := f( a );
        fb := f( b );
        t  := quad_simpsons_mem( f, a, fa, b, fb );
        quad_asr_rec( f, a, fa, b, fb, eps, simpOf(t), mOf(t), fmOf(t) )
    end quad_asr ;

    real a, b, sinx;
    a := 0.0; b := 1.0;
    sinx := quad_asr( sin, a, b, 1'-09 );
    write( s_w := 0, r_format := "A", r_w := 6, r_d := 3
         , "Simpson's integration of sine from ", a, " to ", b, " = "
         , r_w := 12, r_d := 8
         , sinx
         )
end.
Output:
Simpson's integration of sine from  0.000 to  1.000 =   0.45969769

ATS

(* Adaptive Simpson quadrature was, I believe, one of the subjects in
   my undergraduate numerical analysis coursework (in the middle
   1980s). It is most likely the first recursive quadrature method to
   have been published, and it is very easy to write a program that
   does it. *)

#include "share/atspre_staload.hats"

(* The interface is a template function, for any of the floating-point
   types. *)
extern fn {tk : tkind}
quadrature_by_adaptive_simpson
          (* f is an ordinary function (compiled to a C function, as
             opposed to a closure with an environment), with no effect
             restrictions. I make it a function rather than a closure
             so we can use ordinary C functions for "f". I could have
             written "-<1>" or just "->" instead of "-<fun1>" *)
          (f     : g0float tk -<fun1> g0float tk,
           a     : g0float tk,
           b     : g0float tk,
           tol   : g0float tk,  (* Tolerance. *)
           depth : intGte 0) (* The maximum depth of recursion. *)
    (* The type of the return value is "g0float tk", or in other words
       (for purposes of this example program) any of the
       floating-point types "float", "double", or "ldouble" (long
       double), corresponding to the floating-point types of standard
       C. A plain colon means there are no effect restrictions. Had I
       written ":<1>" or ":<fun1>" instead, it would have meant the
       same thing. *)
    : g0float tk

(* The implementation will be for any floating-point type (although it
   can be overridden, either generally or for a specific type). *)
implement {tk}
quadrature_by_adaptive_simpson (f, a, b, tol, depth) =
  let
    (* A shorthand. Sometimes you will see typedefs written instead as
       "stadef", which is a more general notation. *)
    typedef real = g0float tk

    (* In ATS, you can nest functions inside functions to whatever
       depth you like. *)
    fn
    _quad_asr_simpsons (a  : real,
                        fa : real,
                        b  : real,
                        fb : real)
        (* The return type is an unboxed 3-tuple of floating point
           numbers. It will be passed around as a C struct. Usually
           one can leave out the "@" sign, but I like to include it. A
           BOXED tuple (allocated as a struct in the heap, and passed
           around as a pointer) would be written '() instead of @() *)
        : @(real, real, real) =
      let
        (* The type of 0.5 must be casted to type "real" (= "g0float
           tk"). Here I will do that rather explicitly, using a
           macro. Typically it is possible to let the ATS compiler
           infer which type-to-type cast is necessary. See further
           below where I write "g0i2f" to do the casts from integer to
           floating-point type. *)
        macdef one_half = g0float2float<double_kind,tk> 0.5

        val m = one_half * (a + b)
        val fm = f(m)
        val h = b - a
      in
        @(m, fm, (h / g0i2f 6) * (fa + (g0i2f 4 * fm) + fb))
      end

    (* _quad_asr is a recursive function, so I must write "fun"
       instead of "fn". *)
    fun
    _quad_asr {depth : nat}
              (* "depth" is a hypothetical natural number, used only
                 during typechecking. The notations below mean that,
                 on each call to _quad_asr, the value that is of type
                 "int depth" will have to be proven to decrease
                 without passing zero. If this is indeed so, then the
                 function is guaranteed to terminate (unless it would
                 run out of stack space, etc.). *)
              .<depth>.
              (a     : real,
               fa    : real,
               b     : real,
               fb    : real,
               tol   : real,
               whole : real,
               m     : real,
               fm    : real,
               depth : int depth) : g0float tk =
      let
        (* Saying "and" here means that the computations should be
           performed as if at the same time, and perhaps even at the
           same time. In this case that is not necessary, although
           perhaps it is good as documentation. In some instances it
           is REALLY what you want. *)
        val @(lm, flm, left) = _quad_asr_simpsons (a, fa, m, fm)
        and @(rm, frm, right) = _quad_asr_simpsons (m, fm, b, fb)

        val delta = left + right - whole
        and tol_ = g0f2f 0.5 * tol
      in
        (* I test the value of "depth" first, so the typechecker can
           know that, in later branches, the value of depth is greater
           than zero, and thus can be decremented. (The ATS compiler
           will otherwise refuse to let you subtract 1.) I could have
           tested "depth" at the same time as the later tests, by
           using "+" instead of "||" to mean OR, but why that is so is
           getting too technical for this example. *)
        if depth = 0 then
          left + right + (delta / g0i2f 15)
        else if tol_ = tol || abs (delta) <= g0i2f 15 * tol then
          left + right + (delta / g0i2f 15)
        else
          let
            val left_val = _quad_asr (a, fa, m, fm, tol_, left,
                                      lm, flm, depth - 1)
            and right_val = _quad_asr (m, fm, b, fb, tol_, right,
                                       rm, frm, depth - 1)
          in
            left_val + right_val
          end
      end

    (* Here is an alternative to writing "and". Programmers in
       languages such as Python will recognize the technique: *)
    val @(fa, fb) = @(f(a), f(b))

    val @(m, fm, whole) = _quad_asr_simpsons (a, fa, b, fb)
  in
    _quad_asr (a, fa, b, fb, tol, whole, m, fm, depth)
  end

(* We can compile an instance of the template, AND make that instance
   available to BOTH ATS and C. The
   ="ext#quadrature_by_adaptive_simpson_double" tells the compiler to
   use "quadrature_by_adaptive_simpson_double" as the name of the C
   function, rather than some mangled version of that name. *)
extern fn
quadrature_by_adaptive_simpson_double :
  (* The following means "of the same type as
     'quadrature_by_adaptive_simpson' instantiated for the 'double'
     floating-point type". Here it serves simply as a convenient
     shorthand.  *)
  $d2ctype (quadrature_by_adaptive_simpson<double_kind>)
    = "ext#quadrature_by_adaptive_simpson_double"
implement
quadrature_by_adaptive_simpson_double (f, a, b, tol, depth) =
  quadrature_by_adaptive_simpson<double_kind> (f, a, b, tol, depth)

(* Let us call quadrature_by_adaptive_simpson_double from C rather
   than from ATS, simply to demonstrate the process. *)
%{

#include <math.h>

double quadrature_by_adaptive_simpson_double (void *f,
                                              double a, double b,
                                              double tol, int depth);
double integral_of_sine_from_0_to_1 (double tol, int depth);

double
integral_of_sine_from_0_to_1 (double tol, int depth)
{
  return
    quadrature_by_adaptive_simpson_double
      ((void *) sin, 0, 1, tol, depth);
}

%}

(* A "main" that reads no command-line arguments. It must return an
   exit status. *)
implement
main () =
  let
    val tol = 1e-9
    val depth = 1000

    (* One of the ways to make the foreign function call to C. *)
    val result1 = $extfcall (double, "integral_of_sine_from_0_to_1",
                             tol, depth)
    val () = println! ("estimate of ∫ sin x dx from 0 to 1, ",
                       "by call to the C interface: ", result1)

    (* Another way to make a foreign function call. *)
    extern fn sin : double -> double = "mac#sin"

    (* Calling the same compiled function, but this time from ATS
       instead of C. *)
    val result2 =
      quadrature_by_adaptive_simpson_double
        (sin, 0.0, 1.0, tol, depth)
    val () = println! ("estimate of ∫ sin x dx from 0 to 1, ",
                       "by call to a compiled function: ", result2)

    (* Using the template function directly, with single precision
       numbers. *)
    extern fn sinf : float -> float = "mac#sinf"
    val result3 =
      quadrature_by_adaptive_simpson<float_kind>
        (sinf, 0.0f, 1.0f, 1e-5f, depth)
    val () = println! ("estimate of ∫ sin x dx from 0 to 1, ",
                       "by call to the template: ", result3)
  in
    0                           (* Exit status. *)
  end
Output:
$ patscc -std=gnu2x -g -O2 adaptive_simpson_draft_task.dats -lm && ./a.out
estimate of ∫ sin x dx from 0 to 1, by call to the C interface: 0.459698
estimate of ∫ sin x dx from 0 to 1, by call to a compiled function: 0.459698
estimate of ∫ sin x dx from 0 to 1, by call to the template: 0.459698

AWK

BEGIN {
  printf ("estimate of ∫ sin x dx from 0 to 1: %f\n",
          quad_asr(0.0, 1.0, 1.0e-9, 100))
}

function f(x)
{
  return sin(x)
}

function midpoint(a, b)
{
  return 0.5 * (a + b)
}

function simpson_rule(a, fa, b, fb, fm)
{
  return ((b - a) / 6.0) * (fa + (4.0 * fm) + fb)
}

function recursive_simpson(a, fa, b, fb, tol, whole, m, fm, depth,
                           # Local variables:
                           lm, flm, left,
                           rm, frm, right,
                           delta, tol_,
                           leftval, rightval, quadval)
{
  lm = midpoint(a, m)
  flm = f(lm)
  left = simpson_rule(a, fa, m, fm, flm)

  rm = midpoint(m, b)
  frm = f(rm)
  right = simpson_rule(m, fm, b, fb, frm)

  delta = left + right - whole
  tol_ = 0.5 * tol

  if (depth <= 0 || tol_ == tol || abs(delta) <= 15.0 * tol)
    quadval = left + right + (delta / 15.0)
  else
    {
      leftval = recursive_simpson(a, fa, m, fm, tol_,
                                  left, lm, flm, depth - 1)
      rightval = recursive_simpson(m, fm, b, fb, tol_,
                                   right, rm, frm, depth - 1)
      quadval = leftval + rightval
    }
  return quadval
}

function quad_asr(a, b, tol, depth,
                  # Local variables:
                  fa, fb, whole, m, fm)
{
  fa = f(a)
  fb = f(b)

  m = midpoint(a, b)
  fm = f(m)
  whole = simpson_rule(a, fa, b, fb, fm)

  return recursive_simpson(a, fa, b, fb, tol,
                           whole, m, fm, depth)
}

function abs(x)
{
  return (x < 0 ? -x : x)
}
Output:

Just about any modern Awk should work: nawk, mawk, gawk, busybox awk, goawk, etc. But the ancient Old Awk does not work.

$ nawk -f adaptive_simpson_task.awk
estimate of ∫ sin x dx from 0 to 1: 0.459698

C

Translation of: zkl
#include <stdio.h>
#include <math.h>

typedef struct { double m; double fm; double simp; } triple;

/* "structured" adaptive version, translated from Racket */
triple _quad_simpsons_mem(double (*f)(double), double a, double fa, double b, double fb) {
    // Evaluates Simpson's Rule, also returning m and f(m) to reuse.
    double m = (a + b) / 2;
    double fm = f(m);
    double simp = fabs(b - a) / 6 * (fa + 4*fm + fb);
    triple t = {m, fm, simp};
    return t;
}

double _quad_asr(double (*f)(double), double a, double fa, double b, double fb, double eps, double whole, double m, double fm) {
    // Efficient recursive implementation of adaptive Simpson's rule.
    // Function values at the start, middle, end of the intervals are retained.
    triple lt = _quad_simpsons_mem(f, a, fa, m, fm);
    triple rt = _quad_simpsons_mem(f, m, fm, b, fb);
    double delta = lt.simp + rt.simp - whole;
    if (fabs(delta) <= eps * 15) return lt.simp + rt.simp + delta/15;
    return _quad_asr(f, a, fa, m, fm, eps/2, lt.simp, lt.m, lt.fm) +
           _quad_asr(f, m, fm, b, fb, eps/2, rt.simp, rt.m, rt.fm);
}

double quad_asr(double (*f)(double), double a, double b, double eps) {
    // Integrate f from a to b using ASR with max error of eps.
    double fa = f(a);
    double fb = f(b);
    triple t = _quad_simpsons_mem(f, a, fa, b, fb);
    return _quad_asr(f, a, fa, b, fb, eps, t.simp, t.m, t.fm);
}

int main(){
    double a = 0.0, b = 1.0;
    double sinx = quad_asr(sin, a, b, 1e-09);
    printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx);
    return 0;
}
Output:
Simpson's integration of sine from 0 to 1 = 0.459698

COBOL

Works with: GnuCOBOL version 3.1.2.0

Rather than do recursions, I compute the x-intervals iteratively. No attempt is made to avoid recomputing values of sin(x).

       *> -*- mode: cobol -*-  
       *> 
       *> I do the computations (except perhaps for the call to the
       *> "sin" function) in decimal fixed point.
       *>  

       identification division.
       program-id. adaptive_simpson_task.

       data division.
       working-storage section.
       01  numer picture s9(25).  *> 25 decimal digits.
       01  denom picture s9(25).
       01  numer2 picture s9(25).
       01  a picture s9(5)V9(20). *> 5+20 decimal digits.
       01  b picture s9(5)V9(20).
       01  tol picture s9(5)V9(20).
       01  x picture s9(5)V9(20).
       01  y picture s9(5)V9(20).
       01  x0 picture s9(5)V9(20).
       01  y0 picture s9(5)V9(20).
       01  x1 picture s9(5)V9(20).
       01  y1 picture s9(5)V9(20).
       01  x2 picture s9(5)V9(20).
       01  y2 picture s9(5)V9(20).
       01  x3 picture s9(5)V9(20).
       01  y3 picture s9(5)V9(20).
       01  x4 picture s9(5)V9(20).
       01  y4 picture s9(5)V9(20).
       01  ruleval picture s9(5)V9(20).
       01  ruleval0 picture s9(5)V9(20).
       01  ruleval1 picture s9(5)V9(20).
       01  delta picture s9(5)V9(20).
       01  abs-delta picture 9(5)V9(20). *> unsigned.
       01  tol0 picture s9(5)V9(20).
       01  tol1 picture s9(5)V9(20).
       01  delta15 picture s9(5)V9(20).
       01  stepval picture s9(5)V9(20).
       01  quadval picture s9(5)V9(20).
       01  result picture -9.9(9).
 
       procedure division.
           move 0.0 to a
           move 1.0 to b
           move 0.000000001 to tol
           perform adaptive-quadrature
           move quadval to result
           display result
           stop run.

       adaptive-quadrature.
           move 0 to numer
           move 1 to denom
           perform until (numer = denom) and (denom = 1)
               perform get-interval
               perform simpson-rule-thrice
               compute delta = ruleval0 + ruleval1 - ruleval
               compute abs-delta = delta
               compute tol0 = tol * (x4 - x0)
               compute tol1 = tol * (x2 - x0)
               if (tol1 = tol0) or (abs-delta <= 15 * tol0) then
                   compute delta15 = delta / 15
                   compute stepval = ruleval0 + ruleval1 + delta15
                   compute quadval = quadval + stepval
                   perform go-to-next-interval
               else
                   perform bisect-current-interval
               end-if
           end-perform.

       simpson-rule-thrice.
           *> There is no attempt here to minimize the number of
           *> recomputations of sin(x).

           compute x2 = (x0 + x4) / 2
           compute x1 = (x0 + x2) / 2
           compute x3 = (x2 + x4) / 2

           compute x = x0
           perform compute-function
           compute y0 = y

           compute x = x1
           perform compute-function
           compute y1 = y

           compute x = x2
           perform compute-function
           compute y2 = y

           compute x = x3
           perform compute-function
           compute y3 = y

           compute x = x4
           perform compute-function
           compute y4 = y

           compute ruleval = ((x4 - x0) / 6) * (y0 + (4 * y2) + y4)
           compute ruleval0 = ((x2 - x0) / 6) * (y0 + (4 * y1) + y2)
           compute ruleval1 = ((x4 - x2) / 6) * (y2 + (4 * y3) + y4).

       bisect-current-interval.
           compute numer = numer + numer
           compute denom = denom + denom.

       go-to-next-interval.
           compute numer = numer + 1
           divide numer by 2 giving numer2
           if numer2 + numer2 = numer then
               perform until not (numer2 + numer2 = numer)
                   move numer2 to numer
                   divide denom by 2 giving denom
                   divide numer by 2 giving numer2
               end-perform
           end-if.

       get-interval.
           compute x0 = numer / denom
           compute x0 = (a * (1 - x0)) + (b * x0)
           compute x4 = (numer + 1) / denom
           compute x4 = (a * (1 - x4)) + (b * x4).

       compute-function.
           compute y = function sin (x).
Output:
$ cobc -std=cobol2014 -x adaptive_simpson_task.cob && ./adaptive_simpson_task
 0.459697694

Common Lisp

(defun %%quad-asr-simpsons (f a fa b fb)
  (let* ((m (/ (+ a b) 2))
         (fm (funcall f m))
         (h (- b a)))
    ;; Common Lisp supports returning multiple values at once. There
    ;; is no need to return them explicitly as a data structure
    ;; (though that also could be done).
    (values m fm (* (/ h 6) (+ fa (* 4 fm) fb)))))

(defun %%quad-asr (f a fa b fb tol whole m fm depth)
  (multiple-value-bind (lm flm left)
      (%%quad-asr-simpsons f a fa m fm)
    (multiple-value-bind (rm frm right)
        (%%quad-asr-simpsons f m fm b fb)
      (let ((delta (- (+ left right) whole))
            (tol^ (/ tol 2)))
        (if (or (<= depth 0)
                (= tol^ tol)
                (<= (abs delta) (* 15 tol)))
            (+ left right (/ delta 15))
            (+ (%%quad-asr f a fa m fm tol^ left
                           lm flm (1- depth))
               (%%quad-asr f m fm b fb tol^ right
                           rm frm (1- depth))))))))

(defun quad-asr (f a b tol depth)
  (let ((fa (funcall f a))
        (fb (funcall f b)))
    (multiple-value-bind (m fm whole)
        (%%quad-asr-simpsons f a fa b fb)
      (%%quad-asr f a fa b fb tol whole m fm depth))))

(princ "estimated definite integral of sin(x) for x from 0 to 1: ")
(princ (quad-asr #'sin 0 1 1e-9 1000))
(terpri)
Output:
$ clisp adaptive_simpson_task.lisp
estimated definite integral of sin(x) for x from 0 to 1: 0.45969772

D

import std.math;
import std.stdio;
import std.typecons;

double
quad_asr (double function (double) f, double a, double b,
          double tol, int depth)
{
   auto quad_asr_simpsons_ (double a, double fa, double b, double fb)
   {
     auto m = 0.5 * (a + b);
     auto fm = f(m);
     auto h = b - a;
     return tuple!("m", "fm", "quadval")
       (m, fm, (h / 6) * (fa + 4*fm + fb));
   }

   double quad_asr_ (double a, double fa, double b, double fb,
                     double tol, double whole, double m, double fm,
                     int depth)
   {
     //
     // Please do not ask why I do not use multiple return
     // statements. Finding reasons why is left as an exercise.
     //
     auto retval = NaN (0);
     auto left = quad_asr_simpsons_ (a, fa, m, fm);
     auto right = quad_asr_simpsons_ (m, fm, b, fb);
     auto delta = left.quadval + right.quadval - whole;
     auto tol_ = 0.5 * tol;
     if (depth <= 0 || tol_ == tol || abs (delta) <= 15 * tol)
       retval = left.quadval + right.quadval + (delta / 15);
     else
       retval = (quad_asr_ (a, fa, m, fm, tol_, left.quadval,
                           left.m, left.fm, depth - 1)
                 + quad_asr_ (m, fm, b, fb, tol_, right.quadval,
                              right.m, right.fm, depth - 1));
     return retval;
   }

   auto fa = f(a);
   auto fb = f(b);
   auto whole = quad_asr_simpsons_ (a, fa, b, fb);
   return quad_asr_ (a, fa, b, fb, tol, whole.quadval,
                     whole.m, whole.fm, depth);
}

int
main ()
{
  auto result = quad_asr (&sin, 0.0, 1.0, 1e-9, 1000);
  printf ("estimate of ∫ sin x dx from 0 to 1: %lf\n", result);
  return 0;
}
Output:
$ gdc -g adaptive_simpson_task.d && ./a.out
estimate of ∫ sin x dx from 0 to 1: 0.459698

Factor

Translation of: Julia
USING: formatting kernel locals math math.functions math.ranges
sequences ;
IN: rosetta-code.simpsons

:: simps ( f a b n -- x )
    n even?
    [ n "n must be even; %d was given" sprintf throw ] unless
    b a - n / :> h
    1 n 2 <range> 2 n 1 - 2 <range>
    [ [ a + h * f call ] map-sum ] bi@ [ 4 ] [ 2 ] bi*
    [ * ] 2bi@ a b [ f call ] bi@ + + + h 3 / * ; inline

[ sin ] 0 1 100 simps
"Simpson's rule integration of sin from 0 to 1 is: %u\n" printf
Output:
Simpson's rule integration of sin from 0 to 1 is: 0.4596976941573994

Forth

Works with: Gforth version 0.7.3
\ -*- mode: forth -*-
\
\ The following was written for Gforth, using its implementation of
\ local variables and of f=
\

Defer f

: simpson-rule { F: a F: fa F: b F: fb -- m fm quadval }
    a b f+ 0.5e0 f*                     \ -- m
    fdup f                              \ -- m fm
    fdup 4.0e0 f* fa f+ fb f+
    b a f- 6.0e0 f/ f*                  \ -- m fm quadval
;    

: recursive-simpson
    { F: a F: fa F: b F: fb F: tol F: whole F: m F: fm depth -- quadval }
    a fa m fm simpson-rule { F: lm F: flm F: left }
    m fm b fb simpson-rule { F: rm F: frm F: right }
    left right f+ whole f- { F: delta }
    tol 0.5e0 f*           { F: tol_ }
    depth 0 <=
    tol_ tol f= or
    delta fabs 15.0e0 tol f* f<= or
    if
        left right f+ delta 15.0e0 f/ f+
    else
        a fa m fm tol_ left lm flm depth 1 - recurse
        m fm b fb tol_ right rm frm depth 1 - recurse
        f+
    then
;

: quad-asr { F: a F: b F: tol depth -- quadval }
    a f { F: fa }
    b f { F: fb }
    a fa b fb simpson-rule { F: m F: fm F: whole }
    a fa b fb tol whole m fm depth recursive-simpson
;

:noname fsin ; IS f

." estimate of ∫ sin x dx from 0 to 1: "
0.0e0 1.0e0 1.0e-12 100 quad-asr f.
cr
Output:
$ gforth adaptive_simpson_task.fs -e bye
estimate of ∫ sin x dx from 0 to 1: 0.45969769413186

Fortran

program adaptive_simpson_task
  implicit none

  integer, parameter :: dbl = kind (1.0d0)

  interface
     function function_dbl_to_dbl (x) result (y)
       import dbl
       real(dbl), intent(in) :: x
       real(dbl) :: y
     end function function_dbl_to_dbl
  end interface

  write (*, 1000) quad_asr (sine, 0.0_dbl, 1.0_dbl, 1E-9_dbl, 1000)
1000 format ("estimated definite integral of sin(x) ", &
       &  "for x from 0 to 1: ", F15.13)

contains

  function sine (x) result (y)
    real(dbl), intent(in) :: x
    real(dbl) :: y

    y = sin (x)
  end function sine

  function quad_asr (f, a, b, tol, depth) result (quadval)
    procedure(function_dbl_to_dbl) :: f
    real(dbl), intent(in) :: a, b, tol
    integer, intent(in) :: depth
    real(dbl) :: quadval

    real(dbl) :: fa, fb, m, fm, whole

    fa = f(a)
    fb = f(b)
    call quad_asr_simpsons_ (f, a, fa, b, fb, m, fm, whole)
    quadval = quad_asr_ (f, a, fa, b, fb, tol, whole, m, fm, depth)
  end function quad_asr

  recursive function quad_asr_ (f, a, fa, b, fb, tol, whole, &
       &                        m, fm, depth) result (quadval)
    procedure(function_dbl_to_dbl) :: f
    real(dbl), intent(in) :: a, fa, b, fb, tol, whole, m, fm
    integer, intent(in) :: depth
    real(dbl) :: quadval

    real(dbl) :: lm, flm, left
    real(dbl) :: rm, frm, right
    real(dbl) :: delta, tol_

    call quad_asr_simpsons_ (f, a, fa, m, fm, lm, flm, left)
    call quad_asr_simpsons_ (f, m, fm, b, fb, rm, frm, right)
    delta = left + right - whole
    tol_ = tol / 2
    if (depth <= 0 .or. tol_ == tol .or. abs (delta) <= 15 * tol) then
       quadval = left + right + (delta / 15)
    else
       quadval = quad_asr_ (f, a, fa, m, fm, tol_, left, &
            &               lm, flm, depth - 1)
       quadval = quadval + quad_asr_ (f, m, fm, b, fb, tol_, &
            &                         right, rm, frm, depth - 1)
    end if
  end function quad_asr_

  subroutine quad_asr_simpsons_ (f, a, fa, b, fb, m, fm, quadval)
    procedure(function_dbl_to_dbl) :: f
    real(dbl), intent(in) :: a, fa, b, fb
    real(dbl), intent(out) :: m, fm, quadval

    m = (a + b) / 2
    fm = f(m)
    quadval = ((b - a) / 6) * (fa + (4 * fm) + fb)
  end subroutine quad_asr_simpsons_

end program adaptive_simpson_task
Output:
$ gfortran -O -std=f2018 -fbounds-check -g adaptive_simpson_task.f90 && ./a.out
estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941318

On x86-64, if you leave out the optimizer you get a trampoline:

$ gfortran -std=f2018 -fbounds-check -g adaptive_simpson_task.f90 && ./a.out
adaptive_simpson_task.f90:20:2:

   20 |   function sine (x) result (y)
      |  ^
Warning: trampoline generated for nested function ‘sine’ [-Wtrampolines]
/usr/lib/gcc/x86_64-pc-linux-gnu/12/../../../../x86_64-pc-linux-gnu/bin/ld: warning: /tmp/cc2mQOD2.o: requires executable stack (because the .note.GNU-stack section is executable)
estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941318

Such a bit of executable code on the stack often is not allowed on systems hardened against hackers. Otherwise you probably needn't be concerned, even though warnings are being generated both by the compiler and the linker. :)

Go

Like the zkl entry, this is also a translation of the Python code in the Wikipedia article.

package main

import (
    "fmt"
    "math"
)

type F = func(float64) float64

/* "structured" adaptive version, translated from Racket */
func quadSimpsonsMem(f F, a, fa, b, fb float64) (m, fm, simp float64) {
    // Evaluates Simpson's Rule, also returning m and f(m) to reuse.
    m = (a + b) / 2
    fm = f(m)
    simp = math.Abs(b-a) / 6 * (fa + 4*fm + fb)
    return
}

func quadAsrRec(f F, a, fa, b, fb, eps, whole, m, fm float64) float64 {
    // Efficient recursive implementation of adaptive Simpson's rule.
    // Function values at the start, middle, end of the intervals are retained.
    lm, flm, left := quadSimpsonsMem(f, a, fa, m, fm)
    rm, frm, right := quadSimpsonsMem(f, m, fm, b, fb)
    delta := left + right - whole
    if math.Abs(delta) <= eps*15 {
        return left + right + delta/15
    }
    return quadAsrRec(f, a, fa, m, fm, eps/2, left, lm, flm) +
        quadAsrRec(f, m, fm, b, fb, eps/2, right, rm, frm)
}

func quadAsr(f F, a, b, eps float64) float64 {
    // Integrate f from a to b using ASR with max error of eps.
    fa, fb := f(a), f(b)
    m, fm, whole := quadSimpsonsMem(f, a, fa, b, fb)
    return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm)
}

func main() {
    a, b := 0.0, 1.0
    sinx := quadAsr(math.Sin, a, b, 1e-09)
    fmt.Printf("Simpson's integration of sine from %g to %g = %f\n", a, b, sinx)   
}
Output:
Simpson's integration of sine from 0 to 1 = 0.459698

Icon

record simpson_result (m, fm, quadval)

procedure main ()
  write ("estimate of ∫ sin x dx from 0 to 1: ",
         quad_asr (sin, 0.0, 1.0, 1E-9, 100))
end

procedure quad_asr(f, a, b, tol, depth)
  local fa, fb, whole

  fa := f(a)
  fb := f(b)
  whole := simpson_rule (f, a, fa, b, fb)
  return recursive_simpson(f, a, fa, b, fb, tol, whole.quadval,
                           whole.m, whole.fm, depth)
end

procedure recursive_simpson (f, a, fa, b, fb, tol, whole,
                             m, fm, depth)
  local left, right, delta, tol2
  local leftval, rightval, quadval

  left := simpson_rule (f, a, fa, m, fm)
  right := simpson_rule (f, m, fm, b, fb)
  delta := left.quadval + right.quadval - whole
  tol2 := 0.5 * tol
  if depth <= 0 | tol2 = tol | abs (delta) <= 15 * tol then
    quadval := left.quadval + right.quadval + (delta / 15)
  else
  {
    leftval := recursive_simpson (f, a, fa, m, fm, tol2,
                                  left.quadval, left.m, left.fm,
                                  depth - 1)
    rightval := recursive_simpson (f, m, fm, b, fb, tol2,
                                   right.quadval, right.m, right.fm,
                                   depth - 1)
    quadval := leftval + rightval
  }
  return quadval
end

procedure simpson_rule (f, a, fa, b, fb)
  local m, fm, quadval

  m := 0.5 * (a + b)
  fm := f(m)
  quadval := ((b - a) / 6) * (fa + (4 * fm) + fb)
  return simpson_result (m, fm, quadval)
end
Output:
$ icon adaptive_simpson_task_arizona.icn
estimate of ∫ sin x dx from 0 to 1: 0.4596976941

J

Typically one would choose the library implementation:

   load'~addons/math/misc/integrat.ijs'

   NB. integrate returns definite integral and estimated digits of accuracy
   1&o. integrate 0 1
0.459698 9

   NB. adapt implements adaptive Simpson's rule, however recomputes some integrands
   1&o. adapt 0 1 1e_9
0.459698
   


Note'expected answer computed by j www.jsoftware.com'

       1-&:(1&o.d._1)0
    0.459698

    translated from c
)

mp=: +/ .*  NB. matrix product

NB. Evaluates Simpson's Rule, also returning m and f(m) to reuse.
uquad_simpsons_mem=: adverb define
 'a fa b fb'=. y
 em=. a ([ + [: -: -~) b
 fm=. u em
 simp=. ((| b - a) % 6) * 1 4 1 mp fa , fm , fb
 em, fm, simp
)

Simp=: 1 :'2{m'
Fm=: 1 :'1{m'
M=: 1 :'0{m'

NB. Efficient recursive implementation of adaptive Simpson's rule.	    
NB. Function values at the start, middle, end of the intervals are retained.
uquad_asr=: adverb define
 'a fa b fb eps whole em fm'=. y
 lt=. u uquad_simpsons_mem(a, fa, em, fm)
 rt=. u uquad_simpsons_mem(em, fm, b, fb)
 delta=. lt Simp + rt Simp - whole
 if. (| delta) <: eps * 15 do.
  lt Simp + rt Simp + delta % 15
 else.
  (a, fa, em, fm, (-: eps), lt Simp, lt M, lt Fm) +&(u uquad_asr) (em, fm, b, fb, (-: eps), rt Simp, rt M, rt Fm)
 end.
)

NB. Integrate u from a to b using ASR with max error of eps.
quad_asr=: adverb define
 'a b eps'=. y
 fa=. u a
 fb=. u b
 t=. u uquad_simpsons_mem a, fa, b, fb
 u uquad_asr a, fa, b, fb, eps, t Simp, t M, t Fm
)
   echo 'Simpson''s integration of sine from 0 to 1 = ' , ": 1&o. quad_asr 0 1 1e_9
Simpson's integration of sine from 0 to 1 = 0.459698

Java

import java.util.function.Function;

public class NumericalIntegrationAdaptiveSimpsons {

    public static void main(String[] args) {
        Function<Double,Double> f = x -> sin(x);
        System.out.printf("integrate sin(x), x = 0 .. Pi = %2.12f.  Function calls = %d%n", quadratureAdaptiveSimpsons(f, 0, Math.PI, 1e-8), functionCount);
        functionCount = 0;
        System.out.printf("integrate sin(x), x = 0 .. 1 = %2.12f.  Function calls = %d%n", quadratureAdaptiveSimpsons(f, 0, 1, 1e-8), functionCount);
    }
    
    private static double quadratureAdaptiveSimpsons(Function<Double,Double> function, double a, double b, double error) {
        double fa = function.apply(a);
        double fb = function.apply(b);
        Triple t =  quadratureAdaptiveSimpsonsOne(function, a, fa, b ,fb);
        return quadratureAdaptiveSimpsonsRecursive(function, a, fa, b, fb, error, t.s, t.x, t.fx);
    }
    
    private static double quadratureAdaptiveSimpsonsRecursive(Function<Double,Double> function, double a, double fa, double b, double fb, double error, double whole, double m, double fm) {
        Triple left  = quadratureAdaptiveSimpsonsOne(function, a, fa, m, fm);
        Triple right = quadratureAdaptiveSimpsonsOne(function, m, fm, b, fb);
        double delta = left.s + right.s - whole;
        if ( Math.abs(delta) <= 15*error ) {
            return left.s + right.s + delta / 15;
        }
        return quadratureAdaptiveSimpsonsRecursive(function, a, fa, m, fm, error/2, left.s, left.x, left.fx) +
               quadratureAdaptiveSimpsonsRecursive(function, m, fm, b, fb, error/2, right.s, right.x, right.fx);
    }
    
    private static Triple quadratureAdaptiveSimpsonsOne(Function<Double,Double> function, double a, double fa, double b, double fb) {
        double m = (a + b) / 2;
        double fm = function.apply(m);
        return new Triple(m, fm, Math.abs(b-a) / 6 * (fa + 4*fm + fb));
    }
    
    private static class Triple {
        double x, fx, s;
        private Triple(double m, double fm, double s) {
            this.x = m;
            this.fx = fm;
            this.s = s;
        }
    }
    
    private static int functionCount = 0;
    
    private static double sin(double x) {
        functionCount++;
        return Math.sin(x);
    }
    
}
Output:
integrate sin(x), x = 0 .. Pi = 1.999999999998.  Function calls = 121
integrate sin(x), x = 0 .. 1 = 0.459697694131.  Function calls = 33

Julia

Originally from Modesto Mas, https://mmas.github.io/simpson-integration-julia

function simps(f::Function, a::Number, b::Number, n::Number)
    iseven(n) || throw("n must be even, and $n was given")
    h = (b-a)/n
    s = f(a) + f(b)
    s += 4 * sum(f.(a .+ collect(1:2:n) .* h))
    s += 2 * sum(f.(a .+ collect(2:2:n-1) .* h))
    h/3 * s
end

println("Simpson's rule integration of sin from 0 to 1 is: ",  simps(sin, 0.0, 1.0, 100))
Output:
Simpson's rule integration of sin from 0 to 1 is: 0.45969769415739936

K

Works with: ngn/k
simpsonsRule:{[f;a;b] m: (a+b)%2; h: b-a; m, (h%6)*f[a]+f[b]+4*f[m]}

recursiveSimpson:{[f;a;b;tol;whole;m;depth]
 (lm;left): simpsonsRule[f;a;m]
 (rm;right): simpsonsRule[f;m;b]
 delta: left+right-whole
 tol2: tol%2
 $[(1>depth)|(tol2=tol)|~(15*tol)<delta|-delta
  left+right+delta%15
  recursiveSimpson[f;a;m;tol2;left;lm;depth-1]+recursiveSimpson[f;m;b;tol2;right;rm;depth-1]]}

quadAsr:{[f;a;b;tol;depth]; (m;whole): simpsonsRule[f;a;b]; recursiveSimpson[f;a;b;tol;whole;m;depth]}
quadAsr[`sin;0;1;1e-12;100]
0.45969769413186023

Kotlin

Translation of: Go
// Version 1.2.71

import kotlin.math.abs
import kotlin.math.sin

typealias F = (Double) -> Double
typealias T = Triple<Double, Double, Double>

/* "structured" adaptive version, translated from Racket */
fun quadSimpsonsMem(f: F, a: Double, fa: Double, b: Double, fb: Double): T {
    // Evaluates Simpson's Rule, also returning m and f(m) to reuse
    val m = (a + b) / 2
    val fm = f(m)
    val simp = abs(b - a) / 6 * (fa + 4 * fm + fb)
    return T(m, fm, simp)
}

fun quadAsrRec(f: F, a: Double, fa: Double, b: Double, fb: Double,
    eps: Double, whole: Double, m: Double, fm: Double): Double {
    // Efficient recursive implementation of adaptive Simpson's rule.
    // Function values at the start, middle, end of the intervals are retained. 
    val (lm, flm, left) = quadSimpsonsMem(f, a, fa, m, fm)
    val (rm, frm, right) = quadSimpsonsMem(f, m, fm, b, fb)
    val delta = left + right - whole
    if (abs(delta) <= eps * 15)  return left + right + delta / 15
    return quadAsrRec(f, a, fa, m, fm, eps / 2, left, lm, flm) +
        quadAsrRec(f, m, fm, b, fb, eps / 2, right, rm, frm)
}

fun quadAsr(f: F, a: Double, b: Double, eps: Double): Double {
    // Integrate f from a to b using ASR with max error of eps.
    val fa = f(a)
    val fb = f(b)
    val (m, fm, whole) = quadSimpsonsMem(f, a, fa, b, fb)
    return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm)
}

fun main(args: Array<String>) {
    val a = 0.0
    val b = 1.0
    val sinx = quadAsr(::sin, a, b, 1.0e-09)
    println("Simpson's integration of sine from $a to $b = ${"%6f".format(sinx)}")
}
Output:
Simpson's integration of sine from 0.0 to 1.0 = 0.459698

Lambdatalk

1) using the sin's primitive -cos(x)

01 sin(x) = {- {cos 0} {cos 1}} 
-> 0.45969769413186023   // will be the reference value

2) using the adaptative simpson method

{def adasimp

 {def adasimp.calc
  {lambda {:fun :a :fa :b :fb}
   {let { {:delta {- :b :a}}
          {:m {/ {+ :a :b} 2}}
          {:fa :fa} {:fb :fb}
          {:fm {:fun {/ {+ :a :b} 2}}}
        } {A.new :m :fm {/ {* {abs :delta} {+ :fa {* 4 :fm} :fb}} 6}} }}}

 {def adasimp.rec
  {lambda {:fun :a :fa :b :fb :eps :mfs}
   {let { {:fun :fun}
          {:a :a} {:fa :fa} {:b :b} {:fb :fb}
          {:eps :eps} {:mfs :mfs}
          {:left  {adasimp.calc :fun :a :fa {A.get 0 :mfs} {A.get 1 :mfs}}}
          {:right {adasimp.calc :fun {A.get 0 :mfs} {A.get 1 :mfs} :b :fb}}
        } {let { {:fun :fun}
                 {:a :a} {:fa :fa} {:b :b} {:fb :fb}
                 {:eps :eps} {:mfs :mfs}
                 {:left :left} {:right :right}
                 {:delta {+ {A.get 2 :left}
                            {A.get 2 :right}
                         {- {A.get 2 :mfs}}}}
               } {if {<= {abs :delta} {* :eps 15}}
                  then {+ {A.get 2 :left}
                          {A.get 2 :right}
                          {/ :delta 15}}
                  else {+ {adasimp.rec :fun :a :fa 
                                       {A.get 0 :mfs} {A.get 1 :mfs}
                                       {/ :eps 2} :left }
                          {adasimp.rec :fun 
                                       {A.get 0 :mfs} {A.get 1 :mfs}
                                       :b :fb {/ :eps 2} :right }}} }}}}
 
 {lambda {:fun :a :b :eps}
  {adasimp.rec :fun :a {:fun :a} :b {:fun :b} :eps 
               {adasimp.calc :fun :a {:fun :a} :b {:fun :b}}} }}
-> adasimp

{adasimp sin 0 1 1.0e-11}
-> 0.45969769413186023  
  (0.45969769413186023)  // reference value

Mercury

%%% -*- mode: mercury; prolog-indent-width: 2; -*-

:- module adaptive_simpson_task_mercury.

:- interface.
:- import_module io.
:- pred main(io::di, io::uo) is det.

:- implementation.
:- import_module float.
:- import_module int.
:- import_module math.
:- import_module string.

:- pred simpson_rule((func(float) = float)::in,
                     float::in, float::in, float::in, float::in,
                     float::out, float::out, float::out) is det.
simpson_rule(F, A, FA, B, FB, M, FM, QuadVal) :-
  FM = F(M), (M = 0.5 * (A + B)),
  (QuadVal = ((B - A) / 6.0) * (FA + (4.0 * FM) + FB)).

:- func recursive_simpson(func(float) = float,
                          float, float, float, float,
                          float, float, float, float, int) = float.
recursive_simpson(F, A, FA, B, FB, Tol,
                  Whole, M, FM, Depth) = QuadVal :-
  simpson_rule(F, A, FA, M, FM, LM, FLM, Left),
  simpson_rule(F, M, FM, B, FB, RM, FRM, Right),
  (Left + Right - Whole = Delta),
  (0.5 * Tol = Tol_),
  (if (Depth =< 0; Tol_ = Tol; abs(Delta) =< 15.0 * Tol)
   then (QuadVal = Left + Right + (Delta / 15.0))
   else (QuadVal = recursive_simpson(F, A, FA, M, FM, Tol_,
                                     Left, LM, FLM, Depth - 1)
                   + recursive_simpson(F, M, FM, B, FB, Tol_,
                                       Right, RM, FRM, Depth - 1))).

:- func quad_asr(func(float) = float,
                 float, float, float, int) = float.
quad_asr(F, A, B, Tol, Depth) = QuadVal :-
  F(A) = FA, F(B) = FB,
  simpson_rule(F, A, FA, B, FB, M, FM, Whole),
  (QuadVal = recursive_simpson(F, A, FA, B, FB, Tol,
                               Whole, M, FM, Depth)).

main(!IO) :-
  print("estimate of ∫ sin x dx from 0 to 1: ", !IO),
  (QuadVal = quad_asr(sin, 0.0, 1.0, 1.0e-9, 100)),
  print(from_float(QuadVal), !IO),
  nl(!IO).

:- end_module adaptive_simpson_task_mercury.
Output:
$ mmc --no-verbose-make --make --use-subdirs adaptive_simpson_task_mercury && ./adaptive_simpson_task_mercury
estimate of ∫ sin x dx from 0 to 1: 0.4596976941317858

Modula-2

Works with: GCC version 13.1.0
MODULE adaptive_simpson_task_Modula2;

(* ISO Modula-2 libraries. *)
IMPORT RealMath;
IMPORT SRealIO;
IMPORT STextIO;

TYPE functionReal2Real = PROCEDURE (REAL) : REAL;

PROCEDURE simpsonRule (f : functionReal2Real;
                       a, fa, b, fb : REAL;
                       VAR m, fm, quadVal : REAL);
BEGIN
  m := 0.5 * (a + b);
  fm := f(m);
  quadVal := ((b - a) / 6.0) * (fa + (4.0 * fm) + fb)
END simpsonRule;

PROCEDURE recursiveSimpson (f : functionReal2Real;
                            a, fa, b, fb, tol, whole, m, fm : REAL;
                            depth : INTEGER) : REAL;
  VAR lm, flm, left : REAL;
      rm, frm, right : REAL;
      delta, tol2, quadVal : REAL;
BEGIN
  simpsonRule (f, a, fa, m, fm, lm, flm, left);
  simpsonRule (f, m, fm, b, fb, rm, frm, right);
  delta := left + right - whole;
  tol2 := 0.5 * tol;
  IF (depth <= 0) OR (tol2 = tol) OR (ABS (delta) <= 15.0 * tol) THEN
    quadVal := left + right + (delta / 15.0)
  ELSE
    quadVal := (recursiveSimpson (f, a, fa, m, fm, tol2,
                                  left, lm, flm, depth - 1)
                + recursiveSimpson (f, m, fm, b, fb, tol2,
                                    right, rm, frm, depth - 1))
  END;
  RETURN quadVal;
END recursiveSimpson;

PROCEDURE quadASR (f : functionReal2Real;
                   a, b, tol : REAL;
                   depth : INTEGER) : REAL;
  VAR fa, fb, m, fm, whole : REAL;
BEGIN
  fa := f(a);  fb := f(b);
  simpsonRule (f, a, fa, b, fb, m, fm, whole);
  RETURN recursiveSimpson (f, a, fa, b, fb, tol,
                           whole, m, fm, depth)
END quadASR;

BEGIN
  STextIO.WriteString ('estimate of ∫ sin x dx from 0 to 1: ');
  SRealIO.WriteReal (quadASR (RealMath.sin, 0.0, 1.0,
                              0.000000001, 100),
                     10);
  STextIO.WriteLn;
END adaptive_simpson_task_Modula2.
Output:
$ gm2 -fiso adaptive_simpson_task_Modula2.mod && ./a.out
estimate of ∫ sin x dx from 0 to 1: 0.45969769

Nim

Direct translation of Python code from Wikipedia entry.

import math, sugar

type Func = float -> float

proc quadSimpsonsMem(f: Func; a, fa, b, fb: float): tuple[m, fm, val: float] =
  ## Evaluates the Simpson's Rule, also returning m and f(m) to reuse
  result.m = (a + b) / 2
  result.fm = f(result.m)
  result.val = abs(b - a) * (fa + 4 * result.fm + fb) / 6

proc quadAsr(f: Func; a, fa, b, fb, eps, whole, m, fm: float): float =
  ## Efficient recursive implementation of adaptive Simpson's rule.
  ## Function values at the start, middle, end of the intervals are retained.
  let (lm, flm, left) = f.quadSimpsonsMem(a, fa, m, fm)
  let (rm, frm, right) = f.quadSimpsonsMem(m, fm, b, fb)
  let delta = left + right - whole
  result = if abs(delta) <= 15 * eps:
             left + right + delta / 15
           else:
             f.quadAsr(a, fa, m, fm, eps / 2, left, lm, flm) +
             f.quadAsr(m, fm, b, fb, eps / 2, right, rm, frm)

proc quadAsr(f: Func; a, b, eps: float): float =
  ## Integrate f from a to b using Adaptive Simpson's Rule with max error of eps.
  let fa = f(a)
  let fb = f(b)
  let (m, fm, whole) = f.quadSimpsonsMem(a, fa, b, fb)
  result = f.quadAsr(a, fa, b, fb, eps, whole, m, fm)

echo "Simpson's integration of sine from 0 to 1 = ", sin.quadAsr(0, 1, 1e-9)
Output:
Simpson's integration of sine from 0 to 1 = 0.4596976941317859

Oberon-2

Translation of: Modula-2
Works with: Oxford Oberon-2 Compiler version 3.3.0
(* -*- mode: Modula2 -*- *)

MODULE adaptiveSimpsonTask;

IMPORT Math, Out;

TYPE functionReal2Real = PROCEDURE (x : REAL) : REAL;

PROCEDURE simpsonRule (f : functionReal2Real;
                       a, fa, b, fb : REAL;
                       VAR m, fm, quadVal : REAL);
BEGIN
  m := 0.5 * (a + b);
  fm := f(m);
  quadVal := ((b - a) / 6.0) * (fa + (4.0 * fm) + fb)
END simpsonRule;

PROCEDURE recursiveSimpson (f : functionReal2Real;
                            a, fa, b, fb, tol, whole, m, fm : REAL;
                            depth : INTEGER) : REAL;
  VAR lm, flm, left : REAL;
      rm, frm, right : REAL;
      delta, tol2, quadVal : REAL;
BEGIN
  simpsonRule (f, a, fa, m, fm, lm, flm, left);
  simpsonRule (f, m, fm, b, fb, rm, frm, right);
  delta := left + right - whole;
  tol2 := 0.5 * tol;
  IF (depth <= 0) OR (tol2 = tol) OR (ABS (delta) <= 15.0 * tol) THEN
    quadVal := left + right + (delta / 15.0)
  ELSE
    quadVal := (recursiveSimpson (f, a, fa, m, fm, tol2,
                                  left, lm, flm, depth - 1)
                + recursiveSimpson (f, m, fm, b, fb, tol2,
                                    right, rm, frm, depth - 1))
  END;
  RETURN quadVal
END recursiveSimpson;

PROCEDURE quadASR (f : functionReal2Real;
                   a, b, tol : REAL;
                   depth : INTEGER) : REAL;
  VAR fa, fb, m, fm, whole : REAL;
BEGIN
  fa := f(a);  fb := f(b);
  simpsonRule (f, a, fa, b, fb, m, fm, whole);
  RETURN recursiveSimpson (f, a, fa, b, fb, tol,
                           whole, m, fm, depth)
END quadASR;

BEGIN
  Out.String ('estimate of ∫ sin x dx from 0 to 1: ');
  Out.Real (quadASR (Math.Sin, 0.0, 1.0, 0.000001, 100), 7);
  Out.Ln
END adaptiveSimpsonTask.
Output:
$ obc adaptiveSimpsonTask.Mod && ./a.out
estimate of ∫ sin x dx from 0 to 1: 0.459698

ObjectIcon

$encoding UTF-8

import
  io(write),
  util(Math)

procedure main ()
  write (u"estimate of ∫ sin x dx from 0 to 1: ",
         ASR.quadrature (Math.sin, 0.0, 1.0, 1E-9, 100))
end

class Simpson_rule ()
  public m, fm, quadval

  public new (f, a, fa, b, fb)
    m := 0.5 * (a + b)
    fm := f(m)
    quadval := ((b - a) / 6) * (fa + (4 * fm) + fb)
    return
  end
end

final abstract class ASR ()

  public static quadrature (f, a, b, tol, depth)
    local fa, fb, whole

    fa := f(a)
    fb := f(b)
    whole := Simpson_rule (f, a, fa, b, fb)
    return recursive_simpson (f, a, fa, b, fb, tol, whole.quadval,
                              whole.m, whole.fm, depth)
  end

  private static recursive_simpson (f, a, fa, b, fb, tol, whole,
                                    m, fm, depth)
    local left, right, delta, tol2
    local leftval, rightval, quadval

    left := Simpson_rule (f, a, fa, m, fm)
    right := Simpson_rule (f, m, fm, b, fb)
    delta := left.quadval + right.quadval - whole
    tol2 := 0.5 * tol
    if depth <= 0 | tol2 = tol | abs (delta) <= 15 * tol then
      quadval := left.quadval + right.quadval + (delta / 15)
    else
    {
      leftval := recursive_simpson (f, a, fa, m, fm, tol2,
                                    left.quadval, left.m, left.fm,
                                    depth - 1)
      rightval := recursive_simpson (f, m, fm, b, fb, tol2,
                                     right.quadval, right.m, right.fm,
                                     depth - 1)
      quadval := leftval + rightval
    }
    return quadval
  end

end
Output:
$ oiscript ./adaptive_simpson_task_oi.icn
estimate of ∫ sin x dx from 0 to 1: 0.4596976941

OCaml

let quad_asr f a b tol depth =
  let _quad_asr_simpsons a fa b fb =
    let m = 0.5 *. (a +. b) in
    let fm = f m in
    let h = b -. a in
    (m, fm, (h /. 6.0) *. (fa +. (4.0 *. fm) +. fb))
  in
  let rec _quad_asr a fa b fb tol whole m fm depth =
    let (lm, flm, left) = _quad_asr_simpsons a fa m fm in
    let (rm, frm, right) = _quad_asr_simpsons m fm b fb in
    let delta = left +. right -. whole in
    let tol_ = 0.5 *. tol in
    if depth <= 0 || tol_ = tol || Float.abs delta <= 15.0 *. tol then
      left +. right +. (delta /. 15.0)
    else
      _quad_asr a fa m fm tol_ left lm flm (depth - 1)
      +. _quad_asr m fm b fb tol_ right rm frm (depth - 1)
  in
  let fa = f a and fb = f b in
  let (m, fm, whole) = _quad_asr_simpsons a fa b fb in
  _quad_asr a fa b fb tol whole m fm depth
;;

print_string
  "estimated definite integral of sin(x) for x from 0 to 1: ";
print_float (quad_asr sin 0.0 1.0 1.0e-9 1000);
print_newline ()
;;
Output:
$ ocaml adaptive_simpson_task.ml
estimated definite integral of sin(x) for x from 0 to 1: 0.459697694132

Ol

This program also works with R7RS Scheme implementations that cannot run the #Scheme implementation.

Translation of: Scheme
Works with: Ol version 2.4
Works with: Chibi Scheme version 0.10.0
(import (scheme base)
        (scheme inexact)
        (scheme write))

(define (quad-asr f a b tol depth)
  (letrec ;; Recursive let, so %%quad-asr can call itself.
      ((%%quad-asr-simpsons
        (lambda (a fa b fb)
          (let* ((m (/ (+ a b) 2))
                 (fm (f m))
                 (h (- b a)))
            ;; Scheme supports returning multiple values at
            ;; once. There is no need to return them explicitly as
            ;; a data structure (though that also could be done).
            (values m fm (* (/ h 6) (+ fa (* 4 fm) fb))))))
       (%%quad-asr
        (lambda (a fa b fb tol whole m fm depth)
          ;; The R7RS standard specifies there must be
          ;; "let-values" and "let*-values" for receiving multiple
          ;; values. However, I do not want to assume its
          ;; presence. I will use the most widely supported
          ;; method, which is "call-with-values". (Typically
          ;; "let*-values" would be implemented as syntactic sugar
          ;; for the following.)
          (call-with-values
              (lambda () (%%quad-asr-simpsons a fa m fm))
            (lambda (lm flm left)
              (call-with-values
                  (lambda () (%%quad-asr-simpsons m fm b fb))
                (lambda (rm frm right)
                  (let ((delta (- (+ left right) whole))
                        (tol^ (/ tol 2)))
                    (if (or (<= depth 0)
                            (= tol^ tol)
                            (<= (abs delta) (* 15 tol)))
                        (+ left right (/ delta 15))
                        (+ (%%quad-asr a fa m fm tol^ left
                                       lm flm (- depth 1))
                           (%%quad-asr m fm b fb tol^ right
                                       rm frm (- depth 1))))))))))))
    (let ((fa (f a))
          (fb (f b)))
      (call-with-values (lambda () (%%quad-asr-simpsons a fa b fb))
        (lambda (m fm whole)
          (%%quad-asr a fa b fb tol whole m fm depth))))))

(display "estimated definite integral of sin(x) for x from 0 to 1: ")
(display (quad-asr sin 0 1 1e-9 1000))
(newline)
Output:
$ ol adaptive_simpson_task_ol.scm
estimated definite integral of sin(x) for x from 0 to 1: 0.459697694
$ chibi adaptive_simpson_task_ol.scm
estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941317858

Gauche Scheme can be told to start up with the R7RS default environment, and thus to need this version of the Scheme program:

$ gosh -r7 adaptive_simpson_task_ol.scm
estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941317858

Pascal

Works with: Free Pascal Compiler version 3.2.2
program adaptive_simpson_task;

type function_real_to_real = function(value : real) : real;
var fa, fb, m, fm, whole : real;

  function quad_asr (f         : function_real_to_real;
                     a, b, tol : real;
                     depth     : integer) : real;

    procedure quad_asr_simpsons_ (    a, fa, b, fb   : real;
                                  var m, fm, quadval : real);
    begin
      m := (a + b) / 2;
      fm := f(m);
      quadval := ((b - a) / 6) * (fa + (4 * fm) + fb)
    end;

    function quad_asr_ (a, fa, b, fb      : real;
                        tol, whole, m, fm : real;
                        depth             : integer) : real;
    var
      lm, flm, left  : real;
      rm, frm, right : real;
      delta, tol_    : real;
    begin
      quad_asr_simpsons_ (a, fa, m, fm, lm, flm, left);
      quad_asr_simpsons_ (m, fm, b, fb, rm, frm, right);
      delta := left + right - whole;
      tol_ := tol / 2;
      if (depth <= 0) or (tol_ = tol) or (abs(delta) <= 15 * tol) then
        quad_asr_ := left + right + (delta / 15)
      else
        quad_asr_ := (quad_asr_ (a, fa, m, fm, tol_,
                      left , lm, flm, depth - 1)
                      + quad_asr_ (m, fm, b, fb, tol_,
                      right, rm, frm, depth - 1))
    end;

  begin
    fa := f(a);
    fb := f(b);
    quad_asr_simpsons_ (a, fa, b, fb, m, fm, whole);
    quad_asr := quad_asr_ (a, fa, b, fb, tol, whole, m, fm, depth)
  end;

  function sine (x : real) : real;
  begin
    sine := sin (x);
  end;

begin
  writeln ('estimated definite integral of sin(x) ',
           'for x from 0 to 1: ', quad_asr (@sine, 0, 1, 1e-9, 1000))
end.
Output:
estimated definite integral of sin(x) for x from 0 to 1:  4.5969769413178579E-001

Perl

Translation of: Raku
use strict;
use warnings;

sub adaptive_Simpson_quadrature {
    my($f, $left, $right, $eps) = @_;
    my $lf = eval "$f($left)";
    my $rf = eval "$f($right)";
    my ($mid, $midf, $whole) = Simpson_quadrature_mid($f, $left, $lf, $right, $rf);
    return recursive_Simpsons_asr($f, $left, $lf, $right, $rf, $eps, $whole, $mid, $midf);

    sub Simpson_quadrature_mid {
        my($g, $l, $lf, $r, $rf) = @_;
        my $mid = ($l + $r) / 2;
        my $midf = eval "$g($mid)";
        ($mid, $midf, abs($r - $l) / 6 * ($lf + 4 * $midf + $rf))
    }

    sub recursive_Simpsons_asr {
        my($h, $a, $fa, $b, $fb, $eps, $whole, $m, $fm) = @_;
        my ($lm, $flm, $left)  = Simpson_quadrature_mid($h, $a, $fa, $m, $fm);
        my ($rm, $frm, $right) = Simpson_quadrature_mid($h, $m, $fm, $b, $fb);
        my $delta = $left + $right - $whole;
        abs($delta) <= 15 * $eps
            ? $left + $right + $delta / 15
            : recursive_Simpsons_asr($h, $a, $fa, $m, $fm, $eps/2, $left,  $lm, $flm) +
              recursive_Simpsons_asr($h, $m, $fm, $b, $fb, $eps/2, $right, $rm, $frm)
    }
}

my ($a, $b) = (0, 1);
my $sin = adaptive_Simpson_quadrature('sin', $a, $b, 1e-9);
printf "Simpson's integration of sine from $a to $b = %.9f", $sin
Output:
Simpson's integration of sine from 0 to 1 = 0.459697694

Phix

Translation of: Go
with javascript_semantics 
function quadSimpsonsMem(integer f, atom a, fa, b, fb)
    -- Evaluates Simpson's Rule, also returning m and f(m) to reuse.
    atom m = (a + b) / 2,
         fm = f(m),
         simp = abs(b-a) / 6 * (fa + 4*fm + fb)
    return {m, fm, simp}
end function
 
function quadAsrRec(integer f, atom a, fa, b, fb, eps, whole, m, fm)
    -- Efficient recursive implementation of adaptive Simpson's rule.
    -- Function values at the start, middle, end of the intervals are retained.
    atom {lm, flm, left} := quadSimpsonsMem(f, a, fa, m, fm),
         {rm, frm, rght} := quadSimpsonsMem(f, m, fm, b, fb),
         delta := left + rght - whole
    if abs(delta) <= eps*15 then
        return left + rght + delta/15
    end if
    return quadAsrRec(f, a, fa, m, fm, eps/2, left, lm, flm) +
           quadAsrRec(f, m, fm, b, fb, eps/2, rght, rm, frm)
end function
 
function quadAsr(integer f, atom a, b, eps)
    -- Integrate f from a to b using ASR with max error of eps.
    atom fa := f(a),
         fb := f(b),
         {m, fm, whole} := quadSimpsonsMem(f, a, fa, b, fb)
    return quadAsrRec(f, a, fa, b, fb, eps, whole, m, fm)
end function
 
-- we need a mini wrapper to get a routine_id for sin()
-- (because sin() is implemented in low-level assembly)
function _sin(atom a)
    return sin(a)
end function 
 
atom a := 0.0, b := 1.0,
     sinx := quadAsr(_sin, a, b, 1e-09)
printf(1,"Simpson's integration of sine from %g to %g = %f\n", {a, b, sinx})
Output:
Simpson's integration of sine from 0 to 1 = 0.459698

PostScript

Translation of: COBOL

The following program makes little effort to be efficient. It iteratively computes the intervals, rather than do recursions.

%!PS-Adobe-3.0 EPSF-3.0
%%BoundingBox: 0 0 280 30

/origstate save def

% We have to convert radians to degrees.
/f { 57.29577951308232 mul sin } def

/a 0.0 def
/b 1.0 def
/tol 0.000000001 def

/bisect-current-interval {
    numer 2 mul /numer exch def
    denom 2 mul /denom exch def
} def

/go-to-next-interval {
    numer 1 add /numer exch def
    {
        numer 2 mod 1 eq { exit } if
        numer 2 idiv /numer exch def
        denom 2 idiv /denom exch def
    } loop
} def

/there-are-no-more-intervals {
    numer denom eq
    denom 1 eq
    and
} def

/x0 {
    numer denom div
    dup 1 exch sub
    a mul
    exch b mul
    add
} def

/x1 {
    x0 0.75 mul
    x4 0.25 mul
    add
} def

/x2 {
    x0 0.5 mul
    x4 0.5 mul
    add
} def

/x3 {
    x0 0.25 mul
    x4 0.75 mul
    add
} def

/x4 {
    numer 1 add denom div
    dup 1 exch sub
    a mul
    exch b mul
    add
} def

/simpson-rule-thrice {
    x0 f /y0 exch def
    x1 f /y1 exch def
    x2 f /y2 exch def
    x3 f /y3 exch def
    x4 f /y4 exch def

    x4 x0 sub 6.0 div /h04 exch def
    x2 x0 sub 6.0 div /h02 exch def
    x4 x2 sub 6.0 div /h24 exch def

    y2 4.0 mul y0 add y4 add h04 mul /whole exch def
    y1 4.0 mul y0 add y2 add h02 mul /left exch def
    y3 4.0 mul y2 add y4 add h24 mul /right exch def
} def

/adaptive-simpson {
    /numer 0 def
    /denom 1 def
    /sum 0.0 def
    {
        there-are-no-more-intervals { exit } if
        simpson-rule-thrice
        left right add whole sub /delta exch def
        x4 x0 sub tol mul /tol0 exch def
        x2 x0 sub tol mul /tol1 exch def
        tol1 tol0 eq delta abs 15.0 tol0 mul le or
        {
            left right add delta 15.0 div add
            sum add /sum exch def
            go-to-next-interval
        }
        {
            bisect-current-interval
        } ifelse
    } loop
} def

/Times-Roman findfont 12 scalefont setfont
10 10 moveto (estimate of integral of sin x dx from 0 to 1: ) show
adaptive-simpson
sum 20 string cvs show

showpage
origstate restore
%%EOF
Output:

I ran the program thus:

magick convert -flatten adaptive_simpson_task.eps adaptive_simpson_task.png

And obtained:

Output from the PostScript for Adaptive Simpson. It shows 0.459698

Prolog

Works with: SWI-Prolog version 9.1.2
Works with: GNU Prolog version 1.5.0
%%% -*- mode: prolog; prolog-indent-width: 2; -*-

main :-
  quad_asr(sine, 0.0, 1.0, 0.000000001, 1000, QuadVal),
  write('estimate of ∫ sin x dx from 0 to 1: '),
  write(QuadVal),
  write('\n'),
  halt.

sine(X, Y) :- Y is sin(X).

quad_asr(F, A, B, Tol, Depth, QuadVal) :-
  call(F, A, FA),
  call(F, B, FB),
  simpson_rule(F, A, FA, B, FB, M, FM, Whole),
  recursive_simpson(F, A, FA, B, FB, Tol, Whole, M, FM, Depth,
                    QuadVal).

recursive_simpson(F, A, FA, B, FB, Tol, Whole, M, FM, Depth,
                  QuadVal) :-
  simpson_rule(F, A, FA, M, FM, LM, FLM, Left),
  simpson_rule(F, M, FM, B, FB, RM, FRM, Right),
  Delta is (Left + Right - Whole),
  Tol_ is (0.5 * Tol),
  ((Depth > 0,
    Tol_ =\= Tol,
    AbsDelta is abs(Delta),
    Tol15 is (15.0 * Tol),
    AbsDelta > Tol15)
  -> (Depth_ is Depth - 1,
      recursive_simpson(F, A, FA, M, FM, Tol_, Left, LM, FLM,
                        Depth_, QuadValLeft),
      recursive_simpson(F, M, FM, B, FB, Tol_, Right, RM, FRM,
                        Depth_, QuadValRight),
      QuadVal is QuadValLeft + QuadValRight)
  ;  left_right_estimate(Left, Right, Delta, QuadVal)).

left_right_estimate(Left, Right, Delta, Estimate) :-
  Estimate is Left + Right + (Delta / 15.0).

simpson_rule(F, A, FA, B, FB, M, FM, QuadVal) :-
  M is (0.5 * (A + B)),
  call(F, M, FM),
  QuadVal is ((B - A) / 6.0) * (FA + (4.0 * FM) + FB).

:- initialization(main).
Output:

The number is printed differently by different Prolog implementations:

$ swipl adaptive_simpson_task_prolog.pl
estimate of ∫ sin x dx from 0 to 1: 0.4596976941317858
$ gplc adaptive_simpson_task_prolog.pl && ./adaptive_simpson_task_prolog
estimate of ∫ sin x dx from 0 to 1: 0.45969769413178579

Python

#! python3

'''
    example

    $ python3 /tmp/integrate.py
    Simpson's integration of sine from 0.0 to 1.0 = 0.4596976941317858

    expected answer computed by j www.jsoftware.com  

       1-&:(1&o.d._1)0
    0.459698


    translated from c
'''

import math

import collections
triple = collections.namedtuple('triple', 'm fm simp')

def _quad_simpsons_mem(f: callable, a: float , fa: float, b: float, fb: float)->tuple:
    '''
        Evaluates Simpson's Rule, also returning m and f(m) to reuse.
    '''
    m = a + (b - a) / 2
    fm = f(m)
    simp = abs(b - a) / 6 * (fa + 4*fm + fb)
    return triple(m, fm, simp,)
 
def _quad_asr(f: callable, a: float, fa: float, b: float, fb: float, eps: float, whole: float, m: float, fm: float)->float:
    '''
    	Efficient recursive implementation of adaptive Simpson's rule.
    	Function values at the start, middle, end of the intervals are retained.
    '''
    lt = _quad_simpsons_mem(f, a, fa, m, fm)
    rt = _quad_simpsons_mem(f, m, fm, b, fb)
    delta = lt.simp + rt.simp - whole
    return (lt.simp + rt.simp + delta/15
        if (abs(delta) <= eps * 15) else
            _quad_asr(f, a, fa, m, fm, eps/2, lt.simp, lt.m, lt.fm) +
            _quad_asr(f, m, fm, b, fb, eps/2, rt.simp, rt.m, rt.fm)
    )

def quad_asr(f: callable, a: float, b: float, eps: float)->float:
    '''
        Integrate f from a to b using ASR with max error of eps.
    '''
    fa = f(a)
    fb = f(b)
    t = _quad_simpsons_mem(f, a, fa, b, fb)
    return _quad_asr(f, a, fa, b, fb, eps, t.simp, t.m, t.fm)

def main():
    (a, b,) = (0.0, 1.0,)
    sinx = quad_asr(math.sin, a, b, 1e-09);
    print("Simpson's integration of sine from {} to {} = {}\n".format(a, b, sinx))

main()

Raku

(formerly Perl 6)

Works with: Rakudo version 2018.10

Fairly direct translation of the Python code.

sub adaptive-Simpson-quadrature(&f, $left, $right, \ε = 1e-9) {
    my $lf = f($left);
    my $rf = f($right);
    my ($mid, $midf, $whole) = Simpson-quadrature-mid(&f, $left, $lf, $right, $rf);
    return recursive-Simpsons-asr(&f, $left, $lf, $right, $rf, ε, $whole, $mid, $midf);

    sub Simpson-quadrature-mid(&g, $l, $lf, $r, $rf){
        my $mid = ($l + $r) / 2;
        my $midf = g($mid);
        ($mid, $midf, ($r - $l).abs / 6 * ($lf + 4 * $midf + $rf))
    }

    sub recursive-Simpsons-asr(&h, $a, $fa, $b, $fb, $eps, $whole, $m, $fm){
        my ($lm, $flm, $left)  = Simpson-quadrature-mid(&h, $a, $fa, $m, $fm);
        my ($rm, $frm, $right) = Simpson-quadrature-mid(&h, $m, $fm, $b, $fb);
        my $delta = $left + $right - $whole;
        $delta.abs <= 15 * $eps
            ?? $left + $right + $delta / 15
            !! recursive-Simpsons-asr(&h, $a, $fa, $m, $fm, $eps/2, $left,  $lm, $flm) +
               recursive-Simpsons-asr(&h, $m, $fm, $b, $fb, $eps/2, $right, $rm, $frm)
    }
}

my ($a, $b) = 0e0, 1e0;
my $sin = adaptive-Simpson-quadrature(&sin, $a, $b, 1e-9).round(10**-9);;
say "Simpson's integration of sine from $a to $b = $sin";
Output:
Simpson's integration of sine from 0 to 1 = 0.459697694

RATFOR

I make available a copy of the public domain ratfor77, with some bug fixes, at https://sourceforge.net/p/chemoelectric/ratfor77/ It is written in C89.

#
# This is RATFOR 77, which is a preprocessor for FORTRAN 77. In
# neither language can you do anything that might look like it would
# be a recursive call, and have it actually be a recursive call. They
# are languages that assume all variables are equivalent to "extern"
# or "static" variables in C.
#
# There are other options, however.
#
# One method would be to keep track of interval size as one goes
# along, and perform the whole computation "flat": without anything
# resembling a stack. This could result, in this case, in the
# recalculation of intermediate results.  Usually, however, it is my
# favorite way to do adaptive bisection, because it is safest
# (especially if you have big integers). There is no possibility of
# running out of stack.
#
# Another method is to build up a queue or stack of subtasks to
# perform, as if these subtasks were on the system stack. This also
# can be safer than recursion, because it lets you put the subtasks in
# the heap instead of the stack. (A stack of subtasks to perform is,
# in fact, how most compilers for other languages would implement
# recursive calls.)
#
# I will use the latter method, but without a heap. (An external heap
# would require either a language extension or a foreign function
# call.)
#

define(STKSZ,1000)

subroutine simpsn (f, a, fa, b, fb, m, fm, quad)
  implicit none

  external f
  double precision f
  double precision a, fa, b, fb, m, fm, quad

  m = (a + b) / 2
  fm = f(m)
  quad = ((b - a) / 6) * (fa + (4 * fm) + fb)
end

subroutine adapt (sum, f, a, fa, b, fb, tol, quad, m, fm)
  implicit none

  external f, simpsn
  double precision f
  double precision a(STKSZ), fa(STKSZ), b(STKSZ), fb(STKSZ)
  double precision tol(STKSZ), quad(STKSZ), m(STKSZ), fm(STKSZ)
  double precision sum

  integer i
  double precision lm, flm, left
  double precision rm, frm, right
  double precision delta, tol2

  i = 1
  while (i != 0)
    {
      if (i == STKSZ)
        {                      # It is not possible to bisect further.
          sum = sum + quad(i)
          i = i - 1
        }
      else
        {
          call simpsn (f, a(i), fa(i), m(i), fm(i), lm, flm, left)
          call simpsn (f, m(i), fm(i), b(i), fb(i), rm, frm, right)
          delta = left + right - quad(i)
          tol2 = tol(i) / 2
          if (tol2 == tol(i) .or. abs (delta) <= 15 * tol(i))
            {                   # The tolerance is satisfied.
              sum = sum + left + right + (delta / 15)
              i = i - 1
            }
          else
            {        # The tolerance is not satisfied. Bisect further.
              b(i + 1) = b(i);      fb(i + 1) = fb(i)
              a(i + 1) = m(i);      fa(i + 1) = fm(i)
              b(i) = m(i);          fb(i) = fm(i)
              tol(i + 1) = tol2;    tol(i) = tol2
              quad(i + 1) = right;  quad(i) = left
              m(i + 1) = rm;        m(i) = lm
              fm(i + 1) = frm;      fm(i) = flm
              i = i + 1
            }
        }
    }
end

function asimps (f, a0, b0, tol0)
  implicit none
  double precision asimps

  external f, simpsn, adapt
  double precision f, a0, b0, tol0

  double precision sum
  double precision a(STKSZ), fa(STKSZ)
  double precision b(STKSZ), fb(STKSZ)
  double precision m(STKSZ), fm(STKSZ)
  double precision quad(STKSZ), tol(STKSZ)

  sum = 0
  tol(1) = tol0
  a(1) = a0;  fa(1) = f(a0)
  b(1) = b0;  fb(1) = f(b0)
  call simpsn (f, a(1), fa(1), b(1), fb(1), m(1), fm(1), quad(1))
  call adapt (sum, f, a, fa, b, fb, tol, quad, m, fm)
  asimps = sum
end

function sine (x)
  implicit none
  double precision sine, x
  sine = dsin (x)
end

program asrtsk
  implicit none

  external asimps, sine
  double precision asimps, sine
  double precision quad

  quad = asimps (sine, 0.0D0, 1.0D0, 1D-9)
  write (*, 1000) quad
1000 format ("estimated definite integral of sin(x) ", _
             "for x from 0 to 1: ", F15.13)
end
Output:

The output looks different, depending on which of my Fortran compilers I use:

$ ratfor77 adaptive_simpson_task_ratfor.r > asr.f && gfortran -fbounds-check -g asr.f && ./a.out
estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941318
$ ratfor77 adaptive_simpson_task_ratfor.r > asr.f && f2c asr.f 2>/dev/null > asr.c && gcc -g asr.c -lf2c -lm && ./a.out
estimated definite integral of sin(x) for x from 0 to 1:  .4596976941318

REXX

Translation of: Go
/*REXX program performs numerical integration using adaptive Simpson's method.          */
numeric digits length( pi() )  -  length(.)      /*use # of digits in pi for precision. */
a= 0;     b= 1;       f= 'SIN'                   /*define values for  A,  B,  and  F.   */
sinx= quadAsr('SIN',a,b,"1e" || (-digits() + 1) )
say "Simpson's integration of sine from "      a     ' to '      b      ' = '       sinx
exit                                             /*stick a fork in it,  we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
pi:       pi= 3.14159265358979323846;  return pi /*pi  has  twenty-one  decimal digits. */
r2r:      return arg(1)  //  (pi() *2)           /*normalize radians ──► a unit circle, */
/*──────────────────────────────────────────────────────────────────────────────────────*/
quadSimp: procedure; parse arg f,a,fa,b,fb;   m= (a+b) / 2;     interpret  'fm='   f"(m)"
          simp= abs(b-a) / 6 * (fa + 4*fm + fb);    return  m  fm  simp
/*──────────────────────────────────────────────────────────────────────────────────────*/
quadAsr:  procedure; parse arg f,a,b,eps;                       interpret  'fa='   f"(a)"
                                                                interpret  'fb='   f"(b)"
          parse value  quadSimp(f,a,fa,b,fb)   with   m  fm  whole
          return quadAsrR(f,a,fa,b,fb,eps,whole,m,fm)
/*──────────────────────────────────────────────────────────────────────────────────────*/
quadAsrR: procedure; parse arg f,a,fa,b,fb,eps,whole,m,fm;      frac= digits() * 3 / 4
          parse value   quadSimp(f,a,fa,m,fm)   with  lm  flm  left
          parse value   quadSimp(f,m,fm,b,fb)   with  rm  frm  right
          $= left + right - whole                                     /*calculate delta.*/
          if abs($)<=eps*frac  then return left + right + $/frac
          return quadAsrR(f,a,fa,m,fm,eps/2, left,lm,flm) + ,
                 quadAsrR(f,m,fm,b,fb,eps/2,right,rm,frm)
/*──────────────────────────────────────────────────────────────────────────────────────*/
sin:      procedure; parse arg x;   x= r2r(x);   numeric fuzz min(5, max(1, digits() -3) )
          if x=pi*.5           then return 1;       if x==pi * 1.5  then return -1
          if abs(x)=pi | x=0   then return 0
                                           #= x;    _= x;    q= x*x
            do k=2  by 2  until p=#;    p= #;       _= - _ * q / (k * (k+1) );    #= # + _
            end   /*k*/
          return #
output   when using the default inputs:
Simpson's integration of sine from  0  to  1  =  0.459697694131860282602

Scheme

For an R7RS-specific version, see #Ol.

(define (quad-asr f a b tol depth)
  (letrec ;; Recursive let, so %%quad-asr can call itself.
      ((%%quad-asr-simpsons
        (lambda (a fa b fb)
          (let* ((m (/ (+ a b) 2))
                 (fm (f m))
                 (h (- b a)))
            ;; Scheme supports returning multiple values at
            ;; once. There is no need to return them explicitly as
            ;; a data structure (though that also could be done).
            (values m fm (* (/ h 6) (+ fa (* 4 fm) fb))))))
       (%%quad-asr
        (lambda (a fa b fb tol whole m fm depth)
          ;; The R7RS standard specifies there must be
          ;; "let-values" and "let*-values" for receiving multiple
          ;; values. However, I do not want to assume its
          ;; presence. I will use the most widely supported
          ;; method, which is "call-with-values". (Typically
          ;; "let*-values" would be implemented as syntactic sugar
          ;; for the following.)
          (call-with-values
              (lambda () (%%quad-asr-simpsons a fa m fm))
            (lambda (lm flm left)
              (call-with-values
                  (lambda () (%%quad-asr-simpsons m fm b fb))
                (lambda (rm frm right)
                  (let ((delta (- (+ left right) whole))
                        (tol^ (/ tol 2)))
                    (if (or (<= depth 0)
                            (= tol^ tol)
                            (<= (abs delta) (* 15 tol)))
                        (+ left right (/ delta 15))
                        (+ (%%quad-asr a fa m fm tol^ left
                                       lm flm (- depth 1))
                           (%%quad-asr m fm b fb tol^ right
                                       rm frm (- depth 1))))))))))))
    (let ((fa (f a))
          (fb (f b)))
      (call-with-values (lambda () (%%quad-asr-simpsons a fa b fb))
        (lambda (m fm whole)
          (%%quad-asr a fa b fb tol whole m fm depth))))))

(display "estimated definite integral of sin(x) for x from 0 to 1: ")
(display (quad-asr sin 0 1 1e-9 1000))
(newline)
Output:

Running it with Racket's R5RS implementation:

$ plt-r5rs adaptive_simpson_task.scm
estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941317858

Running it with SCM:

$ scm -f adaptive_simpson_task.scm
estimated definite integral of sin(x) for x from 0 to 1: 459.6976941317858e-3

Running it with Gambit Scheme:

$ gsi adaptive_simpson_task.scm
estimated definite integral of sin(x) for x from 0 to 1: .4596976941317858

Shen

(define simpson-rule
  F A FA B FB ->
    (let M (* 0.5 (+ A B))
         FM (F M)
         QuadVal (* (/ (- B A) 6.0) (+ FA (* 4.0 FM) FB))
         [M FM QuadVal]))

(define recursive-simpson
  F A FA B FB Tol Whole M FM Depth ->
    (let TheLeftStuff (simpson-rule F A FA M FM)
         TheRightStuff (simpson-rule F M FM B FB)
         (recursive-simpson-united F A FA B FB Tol Whole M FM Depth
                                   TheLeftStuff TheRightStuff)))

(define recursive-simpson-united
  F A FA B FB Tol Whole M FM Depth
  [LM FLM Left] [RM FRM Right] ->
    (let Delta (- (+ Left Right) Whole)
         Tol_ (* 0.5 Tol)
         (if (or (<= Depth 0)
                 (= Tol_ Tol)
                 (<= (abs Delta) (* 15.0 Tol)))
            (+ Left Right (/ Delta 15.0))
            (+ (recursive-simpson F A FA M FM Tol_
                                  Left LM FLM (- Depth 1))
               (recursive-simpson F M FM B FB Tol_
                                  Right RM FRM (- Depth 1))))))

(define quad-asr
  F A B Tol Depth ->
    (let FA (F A)
         FB (F B)
         TheWholeStuff (simpson-rule F A FA B FB)
         (quad-asr-part-two F A FA B FB Tol Depth TheWholeStuff)))

(define quad-asr-part-two
  F A FA B FB Tol Depth [M FM Whole] ->
    (recursive-simpson F A FA B FB Tol Whole M FM Depth))

\\ The Shen standard library contains only an unserious
\\ implementation of sin(x), so I might as well toss together
\\ my own! Thus: an uneconomized Maclaurin expansion. Gotten
\\ via "taylor(sin(x),x,0,20);" in Maxima. Using up to x**13
\\ should give an error bound on the order of 1e-12, for
\\ 0 <= x <= 1. That is much better than we need.
(define sin
  X -> (compute-sin X (* X X)))
(define compute-sin
  X XX ->
    (* X (+ 1.0
            (* XX (+ (/ -1.0 6.0)
                     (* XX (+ (/ 1.0 120.0)
                              (* XX (+ (/ -1.0 5040.0)
                                       (* XX (+ (/ 1.0 362880.0)
                                                (* XX (+ (/ -1.0 39916800.0)
                                                         (* XX (/ 1.0 6227020800.0)))))))))))))))

(define abs
  \\ Like sin, abs is in the standard library, but it is very
  \\ easy to write.
  X -> (if (< X 0) (- 0 X) X))

(output "estimate of the definite integral of sin x dx from 0 to 1: ")
(print (quad-asr (lambda x (sin x)) 0.0 1.0 1.0e-9 100))
(nl)

\\ local variables:
\\ mode: indented-text
\\ tab-width: 2
\\ end:
Output:
$ shen-scheme script adaptive_simpson_task.shen
estimate of the definite integral of sin x dx from 0 to 1: 0.4596976941318334

Sidef

Translation of: Raku
func adaptive_Simpson_quadrature(f, left, right, ε = 1e-9) {

    func quadrature_mid(l, lf, r, rf) {
        var mid = (l+r)/2
        var midf = f(mid)
        (mid, midf, abs(r-l)/6 * (lf + 4*midf + rf))
    }

    func recursive_asr(a, fa, b, fb, ε, whole, m, fm) {
        var (lm, flm, left)  = quadrature_mid(a, fa, m, fm)
        var (rm, frm, right) = quadrature_mid(m, fm, b, fb)
        var Δ = (left + right - whole)
        abs(Δ) <= 15*ε
            ? (left + right + Δ/15)
            : (__FUNC__(a, fa, m, fm, ε/2, left,  lm, flm) +
               __FUNC__(m, fm, b, fb, ε/2, right, rm, frm))
    }

    var (lf = f(left), rf = f(right))
    var (mid, midf, whole) = quadrature_mid(left, lf, right, rf)
    recursive_asr(left, lf, right, rf, ε, whole, mid, midf)
}

var (a = 0, b = 1)
var area = adaptive_Simpson_quadrature({ .sin }, a, b, 1e-15).round(-15)
say "Simpson's integration of sine from #{a} to #{b}#{area}"
Output:
Simpson's integration of sine from 0 to 1 ≈ 0.45969769413186

Standard ML

fun quad_asr (tol, depth) (f, a, b) =
    let
      fun quad_asr_simpsons_ (a, fa, b, fb) =
          let
            val m = 0.5 * (a + b)
            val fm = f (m)
            val h = b - a
          in
            (m, fm, (h / 6.0) * (fa + (4.0 * fm) + fb))
          end

      fun quad_asr_ (a, fa, b, fb, tol, whole, m, fm, depth) =
          let
            val (lm, flm, left) = quad_asr_simpsons_ (a, fa, m, fm)
            and (rm, frm, right) = quad_asr_simpsons_ (m, fm, b, fb)
            val delta = left + right - whole
            and tol_ = 0.5 * tol
          in
            if depth <= 0
               orelse abs (tol_ - tol) <= 0.0 (* <-- SML silliness *)
               orelse Real.abs (delta) <= 15.0 * tol then
              left + right + (delta / 15.0)
            else
              quad_asr_ (a, fa, m, fm, tol_,
                         left, lm, flm, depth - 1)
              + quad_asr_ (m, fm, b, fb, tol_,
                           right, rm, frm, depth - 1)
          end 

      val fa = f (a) and fb = f (b)
      val (m, fm, whole) = quad_asr_simpsons_ (a, fa, b, fb)
    in
      quad_asr_ (a, fa, b, fb, tol, whole, m, fm, depth)
    end
;

val quadrature = quad_asr (0.000000001, 1000);

print "estimated definite integral of sin(x) for x from 0 to 1: ";
print (Real.toString (quadrature (Math.sin, 0.0, 1.0)));
print "\n";
Output:
$ mlton adaptive_simpson_task.sml && ./adaptive_simpson_task
estimated definite integral of sin(x) for x from 0 to 1: 0.459697694132

Wren

Translation of: Go
Library: Wren-fmt
import "./fmt" for Fmt

/* "structured" adaptive version, translated from Racket */
var quadSimpsonMem = Fn.new { |f, a, fa, b, fb|
    // Evaluates Simpson's Rule, also returning m and f.call(m) to reuse.
    var m = (a + b) / 2
    var fm = f.call(m)
    var simp = (b - a).abs / 6 * (fa + 4*fm + fb)
    return [m, fm, simp]
}

var quadAsrRec // recursive
quadAsrRec = Fn.new { |f, a, fa, b, fb, eps, whole, m, fm|
    // Efficient recursive implementation of adaptive Simpson's rule.
    // Function values at the start, middle, end of the intervals are retained.
    var r1 = quadSimpsonMem.call(f, a, fa, m, fm)
    var r2 = quadSimpsonMem.call(f, m, fm, b, fb)
    var lm    = r1[0]
    var flm   = r1[1]
    var left  = r1[2]
    var rm    = r2[0]
    var frm   = r2[1]
    var right = r2[2]
    var delta = left + right - whole
    if (delta.abs < eps * 15) return left + right + delta/15
    return quadAsrRec.call(f, a, fa, m, fm, eps/2, left, lm, flm) +
           quadAsrRec.call(f, m, fm, b, fb, eps/2, right, rm, frm)
}

var quadAsr = Fn.new { |f, a, b, eps|
    // Integrate f from a to b using ASR with max error of eps.
    var fa = f.call(a)
    var fb = f.call(b)
    var r = quadSimpsonMem.call(f, a, fa, b, fb)
    var m     = r[0]
    var fm    = r[1]
    var whole = r[2]
    return quadAsrRec.call(f, a, fa, b, fb, eps, whole, m, fm)
}

var a = 0
var b = 1
var sinx = quadAsr.call(Fn.new { |x| x.sin }, a, b, 1e-09)
Fmt.print("Simpson's integration of sine from $d to $d = $f", a, b, sinx)
Output:
Simpson's integration of sine from 0 to 1 = 0.459698

zkl

Translation of: Python
# "structured" adaptive version, translated from Racket
fcn _quad_simpsons_mem(f, a,fa, b,fb){
   #Evaluates the Simpson's Rule, also returning m and f(m) to reuse"""
   m,fm := (a + b)/2, f(m);
   return(m,fm, (b - a).abs()/6*(fa + fm*4 + fb));
}
fcn _quad_asr(f, a,fa, b,fb, eps, whole, m,fm){
  # Efficient recursive implementation of adaptive Simpson's rule.
  # Function values at the start, middle, end of the intervals are retained.

   lm,flm,left  := _quad_simpsons_mem(f, a,fa, m,fm);
   rm,frm,right := _quad_simpsons_mem(f, m,fm, b,fb);
   delta:=left + right - whole;
   if(delta.abs() <= eps*15) return(left + right + delta/15);
   _quad_asr(f, a,fa, m,fm, eps/2, left , lm,flm) +
   _quad_asr(f, m,fm, b,fb, eps/2, right, rm,frm)
}
fcn quad_asr(f,a,b,eps){
   #Integrate f from a to b using Adaptive Simpson's Rule with max error of eps
   fa,fb      := f(a),f(b);
   m,fm,whole := _quad_simpsons_mem(f, a,fa, b,fb);
   _quad_asr(f, a,fa, b,fb, eps,whole,m,fm);
}
sinx:=quad_asr((1.0).sin.unbind(), 0.0, 1.0, 1e-09);
println("Simpson's integration of sine from 1 to 2 = ",sinx);
Output:
Simpson's integration of sine from 1 to 2 = 0.459698