Efficently calculating the "integer hypotenuse"

Question

Background: the integer square root function is defined as $\DeclareMathOperator{\isqrt}{isqrt}\isqrt{\left(x\right)}=\lfloor\sqrt{x}\rfloor$. There are methods to compute the integer square root using only integers (see previous link), which makes it useful for approximating the real square root on computers with limited processing power.

I'm interested in calculating what could be called the "integer hypotenuse" function: $\DeclareMathOperator{\ihypot}{ihypot}\ihypot{\left(x,y\right)}=\lfloor\sqrt{x^2 + y^2}\rfloor$.* This can be calculated using the same methods as the integer square root, but are there "better" approaches that somehow exploit the form of the argument (the sum of two squares) while still only requiring integers for computation?

I realize that is somewhat vague, so here's an example: One of the ways to find the integer square root is using binary search. Binary search requires an upper and lower bound. Since

$$ 0 \le \lfloor\sqrt{a}\rfloor \le a$$

$0$ and $a$ are often (always?) used as the bounds. We could naively do the same thing with the integer hypotenuse:

$$ 0 \le \lfloor\sqrt{a^2 + b^2}\rfloor \le a^2 + b^2$$

but we can find tighter bounds by realizing that the argument is the sum of squares and using the triangle inequality:

$$ \max\left(\left|a\right|,\left|b\right|\right) \le \lfloor\sqrt{a^2 + b^2}\rfloor \le \left|a\right| + \left|b\right|$$

Using those bounds would mean fewer iterations would be required to find the solution using binary search, so that's an example of a "better" approach that exploits the form of the argument.

I can't rigorously define "better", but I'm using it in the sense of "requiring less computational resources (including time and memory)" to get the answer. Things that might be better include:

• Approaches that require fewer iterations (like the example above).
• Not having to explicitly calculate $a^2 + b^2$. (Calculating $a^2 + b^2$ could lead to integer overflow issues on a computer.)
• Not having to divide by anything other than a power of 2. (Arbitrary division is expensive on a computer.)

Here's an example of what the second bullet might look like --- except that it doesn't work:

Heron's method for calculating the square root makes use of the inequality of arithmetic and geometric means as follows:

$$\sqrt{a} = \sqrt{x \frac{a}{x}} \leq \frac{1}{2}\left(x + \frac{a}{x}\right)$$

Then you can make an initial guess at the solution $x_0$, iterate $x_{n+1} = \frac{1}{2}\left(x_n + \frac{a}{x_n}\right)$, and $x_n$ will approach $\sqrt{a}$.

I was tempted to do the following:

$$\sqrt{a^2+b^2} = \sqrt{\left(a+bi\right)\left(a-bi\right)} \leq \frac{1}{2}\left(\left(a+bi\right)x + \frac{a-bi}{x}\right)$$

and then iterate in the same way as above. This looks promising and avoids ever having to calculate $a^2$ or $b^2$. However the inequality of arithmetic and geometric means only works for real numbers, so the above line is wrong. (There is a generalization to complex numbers --- eq 1.1 here, but it added enough complication that it didn't look promising to me.)

Also, just to be explicit: I'm not interested in approximations or methods that don't use integers for intermediary steps.

* I realize that this name might sound weird: The hypotenuse is the side of a triangle --- not a function! The motivation is that programming languages typically have separate built-in functions for calculating the square root and the length of the hypotenuse (e.g. math.sqrt and math.hypot in Python). Python even has an integer square root function: math.isqrt. I'm interested in an integer version of the hypotenuse function --- hence the name.

Answer 1

For $0<b<a$ you know from the Taylor expansion $$ \sqrt{a^2+b^2}=a\sqrt{1+\frac{b^2}{a^2}}=a+\frac{b^2}{2a}+... $$ This should give a better first approximation. One can of course use more terms of $\sqrt{1+x}=1+\frac{x}2-\frac{x^2}8+\frac{x^3}{16}-\frac{5x^4}{128}+...$ to get closer results.

Answer 2

In the same spirit as @Lutz Lehmann, using the $[n,n]$ Padé approximant $P_n$ of $\sqrt{1+x}$ (where $x=\left(\frac{b}{a}\right)^2)$ $$P_1=\frac{3 x+4}{x+4}$$ $$P_2=\frac{5 x^2+20 x+16}{x^2+12 x+16}$$ $$P_3=\frac{7 x^3+56 x^2+112 x+64}{x^3+24 x^2+80 x+64}$$ $$P_4=\frac{9 x^4+120 x^3+432 x^2+576 x+256}{x^4+40 x^3+240 x^2+448 x+256}$$

Computing the norm $$\Phi_n=\int_0^1\left(\sqrt{1+x} -P_n\right)^2\,dx$$ $$\Phi_1=3.5 \times 10^{-5} \quad \Phi_2=2.0 \times 10^{-8}\quad \Phi_3=1.3 \times 10^{-11}\quad \Phi_4=9.2 \times 10^{-15}$$

Answer 3

In general, you can compute square roots by applying Newton's method to the equation $$ f(x) = x^2 - \alpha.$$ Newton's method takes the form $$ x_{n+1} = x_n - \frac{x_n^2 -\alpha}{2x_n} $$ where $x_0 = x_0(\alpha)$ must be selected by the user. It is straightforward to verify that the relative error $$ r_n = \frac{\sqrt{\alpha} - x_n}{\sqrt{\alpha}} $$ satisfies the inequality $$ |r_{n}| \leq 2(|r_0|/2)^{2^n}.$$ Hence, Newton's method will converge at least quadratically, provided that the initial relative error is strictly less than $2$.

In the current situation, we have $\alpha = a^2 + b^2$ where $a$ and $b$ are integers and it is desirable to use integer calculations to compute $\lfloor \sqrt{\alpha} \rfloor$.

We shall first consider the selection of the initial guess $x_0$. We can without loss of generality assume that $0 \leq a \leq b$. It follows that $$b \leq \sqrt{\alpha} \leq \sqrt{2} b$$ or equivalently $$0 \leq \sqrt{\alpha} - b \leq (\sqrt{2} - 1) b.$$ It follows that if we choose $x_0 = b$, then the initial relative error satisfies $$ r_0 = \frac{\sqrt{\alpha} - b}{\sqrt{\alpha} } \leq (\sqrt{2} - 1) \frac{b}{\sqrt{\alpha}} \leq \sqrt{2} - 1 < 2.$$ Hence quadratic convergence is assured from the start.

Next we consider how this instance of Newton's iteration can be approximated using integer calculations only. There are at two points to consider, the potential for nonconvergence and integer overflow. The wiki page on the integer square root claims that convergence is assured and I have no reason to doubt that this is true as Newton's method is quite forgiving when it comes to rounding errors. It remains to describe how to calculate $$ \left \lfloor \frac{x^2 - a^2 - b^2}{2x}\right \rfloor $$ In reality this is trivial. Using integer division with remainder, we have $$ x = 2p + r_x, \quad r_x \in \{0,1\}$$ and $$x^2 = 2xp + x r_x.$$ Similarly, we have $$ a = 2x p_a + r_a, \quad 0 \leq r_a < 2x$$ It follows that $$ a^2 = 4x^2 p_a^2 + 4xp_a r_a + r_a^2 = 2x (2x p_a^2 + 2 p_a r_a) + r_a^2.$$ Similarly, we have $$ b^2 = 2x ( 2x p_b^2 + 2 p_b r_b) + r_b^.2$$ It follows that $$ x^2 - a^2 - b^2 = 2x (p - 2x p_a^2 - 2 p_a r_a - 2xp_b^2 - 2 p_b r_b) + x r_x - r_a^2 - r_b^2$$ It remains to write $$ xr_x - r_a^2 - r_b^2 = 2x s + t, \quad 0 \leq t < 2x$$ so that $$ x^2 - a^2 - b^2 = 2x z + t.$$ where $$ z = p - 2x p_a^2 - 2 p_a r_a - 2xp_b^2 - 2 p_b r_b + s. $$ It remains to derive conditions on $a$ and $b$ than will ensure that this formula for $z$ cannot trigger an overflow, but that is a topic for another day.