\ 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] = 10^(-k-1} where 10^k is the most significant
\   digit in a > 1. If a < 1, choose 10^m where 10^{-m} is the most
\   significant digit in a. Use a lookup table.


MARKER -1/a

include ftran111.f

FVARIABLE a
FVARIABLE x
FVARIABLE xp


\ With this table we can do numbers in the range from 10^-4 to 10^4.
create  powers  1e5 f, 1e4 f, 1e3 f, 1e2 f, 1e1  f,  1e0 f,
                1e-1 f,  1e-2  f,  1e-3  f,  1e-4  f,  1e-5  f, 1e-6 f,


: x0   ( f: a -- x0 )   \ initial guess
    pad 1 represent     ( -- exponent f1 f2)
    IF     SWAP  6 +  FLOATS   powers  +   F@   \ look up in table
    ELSE   ." Number not in range"  ABORT  THEN  ;

: 1/a  ( f: a -- 1/a)
    a F!
    f" xp = x0(a)"
    BEGIN   f" x = xp"
            f" xp = x * (2 - abs(a)*x)"
            f" abs(xp - x)"  1e-8  F<
    UNTIL
    f" (xp + x) /2"
    IF  FNEGATE  THEN  ;