function inverse = NewtonReciprocation1(a)
% Compute the reciprocal of a given nonzero number w/o division
% Uses Newton-Raphson method
% NewtonReciprocation1.m
% Usage: inverse = NewtonReciprocation1(a)
% K. Ming Leung, 02/24/03

pt5=0.5; one=1; tol=eps; maxIt=50; it=0;
% Take out the sign of a so that it is always positive
sgn=sign(a); a=abs(a); 

% Scale a until it lies between 0.5 & 1
% Should really be done by bit shifting
s=1;
if a < pt5
    while a < pt5
        s=2*s; a=2*a;
    end
elseif a > one 
    while a > one
        s=pt5*s; a=pt5*a; 
    end
end
a,
% Start Newton-Raphson iteration with initial guess:
%x=1.5  % pick the arith. mean of the interval
%x=pt5*a,   % Crandall's initial guess (10 iterations for the worst case
if a >= 0.66666
    x=2-a;      % My initial guess (6 iterations for the worst case)
else
    x=4*(1-a)
end

x=1.333333
%x=(a-3)*a+3, 
dx=x*(1-a*x);
while abs(dx) > tol & it < maxIt
    x=x+dx;
    dx=x*(1-a*x);
    it=it+1; [it x],
end
% Put back the sign and the overall scale factor
inverse=sgn*s*x;
