function y = edsqrt(x)
% EDSQRT  An awkward Taylor polynomial method to compute sqrt(x)
%         using only + and *.

if x < 0, error('EDSQRT(x) only works for x >= 0'), end
if x == 0, y = 0; return, end          % done with x = 0 case

% extract mantissa and exponent in base 2 scientific notation;
%   (in a "real" implementation this would be done by hardware)
[mant,expo] = log2(x);
mant = mant * 2;                       % adjust to match in-class:
expo = expo - 1;                       %    mantissa is in [1,2]
if mod(expo,2)==1                      % is exponent odd?
  expo = expo - 1;                     % fix exponent to be even
  mant = mant * 2;
end
x = mant;                              % now 1 <= x <= 4

% stored coefficients; computed by a PLOTPOLYSQRT run
c = [2.000000000000e+00  2.500000000000e-01  -1.562500000000e-02 ...
  1.953125000000e-03  -3.051757812500e-04  5.340576171875e-05 ...
  -1.001358032227e-05  1.966953277588e-06  -3.995373845100e-07   ...
  8.323695510626e-08  -1.768785296008e-08  3.818968252745e-09   ...
  -8.353993052879e-10  1.847517694387e-10  -4.123923424970e-11   ...
  9.278827706183e-12  -2.102234402182e-12  4.791857828503e-13   ...
  -1.098134085699e-13  2.528598223648e-14  -5.847383392186e-15   ...
  1.357428287472e-15  -3.162190896951e-16  7.389902639615e-17   ...
  -1.732008431160e-17  4.070219813225e-18  -9.588498598464e-19   ...
  2.263951057971e-19  -5.356669913948e-20  1.269900195117e-20   ...
  -3.016012963403e-21  7.175192130676e-22  -1.709713749888e-22   ...
  4.079998721323e-23  -9.749996944338e-24  2.333034983110e-24   ...
  -5.589562980367e-25  1.340739768939e-25  -3.219539576729e-26   ...
  7.739277828676e-27  -1.862263727525e-27  4.485330319344e-28   ...
  -1.081284987699e-28  2.608914359855e-29  -6.299935243968e-30   ...
  1.522484350626e-30  -3.682095304502e-31  8.911454061429e-32   ...
  -2.158242780502e-32  5.230435309891e-33  -1.268380562649e-33   ...
  3.077688129956e-34  -7.472271661672e-35  1.815198068755e-35   ...
  -4.411939750446e-36  1.072903530222e-36  -2.610412607014e-37   ...
  6.354293846021e-38  -1.547489665087e-38  3.770366768750e-39   ...
  -9.190268998828e-40  2.241069694386e-40  -5.467125665741e-41   ... 
  1.334239001758e-41  -3.257419437886e-42  7.955620550223e-43   ...
  -1.943702748066e-43  4.750467537250e-44  -1.161419453041e-44   ...
  2.840428010154e-45  -6.948904239127e-46];
y = c(1);
for n = 1:70
  y = y + c(n+1) * (x-4)^n;            % "^n" could be multiplication
end

y = y * 2^(expo/2);                    % divide exponent by 2 to sqrt

