\ Inverse of a real number without division
\   version of September 2nd, 1999 - 20:54

\ ---------------------------------------------------
\     (c) Copyright 1999  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

\ Algorithm:
\
\   Newton's method applied to  a - 1/x :
\
\       x[n+1] = x[n] * ( 2 - a * x[n] )
\
\   converges if for a > 1, x[0] < 1 and vice-versa.
\
\   Choose x[0] = 2^(-k-1} where 2^k is the most significant
\   bit in a > 1. If a < 1, choose 2^m where 2^{-m} is the most
\   significant bit in a. Use a machine code (assembly language)
\   subroutine to do this efficiently.


MARKER -1/a

include ftran111.f

FVARIABLE a
FVARIABLE x
FVARIABLE xp


CODE   guess    ( 87: a -- x0)
    fxtract             \ get base-2 exponent and mantissa
    fstp st(0)          \ discard mantissa
    fchs                \ negate exponent
    fld1                \ load 1e0
    fsub st(1), st(0)   ( 87: -k 1 --  -k-1 1)
    fscale              ( 87: --  -k-1 2^[-k-1] )
    fxch                \ swap st1 and st0
    fstp  st(0)         \ discard -k-1 ( 87: -- 2^[-k-1])
    next,
END-CODE

: x0   fpop     \ fp# from software stack to fpu stack
       guess    \ do machine code subroutine
       fpush    \ fp# from fpu stack to software stack
;

: 1/a  ( f: a -- 1/a)
    FDUP  F0<       \ set sign flag
    FABS  a F!      \ save |fp#|
    f" xp = x0(a)"  \ compute first guess
    BEGIN   f" x = xp"                  \ begin loop
            f" xp = x * (2 - a * x)"    \ Newton
            f" abs(xp - x)"  1e-8  F<   \ stopping criterion
    UNTIL                               \ end loop
    f" xp + x * (1 - a * x)^2"          \ polish answer
    IF  FNEGATE  THEN                   \ restore sign
;