function z=mysqrt(alpha);
% MYSQRT  A relatively robust implementation of Newton's method,
% and scaling by factors of 4, to get square roots.  Note that the maximum
% number of nontrivial multiplications or divisions is 13; the division or
% multiplication by powers of 2 can be implemented quickly in binary.

if (alpha<0)|(~isreal(alpha))
    error('mysqrt does not calculate complex results')
elseif alpha<10/inf, z=0; return
elseif alpha<1, z=1/mysqrt(1/alpha); return
else % at this point alpha>=1
    b=alpha; j=0;
    while b>1, b=b/4; j=j+1; end
    % now: alpha = b * 4^j  so sqrt(alpha) = sqrt(b) * 2^j
    % furthermore 1/4 < b <= 1; Taylor series gives start for Newton
    y=(1-b)/2;
    x=1-y*(1+.5*y*(1+y));  % x=1-(1/2)(1-b)-(1/8)(1-b)^2 -(1/16)(1-b)^3
    for k=1:10
        xnew=.5*(x+b/x);
        if abs(x-xnew)<eps, break, end
        x=xnew;
    end
    if k==10, error('max iterations exceeded in mysqrt'), end
    z=xnew*(2^j);
end