```\ Integration by recursion -- an illustration of the concept.
\   Uses 3 point Gaussian integration
\   version of January 21st, 2003 - 10:51
\
\   Has a Fortran-ish interface for clarity
\
\ ---------------------------------------------------
\     (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.   \
\ ---------------------------------------------------

FALSE [IF]

Algorithm:

c = (a+b)/2
integral(a,b) = integral(a,c) + integral(c,b)

uses 3 point Gauss-Legendre with Richardson extrapolation

Environmental dependencies: FLOAT, separate floating point stack

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

use( fsqrt 0e 1e 1e-4 )integral fs. 6.66666744641E-1  ok
use( fsqrt 0e 1e 1e-6 )integral fs. 6.66666670150E-1  ok

: f1     FDUP FSQRT F*  ;  ok
use( f1 0e 1e 1e-4 )integral fs. 3.99994557189E-1  ok
use( f1 0e 2e 1e-6 )integral fs. 2.26274169632E0  ok

[THEN]

MARKER  -int

include ftran201.f      \ be sure to review ftrandoc.txt !!

\ points and weights for 3 point Gauss-Legendre integration
FVARIABLE  w0   f\$" w0 = 8 / 9 "
FVARIABLE  w1   f\$" w1 = 5 / 9 "
FVARIABLE  x1   f\$" x1 = sqrt( 9 / 15 )"

v: fdummy

FVARIABLE  xa       FVARIABLE  xb
FVARIABLE  x0       FVARIABLE  dx

: scale    ( f: xa xb -- )
f" xb ="   f" xa ="
f" dx = 0.5 * ( xb - xa )"
f" x0 = xa + dx"
;

VARIABLE Ntimes             \ count of function evaluations

: gauss3  ( f: xa xb -- gauss_int)    \ 3 point Gauss-Legendre
scale
f" xa = x0 - dx * x1"
f" xb = x0 + dx * x1"
f" dx * ( w0 * fdummy(x0) + w1 * ( fdummy( xa ) + fdummy( xb ) ) )"
3 Ntimes +!
;

FVARIABLE  oldI
FVARIABLE  newI
FVARIABLE  deltaI

FVARIABLE  xa
FVARIABLE  xb
FVARIABLE  xc
FVARIABLE  err

: juggle   ( f: -- xa xc err/2  xc xb err/2 )
f" xa"
f" xc"
f" 0.5 * err"
FOVER  FOVER
f" xb"  FSWAP  ;

: LessThan   F<     ( f: x y --)  ( -- flag)  ;

FVARIABLE  finI

: adaptive      ( f: xa xb err -- result)
f" err ="   f" xb ="   f" xa ="
f" xc = 0.5 * ( xa + xb )"
f" oldI = gauss3( xa, xb )"
f"  newI = gauss3( xa, xc ) + gauss3( xc, xb )"
f" deltaI = newI - oldI"
f" LessThan( abs(deltaI), err)"
IF   f" finI = finI + newI + deltaI / 63"
ELSE juggle  RECURSE  RECURSE  THEN  ;

: )integral     ( xt --)    ( f: xa xb err -- integral)
defines fdummy          \ pass function xt to fdummy
f" finI = 0"            \ initialize integral
0 Ntimes !              \ initialize count