x. Quadratic Residues and Modular Roots

$a$ is a quadratic residue modulo an odd prime $p$ if $x^2\equiv a$ is solvable. Euler’s criterion: $\left(\frac{a}{p}\right)\equiv a^{(p-1)/2}(\mathrm{mod}p)$ ($+1$ residue, $-1$ non-residue, $0$ if $p\mid a$); exactly $\frac{p-1}{2}$ nonzero residues exist. Roots come from Tonelli–Shanks or Cipolla, lift to $p^k$ by Hensel, combine over composite moduli by CRT. Related Diophantine equations: $x^2+d{y}^2=m$ (Cornacchia), sums of two squares.

Task	Method	Cost
is $a$ a residue mod prime $p$	Euler’s criterion / Jacobi symbol	$O(\log p)$
Jacobi $\left(\frac{a}{n}\right)$ for odd $n$	reciprocity, like Euclid	$O(\log n)$
$\sqrt{a}\mathrm{mod}p$	Tonelli–Shanks ($O({\log }^{2}p)$) or Cipolla ($O(\log p)$)	—
$\sqrt{a}\mathrm{mod}{p}^k$	Hensel lifting	$O(k\log p)$
$\sqrt{a}\mathrm{mod}m$	roots mod each prime power + CRT ($2^{\omega }$ combinations)	—
$x^2+d{y}^2=p$	Cornacchia (Euclid on $(p,\sqrt{-d})$)	$O(\log p)$
$x^2+y^2=n$, count $r_2(n)$	factor; $r_2(n)=4(d_1(n)-d_3(n))$	factoring

int jacobi(ll a, ll n) {                                   // n odd > 0; Legendre symbol when n is prime
    a %= n; if (a < 0) a += n; int t = 1;
    while (a) {
        while (a % 2 == 0) { a /= 2; if (n % 8 == 3 || n % 8 == 5) t = -t; }
        swap(a, n); if (a % 4 == 3 && n % 4 == 3) t = -t;
        a %= n;
    }
    return n == 1 ? t : 0;
}
ll sqrtMod(ll a, ll p) {                                   // Tonelli-Shanks: some x with x^2 = a (mod p), -1 if none; p prime
    a %= p; if (a < 0) a += p;
    if (a == 0 || p == 2) return a;
    if (power(a, (p - 1) / 2, p) != 1) return -1;
    if (p % 4 == 3) return power(a, (p + 1) / 4, p);
    ll q = p - 1; int s = 0; while (q % 2 == 0) q /= 2, s++;
    ll z = 2; while (power(z, (p - 1) / 2, p) != p - 1) z++;
    ll c = power(z, q, p), x = power(a, (q + 1) / 2, p), t = power(a, q, p); int m = s;
    while (t != 1) {
        int i = 0; ll tt = t; while (tt != 1) tt = (i128)tt * tt % p, i++;
        ll b = c; for (int j = 0; j < m - i - 1; j++) b = (i128)b * b % p;
        x = (i128)x * b % p; c = (i128)b * b % p; t = (i128)t * c % p; m = i;
    }
    return x;
}
ll cipolla(ll a, ll p) {                                   // same contract, O(log p) multiplications in F_{p^2}
    a %= p; if (a < 0) a += p;
    if (a == 0 || p == 2) return a;
    if (power(a, (p - 1) / 2, p) != 1) return -1;
    ll r = 0, w = 0; mt19937_64 rng(7);
    do { r = rng() % p; w = ((i128)r * r - a + p) % p; } while (power(w, (p - 1) / 2, p) != p - 1);
    auto mul = [&](pair<ll, ll> x, pair<ll, ll> y) { return pair<ll, ll>{(ll)(((i128)x.first * y.first + (i128)x.second * y.second % p * w) % p), (ll)(((i128)x.first * y.second + (i128)x.second * y.first) % p)}; };
    pair<ll, ll> res{1, 0}, b{r, 1};
    for (ll e = (p + 1) / 2; e; e >>= 1, b = mul(b, b)) if (e & 1) res = mul(res, b);
    return res.first;                                      // (r + sqrt(w))^{(p+1)/2} lies in F_p
}
ll sqrtModPrimePower(ll a, ll p, int k) {                  // odd p, gcd(a, p) = 1: x^2 = a (mod p^k), -1 if none
    ll x = sqrtMod(a % p, p); if (x < 0) return -1;
    ll pk = p;
    for (int i = 1; i < k; i++) {                          // x <- x - (x^2 - a) / (2x) mod p^{i+1}
        pk *= p; ll am = ((a % pk) + pk) % pk;
        ll f = ((i128)x * x - am) % pk; if (f < 0) f += pk;
        x = (x - (i128)f * invMod(2 * x % pk, pk)) % pk; if (x < 0) x += pk;
    }
    return x;
}
pair<ll, ll> cornacchia(ll d, ll p) {                      // x^2 + d y^2 = p, p prime, 0 < d < p; {-1, -1} if none
    ll r = sqrtMod(p - d % p, p); if (r < 0) return {-1, -1};
    if (2 * r < p) r = p - r;                              // take the root > p / 2
    ll a = p, b = r, lim = (ll)sqrtl(p); while (lim * lim > p) lim--;
    while (b > lim) { ll t = a % b; a = b; b = t; }
    ll rest = p - b * b; if (rest % d) return {-1, -1};
    ll y = (ll)sqrtl(rest / d); while (y * y > rest / d) y--; while ((y + 1) * (y + 1) <= rest / d) y++;
    if (y * y * d != rest) return {-1, -1};
    return {b, y};
}

Sums of two squares in full (all representations of $n$, and $r_2(n)$), from the library:

pair<i64, i64> primeTwoSquares(i64 p) {
    if (p == 2) return {1, 1};
    auto mulmod = [&](i64 a, i64 b) { return (i64)((i128)a * b % p); };
    auto pw = [&](i64 b, i64 e) {
        i64 r = 1;
        for (b %= p; e; e >>= 1, b = mulmod(b, b)) {
            if (e & 1) r = mulmod(r, b);
        }
        return r;
    };
    i64 x = 0;
    for (i64 c = 2;; c++) {
        if (pw(c, (p - 1) / 2) == p - 1) {
            x = pw(c, (p - 1) / 4);
            break;
        }
    }
    i64 a = p, b = x;
    i64 lim = sqrtl((long double)p);
    while (b > lim) {
        i64 t = a % b;
        a = b, b = t;
    }
    i64 c = (i64)sqrtl((long double)(p - b * b));
    while (c * c > p - b * b) c--;
    while ((c + 1) * (c + 1) <= p - b * b) c++;
    return {b, c};
}
 
vector<pair<i64, i64>> twoSquares(i64 n) {
    vector<pair<i64, i64>> res;
    if (n < 0) return res;
    for (i64 x = 0; x * x <= n; x++) {
        i64 rest = n - x * x;
        i64 y = sqrtl((long double)rest);
        while (y * y > rest) y--;
        while ((y + 1) * (y + 1) <= rest) y++;
        if (y * y == rest) res.push_back({x, y});
    }
    return res;
}
 
i64 r2(i64 n) {
    if (n == 0) return 1;
    i64 d1 = 0, d3 = 0;
    auto count = [&](i64 e) {
        if (e % 4 == 1) d1++;
        else if (e % 4 == 3) d3++;
    };
    for (i64 d = 1; d * d <= n; d++) {
        if (n % d) continue;
        count(d);
        if (d != n / d) count(n / d);
    }
    return 4 * (d1 - d3);
}

Fact	Statement
Euler’s criterion	$a^{(p-1)/2}\equiv \left(\frac{a}{p}\right)$
$-1$, $2$ as residues	$-1$ is a QR iff $p\equiv 1(\mathrm{mod}4)$; $2$ iff $p\equiv \pm 1(\mathrm{mod}8)$
quadratic reciprocity	$\left(\frac{p}{q}\right)\left(\frac{q}{p}\right)=(-1{)}^{\frac{p-1}{2}\frac{q-1}{2}}$
number of roots of $x^2\equiv a(\mathrm{mod}{p}^k)$ (odd $p\nmid a$)	0 or 2; modulo $m$: product over prime powers ($2^k$ has up to 4)
sum of two squares	$n$ is one iff every prime $\equiv 3(\mathrm{mod}4)$ has even exponent; $r_2(n)=4{\sum }_{d\mid n}{\chi }_4(d)$
Fermat	prime $p=x^2+y^2$ iff $p=2$ or $p\equiv 1(\mathrm{mod}4)$

Proof: Euler's criterion and the Jacobi algorithm

With a primitive root $g$, $a=g^t$ is a square iff $t$ is even, and $a^{(p-1)/2}=g^{t(p-1)/2}$ is $1$ iff $t$ is even, else the unique element of order 2, $-1$; half of the $p-1$ exponents are even. The Jacobi symbol is multiplicative in both arguments; the loop uses $\left(\frac{2}{n}\right)=(-1{)}^{(n^2-1)/8}$ (the cases $n\equiv 3,5(\mathrm{mod}8)$) and reciprocity for odd coprime $a,n$ to swap the arguments, reducing like Euclid; a common factor ends with $n\ne 1$ and symbol 0.

Proof: Tonelli–Shanks and Cipolla

Tonelli–Shanks: write $p-1=q{2}^s$. Invariants: $x^2\equiv at$, $c$ has order $2^m$, $t$ has order $2^i$ with $i<m$. Initially $x=a^{(q+1)/2}$, $t=a^q$ (order dividing $2^{s-1}$ since $a$ is a residue), $c=z^q$ (order exactly $2^s$ for a non-residue $z$). Multiplying $x$ by $b=c^{2^{m-i-1}}$ and $t$ by $b^2$ (order $2^i$) keeps $x^2=at$ and makes the order of $t$ drop below $2^i$; with $c\leftarrow b^2$, $m\leftarrow i$ the invariants hold. When $t=1$, $x^2=a$. At most $s$ rounds of $O(s)$ squarings. $p\equiv 3(\mathrm{mod}4)$: $(a^{(p+1)/4}{)}^2=a^{(p+1)/2}=a\cdot a^{(p-1)/2}=a$. Cipolla: with $w=r^2-a$ a non-residue, ${\mathbb{F}}_p[\sqrt{w}]$ is a field ${\mathbb{F}}_{p^2}$ where $(r+\sqrt{w}{)}^p=r-\sqrt{w}$ (Frobenius, ${\sqrt{w}}^p=w^{(p-1)/2}\sqrt{w}=-\sqrt{w}$). So $(r+\sqrt{w}{)}^{p+1}=(r+\sqrt{w})(r-\sqrt{w})=r^2-w=a$, and $x=(r+\sqrt{w}{)}^{(p+1)/2}$ squares to $a$; since $a$ has a root in ${\mathbb{F}}_p$ and roots are unique up to sign in the field, $x\in {\mathbb{F}}_p$. A random $r$ gives a non-residue $w$ with probability about $\frac{1}{2}$.

Proof: Hensel lifting and Cornacchia

Hensel: if $x^2\equiv a(\mathrm{mod}{p}^i)$ and $p\nmid 2x$, set $x^{\prime }=x-\frac{x^2-a}{2x}$; then $x^{\prime 2}=x^2-(x^2-a)+\frac{(x^2-a{)}^2}{4x^2}\equiv a(\mathrm{mod}{p}^{2i})$, in particular modulo $p^{i+1}$, and the lift is unique (the derivative $2x$ is invertible). Cornacchia: if $p=x^2+d{y}^2$ with $\gcd (x,y)=1$, then $x\equiv \pm \sqrt{-d}y$, and running Euclid on $(p,r)$ with $r^2\equiv -d$, $r>p/2$, stops at the first remainder $<\sqrt{p}$, which is $x$; this follows from the continued-fraction characterisation of the convergents of $r/p$ (the pair $(x,y)$ appears as a convergent because $|r/p-k/y|<1/(2y^2)$ for the corresponding $k$). The final check verifies that $(p-x^2)/d$ is a perfect square; if not, no primitive representation exists.