View source for Numerical integration/Adaptive Simpson's method

{{draft task|Arithmetic operations}}
Lyness's (1969) Modified [[wp:Adaptive_Simpson's_method|Adaptive Simpson's method]] ([https://doi.org/10.1145/321526.321537 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.

{| class="mw-collapsible"
|+ 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)
|}

=={{header|11l}}==
{{trans|Python}}

<syntaxhighlight lang="11l">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))</syntaxhighlight>

{{out}}
<pre>
Simpson's integration of sine from 0 to 1 = 0.459697694
</pre>

=={{header|ALGOL 60}}==
{{works with|GNU MARST|2.7}}

<syntaxhighlight lang="algol60">
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
</syntaxhighlight>

{{out}}
<pre>$ 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</pre>

=={{header|ALGOL 68}}==
{{Trans|C}}...with the helper routines nested inside the main quadrature routine.
<syntaxhighlight lang="algol68">
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
</syntaxhighlight>
{{out}}
<pre>
Simpson's integration of sine from  0.000 to  1.000 =  0.459698
</pre>

=={{header|ALGOL W}}==
{{Trans|ALGOL 68}} which is
{{Trans|C}}
<syntaxhighlight lang="algolw">
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.
</syntaxhighlight>
{{out}}
<pre>
Simpson's integration of sine from  0.000 to  1.000 =   0.45969769
</pre>

=={{header|ATS}}==
<syntaxhighlight lang="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
</syntaxhighlight>

{{out}}
<pre>$ 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</pre>

=={{header|C}}==
{{trans|zkl}}
<syntaxhighlight lang="c">#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;
}</syntaxhighlight>

{{out}}
<pre>
Simpson's integration of sine from 0 to 1 = 0.459698
</pre>

=={{header|COBOL}}==
{{works with|GnuCOBOL|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).

<syntaxhighlight lang="cobol">
       *> -*- 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 9(25).  *> 25 decimal digits.
       01  denom picture 9(25).
       01  numer2 picture 9(25).
       01  a picture 9(5)V9(20). *> 5+20 decimal digits.
       01  b picture 9(5)V9(20).
       01  tol picture 9(5)V9(20).
       01  x picture 9(5)V9(20).
       01  y picture 9(5)V9(20).
       01  x0 picture 9(5)V9(20).
       01  y0 picture 9(5)V9(20).
       01  x1 picture 9(5)V9(20).
       01  y1 picture 9(5)V9(20).
       01  x2 picture 9(5)V9(20).
       01  y2 picture 9(5)V9(20).
       01  x3 picture 9(5)V9(20).
       01  y3 picture 9(5)V9(20).
       01  x4 picture 9(5)V9(20).
       01  y4 picture 9(5)V9(20).
       01  ruleval picture 9(5)V9(20).
       01  ruleval0 picture 9(5)V9(20).
       01  ruleval1 picture 9(5)V9(20).
       01  delta picture 9(5)V9(20).
       01  tol0 picture 9(5)V9(20).
       01  tol1 picture 9(5)V9(20).
       01  delta15 picture 9(5)V9(20).
       01  stepval picture 9(5)V9(20).
       01  quadval picture 9(5)V9(20).
       01  result picture 9V9(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
               if delta < 0 then
                   compute delta = -delta
               end-if
               compute tol0 = tol * (x4 - x0)
               compute tol1 = tol * (x2 - x0)
               if (tol1 = tol0) or (delta <= 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).
</syntaxhighlight>

{{out}}
<pre>$ cobc -std=cobol2014 -x adaptive_simpson_task.cob && ./adaptive_simpson_task
0.459697694</pre>

=={{header|Common Lisp}}==
<syntaxhighlight lang="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)
</syntaxhighlight>

{{out}}
<pre>$ clisp adaptive_simpson_task.lisp
estimated definite integral of sin(x) for x from 0 to 1: 0.45969772</pre>

=={{header|D}}==
<syntaxhighlight lang="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;
}
</syntaxhighlight>

{{out}}
<pre>$ gdc -g adaptive_simpson_task.d && ./a.out
estimate of ∫ sin x dx from 0 to 1: 0.459698</pre>

=={{header|Factor}}==
{{trans|Julia}}
<syntaxhighlight lang="factor">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</syntaxhighlight>
{{out}}
<pre>
Simpson's rule integration of sin from 0 to 1 is: 0.4596976941573994
</pre>

=={{header|Fortran}}==
<syntaxhighlight lang="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
</syntaxhighlight>

{{out}}
<pre>$ 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</pre>

On x86-64, if you leave out the optimizer you get a trampoline:
<pre>$ 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</pre>

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. :)

=={{header|Go}}==
Like the zkl entry, this is also a translation of the Python code in the Wikipedia article.
<syntaxhighlight lang="go">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)   
}</syntaxhighlight>

{{out}}
<pre>
Simpson's integration of sine from 0 to 1 = 0.459698
</pre>

=={{header|J}}==
Typically one would choose the library implementation:
<pre>
   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
   
</pre>


<syntaxhighlight lang="j">
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
)
</syntaxhighlight>

<pre>
   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
</pre>

=={{header|Java}}==
<syntaxhighlight lang="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);
    }
    
}
</syntaxhighlight>
{{out}}
<pre>
integrate sin(x), x = 0 .. Pi = 1.999999999998.  Function calls = 121
integrate sin(x), x = 0 .. 1 = 0.459697694131.  Function calls = 33
</pre>

=={{header|Julia}}==
Originally from Modesto Mas, https://mmas.github.io/simpson-integration-julia
<syntaxhighlight lang="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))
</syntaxhighlight>{{out}}
<pre>
Simpson's rule integration of sin from 0 to 1 is: 0.45969769415739936
</pre>

=={{header|Kotlin}}==
{{trans|Go}}
<syntaxhighlight lang="scala">// 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)}")
}</syntaxhighlight>

{{output}}
<pre>
Simpson's integration of sine from 0.0 to 1.0 = 0.459698
</pre>

=={{header|Lambdatalk}}==
<syntaxhighlight lang="scheme">
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
</syntaxhighlight>

=={{header|Nim}}==
Direct translation of Python code from Wikipedia entry.
<syntaxhighlight lang="nim">import math, sugar

type Func = float -> float

func 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

func 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)

func 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)</syntaxhighlight>

{{out}}
<pre>Simpson's integration of sine from 0 to 1 = 0.4596976941317859</pre>

=={{header|OCaml}}==
<syntaxhighlight lang="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 ()
;;
</syntaxhighlight>

{{out}}
<pre>$ ocaml adaptive_simpson_task.ml
estimated definite integral of sin(x) for x from 0 to 1: 0.459697694132</pre>

=={{header|Ol}}==
This program also works with R7RS Scheme implementations that cannot run the [[#Scheme]] implementation.
{{trans|Scheme}}
{{works with|Ol|2.4}}
{{works with|Chibi Scheme|0.10.0}}

<syntaxhighlight lang="scheme">
(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)
</syntaxhighlight>

{{out}}
<pre>$ 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</pre>

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

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

=={{header|Pascal}}==
{{works with|Free Pascal Compiler|3.2.2}}

<syntaxhighlight lang="pascal">
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.
</syntaxhighlight>

{{out}}
<pre>estimated definite integral of sin(x) for x from 0 to 1:  4.5969769413178579E-001</pre>

=={{header|Perl}}==
{{trans|Raku}}
<syntaxhighlight lang="perl">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</syntaxhighlight>
{{out}}
<pre>Simpson's integration of sine from 0 to 1 = 0.459697694</pre>

=={{header|Phix}}==
{{trans|Go}}
<!--<syntaxhighlight lang="phix">(phixonline)-->
 <span style="color: #008080;">with</span> <span style="color: #008080;">javascript_semantics</span> 
 <span style="color: #008080;">function</span> <span style="color: #000000;">quadSimpsonsMem</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">)</span>
     <span style="color: #000080;font-style:italic;">-- Evaluates Simpson's Rule, also returning m and f(m) to reuse.</span>
     <span style="color: #004080;">atom</span> <span style="color: #000000;">m</span> <span style="color: #0000FF;">=</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">a</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">2</span><span style="color: #0000FF;">,</span>
          <span style="color: #000000;">fm</span> <span style="color: #0000FF;">=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">m</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">simp</span> <span style="color: #0000FF;">=</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">-</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">/</span> <span style="color: #000000;">6</span> <span style="color: #0000FF;">*</span> <span style="color: #0000FF;">(</span><span style="color: #000000;">fa</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">4</span><span style="color: #0000FF;">*</span><span style="color: #000000;">fm</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">return</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">simp</span><span style="color: #0000FF;">}</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
  
 <span style="color: #008080;">function</span> <span style="color: #000000;">quadAsrRec</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">eps</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">whole</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fm</span><span style="color: #0000FF;">)</span>
     <span style="color: #000080;font-style:italic;">-- Efficient recursive implementation of adaptive Simpson's rule.
     -- Function values at the start, middle, end of the intervals are retained.</span>
     <span style="color: #004080;">atom</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">lm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">flm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">left</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">quadSimpsonsMem</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fm</span><span style="color: #0000FF;">),</span>
          <span style="color: #0000FF;">{</span><span style="color: #000000;">rm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">frm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rght</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">quadSimpsonsMem</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">delta</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">left</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">rght</span> <span style="color: #0000FF;">-</span> <span style="color: #000000;">whole</span>
     <span style="color: #008080;">if</span> <span style="color: #7060A8;">abs</span><span style="color: #0000FF;">(</span><span style="color: #000000;">delta</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;"><=</span> <span style="color: #000000;">eps</span><span style="color: #0000FF;">*</span><span style="color: #000000;">15</span> <span style="color: #008080;">then</span>
         <span style="color: #008080;">return</span> <span style="color: #000000;">left</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">rght</span> <span style="color: #0000FF;">+</span> <span style="color: #000000;">delta</span><span style="color: #0000FF;">/</span><span style="color: #000000;">15</span>
     <span style="color: #008080;">end</span> <span style="color: #008080;">if</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">quadAsrRec</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">eps</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">left</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">lm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">flm</span><span style="color: #0000FF;">)</span> <span style="color: #0000FF;">+</span>
            <span style="color: #000000;">quadAsrRec</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">eps</span><span style="color: #0000FF;">/</span><span style="color: #000000;">2</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rght</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">rm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">frm</span><span style="color: #0000FF;">)</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
  
 <span style="color: #008080;">function</span> <span style="color: #000000;">quadAsr</span><span style="color: #0000FF;">(</span><span style="color: #004080;">integer</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">eps</span><span style="color: #0000FF;">)</span>
     <span style="color: #000080;font-style:italic;">-- Integrate f from a to b using ASR with max error of eps.</span>
     <span style="color: #004080;">atom</span> <span style="color: #000000;">fa</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">),</span>
          <span style="color: #000000;">fb</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">f</span><span style="color: #0000FF;">(</span><span style="color: #000000;">b</span><span style="color: #0000FF;">),</span>
          <span style="color: #0000FF;">{</span><span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fm</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">whole</span><span style="color: #0000FF;">}</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">quadSimpsonsMem</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">return</span> <span style="color: #000000;">quadAsrRec</span><span style="color: #0000FF;">(</span><span style="color: #000000;">f</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fa</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fb</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">eps</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">whole</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">m</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">fm</span><span style="color: #0000FF;">)</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span>
  
 <span style="color: #000080;font-style:italic;">-- we need a mini wrapper to get a routine_id for sin()
 -- (because sin() is implemented in low-level assembly)</span>
 <span style="color: #008080;">function</span> <span style="color: #000000;">_sin</span><span style="color: #0000FF;">(</span><span style="color: #004080;">atom</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
     <span style="color: #008080;">return</span> <span style="color: #7060A8;">sin</span><span style="color: #0000FF;">(</span><span style="color: #000000;">a</span><span style="color: #0000FF;">)</span>
 <span style="color: #008080;">end</span> <span style="color: #008080;">function</span> 
  
 <span style="color: #004080;">atom</span> <span style="color: #000000;">a</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">0.0</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">1.0</span><span style="color: #0000FF;">,</span>
      <span style="color: #000000;">sinx</span> <span style="color: #0000FF;">:=</span> <span style="color: #000000;">quadAsr</span><span style="color: #0000FF;">(</span><span style="color: #000000;">_sin</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">1e-09</span><span style="color: #0000FF;">)</span>
 <span style="color: #7060A8;">printf</span><span style="color: #0000FF;">(</span><span style="color: #000000;">1</span><span style="color: #0000FF;">,</span><span style="color: #008000;">"Simpson's integration of sine from %g to %g = %f\n"</span><span style="color: #0000FF;">,</span> <span style="color: #0000FF;">{</span><span style="color: #000000;">a</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">b</span><span style="color: #0000FF;">,</span> <span style="color: #000000;">sinx</span><span style="color: #0000FF;">})</span>
<!--</syntaxhighlight>-->
{{out}}
<pre>
Simpson's integration of sine from 0 to 1 = 0.459698
</pre>

=={{header|Prolog}}==
{{works with|SWI-Prolog|9.1.2}}
{{works with|GNU Prolog|1.5.0}}

<syntaxhighlight lang="prolog">
%%% -*- 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).
</syntaxhighlight>

{{out}}
The number is printed differently by different Prolog implementations:
<pre>$ 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
</pre>

=={{header|Python}}==
<syntaxhighlight lang="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()
</syntaxhighlight>

=={{header|Raku}}==
(formerly Perl 6)
{{works with|Rakudo|2018.10}}
Fairly direct translation of the Python code.

<syntaxhighlight lang="raku" line>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";</syntaxhighlight>
{{out}}
<pre>Simpson's integration of sine from 0 to 1 = 0.459697694</pre>

=={{header|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.

<syntaxhighlight lang="ratfor">
#
# 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
</syntaxhighlight>

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

<pre>$ 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</pre>

=={{header|REXX}}==
{{trans|Go}}
<syntaxhighlight lang="rexx">/*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 #</syntaxhighlight>
{{out|output|text=&nbsp; when using the default inputs:}}
<pre>
Simpson's integration of sine from  0  to  1  =  0.459697694131860282602
</pre>

=={{header|Scheme}}==

For an R7RS-specific version, see [[#Ol]].

<syntaxhighlight lang="scheme">
(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)
</syntaxhighlight>

{{out}}
Running it with Racket's R5RS implementation:
<pre>$ plt-r5rs adaptive_simpson_task.scm
estimated definite integral of sin(x) for x from 0 to 1: 0.4596976941317858</pre>

Running it with SCM:
<pre>$ scm -f adaptive_simpson_task.scm
estimated definite integral of sin(x) for x from 0 to 1: 459.6976941317858e-3</pre>

Running it with Gambit Scheme:
<pre>$ gsi adaptive_simpson_task.scm
estimated definite integral of sin(x) for x from 0 to 1: .4596976941317858</pre>

=={{header|Sidef}}==
{{trans|Raku}}
<syntaxhighlight lang="ruby">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}"</syntaxhighlight>
{{out}}
<pre>
Simpson's integration of sine from 0 to 1 ≈ 0.45969769413186
</pre>

=={{header|Standard ML}}==
<syntaxhighlight lang="sml">
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";
</syntaxhighlight>

{{out}}
<pre>$ mlton adaptive_simpson_task.sml && ./adaptive_simpson_task
estimated definite integral of sin(x) for x from 0 to 1: 0.459697694132</pre>

=={{header|Wren}}==
{{trans|Go}}
{{libheader|Wren-fmt}}
<syntaxhighlight lang="ecmascript">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)</syntaxhighlight>

{{out}}
<pre>
Simpson's integration of sine from 0 to 1 = 0.459698
</pre>

=={{header|zkl}}==
{{trans|Python}}
<syntaxhighlight lang="zkl"># "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);
}</syntaxhighlight>
<syntaxhighlight lang="zkl">sinx:=quad_asr((1.0).sin.unbind(), 0.0, 1.0, 1e-09);
println("Simpson's integration of sine from 1 to 2 = ",sinx);</syntaxhighlight>
{{out}}
<pre>
Simpson's integration of sine from 1 to 2 = 0.459698
</pre>