The following Python program implements Newton’s method for computing the square root of a number:

def sqrt(x):
    def sqrt_iter(guess):
        return guess if good_enough(guess) else sqrt_iter(improve(guess))
    def good_enough(guess):
        tolerance = 1e-3
        return abs(guess**2 - x) < tolerance
    def improve(guess):
        return (guess + x/guess)/2
    initial_guess = 1.0
    return sqrt_iter(initial_guess)



Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 12, in sqrt
  File "<stdin>", line 3, in sqrt_iter
  File "<stdin>", line 3, in sqrt_iter
  File "<stdin>", line 3, in sqrt_iter
  [Previous line repeated 993 more times]
  File "<stdin>", line 6, in good_enough
RecursionError: maximum recursion depth exceeded while calling a Python object

As we can see, this naive program does not perform well:

  • for small numbers (from $x = 10^{-4}$ and below), the tolerance $10^{-3}$ is too large so the result is inaccurate;
  • for large numbers (from $x = 10^{13}$ and above), the tolerance $10^{-3}$ is too small so the program enters an infinite recursion.

Both problems can be solved by redefining the good_enough procedure like this:

def good_enough(guess):
    return guess == 0 or improve(guess) == guess

But discussing the various solutions is not the point of my post. Instead, I would like to reliably predict for a given $x$ if the above naive program will return.

I have not read the IEEE Standard for Floating-Point Arithmetic (IEEE 754), but to my understanding, floating-point numbers are not uniformly distributed on the real line. Their spacing is large for small numbers and small for large numbers (this Wikipedia figure seems to confirm this). In other words, small floating-point numbers are dense and large floating-point numbers are sparse. The consequence of this is that the naive program will enter an infinite recursion if guess has not reached the tolerance range yet and the improve procedure cannot improve the guess anymore (meaning that a fixed point of the improve procedure has been reached), because the new guess would be off the old guess by a distance below the spacing of the old guess.

So to guarantee that the naive program will return for a given $x$, it seems to me intuitively that this predicate must hold:

tolerance > spacing($\sqrt{x}$).

If we choose a tolerance of $10^{-3}$ like in the naive program, that means that the spacing of $\sqrt{x}$ should be less than $10^{-3}$. Consequently, according to the Wikipedia figure above for binary64 floating-point numbers, $\sqrt{x}$ should be less than $10^{13}$ and therefore $x$ should be less than $10^{26}$.

Here are my questions:

  1. Is my predicate correct?
  2. Why does the program enter an infinite recursion from $x = 10^{13}$ whereas my predicate guarantees that it cannot happen below $x = 10^{26}$?

Note. — I am using the CPython interpreter on a 64-bit MacBook Pro so the IEEE 754 binary64 format.

  • $\begingroup$ If you formulate your algorithm without using scheme, then you much more likely to get an answer. Your use of the word precision is wrong. Precision refers to the (relative) accuracy of the basic arithmetic operations. I believe you want to write "error" or "absolute error" instead of "absolute precision". Relevant information about the behavior of Newton's method for square roots and the construction of initial guesses can be found here $\endgroup$ Jan 28 '20 at 14:06
  • $\begingroup$ @CarlChristian Which language do you recommend? Actually I took the phrase "absolute precision" from the IEEE 754 Wikipedia article of the figure above (unfortunately it is not defined). How do you understand the quantity plotted in the figure? My qualitative understanding was that small floats are denser than large floats in the IEEE 754 representation. (Thanks for the linked answer by the way.) $\endgroup$
    – Maggyero
    Jan 28 '20 at 15:24
  • 1
    $\begingroup$ I would avoid using any language other than mathematics. If absolutely necessary I would write Matlab as this language is at least understood by everybody in the numerical linear algebra community. The figure plots the upper bound on the absolute error when representing a real number. The label is wrong, but this is a common mistake. You qualitative understanding is correct. The value of the floating point systems is that we have a uniform bound on the relative error for all reals in the representational range. $\endgroup$ Jan 28 '20 at 15:53
  • $\begingroup$ @CarlChristian Scheme is one of the simplest language and I don’t use Matlab. I could rewrite the program in Python though (but I am not convinced that it is a programming language problem). It seems to me that it is more likely that I did not expose the problem very well. For instance, do you understand the 3 questions or it is not clear what I am looking for? $\endgroup$
    – Maggyero
    Jan 28 '20 at 17:29
  • $\begingroup$ Simplicity in the eye of the beholder. You want people to read your questions. As it is currently written your post is a massive wall of code and text. I expect this to reduce the number of readers which make it to the bottom. You could increase this number by simply writing out the iteration you have implemented and stating your convergence criteria. That is about two lines of mathematics. Yes, I understand your questions, but I can't help with your program because I do not understand scheme. $\endgroup$ Jan 28 '20 at 17:57

Your program stalls because of the limitations of floating-point arithmetic. In general, when computing approximations $x$ of $\sqrt{\alpha}$, the computed value $\hat{y}$ of the residual $y = x^2 - \alpha$ satisfies $$|y - \hat{y}| \leq \gamma_2 (x^2 + \alpha),$$ where $\gamma_k = \frac{ku}{1 - ku}$ and $u$ is the unit roundoff.* When we get very close to the target, i.e. when $x \approx \sqrt{\alpha}$ and $y \approx 0$, the computed value of the residual will satisfy $$|\hat{y}| \lesssim 4u\alpha.$$ For $\alpha = 10^{13}$ and IEEE double precision floating-point arithmetic, the right hand side is about $4 \times 10^{-3}$ which is larger than the (absolute) tolerance which we are currently using. We can resolve this issue by terminating the iteration when $|\hat{y}| \leq \tau \alpha$, where $\tau$ is some user defined tolerance level. We should be able to do $\tau \approx 10^{-15}$ on machines which have IEEE double precision floating-point arithmetic. It is worth recognizing that we obtain an accurate bound for the relative error because $$\frac{x - \sqrt{\alpha}}{\sqrt{\alpha}} = \frac{x^2 - \alpha}{\sqrt{\alpha} (x + \sqrt{\alpha})} \approx \frac{1}{2} \frac{x^2 - \alpha}{\alpha}$$ is a good approximation when $x$ is close to $\sqrt{\alpha}$.

Edit. — Let $x$ and $\alpha$ denote nonnegative floating-point numbers. In the absence of floating-point overflow/underflow, the computed value $\hat{y}$ of the residual $y = x^2 - \alpha$ can be written as $$\hat{y} = (x^2(1 + \epsilon) - \alpha) (1 + \delta) = y + x^2 (\epsilon + \delta + \epsilon\delta) - \alpha \delta,$$ where $|\epsilon| \leq u$ and $|\delta| \leq u$. It follows that $$|y - \hat{y}| \leq x^2(2u + u^2) + \alpha u.$$ When we get very close to the target, i.e. when $x \approx \sqrt{\alpha}$ and $y \approx 0$, $$|\hat{y}| \lesssim \alpha(2u + u^2) + \alpha u \approx 3u \alpha.$$ From a purely mathematical point of view, the fact that the right-hand side of this inequality exceeds the threshold does not imply that the left-hand side of the of the inequality will also exceed the threshold. Experience, however, suggests, that once we no longer have a mathematical guarantee that the machine will behave, it will seize the opportunity to misbehave.

A closer analysis is possible using Sterbenz lemma regarding the subtraction of floating-point numbers which are sufficiently close to each other. When $x^2$ is a good approximation of $\alpha$, there is no error in the subtraction when computing the residual. It follows that $$\hat{y} = x^2 (1 + \epsilon) - \alpha = y + x^{2} \epsilon \approx x^2 \epsilon.$$ The very best we can hope for is that $|\epsilon| \leq u$, and while it does happen that $\epsilon = 0$, we have no right to expect that $|\epsilon| \ll u$ and in general this is not the case and $|\hat{y}| \approx \alpha u$, something which we can easily verify.

* The unit roundoff $u$ is $u = 2^{-53}$ in IEEE double precision arithmetic, and $u = 2^{-24}$ in IEEE single precision arithhmetic. Machine epsilon $\text{eps}$ is the distance between $1$ and the next floating-point number $1 + 2u$, i.e. $\text{eps} = 2u$. Some authors will use the terms as if they are identical. They are not. When in doubt, follow the notation of Niclas J. Higham’s textbook Accuracy and Stability of Numerical Algorithms. The $\gamma$-factor which I used earlier is also explained in this textbook. Briefly, it is a convenient tool which allows us to simplify rounding error analysis considerably at the cost of a bit of exactness. It was not strictly necessary here, so I decided not to pursue it in these edit.

  • $\begingroup$ Thanks for this quantitative answer. I am trying to follow your derivation but I have some trouble doing so. How do you derive the first inequation? What you call the unit roundoff is the machine epsilon $2^{-53} \approx 1.11\times10^{-16}$ for the binary64 format, right? (Then you may give its value because you use it in the second inequation to get $4 \times 10^{-3}$.) How do you derive the second inequation from the first? How the fact that the right-hand side of the second inequation is larger than the tolerance makes also the left-hand side larger than the tolerance? $\endgroup$
    – Maggyero
    Jan 28 '20 at 22:04
  • 1
    $\begingroup$ @Maggyero I added paragraphs which should address your comments. $\endgroup$ Jan 29 '20 at 12:45
  • 1
    $\begingroup$ @Maggyero: We can make no assumption about the sign of $\epsilon$ and $\delta$. We only have a bound on the absolute value of both. This limits the class of reductions that we can do. If you have another look, you will see that I am in fact not adding $+\alpha \delta$, but the non-negative number $\alpha u$. $\endgroup$ Jan 30 '20 at 22:18
  • 1
    $\begingroup$ @Maggyero. The best we can hope for is $x = \text{fl}(\sqrt{\alpha})$ the floating point representation of $\sqrt{\alpha}$. We have $x = \sqrt{\alpha}(1+\delta)$ for $|\delta| \leq u$. The very best we have hope for is $\hat{y} = \text{fl}(x^2-\alpha) = \alpha(1+\delta)^2-\alpha$. This assumes that $\alpha$ is floating point number, that the multiplication is exact and it applies Sterbenz's lemma for the subtraction. It follows that $\hat{y} = (2\delta + \delta^2)\alpha$ and $|\hat{y}| \leq 2u \alpha$ is really the best that we can hope for. Realistically, the multiplication is not exact. $\endgroup$ Jan 19 at 15:30
  • 1
    $\begingroup$ @Maggyero I think that there is an issue that should be addressed before you do further development. As it is currently written your program will misbehave when your guess is larger than $\sqrt(\Omega)$ where $\Omega$ is the largest positive floating point number. $\Omega$ is slightly less than $2^{1024}$ is double precision, so a guess of $2^{512}$ will trigger an overflow and a NaN is sure to follow. You can reduce the problem of computing a square root down to the interval $[1,4]$ if you have access to the binary representation of the input. $\endgroup$ Feb 4 at 22:00

First off, it's not nice to your computer to torture it with such a bad initial guess. In python it seems you can provide a guess good enough to converge in $4$ iterations:


                            Online Python Interpreter.
                Code, Compile, Run and Debug python program online.
Write your code in this editor and press "Run" button to execute it.

import struct
import math
x = float(input('Enter a number:> '));
temp = struct.unpack('L',struct.pack('d',x))[0]
magic = 0X3FEED9EBA16132A9;
y = struct.unpack('d',struct.pack('L',(magic+temp) >> 1))[0]
print(y, math.sqrt(x))
print(y, math.sqrt(x))
print(y, math.sqrt(x))
print(y, math.sqrt(x))
print(y, math.sqrt(x))

This interpreter produced the output

Enter a number:> 1.0e49                                                                                                                                  
3.190080382312335e+24 3.162277660168379e+24                                                                                                              
3.1623988156356365e+24 3.162277660168379e+24                                                                                                             
3.162277662489188e+24 3.162277660168379e+24                                                                                                              
3.162277660168379e+24 3.162277660168379e+24                                                                                                              
3.162277660168379e+24 3.162277660168379e+24                                                                                                              

...Program finished with exit code 0                                                                                                                     
Press ENTER to exit console.                                                                                                                             

But going back to your questions it should be clear that your predicate isn't good enough. You are comparing $\text{guess}^2$ to $x$ so the roundoff happens at the level of $x$, not $\sqrt x$ so you want a condition like $\text{tolerance}>\text{spacing}(x)$, that is to say the distance between the floating point representation of $x$ and the next representable number. IEEE-754 double precision floating point numbers have $52$ mantissa bit plus an implicit leading $1$ bit so the least significant bit is $2^{-52}$ times the most significant bit and the number is less than $2$ times the most significant bit so the spacing of floating point numbers is between $2^{-52}$ and $2^{-53}$ times the magnitude of the number itself. To clarify: if $x=1+2^{-52}$ then $\text{spacing}(x)/x\approx2^{-52}$ but if $y=2-2^{52}$ then $y$ has the same exponent as $x$ so $\text{spacing}(x)=\text{spacing}(y)$ but $\text{spacing}(y)/y\approx2^{-53}$. It's worse for bigger radices: both $9.999$ and $1.001$ have $3$ decimal place accuracy but the relative uncertainty in the second number is about $10\times$ the first. $10^{13}/2^{52}=0.00222$ so the spacing is bigger than the tolerance for sure. If the best possible square root squares out to a number than is different by $1$ Unit in the Last Place (ulp) from $10^{13}$ it will not meet your tolerance condition.

The situation is different for even powers of $10$. An IEEE-754 double precision floating point number can represent exactly any number of the for $2^n*q$ where $-1022\le n\le1023$ so it's not an issue but $q$ is represented by a $53$ bit number so $2^{52}\le q\le2^{53}-1$. $5^{23}=2^{53}\times1.32$ so it needs $54$ bits to represent it so $10^{23}=2^{23}*5^{23}$ is not exactly representable, and so $\sqrt{10^{46}}$ must be represented inexactly. It is surprising to me that your program didn't hang for that input and even for $\sqrt{10^{48}}$, but it did for higher powers of $10$. In my own tests shown in edits below it did hang for $\sqrt{10^{46}}$.

Typically you would compare the change in your $\text{guess}$ to $\text{guess}$ to see if the relative error is less than some tolerance between machine epsilon and the square root of machine epsilon. Bigger than machine epsilon because that's where infinite loops start and smaller than its square root because Newton's method is quadratically convergent so when your $\text{guess}$ is updated with that size of change it's about as good as it's gonna get. In the case of the Babylonian algorithm for the square root you know it will get there in $4$ iterations given the canonical initial value so there is no need to check a stopping condition.

Of course it's more usual to perform iterations for $1/\sqrt x$; see for example tihs discussion on an initial estimate. Whoa, my computer is crashing so I will have to stop now.

EDIT: I an going to try to explain my first ever python program today. The IEEE-754 double precision floating point format has bit $63$ as the sign, bits $52:62$ as the biased exponent and bits $0:51$ as the mantissa with an implicit leading $1$ bit. For example $\pi=\color{red}{0}\color{green}{400}\color{blue}{921\text{FB}54442\text{D}18}$. The sign bit is $\color{red}{0}$ so we know the sign of the number is $s=(-1)^{\color{red}{0}}=+1$. The mantissa is $\color{green}{400}$ and the bias is $3\text{FF}$ so the exponent is $e=\color{green}{400}-3\text{FF}=1$ so we know the mantissa will be multiplied by $(+1)2^e=2$ to get the final number. The visible part of the mantissa is $\color{blue}{921\text{FB}54442\text{D}18}$ so the actual mantissa is $1.\color{blue}{921\text{FB}54442\text{D}18}=1\color{blue}{921\text{FB}54442\text{D}18}/2^{52}=1.57079632679489655800$ and so the number is $2\times1.57079632679489655800=3.141592653589793116$.

With this background in IEEE-754, now consider what happens when the bits of the internal representation are shifted right by $1$: the exponent is sort of divided by $2$ but it's a little more complicated that that because it's a biased exponent. What happens to the mantissa is even more complicated than that. We propose to add another number to the internal representation before carrying out the right shift to handle the effect of the bias and also maybe to make the bits of the resulting mantissa a little better approximation to the square root. The approximation will be piecewise linear because changing the input by $2$ ulps changes the output by $1$ ulp but the relative scale of the changes will change when the input or output crosses a power of $2$. Thus we will analyze what happens at the kink points due to input being $1$, $2$, $4$, or $x_0$ where $f(x_0)=2^n$.

If $1<x_0<2$, then $f(x_0)=1$ or else it will be a very poor approximation to $\sqrt x_0$. The exponent will be $0$ so the biased exponent will be $3\text{FF}=1023$ and subtracting the hidden bit the mantissa will be $x_0-1$. The sign bit is always $0$ so we can even take the square root. Packing this into the internal representation we get TRANSFER(x0,0_INT64)=(1023+x_0-1)*2**52. You may be able to see that TRANSFER(1.0_REAL64,0_INT64)=(1023+1.0-1)*2**52=(1023)*2**52. We want $x_0$ to map to REAL(1,REAL64) so $$((1023+x_0-1)\times2^{52}+\text{magic})/2=(1023)\times2^{52}$$ Thus $\text{magic} = (1023+1-x_0)\times2^{52}$. Now we map the other kink points: $x=1$ maps to $$(1023+1023+1-x_0)\times2^{52}/2=\left(1022+\frac{3-x_0}2\right)\times2^{52}$$ Now we have to add the implicit $1$ to the mantissa and then multiply by $2^e=2^{1022-1023}=2^{-1}$ to get $f(1)=\frac{5-x_0}4$. Then mapping $2$ we get $$(1024+1023+1-x_0)\times2^{52}/2=\left(1023+\frac{2-x_0}2\right)\times2^{52}$$ Adding the implicit $1$ to the mantissa and multiplying by $2^e=1$ we get $f(2)=\frac{4-x_0}2$. We don't have to work hard to get $f(4)=2f(1)=\frac{5-x_0}2$ because multiplying the input by $4$ multiplies the output by $2$. That's why we are doing this because it approximates the square roots to the same relative accuracy at any scale. We can use the $2$-point formula for a line to connect the dots between the kink points to get, if $1<x_0<2$, $$f(x)=\begin{cases}\frac14x+\frac{4-x_0}4&1\le x\le x_0\\ \frac12x+\frac{2-x_0}2&x_0\le x\le2\\ \frac14x+\frac{3-x_0}2&2\le x\le4\end{cases}$$ It is left as an exercise for the reader to show that if $2<x_0<4$, $$f(x)=\begin{cases}\frac12x+\frac{6-x_0}4&1\le x\le2\\ \frac14x+\frac{8-x_0}4&2\le x\le x_0\\ \frac12x+\frac{4-x_0}2&x_0\le x\le4\end{cases}$$ Here, let's plot these curves:
fig. 1
It can be seen that the curve for $2<x_0<4$ is a worse approximation at $x=2$ and can't be fixed so we choose $1<x_0<2$. If $r=\sqrt x$ then $$y_{n+1}=r+e_{n+1}=\frac{x+y_n^2}{2y_n}=\frac{r^2+(r+e_n)^2}{2(r+e_n)}=r+\frac{e_n^2}{2(r+e_n)}$$ Allows us to find the absolute error $e_{n+1}$ after iteration $n+1$ given the error $e_n$ before iteration $n+1$. But more interesting is the relative error $$\epsilon_{n+1}=\frac{e_{n+1}}r=\frac{\epsilon_n^2}{2(1+\epsilon_n)}$$ It can be seen that relative error is much more convenient for analysis of this problem because it propagates independently of $x$, hence of the scale of the problem. Since negative errors from the above formula are worse than positive errors we find the relative error after the first iteration, which will always be positive: $$\epsilon_1=\frac{\left(\frac{y_0}{\sqrt x}-1\right)^2}{\left(1+\frac{y_0}{\sqrt x}-1\right)}=\frac{y_0^2-2y_0\sqrt x+x}{2y_0\sqrt x}$$ This is also plotted: fig. 2
As can be seen from fig. 2, the worst $2$ points are at $x=x_0$ and $x=2$. Looking at fig. 1 we can see that moving the kink point left would improve the error at $x=x_0$ but make it worse at $x=2$ while moving it right would improve the error at $x=2$ but make it worse at $x=x_0$. Thus the minimum error happens when the relative error after the first iteration is the same for both inputs: $$\frac{\left(\frac{4-x_0-2\sqrt2}{2\sqrt2}\right)^2}{2\left(\frac{4-x_0}{2\sqrt2}\right)}=\frac{\left(\frac{1-\sqrt{x_0}}{\sqrt{x_0}}\right)^2}{2\left(\frac1{\sqrt{x_0}}\right)}$$ On simplification we have $$\begin{align}g(x_0)&=x_0^{5/2}+2\sqrt2\,x_0^2-8x_0^{3/2}-6\sqrt2\,x_0+24x_0^{1/2}-8\sqrt2\\ &=\left(x_0^2+3\sqrt2\,x_0^{3/2}-2x_0-8\sqrt2\,x_0^{1/2}+8\right)\left(x_0^{1/2}-\sqrt2\right)=0\end{align}$$ Now we should make the substitution $x_0=y^2/2$, work out the resolvent cubic, factor into quadratics and get a hard-won result but we just used this method to spit out $x_0=8-4\sqrt3$. Then we could get $$\text{magic}=\left(1023+1-8+4\sqrt3\right)\times2^{52}=3\text{FEED}9\text{EBA}16132\text{A}9$$ Whoa, I actually had this last bit of arithmetic wrong before this edit. I had added $x_0-1$ rather than subtracting. Now fixed. It's getting late, so I don't have time to fix other typos and stuff right now. Here is the Matlab program that made the above graphs:

% psqrt.m

clear all;
close all;
x0 = 8-4*sqrt(3);
x1 = 3.99;
npts = 200;
% Consruct error curve for 1 < x0 < 2
x = linspace(1,x0,npts); % x in [0,x0]
y = x/4+(4-x0)/4;        % Initial approximation
xa = x;
ya = y;
x = linspace(x0,2,npts); % x in [x0,2]
y = x/2+(2-x0)/2;        % Initial approximation
xa = [xa(1:end-1) x];
ya = [ya(1:end-1) y];
x = linspace(2,4,npts);  % x in [2,4]
y = x/4+(3-x0)/2;        % Initial approximation
xa = [xa(1:end-1) x];
ya = [ya(1:end-1) y];
e0a = ya./sqrt(xa)-1;     % Relative error in initial approximation
e1a = e0a.^2./(2*(1+e0a));  % Relative error after first iteration
% Consruct error curve for 1 < x0 < 2
x = linspace(1,2,npts);  % x in [1,2]
y = x/2+(6-x1)/4;        % Initial approximation
xb = x;
yb = y;
x = linspace(2,x1,npts); % x in [2,x0]
y = x/4+(8-x1)/4;        % Initial approximation
xb = [xb(1:end-1) x];
yb = [yb(1:end-1) y];
x = linspace(x1,4,npts); % x in [x0,4]
y = x/2+(4-x1)/2;        % Initial approximation
xb = [xb(1:end-1) x];
yb = [yb(1:end-1) y];
e0b = yb./sqrt(xb)-1;     % Relative error in initial approximation
e1b = e0b.^2./(2*(1+e0b));  % Relative error after first iteration
x = linspace(1,4,300);
y = sqrt(x);
title('Approximations to $\sqrt{x}$','Interpreter','latex');
legend(['x_0 = ' num2str(x0)],['x_0 = ' num2str(x1)],'Exact', ...
title('Relative error after first iteration');
legend(['x_0 = ' num2str(x0)],['x_0 = ' num2str(x1)]);

EDIT: To go back to paragraph $2$, for odd powers of $10$, $\sqrt x$ has no exact representation as a rational number of which IEEE-754 numbers are a subset. I have written a Fortran program that, for certain odd powers of $10$, write out the number $x$ and then its mantissa in hex. Then $y=\sqrt x$ is computed and its mantissa extracted and printed. Then the mantissa of $y$ is squared and printed, then rounded to its final form and printed in hex.

! ieee754.f90
module funmod
   use ISO_FORTRAN_ENV, only: dp=> REAL64, qp=>REAL128, ik8=>INT64
   use ISO_FORTRAN_ENV, only: integer_kinds
   implicit none
   private integer_kinds
   integer, parameter :: ik16 = integer_kinds(size(integer_kinds))
      function extract_mantissa(x)
         integer(ik8) extract_mantissA
         real(dp), value:: x
         extract_mantissa = ior(iand(transfer(x,0_ik8),maskr(52,ik8)),shiftl(1_ik8,52))
      end function extract_mantissa

      function add_round(x)
         integer(ik16) add_round
         integer(ik16), value :: x
         integer lz
         lz = leadz(x)
         if(lz < 75) then
            add_round = x+maskr(74-lz,ik16)+ibits(x,75-lz,1)
            add_round = x
         end if
      end function add_round

      function round(x)
         integer(ik8) round
         integer(ik16), value :: x
         round = shiftr(x,75-leadz(x))
      end function round
end module funmod

program main
   use funmod
   implicit none
   real(dp) x
   real(dp) y
   real(qp) qx
   integer(ik8) i
   integer(ik16) qi
   integer j
   character(20) number
   do j = 7, 27, 2
      write(number,'(*(g0))') '1.0e',j
      read(number,'(f20.20)') x
      write(*,'("x = ",1p,e22.15)') x
      write(*,'("mantissa = ",Z14.14)') extract_mantissa(x)
      qx = x
      y = sqrt(qx)
      i = extract_mantissa(y)
      write(*,'("sqrt =",Z14.14)') i
      qi = i
      qi = qi**2
      write(*,'("squared = ",Z0)') qi
      i = round(add_round(qi))
      write(*,'("rounded = ",Z14.14)') i
      write(*,'(a)') repeat('*',40)
   end do
end program main

Output with gfortran 8.1.0:

x =  1.000000000000000E+07
mantissa = 1312D000000000
sqrt =18B48E29793D2F
squared = 2625A00000000120727F30D6EA1
rounded = 1312D000000001
x =  1.000000000000000E+09
mantissa = 1DCD6500000000
sqrt =1EE1B1B3D78C7A
squared = 3B9AC9FFFFFFFEDD8A3DE01AA24
rounded = 1DCD64FFFFFFFF
x =  1.000000000000000E+11
mantissa = 174876E8000000
sqrt =134D0F1066B7CC
squared = 174876E7FFFFFEF42187A854A90
rounded = 174876E7FFFFFF
x =  1.000000000000000E+13
mantissa = 12309CE5400000
sqrt =182052D48065C0
squared = 246139CA800001617EBE8711000
rounded = 12309CE5400001
x =  1.000000000000000E+15
mantissa = 1C6BF526340000
sqrt =1E286789A07F2F
squared = 38D7EA4C67FFFE6349187EFAAA1
rounded = 1C6BF52633FFFF
x =  1.000000000000000E+17
mantissa = 16345785D8A000
sqrt =12D940B6044F7E
squared = 16345785D8A000D7C19BC9F0204
rounded = 16345785D8A001
x =  1.000000000000000E+19
mantissa = 1158E460913D00
sqrt =178F90E385635D
squared = 22B1C8C12279FFD825753330FC9
rounded = 1158E460913D00
x =  1.000000000000000E+21
mantissa = 1B1AE4D6E2EF50
sqrt =1D73751C66BC34
squared = 3635C9ADC5DE9ED61EDE3CC6A90
rounded = 1B1AE4D6E2EF4F
x =  9.999999999999999E+22
mantissa = 152D02C7E14AF6
sqrt =12682931C035A0
squared = 152D02C7E14AF4E5217BB3BA400
rounded = 152D02C7E14AF5
x =  1.000000000000000E+25
mantissa = 108B2A2C280291
sqrt =1702337E304309
squared = 21165458500521864AC10EDB651
rounded = 108B2A2C280291
x =  1.000000000000000E+27
mantissa = 19D971E4FE8402
sqrt =1CC2C05DBC53CB
squared = 33B2E3C9FD08037BBECAB9542F9
rounded = 19D971E4FE8402

As can be seen, often the best square root you can get is off by $1$ ulp after squaring. But sometimes it isn't. It seems to me that your program should terminate for $x\in\{1.0e19,1.0e25,1.0e27\}$ at least. Can you double check for these inputs?

It's getting late now so not much time for the third paragraph. However there is a proof that $|x|=\sqrt{x^2}$ in IEEE-754 arithmetic provided $x^2$ isn't infinite or denormal.

EDIT: I tried a Fortran program that test some power of $10$ inputs:

program test1
   use ISO_FORTRAN_ENV, only: dp=>REAL64, INT64
   integer i, j
   character(20) number
   real(dp) x, guess
   integer(INT64), parameter :: magic = int(Z'3FEED9EBA16132A9',INT64)
   real(dp), parameter :: tolerance = 1.0e-3_dp
   do i = -12, 100
      write(number,'(*(g0))') '1.0e',i
      read(number,'(f20.20)') x
      guess = transfer(shiftr(magic+transfer(x,magic),1),guess)
      do j = 1, 10
         if(abs(guess**2-x) < tolerance) exit
         guess = (guess+x/guess)/2
      end do
      if(j <= 10) then
         write(*,'(a,ES10.2,a,i0,a)') 'x = ',x,' terminated after ',j-1,' iterations'
         write(10,'(i0,",")',advance='no') i
         write(*,'(a,ES10.2,a)') 'x = ',x,' did not terminate'
         write(20,'(i0,",")',advance='no') i
      end if
   end do
end program test1

The loop terminated for $\log x\in\{-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,14,16,18,19,20,22,24,25,26,27,28,30,32,33,34,36,38,39,40,41,42,43,44,45,47,48,54,55,57,63,66,74,77,79,80,82,83,84,85,86,87,89,91,92,94,95,96,99\}$ and did not for $\log x\in\{13,15,17,21,23,29,31,35,37,46,49,50,51,52,53,56,58,59,60,61,62,64,65,67,68,69,70,71,72,73,75,76,78,81,88,90,93,97,98,100\}$ so it confirms my suspicion in paragraph $3$ that it could hang for $\log x=46$. Also it did indeed terminate for $\log x\in\{19,25,27\}$ as my previous program would have led one to assume. As can be seen, when the input is not guaranteed to hang the program it's a bit of a toss-up whether squaring and rounding will exactly reproduce the input thus terminating the loop or not.

BTW, I tried $10^9$ random inputs and in no case was the oscillation alluded to in paragraph $4$ observed. Maybe it's provable that it can't happen in the case of this algorithm, but in general it's considered risky to to check floating point numbers for equality. Normally one checks whether their difference is within some tolerance. Like you could have set $\text{tolerance}=1.0e-12$ and then compared $\lvert x-\text{guess}^2\rvert<x*\text{tolerance}$ and given quadratic convergence would have at most one more iteration to go when the condition was met and would have terminated easily because this is thousands of ulps away from an exact match.

EDIT: Also some code showing a possible assembly language sequence for computing reciprocal square root and then square root avoiding divisions.
Assemble this via fasm sqrt1.asm

; sqrt1.asm
format MS64 COFF
section '.text' code readable executable align 16
public sqrt1
sqrt1:                           ; xmm0 = D
movq xmm1, [magic]               ; xmm1 = magic
movq xmm2, [three_halves]        ; xmm2 = 1.5
psubq xmm1, xmm0                 ; xmm1 = magic-D
psrlq xmm1, 1                    ; xmm1 = x (close to 1/sqrt(D))
vmulsd xmm3, xmm1, xmm2          ; xmm3 = 1.5*x
vmulsd xmm4, xmm1, xmm0          ; xmm4 = D*x
vpsubq xmm5, xmm1, [half_bit]    ; xmm5 = 0.5*x
mulsd xmm1, xmm5                 ; xmm1 = 0.5*x**2
vfnmadd213sd xmm1, xmm4, xmm3    ; xmm1 = 1.5*x-0.5*D*x**3 = x
vmulsd xmm3, xmm1, xmm2
vmulsd xmm4, xmm1, xmm0
vpsubq xmm5, xmm1, [half_bit]
mulsd xmm1, xmm5
vfnmadd213sd xmm1, xmm4, xmm3
vmulsd xmm3, xmm1, xmm2
vmulsd xmm4, xmm1, xmm0
vpsubq xmm5, xmm1, [half_bit]
mulsd xmm1, xmm5
vfnmadd213sd xmm1, xmm4, xmm3
vmulsd xmm3, xmm1, xmm2
vmulsd xmm4, xmm1, xmm0
vpsubq xmm5, xmm1, [half_bit]
mulsd xmm1, xmm5
vfnmadd213sd xmm1, xmm4, xmm3    ; xmm1 = 1.5*x-0.5*D*x**3 = x
mulsd xmm0, xmm1                 ; xmm0 = D*x = sqrt(D)
section '.data' data readable align 16
magic: dq 0xBFCDD6A18F6A6F53
three_halves: dq 0x3FF8000000000000
half_bit: dq 0x0010000000000000

And then compile this via gfortran sqrt2.f90 sqrt1.obj -osqrt2

! sqrt2.f90
program sqrt2
   use ISO_FORTRAN_ENV, only: dp=>REAL64
   implicit none
      function sqrt1(x) bind(C,name='sqrt1')
         implicit none
         real(dp) sqrt1
         real(dp), value :: x
      end function sqrt1
   end interface
   real(dp) x
   x = 2
   write(*,*) sqrt1(x), sqrt(x), sqrt1(x)-sqrt(x)
end program sqrt2

Output was:

   1.4142135623730951        1.4142135623730951        0.0000000000000000
  • $\begingroup$ Thanks for this sharp answer. About the 1st paragraph, are you computing the initial guess by extracting the exponent of the floating-point representation of $x$ and dividing it by 2? About the 2nd paragraph, I guess that you give two spacings between consecutive floating-point numbers ($2^{-53+e}$ and $2^{-52+e}$) because it depends if they have the same exponent $e$ (in which case it is the latter) or different exponents $e$ and $e + 1$ respectively (in which case it is the former), right? About the 3rd paragraph, I did not understand. $\endgroup$
    – Maggyero
    Jan 30 '20 at 13:02
  • $\begingroup$ About the 4th paragraph, when comparing the new guess with the old guess, why using a tolerance (like you suggested) instead of 0 (like I suggested: def good_enough(guess): return improve(guess) == guess)? Using 0 seems to always succeed (which is expected since at some point you will not be able to improve you guess anymore because of the machine epsilon so the equation will be verified), with the best approximation and without any need for a tolerance. $\endgroup$
    – Maggyero
    Jan 30 '20 at 13:17
  • $\begingroup$ In the first paragraph the only actual computation is (magic+temp)>>1, the pack/unpack stuff is actually a noop that changes how we are thinking about the same data. In the second paragraph, number with the same exponent of $2^{52}$ can range from $2^{52}$ up to almost $2^{53}$; it's the range of the mantissa I'm taking into account. When I have more time I'll fix typos in the third paragraph. In the fourth paragraph, Newton's method is notorious for getting caught in loops and your newer test is not safe. You could oscillate between two values differing by one ulp. $\endgroup$ Jan 30 '20 at 16:23
  • $\begingroup$ In the 1st paragraph, coud you elaborate on that calculation (magic and temp)? In the 2nd paragraph, I still don’t get where your $2e^{-53}$ value comes from. To me the spacing between two consecutive floating-point numbers is always $\beta^{-(p - 1)}\beta^e = 2^{e - 52}$ where $e$ is the exponent of the smaller number of the two. I got this value by computing the difference between two consecutive floating-point numbers with the same exponent $e$: $d_0.d_1d_2 \cdots d_{p - 2}(d_{p - 1} + 1) \times \beta^e - d_0.d_1d_2 \cdots d_{p - 2}d_{p - 1} \times \beta^e = \beta^{-(p - 1)}\beta^e$. $\endgroup$
    – Maggyero
    Jan 31 '20 at 8:55
  • $\begingroup$ And by computing the difference between two consecutive floating-point numbers with different exponents ($e$ for the smaller, $e + 1$ for the larger), which also yields the same result: $1.00 \cdots 0 \times \beta^{e + 1} - (\beta - 1).(\beta - 1)(\beta - 1) \cdots (\beta - 1) \times \beta^e = \beta^{-(p - 1)}\beta^e$. About the 3rd paragraph, okay thanks. About the 4th paragraph, that is very interesting, I did not know that the guess could oscillate that way. Do you have a value for $x$ which triggers that behaviour? $\endgroup$
    – Maggyero
    Jan 31 '20 at 9:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.