```\ Adaptive integration using Simpson's rule
\   with Richardson extrapolation

\ Integrate a real function from xa to xb

\ ---------------------------------------------------
\     (c) Copyright 2003  Julian V. Noble.          \
\       Permission is granted by the author to      \
\       use this software for any application pro-  \
\       vided this copyright notice is preserved.   \
\ ---------------------------------------------------

\ Usage:  use( fn.name xa xb err )integral
\ Examples:

\ 10 set-precision  ok
\ : f1   fdup  fsqrt  f*  ;  ok
\ use( f1 0e0 1e0 1e-6 )integral cr f. 41  function calls
\ .4000000556  ok
\ use( f1 0e0 1e0 1e-8 )integral cr f. 101  function calls
\ .4000000003  ok
\ use( fsqrt 0e0 1e0 1e-8 )integral cr f. 189  function calls
\ .6666666623  ok

\ Programmed by J.V. Noble; version of February 10th, 2003

\ This is an ANS Forth program requiring:
\      The FLOAT and FLOAT EXT word sets
\ Environmental dependencies:
\       Assumes independent floating point stack
\       Uses a FORmula TRANslator for clarity

MARKER  -int

needs ftran201.f    \ FORmula TRANslator

needs arrays.f      \ array lexicon

FALSE [IF]
Non STANDARD words:

f"           treat everything up to the closing "
as a Fortran formula, translate into
Forth, either interactively (display
output) or built into a compiling word.

f\$"          variant of f" except it evaluates the
formula, leaving the result on the fp stack.

1array       create a one-dimensional array
Ex: 20  1 FLOATS  1array A{
}            dereference a one-dimensional array

Vectoring words included in ftran111.f

v:           define a function vector
Ex:   v: dummy
defines      set a vector
Ex:   ' *  defines  dummy  7 2 dummy .  14 ok
: test    ( xt -- ) defines dummy  dummy  ;
3 5 ' * test .  15 ok
use(         get the xt of a word
[THEN]

\ Data structures

FVARIABLE  half         f\$" half = 0.5"
FVARIABLE  third        f\$" third = 1/3"
FVARIABLE  c43          f\$" C43 = 4/3"
FVARIABLE  fifteenth    f\$" fifteenth = 1/15"

40 CONSTANT Nmax
1 FLOATS  CONSTANT  float_len

Nmax 2* long float_len    1array x{
Nmax    long float_len    1array E{
Nmax 2* long float_len    1array f{
Nmax    long float_len    1array I{

0 VALUE  N

FVARIABLE  oldI
FVARIABLE  newI
FVARIABLE  finalI
FVARIABLE  deltaI

\ Begin program

: )int  ( n --)  LOCALS| nn |   \ Simpson's rule
f" I{ nn 2/ 1- } =
( (f{ nn } + f{ nn 2 - }) * third  + f{ nn 1- }*c43 )
* ( x{ nn } - x{ nn 1- } ) "  ;

v: dummy                                  \ dummy function name

VARIABLE Ntimes

: initialize  ( xt --)  ( f: xa xb eps -- integral)
LOCALS| xt |
2 TO n
xt defines  dummy

E{ 0 } F!                   \ initial error
x{ 2 } F!   x{ 0 } F!       \ calculate initial points
f" x{ 1 } = ( x{ 0 } + x{ 2 } ) * half"

f" f{ 0 } = dummy( x{ 0 } ) "   \ f_0, f_1 and f_2
f" f{ 1 } = dummy( x{ 1 } ) "
f" f{ 2 } = dummy( x{ 2 } ) "

2 )int                      \ I_0
f" finalI = 0"
3 Ntimes !  ;

: E/2    f" E{ N 2/ 1- } = half * E{ N 2/ 2 - } "  ;     \ eps -> eps/2

: move_down
f" x{ N } = x{ N 2 - }"      f" f{ N } = f{ N 2 - }"
f" x{ N 2 - } = x{ N 3 - }"   f" f{ N 2 - } = f{ N 3 - }"
;

: new_points
f" x{ N 1- } = half * ( x{ N } + x{ N 2 - } ) "
f" x{ N 3 - } = half * ( x{ N 4 - } + x{ N 2 - } ) "
f" f{ N 1- } = dummy( x{ N 1- } ) "
f" f{ N 3 - } = dummy( x{ N 3 - } ) "
2 Ntimes +!  ;

: N_too_big?       N  [ Nmax 1- ] LITERAL >
ABORT" Too many subdivisions!"  ;

: N=N+2   N 2 +   TO N    N_too_big?  ;

: N=N-4   N 4 -   TO N  ;

: subdivide
N=N+2
f" oldI = I{ N 2/ 2 - } "
E/2     move_down     new_points
N )int   N 2 - )int    ;

: LessThan  F< ;

: converged?   ( f: --)  ( -- f)
f" newI = I{ N 2/ 1- } + I{ N 2/ 2 - } "
f" deltaI = newI - oldI "
f" LessThan( abs(deltaI), E{ N 2/ 2 - })"  ;

: interpolate      \ use Richardson interpolation and store
f" finalI = deltaI * fifteenth + newI + finalI "
N=N-4  ;

: )integral    ( f: xa xb err -- I[xa,xb]) ( xt --)
initialize
BEGIN  N 0>  WHILE
subdivide
converged?  IF   interpolate   THEN
REPEAT
f" finalI"                         ( f: I[xa,xb] )
Ntimes @  .  ."  function calls"   \ optional line
;

```