\ 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
    adaptive
    f" finI"
    Ntimes @  . ."  function calls" ;