Calculate fractional part of square root without taking square root

Let's say I have a number $x>0$ and I need to calculate the fractional part of its square root:

$$f(x) = \sqrt x-\lfloor\sqrt x\rfloor$$

If I have $\lfloor\sqrt x\rfloor$ available, is there a way I can achieve (or approximate) this without having to calculate any other square roots? I would like to avoid square roots for performance reasons as they are usually a rather slow operation to calculate.

As an example of what I am thinking about, here is one option:

$$f(x)\approx\frac{x-\lfloor\sqrt x\rfloor^2}{(\lfloor\sqrt x\rfloor+1)^2-\lfloor\sqrt x\rfloor^2}$$

As $x$ gets large, this approximation becomes better and better (less than 1% error if $\sqrt x\ge40$), but for small $x$, it's quite terrible (unbounded error as $x$ goes to $0$).

Can I do better? Preferably an error of less than 1% for any $\sqrt x<10000$.

The holy grail would be to find a formula that doesn't need any division either.

Aside: In case anyone is interested why I need this: I'm trying to see if I can speed up Xiaolin Wu's anti-aliased circle drawing functions by calculating the needed variables incrementally (Bresenham-style)

• – MJD
Commented May 26, 2015 at 19:55

If you are using IEEE 754 compliant floating point numbers, it may be that square root operations are faster than you might suppose, as the square root operation is directly supported in the standard and any compliant hardware is guaranteed to round correctly (unlike for sine and cosine) and is often using a hardware implementation of the Newton-Raphson method Alijah described.

That said, the algorithm you linked to only uses the fractional part of the square root to calculate pixel opacity, and consequently the final opacity value ranges only from 0 to 255. Because of the small range, floating point numbers may be overkill and a fixed-point integer representation might work better. If the range is truly only a byte and the maximum size of your radii aren't too large, you can use a look-up table with fixed-point input to skip the expensive square root calculation. A 16-bit fixed-point number would give you a 64KB look-up table, which isn't too bad.

You might also be able to avoid a division and square root operation for calculating the $45^\circ$ value (ffd in the algorithm) by using the Fast Inverse Square Root hack.

Now for the question of whether there is a method to calculate the fractional part of a square root knowing only the integer part, there is one approach that iteratively calculates a square root and minimizes divisions: Continued Fractions. The simplest approach I know of for what you are wanting to do (fast convergence) is detailed here and works as follows:

$a_0 = x - \lfloor\sqrt x\rfloor^2$

$b_0 = 2\lfloor\sqrt x\rfloor$

$a_n = a_0b_{n-1}$

$b_n = b_0b_{n-1} + a_{n-1}$

which gives you quadratically better and better approximations of the fractional part of $\sqrt x$, and you divide $\frac{a_n}{b_n}$ when you've done enough iterations to ensure accuracy. If you are ultimately needing only byte sized opacity values, it should only take 3 or 4 iterations or so, and we save the division until the end, which is the only significant difference between this and the Newton-Raphson method other than the fact that it gives you the fractional part directly.

If you really want to pursue incrementally calculating variables as far as it can go, you can use Gosper's continued fraction algorithms (see especially the section on square roots) and calculate all the variables involved as continued fractions one term at a time. This allows you to avoid square root calculations and divisions other than bit shifts, as well as abort as soon as you know the correct pixel opacity (or whatever else you're calculating) without wasting time calculating digits you don't need, but it involves a serious overhaul to the algorithm you linked, and if I went into detail this answer would turn into a book.

So essentially, if you have memory to spare and the max length of your radii isn't huge and your output size is small, go with a look-up table with fixed-point numbers. If you want simplicity of implementation, go with floating point or fixed-point numbers. If you want to absolutely avoid square root calculation without a look-up table go with Newton-Raphson or the continued fraction variant on it. If you want to absolutely minimize wasted computation at the expense of some up-front overhead, go with Gosper's continued fraction algorithms.

• That's an awesome amount of detail! It'll probably take me a little while to dig through your different suggestions to find the best one. The Gosper paper definitely looks very interesting and worth understanding. Thank you! +1 (I do only need 8-bit precision so I could figure out how many iterations I would need in the worst case and if it's a sufficiently small number, I can hard-code them. With only a single division at the end, it might actually be faster than using the general sqrt function of some mobile processor...) Commented May 13, 2015 at 20:34
• The continued fraction variant on Newton-Raphson comes closest to what you asked for in your question (the fractonal part directly), but I mentioned the other possibilities when I saw the nature of the (awesome) anti-aliased circle drawing functions. Check out the web page I linked before the description - "Ted's Math World - Pencil and Paper Square Roots" Commented May 13, 2015 at 20:40
• Another paper by Gosper Commented May 15, 2015 at 15:25
• This continued fractions stuff is amazing! I especially like the simple explanation from Ted's Math World that you linked. I did some tests and if I only do TWO steps (i.e. to 40/69 in Ted's last example), I get it right on in pretty much every case. The worst error is for sqrt(2) where my fractional part ends up being 102 instead of 106 (base 256), the second worst is sqrt(6) with 113 instead of 115 and everything else is within one bit! Awesome! And with only two steps, it should be fairly trivial to keep the numerator and denominator up to date incrementally as well. That'll be next... :) Commented May 16, 2015 at 1:47
• But even without that, this works out to only 3 multiplications, 3 additions, 1 shift, and 1 division! No branches (if-statements) either! Totally awesome! I wonder if this can be used to draw not just circles and ellipses, but arbitrary quadratic and/or cubic bezier curves. That would be insane!!! Commented May 16, 2015 at 1:56

You could try using the Newton-Raphson method, where you calculate the square root of $x$ in its entirety, using the following update (for iteration $n+1$):- $$y(n+1)=y(n)-\frac{y(n)^2-x}{2y(n)}=0.5\left(y(n)+\frac{x}{y(n)}\right)$$ and setting $y(1)=\lfloor\sqrt x\rfloor$ as a starting value. Thus, for iteration $n+1$, the fractional part $f(x)$ would be $y(n+1)-\lfloor\sqrt x\rfloor$.

The advantages of this approach is that the algorithm converges quite rapidly (say in about $4$ to $5$ iterations), and there is no need to calculate any other square roots. Unfortunately, you need to carry out a division per iteration.

• That could work. But I don't know if it would actually be faster than to just have the processor evaluate the square root itself (which is what I was trying to avoid). I just realized, though, that I forgot to specify why I wanted to avoid the square root. So, in computational environments where square roots simply aren't available, your solution can really help out. +1 Commented May 11, 2015 at 20:52

After a bit more reading on Ted's Math World (first pointed out by @hatch22 in his answer), I came across a quite beautiful way to perform this calculation using the Pythagorean theorem:

Given some estimate $y$ of the square root (in the case here we have$\lfloor\sqrt{x}\rfloor$), if we let the hypotenuse of a right triangle be $x+y^2$ and one of its other two sides $x-y^2$, then the remaining side will have length $2y\sqrt{x}$, i.e. it will contain the desired answer.

The angle $\alpha$ formed by the side of interest and the hypotenuse can be calculated as:

$$\alpha=\sin^{-1}\frac{x-y^2}{x+y^2}\approx\frac{x-y^2}{x+y^2}$$

where the approximation is the first term of the Maclaurin Series for $\sin^{-1}$.

The side of interest can then be calculated from:

$$2y\sqrt{x}=(x+y^2)\cos\alpha\approx(x+y^2)\cos\frac{x-y^2}{x+y^2}\approx(x+y^2)\left(1-\frac{1}{2}\left(\frac{x-y^2}{x+y^2}\right)^2\right)$$

Where the second approximation are the first two terms of the Maclaurin Series for $\cos$.

From this, we can now get:

$$\sqrt{x}\approx\frac{x+y^2}{2y}\left(1-\frac{1}{2}\left(\frac{x-y^2}{x+y^2}\right)^2\right)=\frac{x^2+6xy^2+y^4}{4y(x+y^2)}$$

To get the fractional part of $\sqrt{x}$ in the range $0..255$, this can be optimized to:

$$y_{\,\text{square}}=y\times y$$ $$s=x+y_{\,\text{square}}$$ $$r=\frac{(s\times s\ll6) + (x\times y_{\,\text{square}}\ll8)}{s\times y}\,\,\&\,\,255$$

where $\ll$ signifies a bit-wise shift to the left (i.e. $\ll6$ and $\ll8$ are equivalent to $\times\,64$ and $\times\,256$ respectively) and $\&$ signifies a bit-wise and (i.e. $\&\,255$ is equivalent to $\%\,256$ where $\%$ stands for the modulus operator).

The amazing part is that despite the minimal Maclaurin Series used, if we can use the closer of $\lfloor\sqrt{x}\rfloor$ and $\lceil\sqrt{x}\rceil$ as the estimate $y$ (I have both available), the answer in the range $0..255$ is actually EXACT!!! for all values of $x\ge1$ that don't lead to an error due to overflow during the calculation (i.e. $x<134\,223\,232$ for 64-bit signed integers and $x<2\,071$ for 32-bit signed integers).

It is possible to expand the usable range of the approximation to $x<2\,147\,441\,941$ for 64-bit signed integers and $x<41\,324$ for 32-bit signed integers by changing the formula to:

$$r=\left(\frac{s\ll6}{y} + \frac{x\times y\ll8}{s}\right)\,\&\,\,255$$

But due to the earlier rounding, this leads to a reduction in the accuracy such that the value is off by $1$ in many cases.

Now the problem: A little bit of benchmarking and reading indicates that on many processors a division operation is actually not much faster than a square root. So unless I can find a way to get rid of the division as well, this approach isn't actually going to help me much. :(

Update:

If an accuracy of $\pm 1$ is acceptable, the range can be increased significantly with this calculation: $$k = \frac{(x + y \times y) \ll 5}{y}$$ $$r = \left(\left(k + \frac{x \ll 12}{k}\right)\ll 1\right)\,\&\,\,255$$

For 32-bit signed integers, this works for any $x<524\,288$, i.e. it breaks down as soon as $x \ll 12$ overflows. So, it can be used for circles up to radius 723 pixels. Note, that $y$ does not change on every step of the Bresenham algorithm, so $1/y$ can be pre-calculated and therefore does not add a full second division to the algorithm.