\ 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

needs ftran111.f

v: fdummy


WORDLIST  CONSTANT  qmin                \ create separate wordlist
qmin SET-CURRENT                        \ for quad min def'ns
GET-ORDER  qmin   SWAP 1+  SET-ORDER    \ make qmin findable

: fvariables     0 DO  FVARIABLE  LOOP  ;

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



: moments
    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
    xd F@  xa F@  F<  AND                   \ xd < xa  ?
    2
    xa F@  xb F@  xd F@  fwithin AND  +     \ xa <= xd <= xb ?
    3
    xb F@  xc F@  xd F@  fwithin 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   xc F@  xa F@  F-  FABS  1e-6 F>  ;


100 VALUE Nmax

FORTH-WORDLIST   SET-CURRENT            \ definitions to FORTH

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

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

    CR  ." x_min: " xb f@  f.       \ output results
    CR  ." f(x_min): " fb f@  f.
    CR  ." number of iterations needed: " n .
;

: )tabulate     ( xt Npts --)    ( f: a b--)
    LOCALS| Npts |
    defines fdummy
    FOVER FOVER  F-
    Npts  NEGATE  S>D  D>F  F/   ( f: a b dx)
    xc F! xb F! xa F!
    f" xd=xa"
    BEGIN   xd F@  xb F@  F<
    WHILE   f" fdummy(xd)"   f" xd"  cr  f. f.
            f" xd = xd+xc"
    REPEAT
;   \ say use( fn.name a b )tabulate

GET-ORDER  NIP  1-  SET-ORDER           \ hide qmin definitions


FALSE [IF]

: 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+  ;


[THEN]