Legendre and Jacobi Symbols, in depth: quadratic residues, reciprocity, a gcd-speed algorithm and where primality tests use them

Quadratic residues from first principles

Fix an odd prime p. An integer a not divisible by p is a quadratic residue mod p if some x satisfies x squared = a (mod p), and a quadratic non-residue otherwise. Modulo 7 the squares of 1 to 6 are 1, 4, 2, 2, 4, 1, so the residues are {1, 2, 4} and the non-residues are {3, 5, 6}.

The count is not a coincidence. The map x to x squared sends x and p minus x to the same value and is otherwise one-to-one, so exactly (p-1)/2 of the nonzero values are residues and (p-1)/2 are not. Residues also behave like signs under multiplication: residue times residue is a residue, residue times non-residue is a non-residue, and, less obviously, non-residue times non-residue is a residue. That last fact follows from the multiplicative group mod p being cyclic: residues are exactly the even powers of a generator.

The Legendre symbol and Euler's criterion

The Legendre symbol (a/p) packages that sign: it is 1 if a is a nonzero residue mod p, -1 if a is a non-residue, and 0 if p divides a. The sign-like behaviour becomes a formula, (ab/p) = (a/p)(b/p), and the symbol depends only on a mod p.

Euler's criterion gives a direct way to compute it: (a/p) = a^((p-1)/2) mod p, where the result is always 1, p-1 (meaning -1) or 0. The reason: by Fermat's little theorem a^(p-1) = 1, so its square root a^((p-1)/2) must be 1 or -1, and it is 1 precisely for the even powers of a generator. With fast modular exponentiation this costs O(log p) multiplications of log-p-bit numbers, roughly O(log^3 p) bit operations with schoolbook multiplication. We can do better.

Quadratic reciprocity and the two supplements

Gauss's law of quadratic reciprocity relates two different questions: is p a square mod q, and is q a square mod p? For distinct odd primes p and q:

• (p/q)(q/p) = -1 if both p and q are 3 mod 4, and +1 otherwise. So you may swap top and bottom, flipping the sign only when both are 3 mod 4.
• First supplement: (-1/p) = 1 if p = 1 mod 4 and -1 if p = 3 mod 4.
• Second supplement: (2/p) = 1 if p = 1 or 7 mod 8 and -1 if p = 3 or 5 mod 8.

With multiplicativity, these rules let you evaluate any Legendre symbol by hand: reduce the top mod the bottom, factor out 2s and -1, swap using reciprocity, and repeat. The catch is that after a swap the new bottom number is whatever the old top was, which may not be prime, and factoring it to continue is exactly the expensive step we want to avoid. The Jacobi symbol removes that obstacle.

The Jacobi symbol: same rules, any odd modulus

For an odd positive n with prime factorization p1 p2 ... pk (repeats allowed), the Jacobi symbol is defined as the product (a/n) = (a/p1)(a/p2)...(a/pk). When n is prime it is the Legendre symbol. It is 0 exactly when gcd(a, n) > 1.

The remarkable fact is that the Jacobi symbol obeys the same laws as the Legendre symbol: it is multiplicative in both the top and the bottom, depends only on a mod n, and satisfies reciprocity and both supplements with n in place of p, for any odd positive coprime m and n. Those laws never mention primality, so we can apply them to whatever number lands in the bottom position without factoring it. The loop has the same shape as Euclid's gcd algorithm.

The algorithm

Here is the whole algorithm. It accepts any integer a, including negative ones, and any odd positive n:

def jacobi(a: int, n: int) -> int:
    """Jacobi symbol (a/n) for odd n > 0. Returns 1, -1, or 0 (0 when gcd(a, n) > 1)."""
    if n <= 0 or n % 2 == 0:
        raise ValueError("n must be a positive odd integer")
    a %= n                       # also handles negative a
    result = 1
    while a != 0:
        while a % 2 == 0:        # pull out factors of 2: (2/n) = -1 iff n = 3 or 5 mod 8
            a //= 2
            if n % 8 in (3, 5):
                result = -result
        a, n = n, a              # reciprocity: flip iff both are 3 mod 4
        if a % 4 == 3 and n % 4 == 3:
            result = -result
        a %= n
    return result if n == 1 else 0

Why it terminates fast: each swap-and-reduce is one step of the Euclidean algorithm on (a, n), and stripping twos only shrinks a. The number of iterations is therefore O(log n), and each iteration does a division on O(log n)-bit numbers, so the total is O(log^2 n) bit operations with schoolbook arithmetic, a factor of log n better than Euler's criterion. For why that matters once numbers outgrow a machine word, see Big-O analysis, which covers the cost of arithmetic on multi-word integers.

Two details are easy to get wrong. The sign flip for 2 uses the current n, not the original. And the result must be 0 when the loop ends with n other than 1, because the final n is then gcd(a, n).

Worked example: (1001/9907) by hand

9907 is prime, so this is the Legendre symbol from the introduction. Follow the loop:

(1001 / 9907)          9907 is prime, so this is a Legendre symbol
 swap   -> (9907 / 1001)   1001 = 1 mod 4: no flip          sign +1
 reduce -> (898 / 1001)
 strip 2 -> (449 / 1001)   1001 = 1 mod 8: no flip          sign +1
 swap   -> (1001 / 449)    449 = 1 mod 4: no flip
 reduce -> (103 / 449)
 swap   -> (449 / 103)     449 = 1 mod 4: no flip
 reduce -> (37 / 103)
 swap   -> (103 / 37)      37 = 1 mod 4: no flip
 reduce -> (29 / 37)
 swap   -> (37 / 29)       29 = 1 mod 4: no flip
 reduce -> (8 / 29)
 strip 2 three times       29 = 5 mod 8: flip, flip, flip   sign -1
 -> (1 / 29) = 1, loop ends with n = 1                       answer -1
check: pow(1001, 4953, 9907) = 9906 = -1 mod 9907

The answer is -1, so 1001 is not a square mod 9907, and no x satisfies x squared = 1001 (mod 9907). Notice that the numbers that appeared in the bottom position, 1001, 449, 103, 37 and 29, were never checked for primality. 1001 = 7 x 11 x 13 is composite, and the algorithm did not care. Euler's criterion confirms the answer: 1001^4953 mod 9907 is 9906, which is -1.

Test it against a brute-force reference

Number-theory code fails silently: a wrong sign produces a plausible-looking 1 or -1. Test the fast function against the definitions on every small case, which runs in about a second:

def legendre_euler(a: int, p: int) -> int:        # reference: Euler's criterion
    r = pow(a, (p - 1) // 2, p)
    return -1 if r == p - 1 else r

def odd_prime_factors(n: int) -> list:            # with repeats, e.g. 45 -> [3, 3, 5]
    f, d = [], 3
    while d * d <= n:
        while n % d == 0:
            f.append(d)
            n //= d
        d += 2
    return f + ([n] if n > 1 else [])

def check(limit=600):
    for n in range(3, limit, 2):
        f = odd_prime_factors(n)
        for a in range(-50, 2 * n):
            expect = 1                             # product definition; one factor if n is prime
            for q in f:
                expect *= legendre_euler(a % q, q)
            assert jacobi(a, n) == expect, (a, n)

Checking composites matters: a bug in the 2-stripping or swap rule can pass every prime and still be wrong for composite moduli.

The trap: Jacobi = 1 does not mean square

For a prime modulus, (a/p) = 1 means a is a square. For a composite modulus that is false. (2/15) = (2/3)(2/5) = (-1)(-1) = 1, yet the squares mod 15 are 0, 1, 4, 6, 9 and 10, so 2 is not a square mod 15. The Jacobi symbol only guarantees the implication in one direction: if (a/n) = -1, then a is certainly not a square mod n. If it is 1, you know nothing without the factorization.

This is not a bug; it is the property some cryptosystems are built on. Deciding squareness mod n = pq for numbers with Jacobi symbol 1 is believed to be as hard as factoring n, and that assumption underpins the Goldwasser-Micali encryption scheme.

Where software uses the symbols

Solovay-Strassen primality test. For a prime n, Euler's criterion says a^((n-1)/2) = (a/n) mod n for every a. For an odd composite n, at most half of the units satisfy that congruence, so a random a exposes a composite with probability at least 1/2 per round:

import secrets
from math import gcd

def solovay_strassen(n: int, rounds: int = 40) -> bool:
    """False: n is certainly composite. True: n is prime with error below 2**-rounds."""
    if n < 2 or n % 2 == 0:
        return n == 2
    if n == 3:
        return True
    for _ in range(rounds):
        a = 2 + secrets.randbelow(n - 3)        # a in [2, n-2]
        if gcd(a, n) != 1:
            return False
        if pow(a, (n - 1) // 2, n) != jacobi(a, n) % n:   # -1 becomes n-1
            return False
    return True

The Carmichael number 561 = 3 x 11 x 17 fools the Fermat test for every base coprime to it, but only 80 of its 320 units are Euler liars. In practice Miller-Rabin has replaced Solovay-Strassen, since its liars are at most a quarter of the bases and every Miller-Rabin liar is also an Euler liar, but Solovay-Strassen remains the cleanest illustration of the idea.

Baillie-PSW. This widely used combined test runs a strong base-2 Miller-Rabin test and a strong Lucas test. Selfridge's method for choosing the Lucas parameter walks D through 5, -7, 9, -11, 13, ... until (D/n) = -1; for n = 1009 it stops at D = -11. If some (D/n) = 0 with |D| < n, n has a small factor. Note that a perfect square n never yields -1, so implementations check for squares first or cap the search.

Square roots mod p. The Tonelli-Shanks algorithm needs one known non-residue. Trying z = 2, 3, 4, ... and testing (z/p) = -1 finds one quickly, since half the candidates qualify; for p = 41 it is 3. It also uses the symbol up front to reject inputs that have no root.

Factoring and elliptic curves. The quadratic sieve builds its factor base from primes p with (N/p) = 1, since only those can divide the values it sieves. Point decompression and hash-to-curve routines ask whether x^3 + ax + b is a square in the field before taking a root.

Failure modes and operational guidance

In production code, use a library rather than your own loop: GMP provides mpz_jacobi, gmpy2 exposes gmpy2.jacobi, and SymPy has jacobi_symbol. Write your own only to learn, for a language without big-integer support, or when you need an audited constant-time version, and then test it as above. Where the symbol meets cryptography, as with the signature checks in JWT validation or key handling in KMS envelope encryption, the rule is the same: rely on vetted primitives. The same modular arithmetic reappears in universal hashing, where table sizes and multipliers interact with primes, as discussed in hash tables.

What to do next

1. Compute the residues mod 7, 11 and 13 by hand and confirm there are (p-1)/2 of each kind.
2. Implement jacobi() from the code above and run the brute-force test against Euler's criterion and the product definition.
3. Trace (1001/9907) yourself, then pick another pair and check it with pow(a, (p-1)//2, p).
4. Show that (2/15) = 1 while 2 is not a square mod 15, so the trap stays memorable.
5. Implement Solovay-Strassen and count the Euler liars of 561 to see the 1/2 bound in action.
6. Read how your big-integer library implements its primality test; many use Baillie-PSW with Selfridge's D.
7. For anything cryptographic, switch to the library's constant-time routine.