Correctness of FIPS 186-4 square test algorithm

Federal Information Processing Standard 186-4 appendix C.4 gives (without reference) an algorithm intended to test if a positive integer $C$ (which can be thousands bits) is a square:

• Set $n$, such that $2^n>C\ge2^{(n-1)}$.

• $m=\lceil n/2\rceil$.

• $i=0$.

• Select $X_0$, such that $2^m>X_0\ge2^{(m-1)}$.

• Repeat

  5.1   $i=i+1$.

  5.2   $X_i=((X_{i-1})^2+C)/(2X_{i-1})$.

Until $(X_i)^2

• If $C=\lfloor X_i\rfloor^2$, then status = PERFECT SQUARE.

Else status = NOT A PERFECT SQUARE.

• Return status.

At step 5.2, can we compute $X_i$ by Euclidean division (that is rounding down to the nearest integer)? In the affirmative, how do we prove the algorithm terminates, and yields correct status, regardless of the choice made at step 4? Towards this, what invariant condition should we prove is maintained during step 5?

I observe that in step 6 it's used $\lfloor X_i\rfloor$ rather then just $X_i$, and that's strange if $X_i$ is intended to be an integer. Also, the notes after the statement of the algorithm, intended to prove it correct, seem written with real values in mind.

Remark: A practical implementation could suppress $i$ (thus steps 3 and 5.1), and use a single $X$ for all the $X_i$. It could use an arbitrary integer type or class for at least $C$ and $X$.

This is more or less Newton's method applied to the function $f:x\to x^2 - C$. The algorithm finds $r = \left\lfloor \sqrt{C}\right\rfloor$ and check whether $r^2 = C$ or not.

If we define $g(x) = x - \frac{f(x)}{f'(x)} = x - \frac{x^2-C}{2x} = \frac{x^2+C}{2x}$, then the computed sequence is $(x_k)_{k\in\mathbb{N}}$ such that $x_{k+1} = \left\lfloor g(x_k)\right\rfloor$.

Assume $x_k > \sqrt{C}$ for some $k$. Then:
$$x_{k+1}=\left\lfloor g(x_k)\right\rfloor\leqslant g(x_k) = \dfrac{x_k^2+C}{2x_k}0$, let us distinguish:

• if $x_0 0$. That means that $g(x_0)> \sqrt{C}$ and that $x_1 \geqslant \left\lfloor \sqrt{C}\right\rfloor$;

• otherwise $x_0 \geqslant \left\lfloor \sqrt{C}\right\rfloor$.

With the previous properties, we conclude that there exists $k>0$ such that $x_k \leqslant \left\lfloor \sqrt{C}\right\rfloor$ (otherwise $(x_k)_{k>0}$ is a strictly decreasing sequence with $\left\lfloor \sqrt{C}\right\rfloor$ as a lower bound).

Let $k$ be the first integer $>0$ such that $x_k \leqslant \left\lfloor \sqrt{C}\right\rfloor$. That means that $g(x_{k-1}) > \sqrt{C}$. Since $x_k = \left\lfloor g(x_{k-1})\right\rfloor \leqslant \left\lfloor \sqrt{C}\right\rfloor$, we conclude that $x_k = \left\lfloor \sqrt{C}\right\rfloor$.

Since (as stated in your document), $x_k^2 < 2^m +C$ implies $x_k-\sqrt{C}<1$, that means that $x_k \leqslant \left\lfloor \sqrt{C}\right\rfloor$ so we conclude on the termination and the correctness of the algorithm.

StackExchange Computer Science Q#156858, answer score: 10