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:

1. Set $$n$$, such that $$2^n>C\ge2^{(n-1)}$$.
2. $$m=\lceil n/2\rceil$$.
3. $$i=0$$.
4. Select $$X_0$$, such that $$2^m>X_0\ge2^{(m-1)}$$.
5. Repeat
5.1   $$i=i+1$$.
5.2   $$X_i=((X_{i-1})^2+C)/(2X_{i-1})$$.
Until $$(X_i)^2<2^m+C$$.
6. If $$C=\lfloor X_i\rfloor^2$$, then status = PERFECT SQUARE.
Else status = NOT A PERFECT SQUARE.
7. 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$$.

• Conclusion: it directly follows from the accepted answer that yes we can and should use Euclidean division to compute $X_i$ at step 5.2, and use $X_i$ where there is $\lfloor X_i\rfloor$ at step 6. Commented Jan 16, 2023 at 7:36

1 Answer

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}<\frac{2x_k^2}{2x_k} = x_k$$

But since $$x_0 >0$$, let us distinguish:

• if $$x_0 < \sqrt{C}$$, then $$2x_0(g(x_0) - \sqrt{C}) = x_0^2 + C -2\sqrt{C}x_0 = (x_0 -\sqrt{C})^2 > 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.

• I like the answer's approach of viewing the algorithm as finding the zero of $f$ by Newton's method (that is by iterating $g$), except for rounding down to the nearest integer at each step. Commented Jan 15, 2023 at 13:22