
\ ---------------------------------------------------
\     (c) Copyright 2001  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, FLOAT EXT, FILE and TOOLS EXT wordsets.
\
\ Environmental dependences:
\       Assumes independent floating point stack
\       Uses FORmula TRANslator for clarity

\ Algorithm:
\
\   Assume roots are bracketed between xa and xb
\   and that f(xa)*f(xb) < 0   (at least 1 root in interval).
\   Then take next guess to be xp = (xa+xb)/2 and evaluate
\   f(xp}. If f(xa)*f(xp} < 0 then xb = xp, else xa = xp,
\   and continue until done.


MARKER -binroots

\ Data structures

FVARIABLE Ra                      \ f(xa)
FVARIABLE Rb                      \ f(xb)
FVARIABLE Rp                      \ f(xp)
FVARIABLE xa                      \ lower end of interval
FVARIABLE xb                      \ upper end of interval
FVARIABLE xp                      \ new guess
FVARIABLE epsilon                 \ precision

v: dummy                          \ dummy dictionary entry


: MoreThan    POSTPONE  F>  ;  IMMEDIATE

: initialize    ( xt --) ( F: lower upper precision --)
        f" epsilon="    f" xb="   f" xa="   \ store parameters
        defines dummy                       \ xt -> DUMMY
        f" Ra=dummy(xa)"
        f" Rb=dummy(xb)"
        f" MoreThan( Ra*Rb, 0)"             \ same sign?
        ABORT" Even # of roots in interval!"
;


: NotConverged?    ( -- f)
    F" MoreThan( ABS( xa - xb ), epsilon )"   ;

: NewPoint
        f" xp = (xa+xb)/2"      \ new point
        f" Rp = dummy(xp)"
        f" MoreThan( Ra*Rp, 0)" \ xp on same side of root as xa?

        IF      f" xa=xp"  f" Ra=Rp"
        ELSE    f" xb=xp"  f" Rb=Rp"   THEN
;

: )binroot  ( xt --)   ( f: Low High Precision -- root)
    initialize
    BEGIN   NotConverged?   WHILE   NewPoint   REPEAT
    f" (xa+xb)/2"           ( f: -- root)
;

FALSE [IF]
Usage example:

    : f1   fdup  fexp  f*  1e0  f-  ;  ok
    use( f1  0e0 2e0 1e-5 )binroot f. .567142  ok

[THEN]

