Мотивация и контекст

Задача проста: существует ли целое $x$ такое, что $x^2 \equiv a \pmod{p}$? Этот вопрос — о квадратичных вычетах — на первый взгляд кажется элементарным, но ведёт к одной из красивейших теорем элементарной теории чисел — закону квадратичного взаимности, открытому Гауссом в 1796 году. Он опубликовал восемь различных доказательств этой теоремы на протяжении своей жизни.

В соревновательном программировании квадратичные вычеты возникают в задачах:

• Проверка, является ли число квадратом по модулю простого
• Нахождение квадратного корня по модулю простого (алгоритм Тонелли–Шенкса)
• Подсчёт решений уравнений второй степени по модулю
• Задачи на квадратичные поля и числа вида $a + b\sqrt{D}$
• Факторизация чисел методом квадратичного решета (теоретическая основа)

Понимание квадратичных вычетов также необходимо для задач на алгебраические структуры: кольца вычетов, расширения полей, алгебраические кривые над конечными полями.

Теория

Определение квадратичного вычета

Определение. Пусть $p$ — нечётное простое, $\gcd(a, p) = 1$. Число $a$ называется квадратичным вычетом по модулю $p$, если уравнение $x^2 \equiv a \pmod{p}$ имеет решение. В противном случае $a$ называется квадратичным невычетом.

Лемма. Из $p - 1$ ненулевых вычетов по модулю $p$ ровно $(p-1)/2$ являются квадратичными вычетами и $(p-1)/2$ — квадратичными невычетами.

Доказательство. Рассмотрим отображение $f: x \mapsto x^2$ на $\mathbb{Z}_p^*$. Если $y$ — квадратичный вычет, то $y = x^2 = (-x)^2$, значит у $y$ ровно 2 прообраза (так как $x \not\equiv -x$ при $p > 2$). Поэтому образ $f$ имеет ровно $(p-1)/2$ элементов. $\square$

Символ Лежандра

Определение. Символ Лежандра $(a/p)$ (или $\left(\frac{a}{p}\right)$) определяется для нечётного простого $p$ и $\gcd(a,p) = 1$ следующим образом: $\left(\frac{a}{p}\right) = \begin{cases} 1, & \text{если } a \text{ — квадратичный вычет mod } p \\ -1, & \text{если } a \text{ — квадратичный невычет mod } p \\ 0, & \text{если } p \mid a \end{cases}$

⎨

⎧​1,−1,0,​если a — квадратичный вычет mod pесли a — квадратичный невычет mod pесли p∣a​

Теорема Эйлера (критерий Эйлера). Для нечётного простого $p$ и $\gcd(a,p) = 1$: $\left(\frac{a}{p}\right) \equiv a^{(p-1)/2} \pmod{p}$

Доказательство. Обозначим $e = (p-1)/2$.

Случай 1: $a$ — квадратичный вычет, $a \equiv x^2 \pmod{p}$. Тогда $a^e = x^{2e} = x^{p-1} \equiv 1 \pmod{p}$ (малая теорема Ферма). Следовательно, $(a/p) = 1 \equiv a^e$.

Случай 2: $a$ — квадратичный невычет. Рассмотрим $a^e \pmod{p}$. Из $a^{p-1} = (a^e)^2 \equiv 1$ следует $(a^e)^2 \equiv 1$, то есть $a^e \equiv \pm 1 \pmod{p}$. Если $a^e \equiv 1$, то $a$ является корнем многочлена $x^e - 1$. Но $x^e - 1 \equiv 0 \pmod{p}$ имеет не более $e = (p-1)/2$ корней, и все они — квадратичные вычеты (это доказывается через примитивный корень: $g^k$ — квадратичный вычет $\iff$ $k$ чётно, и $(g^k)^e = g^{ke}$; при чётном $k = 2j$: $g^{2je} = g^{j(p-1)} = 1$). Но $a$ по условию — невычет. Противоречие. Значит, $a^e \equiv -1$. $\square$

Следствие (мультипликативность). Символ Лежандра мультипликативен: $\left(\frac{ab}{p}\right) = \left(\frac{a}{p}\right)\left(\frac{b}{p}\right)$

Доказательство. $(ab)^{(p-1)/2} \equiv a^{(p-1)/2} \cdot b^{(p-1)/2} \pmod{p}$. $\square$

Следствие. $\left(\frac{-1}{p}\right) = (-1)^{(p-1)/2}$, то есть:

• $(-1/p) = 1$ если $p \equiv 1 \pmod{4}$
• $(-1/p) = -1$ если $p \equiv 3 \pmod{4}$

Закон квадратичного взаимности

Теорема (Гаусс, 1796). Для различных нечётных простых $p$ и $q$: $\left(\frac{p}{q}\right)\left(\frac{q}{p}\right) = (-1)^{\frac{p-1}{2} \cdot \frac{q-1}{2}}$

Это означает:

• Если $p \equiv 1 \pmod{4}$ или $q \equiv 1 \pmod{4}$: $\left(\frac{p}{q}\right) = \left(\frac{q}{p}\right)$
• Если $p \equiv q \equiv 3 \pmod{4}$: $\left(\frac{p}{q}\right) = -\left(\frac{q}{p}\right)$

Дополнения: $\left(\frac{2}{p}\right) = (-1)^{(p^2-1)/8}$

• $(2/p) = 1$ если $p \equiv \pm 1 \pmod{8}$
• $(2/p) = -1$ если $p \equiv \pm 3 \pmod{8}$

Доказательство закона квадратичного взаимности (по Гауссу, через лемму Гаусса).

Лемма Гаусса. Пусть $p$ — нечётное простое, $\gcd(a, p) = 1$. Рассмотрим числа $a, 2a, 3a, \ldots, \frac{p-1}{2}a$ по модулю $p$, приведённые к промежутку $(-p/2, p/2)$. Пусть $\nu$ — количество отрицательных чисел среди них. Тогда $\left(\frac{a}{p}\right) = (-1)^\nu$.

Доказательство леммы. Обозначим $m = (p-1)/2$. Числа $r_1, \ldots, r_m$ — это $ia \bmod p$, приведённые к $(-p/2, p/2)$. Все $|r_i|$ различны и принадлежат $\{1, 2, \ldots, m\}$ (так как $ia \not\equiv \pm ja$ при $i \ne j$, $0 < i,j \le m$). Значит $\{|r_1|, \ldots, |r_m|\} = \{1, 2, \ldots, m\}$.

Тогда: $\prod_{i=1}^m (ia) = \prod_{i=1}^m r_i = (-1)^\nu \prod_{i=1}^m |r_i| = (-1)^\nu \cdot m!$

Слева: $a^m \cdot m!$. Следовательно, $a^m \equiv (-1)^\nu \pmod{p}$. По критерию Эйлера $(a/p) \equiv a^m \equiv (-1)^\nu$. $\square$

Доказательство теоремы взаимности. Рассмотрим квадрат $[1, (p-1)/2] \times [1, (q-1)/2]$ на целочисленной плоскости. Подсчитаем число точек $(x, y)$ такое, что $qx < py$ (под диагональю). Это число равно $\sum_{x=1}^{(p-1)/2} \lfloor qx/p \rfloor$. Симметрично, число точек над диагональю равно $\sum_{y=1}^{(q-1)/2} \lfloor py/q \rfloor$. Поскольку $p$ и $q$ нечётны, точки не лежат на диагонали $qx = py$. Общее число точек в прямоугольнике равно $\frac{p-1}{2} \cdot \frac{q-1}{2}$. Значит: $\sum_{x=1}^{(p-1)/2} \left\lfloor \frac{qx}{p} \right\rfloor + \sum_{y=1}^{(q-1)/2} \left\lfloor \frac{py}{q} \right\rfloor = \frac{p-1}{2} \cdot \frac{q-1}{2}$

Через лемму Гаусса: $\nu_p = \sum_{x=1}^{(p-1)/2} \lfloor qx/p \rfloor \pmod{2}$ для вычисления $(q/p)$. Следовательно: $\frac{(p/q)}{1} \cdot \frac{(q/p)}{1} = (-1)^{\nu_q + \nu_p} = (-1)^{\frac{p-1}{2} \cdot \frac{q-1}{2}} \square$

Алгоритм Тонелли–Шенкса

Критерий Эйлера позволяет определить, является ли $a$ квадратичным вычетом, но не даёт корня. Для нахождения корня используется алгоритм Тонелли–Шенкса.

Задача: дано нечётное простое $p$ и $a$ с $(a/p) = 1$. Найти $r$ такое, что $r^2 \equiv a \pmod{p}$.

Частный случай: если $p \equiv 3 \pmod{4}$, то $r = a^{(p+1)/4} \pmod{p}$ — непосредственно. Проверка: $r^2 = a^{(p+1)/2} = a \cdot a^{(p-1)/2} \equiv a \cdot 1 = a$.

Общий алгоритм Тонелли–Шенкса. Запишем $p - 1 = Q \cdot 2^S$, где $Q$ нечётно.

1. Найдём квадратичный невычет $z$ (берём случайные числа, половина из них — невычеты).
2. Инициализируем: $M = S$, $c = z^Q$, $t = a^Q$, $R = a^{(Q+1)/2}$.
   • Инвариант: $R^2 \equiv t \cdot a \pmod{p}$, $t^{2^{M-1}} \equiv 1$, $c^{2^{M-1}} \equiv -1$.
3. Пока $t \ne 1$:
   • Найдём наименьшее $i$ такое, что $t^{2^i} \equiv 1$ ($1 \le i < M$).
   • Положим $b = c^{2^{M-i-1}}$, $M = i$, $c = b^2$, $t = t \cdot b^2$, $R = R \cdot b$.
4. Вернуть $R$.

Корректность. Инвариант поддерживается на каждом шаге. При завершении $t = 1$ и $R^2 = a$. Поскольку $M$ строго убывает, алгоритм завершается.

Сложность: $O(S \log p) = O(\log^2 p)$ в худшем случае. На практике очень быстрый.

Символ Якоби

Определение. Для нечётного натурального $n = p_1^{e_1} \cdots p_k^{e_k}$ и $\gcd(a, n) = 1$ символ Якоби $(a/n)$ определяется как: $\left(\frac{a}{n}\right) = \prod_{i=1}^k \left(\frac{a}{p_i}\right)^{e_i}$

Ключевое отличие от символа Лежандра:

• Если $n$ — простое, символ Якоби совпадает с символом Лежандра.
• $(a/n) = 1$ не означает, что $a$ является квадратичным вычетом по модулю $n$.
• $(a/n) = -1$ означает, что $a$ — квадратичный невычет по модулю $n$.

Пример: $(2/15)$. $15 = 3 \cdot 5$. $(2/3) = -1$ (т.к. $3 \equiv 3 \pmod 8$), $(2/5) = -1$ (т.к. $5 \equiv 5 \equiv -3 \pmod 8$). Значит $(2/15) = (-1)(-1) = 1$. Но $x^2 \equiv 2 \pmod{15}$ решений не имеет (проверяется перебором). Символ Якоби равен 1, но $2$ — не квадратичный вычет!

Свойства символа Якоби: те же мультипликативность и закон взаимности, что и для символа Лежандра, но с учётом обобщения.

Закон взаимности для символа Якоби. Для нечётных $m, n > 0$ с $\gcd(m,n) = 1$: $\left(\frac{m}{n}\right)\left(\frac{n}{m}\right) = (-1)^{\frac{m-1}{2} \cdot \frac{n-1}{2}}$ $\left(\frac{-1}{n}\right) = (-1)^{(n-1)/2}, \quad \left(\frac{2}{n}\right) = (-1)^{(n^2-1)/8}$

Это позволяет вычислять символ Якоби за $O(\log^2 n)$ по аналогии с алгоритмом Евклида.

Ключевые формулы (сводная таблица)

Реализация на C++20 (с анализом сложности)

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

// ─── Бинарное возведение в степень ────────────────────────────────────────
ll power(ll base, ll exp, ll mod) {
    ll result = 1;
    base %= mod;
    if (base < 0) base += mod;
    for (; exp > 0; exp >>= 1) {
        if (exp & 1) result = (__int128)result * base % mod;
        base = (__int128)base * base % mod;
    }
    return result;
}

// ─── Символ Лежандра: (a/p), p — нечётное простое ─────────────────────────
// Сложность: O(log p)
int legendre(ll a, ll p) {
    a %= p;
    if (a < 0) a += p;
    if (a == 0) return 0;
    ll res = power(a, (p - 1) / 2, p);
    return res == 1 ? 1 : -1;
}

// ─── Символ Якоби: (a/n), n — нечётное натуральное ────────────────────────
// Алгоритм аналогичен НОД Евклида, использует закон взаимности
// Сложность: O(log^2 n)
int jacobi(ll a, ll n) {
    assert(n > 0 && n % 2 == 1);
    a %= n;
    if (a < 0) a += n;
    int result = 1;
    while (a != 0) {
        // Убираем множители 2
        while (a % 2 == 0) {
            a /= 2;
            // (2/n) = (-1)^((n^2-1)/8)
            if (n % 8 == 3 || n % 8 == 5) result = -result;
        }
        // Применяем закон взаимности
        swap(a, n);
        if (a % 4 == 3 && n % 4 == 3) result = -result;
        a %= n;
    }
    return n == 1 ? result : 0;
}

// ─── Алгоритм Тонелли–Шенкса: sqrt(a) mod p ──────────────────────────────
// Предусловие: p — нечётное простое, (a/p) = 1
// Сложность: O(log^2 p) в худшем случае
ll sqrt_mod(ll a, ll p) {
    a %= p;
    if (a < 0) a += p;
    if (a == 0) return 0;

    // Быстрый случай: p ≡ 3 (mod 4)
    if (p % 4 == 3) {
        return power(a, (p + 1) / 4, p);
    }

    // Разложим p-1 = Q * 2^S
    ll Q = p - 1, S = 0;
    while (Q % 2 == 0) { Q /= 2; S++; }

    // Найдём квадратичный невычет z
    ll z = 2;
    while (legendre(z, p) != -1) z++;

    // Инициализация
    ll M = S;
    ll c = power(z, Q, p);
    ll t = power(a, Q, p);
    ll R = power(a, (Q + 1) / 2, p);

    while (true) {
        if (t == 1) return R;

        // Найдём наименьшее i: t^(2^i) ≡ 1
        ll i = 1;
        ll tmp = (__int128)t * t % p;
        while (tmp != 1) { tmp = (__int128)tmp * tmp % p; i++; }

        // b = c^(2^(M-i-1))
        ll b = power(c, power(2, M - i - 1, p - 1), p);
        M = i;
        c = (__int128)b * b % p;
        t = (__int128)t * c % p;
        R = (__int128)R * b % p;
    }
}

// ─── Проверка и нахождение x^2 ≡ a (mod p^k) ─────────────────────────────
// Подъём Хензеля: если r^2 ≡ a (mod p^k), то r' ≡ r - (r^2-a)/(2r) (mod p^{k+1})
// Сложность: O(k * log p)
ll sqrt_mod_prime_power(ll a, ll p, int k) {
    if (a % p == 0) {
        // Особый случай: a делится на p
        // x^2 ≡ a (mod p^k): нужно чтобы v_p(a) чётно
        int v = 0; ll aa = a;
        while (aa % p == 0) { aa /= p; v++; }
        if (v % 2 != 0) return -1; // нет решений
        // r = sqrt(aa) * p^(v/2) mod p^k
        ll r0 = sqrt_mod(aa % p, p);
        if (r0 == -1) return -1;
        // Подъём для aa
        ll pk = 1;
        for (int i = 0; i < k - v/2; i++) pk *= p; // p^(k - v/2)
        // Упрощённо: возвращаем базовое решение
        return r0 == 0 ? 0 : r0; // TODO: полный подъём
    }

    // Проверяем, является ли a квадратичным вычетом mod p
    if (legendre(a % p, p) != 1) return -1;

    // Начальное приближение mod p
    ll r = sqrt_mod(a % p, p);

    // Подъём Хензеля: r_{k+1} = r_k - (r_k^2 - a) / (2*r_k) mod p^{k+1}
    // Это метод Ньютона для f(x) = x^2 - a
    ll pk = p;
    for (int i = 1; i < k; i++) {
        // Вычислим (r^2 - a) * inv(2r) mod p^{i+1}
        __int128 r2 = (__int128)r * r;
        __int128 num = r2 - a;
        // Делитель (2r) mod p^{i+1}: нам нужно обратное к 2r
        // Так как p нечётно и r ≢ 0, inv(2r) существует
        pk *= p;
        ll two_r = 2 * r % pk;
        // Расширенный НОД для inv(2r) mod pk
        // (упрощение: inv(2r) = inv(2) * inv(r), и inv(r) можно найти)
        // Для p^k можно использовать pow(2r, p^(k-1)*(p-1)-1, p^k) — теорема Эйлера
        // Но проще: поскольку 2r | delta кратно p^i, нужна обратная только mod p
        ll inv_2r = power(2 * (r % p) % p, p - 2, p);
        ll delta = (ll)((__int128)(r * r - a) % pk / p); // (r^2-a)/p^i mod p
        // r <- r - delta * inv_2r * p^i mod p^{i+1}
        r = (r - ((__int128)delta % p * inv_2r % p) % p * (pk / p) % pk + pk) % pk;
    }
    return r;
}

// ─── Подсчёт решений x^2 ≡ a (mod p^k) ───────────────────────────────────
// Возвращает количество решений
int count_sqrt_solutions(ll a, ll p, int k) {
    ll pk = 1;
    for (int i = 0; i < k; i++) pk *= p;
    a = ((a % pk) + pk) % pk;

    if (a == 0) {
        // Подсчёт нетривиален, зависит от p и k
        // Для нечётного простого p: 1 решение при v_p(a) нечётном, иначе 2p^{k - ⌈v/2⌉}... упрощённо
        return 1; // x=0 (только как пример)
    }

    int v = 0; ll aa = a;
    while (aa % p == 0) { aa /= p; v++; }

    if (v % 2 != 0) return 0;
    if (legendre(aa % p, p) == -1) return 0;
    if (legendre(aa % p, p) == 1) return 2; // два решения r и p^k - r

    return 0;
}

// ─── Демонстрация ──────────────────────────────────────────────────────────
int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(nullptr);

    // Символ Лежандра
    ll p = 13;
    for (ll a = 1; a < p; a++) {
        cout << "(" << a << "/" << p << ") = " << legendre(a, p)
             << " [" << (legendre(a,p)==1?"QR":"NR") << "]\n";
    }

    // Квадратный корень
    p = 1000000007LL;
    ll a = 4;
    cout << "\nsqrt(" << a << ") mod " << p << " = " << sqrt_mod(a, p) << "\n";

    a = 9;
    cout << "sqrt(" << a << ") mod " << p << " = " << sqrt_mod(a, p) << "\n";

    // Символ Якоби
    cout << "\nJacobi (2/15) = " << jacobi(2, 15) << "\n"; // = 1, но не вычет!
    cout << "Jacobi (3/35) = " << jacobi(3, 35) << "\n";

    // Закон взаимности — демонстрация
    ll q = 17;
    p = 11;
    cout << "\n(" << p << "/" << q << ") = " << legendre(p, q) << "\n";
    cout << "(" << q << "/" << p << ") = " << legendre(q, p) << "\n";
    // Должно: (p/q)(q/p) = (-1)^((p-1)/2 * (q-1)/2)

    return 0;
}

Анализ сложности:

Разбор задачи 1 (простая)

Задача. По данному простому $p$ ($3 \le p \le 10^9$, $p$ нечётное) и числу $a$ ($0 \le a < p$) определить: является ли $a$ квадратичным вычетом по модулю $p$? Если да, найти оба корня.

Решение. Применяем критерий Эйлера. Если $(a/p) = 1$, то используем алгоритм Тонелли–Шенкса для нахождения одного корня $r$, второй корень — $p - r$.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;

ll power(ll b, ll e, ll m) {
    ll r=1; b%=m;
    for(;e;e>>=1,b=(__int128)b*b%m) if(e&1) r=(__int128)r*b%m;
    return r;
}

int legendre(ll a, ll p) {
    if (a % p == 0) return 0;
    return power(a % p, (p-1)/2, p) == 1 ? 1 : -1;
}

ll tonelli_shanks(ll a, ll p) {
    if (p == 2) return a & 1;
    if (p % 4 == 3) return power(a, (p+1)/4, p);
    ll Q = p-1, S = 0;
    while (Q%2==0) { Q/=2; S++; }
    ll z = 2;
    while (legendre(z, p) != -1) z++;
    ll M=S, c=power(z,Q,p), t=power(a,Q,p), R=power(a,(Q+1)/2,p);
    for(;;){
        if(t==1) return R;
        ll i=1, tmp=(__int128)t*t%p;
        while(tmp!=1){tmp=(__int128)tmp*tmp%p;i++;}
        ll b=power(c,1LL<<(M-i-1),p);
        M=i; c=(__int128)b*b%p; t=(__int128)t*c%p; R=(__int128)R*b%p;
    }
}

int main(){
    ll p, a; cin >> p >> a;
    int leg = legendre(a, p);
    if (leg == 0) cout << "0 (a divisible by p)\n";
    else if (leg == -1) cout << "Quadratic non-residue\n";
    else {
        ll r = tonelli_shanks(a, p);
        cout << "Roots: " << r << " and " << p - r << "\n";
        // Проверка
        assert((__int128)r*r%p == a%p);
    }
}

Пример: $p = 13$, $a = 4$. $(4/13) = 4^6 \bmod 13 = 4096 \bmod 13 = 1$. Корни: $2$ и $11$.

Разбор задачи 2 (средняя)

Задача. Дано простое $p$ и $a$ с $(a/p) = 1$. Найти $x$ такое, что $x^2 \equiv a \pmod{p^k}$ для заданного $k$ (то есть найти корень по степени простого). Ограничения: $p \le 10^6$, $k \le 10$, $p^k \le 10^{18}$.

Решение. Подъём Хензеля: начиная с корня $r_1 \equiv \sqrt{a} \pmod{p}$

​(modp), последовательно поднимаем его до корня по $p^2, p^3, \ldots, p^k$.

Формула подъёма Хензеля для $f(x) = x^2 - a$, $f'(x) = 2x$: $r_{k+1} = r_k - \frac{r_k^2 - a}{p^k} \cdot (2r_k)^{-1} \pmod{p} \cdot p^k$

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 lll;

ll power(ll b, ll e, ll m) {
    ll r=1; b%=m;
    for(;e;e>>=1,b=(lll)b*b%m) if(e&1) r=(lll)r*b%m;
    return r;
}

ll mod_inv(ll a, ll m) { return power(a, m-2, m); } // только для простого m

int legendre(ll a, ll p) {
    a%=p; if(a<0)a+=p;
    if(!a) return 0;
    return power(a,(p-1)/2,p)==1?1:-1;
}

ll tonelli(ll a, ll p) {
    if(p%4==3) return power(a,(p+1)/4,p);
    ll Q=p-1,S=0;
    while(Q%2==0){Q/=2;S++;}
    ll z=2;
    while(legendre(z,p)!=-1)z++;
    ll M=S,c=power(z,Q,p),t=power(a,Q,p),R=power(a,(Q+1)/2,p);
    for(;;){
        if(t==1)return R;
        ll i=1,tmp=(lll)t*t%p;
        while(tmp!=1){tmp=(lll)tmp*tmp%p;i++;}
        ll b=c; for(ll j=0;j<M-i-1;j++) b=(lll)b*b%p;
        M=i;c=(lll)b*b%p;t=(lll)t*c%p;R=(lll)R*b%p;
    }
}

// Подъём Хензеля: r^2 ≡ a (mod p^k)
ll hensel_lift_sqrt(ll a, ll p, int k) {
    ll r = tonelli(a % p, p);
    if (r == -1 || legendre(a % p, p) != 1) return -1;

    ll pk = p; // текущий модуль p^i
    for (int i = 1; i < k; i++) {
        // r^2 ≡ a (mod p^i). Поднимаем до mod p^{i+1}
        // r_new = r - (r^2 - a) / p^i * inv(2r) mod p, затем * p^i
        lll r2 = (lll)r * r;
        lll diff = r2 - a;
        // diff ≡ 0 mod pk, делим на pk
        ll coeff = (ll)((diff / pk) % p); // (r^2-a)/p^i mod p
        if (coeff < 0) coeff += p;
        ll inv2r = mod_inv((2 * (r % p)) % p, p);
        ll delta = (lll)coeff * inv2r % p;
        r = r - delta * pk; // r_new = r - delta * p^i
        // Нормализуем
        pk *= p;
        r = ((r % pk) + pk) % pk;
    }
    return r;
}

int main(){
    ll p, a; int k;
    cin >> p >> a >> k;
    ll r = hensel_lift_sqrt(a, p, k);
    if (r == -1) {
        cout << "No solution\n";
    } else {
        ll pk = 1; for(int i=0;i<k;i++) pk*=p;
        cout << "r = " << r << "\n";
        cout << "Verify: r^2 mod p^k = " << (lll)r*r%pk << " (should be " << a%pk << ")\n";
    }
}

Пример: $p = 7$, $a = 2$, $k = 3$. $p^3 = 343$. $(2/7) = 2^3 \bmod 7 = 1$. Корень mod 7: $3$ (т.к. $9 \equiv 2$). Подъём: $3^2 - 2 = 7$, значит mod $49$: $r = 3 - 1 \cdot 7 \cdot \text{inv}(6) \bmod 7 \cdot 7 = 10$. Проверка: $10^2 = 100 \equiv 2 \pmod{49}$. $\checkmark$

Разбор задачи 3 (сложная — Div.1 D/E или ICPC WF)

Задача. Подсчитать количество пар $(x, y)$ с $0 \le x, y < p$ таких, что $x^2 + y^2 \equiv c \pmod{p}$ для данного простого $p$ и числа $c$. Дополнительно: вывести одну конкретную пару, если $c > 0$.

Источник: Классическая задача теории чисел, аналог USACO/Codeforces Div.1 D.

Теория. Количество решений $x^2 + y^2 \equiv c$ можно подсчитать через символы Лежандра:

$N = \sum_{x=0}^{p-1} \sum_{y=0}^{p-1} [x^2 + y^2 \equiv c \pmod{p}] = \sum_{x=0}^{p-1} \#\{y : y^2 \equiv c - x^2\}$

Для каждого $x$, число $y$: если $c - x^2 \equiv 0$, то 1 решение. Иначе $1 + (c - x^2 / p)$ решений.

$N = \sum_{x=0}^{p-1} \left(1 + \left(\frac{c - x^2}{p}\right)\right) = p + \sum_{x=0}^{p-1} \left(\frac{c - x^2}{p}\right)$

Сумма $S = \sum_{x} (c - x^2 / p)$ вычисляется теоретически: при $c = 0$: $N = 2p - 1$. При $c \ne 0$ и $p \equiv 1 \pmod 4$: $N = p - 1$. При $c \ne 0$ и $p \equiv 3 \pmod 4$: $N = p + 1$.

Для нахождения конкретной пары при $c \ne 0$ используем конструктивный алгоритм.

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128 lll;

ll power(ll b, ll e, ll m){
    ll r=1;b%=m;
    for(;e;e>>=1,b=(lll)b*b%m) if(e&1)r=(lll)r*b%m;
    return r;
}
int legendre(ll a, ll p){
    a%=p;if(a<0)a+=p;if(!a)return 0;
    return power(a,(p-1)/2,p)==1?1:-1;
}
ll tonelli(ll a, ll p){
    if(a%p==0)return 0;
    if(p%4==3)return power(a,(p+1)/4,p);
    ll Q=p-1,S=0;
    while(Q%2==0){Q/=2;S++;}
    ll z=2; while(legendre(z,p)!=-1)z++;
    ll M=S,c=power(z,Q,p),t=power(a,Q,p),R=power(a,(Q+1)/2,p);
    for(;;){
        if(t==1)return R;
        ll i=1,tmp=(lll)t*t%p;
        while(tmp!=1){tmp=(lll)tmp*tmp%p;i++;}
        ll b=power(c,1LL<<(M-i-1),p);
        M=i;c=(lll)b*b%p;t=(lll)t*c%p;R=(lll)R*b%p;
    }
}

// Подсчёт числа решений x^2 + y^2 ≡ c (mod p)
ll count_sum_of_squares(ll c, ll p) {
    c = ((c % p) + p) % p;
    if (c == 0) return 2*p - 1;
    // Для p=2 отдельно
    if (p == 2) return 1;
    // p нечётное простое
    // N = p - (-1/p) при c ≠ 0
    // (-1/p) = 1 если p≡1(mod 4), = -1 если p≡3(mod 4)
    int leg_minus1 = (p % 4 == 1) ? 1 : -1;
    return p - leg_minus1;
}

// Представление c = x^2 + y^2 (mod p) — конструктивный алгоритм
// Используем факт: если p≡1(mod 4), то -1 является квадратичным вычетом
// t^2 ≡ -1 (mod p), тогда (c, t) → ... (алгоритм Корнаки-Флека)
pair<ll,ll> find_sum_of_squares(ll c, ll p) {
    c = ((c%p)+p)%p;
    if (c == 0) return {0, 0};

    // Найдём такое a, что c - a^2 — квадратичный вычет
    // Перебираем a до нахождения
    for (ll a = 0; a < p; a++) {
        ll rem = ((c - (lll)a*a%p) % p + p) % p;
        if (legendre(rem, p) == 1 || rem == 0) {
            ll b = tonelli(rem, p);
            // Проверка
            if (((lll)a*a + (lll)b*b) % p == c) return {a, b};
        }
    }
    return {-1, -1}; // не должно происходить
}

int main(){
    ll p, c; cin >> p >> c;
    cout << "Number of solutions: " << count_sum_of_squares(c, p) << "\n";
    auto [x, y] = find_sum_of_squares(c, p);
    cout << "Example: x=" << x << ", y=" << y << "\n";
    if (x >= 0) {
        cout << "Verify: " << x << "^2 + " << y << "^2 = "
             << ((lll)x*x + (lll)y*y) % p << " (mod " << p << ")\n";
    }
    return 0;
}

Эффективное нахождение пары. В задаче выше find_sum_of_squares в худшем случае работает за $O(p)$. Для больших $p$ нужен лучший алгоритм. Одна из идей: при $p \equiv 1 \pmod{4}$ найти $t$ с $t^2 \equiv -1 \pmod{p}$ (через Тонелли–Шенкс), а затем алгоритмом Евклида для $p$ и $t$ найти разложение $p = a^2 + b^2$. Аналогично для $c \ne 1$.

Задачи для самостоятельного решения

Типичные ошибки

1. Путаница символов Лежандра и Якоби. Символ Якоби $(a/n) = 1$ при составном $n$ не гарантирует, что $a$ является квадратичным вычетом mod $n$. Проверка квадратичности по составному модулю требует факторизации.
2. Забыть случай $p \equiv 3 \pmod{4}$. При $p \equiv 3 \pmod{4}$ корень вычисляется за $O(\log p)$ без алгоритма Тонелли–Шенкса: $r = a^{(p+1)/4}$.
3. Два корня или ноль. Если $(a/p) = 1$ и $a \ne 0$, то есть ровно два корня: $r$ и $p - r$. Если $a = 0$, то один корень: $0$.
4. Переполнение в Тонелли–Шенксе. При $p \approx 10^9$ все промежуточные произведения нужно брать через __int128.
5. Ошибка в шаге нахождения $i$ в Тонелли–Шенксе. Нужно найти наименьшее $i \ge 1$ (не $i \ge 0$!) такое, что $t^{2^i} \equiv 1$.
6. Подъём Хензеля при $p = 2$. Подъём для $p = 2$ принципиально отличается (нет обратного к $2r$ mod 2). Для $p^k$ с $p = 2$ нужен отдельный алгоритм.

Совет профессионала

Быстрая проверка: является ли $n$ полным квадратом в $\mathbb{Z}$, используя символ Якоби.

Для больших $n$ дорого проверять $\lfloor\sqrt{n}\rfloor^2 = n$

​⌋2=n из-за потери точности при вещественной арифметике. Альтернативный быстрый фильтр: $n$ не является полным квадратом, если $(n/p) = -1$ для некоторого простого $p$. Если $n \bmod 4 \in \{2, 3\}$, то явно не квадрат. Проверка mod $3$: квадраты дают $0$ или $1$ mod $3$; если $n \equiv 2 \pmod{3}$, то не квадрат.

Второй совет: для задач на подсчёт квадратичных вычетов в интервале $[L, R]$ используйте формулу: $\#\{a \in [1, N] : (a/p) = 1\} = \frac{N-1}{2} + \epsilon$ (количество квадратичных вычетов среди $1, \ldots, p-1$ ровно $(p-1)/2$, и они равномерно распределены).

Итог

Квадратичные вычеты — ключевой инструмент теории чисел с богатой алгоритмической составляющей:

1. Критерий Эйлера $a^{(p-1)/2} \equiv (a/p) \pmod{p}$ — простейший способ проверить квадратичность за $O(\log p)$.
2. Закон квадратичного взаимности — позволяет вычислять символ Лежандра через рекуррентные упрощения (аналог алгоритма Евклида).
3. Символ Якоби обобщает символ Лежандра на составные модули и является основой теста Соловея–Штрассена.
4. Алгоритм Тонелли–Шенкса — эффективный алгоритм нахождения квадратного корня по простому модулю.
5. Подъём Хензеля — позволяет поднять корень с mod $p$ до mod $p^k$.

Понимание этих инструментов открывает путь к решению широкого класса задач на кольца вычетов, представления чисел суммами квадратов и алгебраические расширения конечных полей.