Editing Integer square root (section)

==Algorithm using Newton's method==
One way of calculating <math>\sqrt{n}</math> and <math>\operatorname{isqrt}(n)</math> is to use [[Methods_of_computing_square_roots#Heron's method|Heron's method]], which is a special case of [[Newton's method]], to find a solution for the equation <math>x^2 - n = 0</math>, giving the iterative formula
<math display="block">x_{k+1} = \frac{1}{2}\!\left(x_k + \frac{n}{x_k}\right), \quad k \ge 0, \quad x_0 > 0.</math>

The [[sequence]] <math>\{x_k\}</math> [[limit of a sequence|converges]] [[Rate of convergence|quadratically]] to <math>\sqrt{n}</math> as <math>k\to\infty</math>. 

===Stopping criterion===
One can prove{{Citation needed|date=August 2021}} that <math>c=1</math> is the largest possible number for which the stopping criterion 
<math display="block">|x_{k+1} - x_{k}| < c</math>
ensures <math>\lfloor x_{k+1} \rfloor=\lfloor \sqrt n \rfloor</math> in the algorithm above.

In implementations which use number formats that cannot represent all [[rational number]]s exactly (for example, floating point), a stopping constant less than 1 should be used to protect against round-off errors.

===Domain of computation===
Although <math>\sqrt{n}</math> is [[irrational number|irrational]] for many <math>n</math>, the sequence <math>\{x_k\}</math> contains only rational terms when <math>x_0</math> is rational. Thus, with this method it is unnecessary to exit the [[field (mathematics)|field]] of rational numbers in order to calculate <math>\operatorname{isqrt}(n)</math>, a fact which has some theoretical advantages.

===Using only integer division===
For computing <math>\lfloor \sqrt n \rfloor</math> for very large integers ''n'', one can use the quotient of [[Euclidean division]] for both of the division operations.  This has the advantage of only using integers for each intermediate value, thus making the use of [[floating point]] representations of large numbers unnecessary.  It is equivalent to using the iterative formula
<math display="block">x_{k+1} = \left\lfloor \frac{1}{2}\!\left(x_k + \left\lfloor \frac{n}{x_k} \right\rfloor \right) \right\rfloor, \quad k \ge 0, \quad x_0 > 0, \quad x_0 \in \mathbb{Z}.</math>

By using the fact that
<math display="block">\left\lfloor \frac{1}{2}\!\left(x_k + \left\lfloor \frac{n}{x_k} \right\rfloor \right) \right\rfloor = \left\lfloor \frac{1}{2}\!\left(x_k + \frac{n}{x_k} \right) \right\rfloor,</math>

one can show that this will reach <math>\lfloor \sqrt n \rfloor</math> within a finite number of iterations.

In the original version, one has <math>x_k \ge \sqrt n</math> for <math>k \ge 1</math>, and <math>x_k > x_{k+1}</math> for <math>x_k > \sqrt n</math>. So in the integer version, one has <math>\lfloor x_k \rfloor \ge \lfloor\sqrt n\rfloor</math> and <math>x_k \ge \lfloor x_k \rfloor > x_{k+1} \ge \lfloor x_{k+1}\rfloor</math> until the final solution <math>x_s</math> is reached. For the final solution <math>x_s</math>, one has <math>\lfloor \sqrt n\rfloor\le\lfloor x_s\rfloor \le \sqrt n</math> and <math>\lfloor x_{s+1} \rfloor \ge \lfloor x_s \rfloor</math>, so the stopping criterion is <math>\lfloor x_{k+1} \rfloor \ge \lfloor x_k \rfloor</math>.

However, <math>\lfloor \sqrt n \rfloor</math> is not necessarily a [[Fixed point (mathematics)|fixed point]] of the above iterative formula. Indeed, it can be shown that <math>\lfloor \sqrt n \rfloor</math> is a fixed point if and only if <math>n + 1</math> is not a perfect square. If <math>n + 1</math> is a perfect square, the sequence ends up in a period-two cycle between <math>\lfloor \sqrt n \rfloor</math> and <math>\lfloor \sqrt n \rfloor + 1</math> instead of converging.

=== Example implementation in C ===
<syntaxhighlight lang="c" line="1">
// Square root of integer
unsigned int int_sqrt(unsigned int s)
{
	// Zero yields zero
    // One yields one
	if (s <= 1) 
		return s;

    // Initial estimate (must be too high)
	unsigned int x0 = s / 2;

	// Update
	unsigned int x1 = (x0 + s / x0) / 2;

	while (x1 < x0)	// Bound check
	{
		x0 = x1;
		x1 = (x0 + s / x0) / 2;
	}		
	return x0;
}
</syntaxhighlight>

=== Numerical example ===

For example, if one computes the integer square root of {{math|2000000}} using the algorithm above, one obtains the sequence
<math display="block">\begin{align}
& 1000000 \rightarrow 500001 \rightarrow 250002 \rightarrow 125004 \rightarrow 62509 \rightarrow 31270 \rightarrow 15666 \rightarrow 7896 \\
& \rightarrow 4074 \rightarrow 2282 \rightarrow 1579 \rightarrow 1422 \rightarrow 1414 \rightarrow 1414
\end{align}</math>
In total 13 iteration steps are needed. Although Heron's method converges quadratically close to the solution, less than one bit precision per iteration is gained at the beginning. This means that the choice of the initial estimate is critical for the performance of the algorithm.

When a fast computation for the integer part of the [[binary logarithm]] or for the [[bit-length]] is available (like e.g. <code>std::bit_width</code> in [[C++20]]), one should better start at 
<math display="block">x_0 = 2^{\lfloor (\log_2 n) /2 \rfloor+1},</math>
which is the least [[power of two]] bigger than <math>\sqrt n</math>. In the example of the integer square root of {{math|2000000}}, <math>\lfloor \log_2 n \rfloor = 20</math>, <math>x_0 = 2^{11} = 2048</math>, and the resulting sequence is
<math display="block">2048 \rightarrow 1512 \rightarrow 1417 \rightarrow 1414 \rightarrow 1414.</math>
In this case only four iteration steps are needed.