I'm working on modernizing Rosetta Code's infrastructure. Starting with communications. Please accept this time-limited open invite to RC's Slack.. --Michael Mol (talk) 20:59, 30 May 2020 (UTC)

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.

Lychee (1969)'s 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.

 ```; Lychee'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) + f(b) + 4*sum1 + 2*sum2)] 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) ```

## 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.rangessequences ;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 ] [email protected] [ 4 ] [ 2 ] bi*    [ * ] [email protected] a b [ f call ] [email protected] + + + 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

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 * send 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.absimport kotlin.math.sin typealias F = (Double) -> Doubletypealias 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
```

## 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
`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 = call_func(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 := call_func(f,{a}),         fb := call_func(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(routine_id("_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 collectionstriple = 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      ' = '       sinxexit                                             /*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 // recursivequadAsrRec = 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    var flm   = r1    var left  = r1    var rm    = r2    var frm   = r2    var right = r2    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    var fm    = r    var whole = r    return quadAsrRec.call(f, a, fa, b, fb, eps, whole, m, fm)} var a = 0var b = 1var 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 Racketfcn _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
```