Numerical integration/Adaptive Simpson's method

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.

Pseudocode: Simpson's method, adaptive

; Lyness's ASR, Modifications 1, 2, 3 procedure _quad_asr_simpsons(f, a, fa, b, fb) m := (a + b) / 2 fm := f(m) h := b - a return multiple [m, fm, (h / 6) * (f(a) + 4*fm + f(b))] procedure _quad_asr(f, a, fa, b, fb, tol, whole, m, fm, depth) lm, flm, left  := _quad_asr_simpsons(f, a, fa, m, fm) rm, frm, right := _quad_asr_simpsons(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 _quad_asr(f, a, fa, m, fm, tol', left , lm, flm, depth - 1) + _quad_asr(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 := _quad_asr_simpsons(f, a, fa, b, fb) return _quad_asr(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

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

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

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

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

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

Nim

Direct translation of Python code from Wikipedia entry.

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)

Output:

Simpson's integration of sine from 0 to 1 = 0.4596976941317859

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

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

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

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

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