\ Adaptive integration using 3 point Gauss-Legendre
\   with Richardson extrapolation

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

\ 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

\ Usage:  use( fn.name xa xb err )integral
\ Examples:
\   12 set-precision  ok
\   use( fsqrt 0e 1e 1e-8 )integral cr fs. 609  function calls
\   6.66666666670E-1  ok
\
\   : x^1.5   fdup  fsqrt  f*  ;  ok
\   use( x^1.5  0e 1e 1e-8 )integral cr fs. 165  function calls
\   3.99999999994E-1  ok

MARKER  -int

needs ftran109.f

needs arrays.f

FALSE [IF]
    Non STANDARD words:

       1array       create a 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 ftran109.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]

\ points and weights for 3 point Gauss-Legendre integration
    8e  9e   F/         FVARIABLE  w0   w0 F!
    5e  9e   F/         FVARIABLE  w1   w1 F!
    9e  15e  F/ FSQRT   FCONSTANT  x1

\ Data structures

64e  63e  F/  FVARIABLE Cinterp   Cinterp F!

20 CONSTANT Nmax
1 FLOATS  CONSTANT  float_len

Nmax long float_len    1array x{
Nmax long float_len    1array E{
Nmax long float_len    1array I{

0 VALUE  N
0 VALUE  nn

FVARIABLE  oldI
FVARIABLE  finI
FVARIABLE  deltaI

FVARIABLE dx
FVARIABLE dxi
FVARIABLE xi

\ Begin program

: scale    ( f: xa xb -- )
    FOVER    F-    F2/      ( f: [xb-xa]/2)
    FSWAP  FOVER  F+        ( f: [xb-xa]/2  [xa+xb]/2 )
    xi F!    FDUP  dx F!   x1 F*  dxi F!  ;

VARIABLE Ntimes             \ count of function evaluations

v: fdummy

: )int  ( n -- )     \ 3 point Gauss-Legendre
    TO nn
    x{ nn 1- }  F@  x{ nn }  F@   scale
    f" I{nn_1-}=dx*(w0*fdummy(xi)+w1*(fdummy(xi+dxi)+fdummy(xi-dxi))) "
    3 Ntimes +! ;

: initialize  ( xt --)  ( f: xa xb eps -- integral)
     defines  fdummy
     1 TO N
     E{ 0 } F!   X{ 1 } F!   X{ 0 } F!
     0 Ntimes !
     1 )int
     0e  finI F!
;

: check.n       N  [ Nmax 1- ] LITERAL >
        ABORT" Too many subdivisions!"  ;

: E/2   E{ N  1- }  DUP   F@   F2/    F! ;

: }down    ( adr n --)
        OVER @  >R   }   DUP   R@ +   R>   MOVE  ;

: move.down    E{ N  1- }down
               x{ N     }down ;

: N+1   N 1+   TO N  ;
: N-2   N 2 -  TO N  ;

: subdivide
        check.n     E/2   move.down
        f" oldI = I{N_1-} "
        f" x{N} + x{N_1-} "  F2/  x{ N } F!
        N )int   N 1+ )int    ;

: converged?   ( f: --)  ( -- f)
        f" I{N} + I{N_1-} - oldI "
        FDUP   deltaI  F!   FABS
        E{ N 1- } F@   F2*    F<  ;

: interpolate    f" finI = deltaI * Cinterp + oldI + finI "  ;


: )integral    ( f: xa xb err -- I[xa,xb]) ( xt --)
     initialize
     BEGIN   N 0>   WHILE
        subdivide
        converged?    N+1
        IF    interpolate  N-2    THEN
     REPEAT   finI  F@
     Ntimes @  .  ."  function calls" ;