\ Brute force Monte Carlo integration in 1 dimension
\
\ ---------------------------------------------------
\     (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.   \
\ ---------------------------------------------------
\
\ Usage:  use( fn.name xa xb ya yb )bfmc
\ Examples:

\ use( fsqrt  10000 0e 2e 0e 2e fsqrt )bfmc fs. 1.88006E0  ok
\ use( fsqrt  10000 0e 2e 0e 2e fsqrt )bfmc fs. 1.90806e0  ok
\ use( fsqrt  10000 0e 2e 0e 2e fsqrt )bfmc fs. 1.88486e0  ok
\ use( fsqrt  10000 0e 2e 0e 2e fsqrt )bfmc fs. 1.89363E0  ok

\ : f1   fdup  fsqrt  f*  ;  ok
\ use( f1  10000 0e 2e 0e 2e f1 )bfmc fs. 2.28933e0  ok
\ use( f1  10000 0e 2e 0e 2e f1 )bfmc fs. 2.25313e0  ok
\ use( f1  10000 0e 2e 0e 2e f1 )bfmc fs. 2.26444e0  ok
\ use( f1  10000 0e 2e 0e 2e f1 )bfmc fs. 2.25935E0  ok

MARKER -bfmc

\ Conditional definition of non-Standard words
: undefined    BL WORD  FIND  NIP  0=  ;

undefined  prng  [IF]  include  prng.f  [THEN]


undefined s>f   [IF] : s>f   S>D  D>F  ;    [THEN]
undefined f^2   [IF] : f^2   FDUP  F*  ;    [THEN]
undefined ftuck [IF] : ftuck  FSWAP  FOVER  ;   [THEN]

undefined  use( [IF]
\ Vectoring: for using function names as arguments
: use(      '       \ state-smart ' for syntactic sugar
    STATE @  IF  POSTPONE LITERAL  THEN  ;  IMMEDIATE

' NOOP  CONSTANT  'noop
: v:   CREATE  'noop  ,  DOES> PERFORM  ;   \ create dummy def'n
: 'dfa   ' >BODY  ;                         ( -- data field address)
: defines    'dfa   STATE @
             IF   POSTPONE  LITERAL    POSTPONE  !
             ELSE   !   THEN  ;  IMMEDIATE
\ end vectoring
[THEN]

\ Data structures
    v: fdummy

    1000 VALUE  Nmax
    0 VALUE Npoints
    0 VALUE Nhits

    0.1 seed 2!

    FVARIABLE xa    FVARIABLE xb-xa
    FVARIABLE ya    FVARIABLE yb-ya

\ Program begins here

: x     ( f: -- x = xa + xi*[xb-xa])    \ guess a new point
    prng   xb-xa F@  F*   xa F@   F+  ;

: y     ( f: -- y = ya + xi*[yb-ya])    \ guess a new point
    prng   yb-ya F@  F*   ya F@   F+  ;


: new_point
    Npoints  1+  TO Npoints
    x fdummy
    y F>    IF  Nhits  1+  TO Nhits  THEN  ;

: initialize       ( xt n --)   ( f:  xa xb ya yb --)
    TO Nmax     0 TO Npoints    0 TO Nhits
    defines fdummy
    FOVER  F-  yb-ya F!     ya F!
    FOVER  F-  xb-xa F!     xa F!   ;


: )bfmc   ( xt --)   ( f:  xa xb ya yb -- integral)
    initialize
    BEGIN   Npoints Nmax <
    WHILE   new_point
    REPEAT
    Nhits  s>f   Npoints  s>f   F/
    xb-xa F@   yb-ya F@  F*  F*  ;