\ Quadratic search minimization in one dimension

FALSE [IF]

---------------------------------------------------
    (C) Copyright 2000  Julian V. Noble.
      Permission is granted by the author to      
      use this software for any application pro-  
      vided this copyright notice is preserved,
      as per GNU Public License agreement.
---------------------------------------------------

This is an ANS Forth compatible program with the following
environmental dependence:

      ANS FLOAT and FLOAT EXT wordsets
      ANS TOOLS EXT wordsets

      Assumes separate floating point stack

Algorithm:

    Given an interval [xa,xc] believed to contain a minimum, let

        xb = (xa+xc)/2

    and calculate fa = f(xa), etc.

    We approximate f(x) by

        f(x) ~ [ fa * (x-xb) * (x-xc) * (xc-xb) + cyc perm ] /
                        [ (xa-xb) * (xb-xc) * (xc-xa) ]


    Calculate

        alpha = fa * (xc-xb) + cyclic perms

        beta  = fa * (xc+xb) * (xc-xb) + cyclic perms

        x' = 0.5*beta/alpha


[THEN]

marker -quad

include ftran201.f

: fvariables     0 DO  FVARIABLE  LOOP  ;

13 fvariables  xa fa xb fb xc fc xd fd  alpha beta ma mb mc

v: fdummy

: ToFloat    ( n -- )  ( f: -- n 0)  \ use as ToFloat(0)
    fdrop  S>D  D>F  ;


: LessThan      ( f: a b --)    ( -- flag)  F<  ;

: MoreThan      ( f: a b --)    ( -- flag)  F>  ;

: report    f" fdummy(xd)"   f" xd"   cr  f. f.  ;

: step      f" xd = xd+xc" ;

: )tabulate     ( xt Npts --)    ( f: a b --)
    LOCALS| Npts |          \ initialize
    defines fdummy
    f" xb="  f" xa="

    Npts   f" xc = (xb - xa)/ToFloat(0)"    \ xc =step
    f" xd=xa"

    BEGIN   f" LessThan(xd,xb)"     \ xd < xb
    WHILE   report  step    REPEAT
;   \ say use( fn.name a b )tabulate

: new_guess
    f" ma = fa * (xc-xb)"
    f" mb = fb * (xa-xc)"
    f" mc = fc * (xb-xa)"

    f" alpha = (ma + mb + mc)*2 "
    f" beta  = ma * (xc+xb) + mb * (xa+xc) + mc * (xb+xa)"

    f" xd = beta / alpha"
    f" fd = fdummy(xd) "
;

: fwithin    ( f: a b x --) ( -- true if a <= x <= b)
    FSWAP  FOVER   FMIN  FROT FROT   FMAX   F=  ;

: shuffle   \ rearrange the points and function values
    \ tests
    1  f" LessThan(xd, xa)"    AND          \ xd < xa  ?
    2  f" fwithin(xa, xb, xd)  AND  +       \ xa <= xd <= xb ?
    3  f" fwithin(xb, xc, xd)  AND  +       \ xb <= xd <= xc ?
    CASE
      1 OF  f" xc = xb"     f" xb = xa"     f" xa = xd"
            f" fc = fb"     f" fb = fa"     f" fa = fd"
        ENDOF


      2 OF  f" xa = xd"     f" xc = (xb + xc) / 2"
            f" fa = fd"     f" fc = fdummy(xc)"
        ENDOF


      3 OF  f" xc = xd"     f" xa = (xa + xb) / 2"
            f" fc = fd"     f" fa = fdummy(xa)"
        ENDOF

        \ otherwise:

            f" xa = xb"     f" xb = xc"     f" xc = xd"
            f" fa = fb"     f" fb = fc"     f" fc = fd"
    ENDCASE
;


: apart     1e-6  f" MoreThan(abs( xc - xa ))"  ;


100 VALUE Nmax

: print_results    ( n --)   ( f: f[x] x -- )
    CR  ." x_min: "       f.
    CR  ." f(x_min): "    f.
    CR  ." number of iterations needed: "  .

: )quadmin   ( xt --)    ( f: xa xc --)
    0  LOCALS| n |               \ initialize
    defines  fdummy
    f" xc="  f" xa="
    f" xb = (xa+xc)/2"
    f" fa = fdummy(xa)"
    f" fb = fdummy(xb)"
    f" fc = fdummy(xc)"

    BEGIN  apart  n Nmax <  AND     \ iterate
    WHILE  new_guess   shuffle
           n 1+ TO n
    REPEAT

    n  f" fb"  f" xb"  print_results
;

: f1   ( f: x -- 1 -  x / e^x )
    FDUP  FEXP  F/  FNEGATE  1e  F+  ;

: f2   ( f: x -- [[x-1]^2-0.75]*[[x-2]^2-2] )
    FDUP  1e F-  FDUP F*  0.75e F-
    FSWAP 2e F-  FDUP F*  2e F-  F*  ;

: f3    ( f: x -- [x^2-2*x+.25]*[x^2-4*x+2]^3 )
    FDUP FDUP  2E0  F-  F*  0.25E0  F+
    FSWAP
    FDUP  4E0 F-  F*  2E0 F+  FDUP FDUP F* F*
    F*  ;

: f4  fdup  1e0 f+  fover f*  fover f*  8e0 f+ f*  8e0 f+  ;