\ 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.


MARKER -1/a

include ftran111.f

FVARIABLE a
FVARIABLE x
FVARIABLE xp

FVARIABLE  dix      10.e0  dix      F!
FVARIABLE  dixieme  0.1e0  dixieme  F!

: x0   ( f: a -- x0 )   \ initial guess
    f" x = 1"
    FDUP  FABS a F!  1e0 F>
    IF      BEGIN   f" x = x * dixieme"
                    f" 2 - a * x"   F0>
            UNTIL   x F@
    ELSE
            BEGIN   f" x = x * dix"
                    f" 2 - a * x"   F0<
            UNTIL   f" x * dixieme"
    THEN   ;


: 1/a  ( f: a -- 1/a)
    FDUP  F0<
    x0  xp  F!
    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  ;