Square roots modulo a prime

nmod_sqrt is the function to use for square roots modulo a word-size prime; _nmod_sqrt_small is the exhaustive search it falls back to for small odd moduli. Of the two square roots r and n - r of a nonzero square modulo a prime n, both functions return the one in [0, n/2], so the result does not depend on the algorithm that found it.

Composite and even moduli#

nmod_sqrt does not check that n is prime. For an arbitrary modulus it behaves as follows.

• For a = 0 and a = 1 it returns 1 and sets *res to a, for every n.
• For even n and a \ge 2 it returns 0, even when a is a square.
• For odd n < 128 it calls _nmod_sqrt_small, which is exact for every odd n: it returns the least root, or 0 if there is none.
• For odd composite n \ge 128, a return value of 1 comes with a genuine root (its square is a modulo n), but it need not be the least root, and a return value of 0 does not prove that a is a non-square.

To find all square roots modulo a composite number, factor it and use n_sqrtmodn from ulong_extras.

Algorithm#

For odd n < 128, the function searches through t = 1, \ldots, (n-1)/2, computing t^2 incrementally. Otherwise all arithmetic is done in Montgomery form, with half-word, full-word lazy or full-word arithmetic depending on the size of n, and the method depends on n:

Here s is the exponent of the largest power of 2 dividing n - 1. Müller's method is used when 6s is at least the bit length of n plus 120, which requires n \ge 2^{25}; it is typical of primes chosen for number-theoretic transforms, such as 469762049 = 7 \cdot 2^{26} + 1.

None of the methods tests whether a is a square first. Each produces a candidate root, which is accepted only if its square is a. For prime n this check is the residuosity test, and it is also what makes the result trustworthy for composite n.

Example#

The four primes below exercise the four main code paths: n \equiv 3 \pmod 4, n \equiv 5 \pmod 8, Tonelli–Shanks with s = 5 (where 10 is not a square), and Müller's method with s = 26.

#include <flint/flint.h>
#include <flint/nmod.h>

int main(void)
{
    ulong primes[4] = { 1000003, 1000037, 1000033, 469762049 };
    nmod_t mod;
    ulong r, a;
    int i;

    for (i = 0; i < 4; i++)
    {
        nmod_init(&mod, primes[i]);
        a = nmod_set_ui(10, mod);

        if (nmod_sqrt(&r, a, mod))
            flint_printf("n = %wu: sqrt(10) = %wu, r^2 = %wu\n",
                         mod.n, r, nmod_mul(r, r, mod));
        else
            flint_printf("n = %wu: 10 is not a square\n", mod.n);
    }

    /* composite n: the root found is genuine */
    nmod_init(&mod, 135);
    if (nmod_sqrt(&r, 81, mod))
        flint_printf("n = 135: sqrt(81) = %wu\n", r);
    return 0;
}

n = 1000003: sqrt(10) = 394215, r^2 = 10
n = 1000037: sqrt(10) = 122489, r^2 = 10
n = 1000033: 10 is not a square
n = 469762049: sqrt(10) = 115267370, r^2 = 10
n = 135: sqrt(81) = 54

Modulo 135, the square roots of 81 in [0, 67] are 9, 36 and 54. The function returns 54, a genuine root but not the least one.