\ 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
                    Ex: A{ I }  ( base_adr -- base_adr + offset )
 

       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
;