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

Задача "найти количество способов представить n в виде суммы двух квадратов a² + b²" встречается на соревнованиях неожиданно часто. Наивный подход — перебрать все a ≤ √n и проверить, является ли n − a² точным квадратом — работает за O(√n). Но если нужно обработать 10^5 запросов с n ≤ 10^12, это слишком медленно.

Ответ скрывается в алгебре: кольцо гауссовых целых чисел Z[i] = {a + bi : a, b ∈ Z} является евклидовой областью, и в ней работает аналог основной теоремы арифметики — единственность разложения на простые. Понимая структуру Z[i], мы можем:

• За O(√n · log n) предобработки отвечать на запросы о представлении n = a² + b² за O(1)
• Подсчитывать решётчатые точки на окружностях и в кругах
• Решать задачи о квадратичных формах

Гауссовы целые числа — первый пример алгебраических целых чисел, и изучение Z[i] даёт интуицию для более общей теории.

---

Теория

Определения: кольцо Z[i]

Определение. Гауссово целое число — это число вида z = a + bi, где a, b ∈ Z, i² = −1. Множество всех таких чисел обозначается Z[i].

Норма. Для z = a + bi норма N(z) = a² + b². Норма мультипликативна: N(zw) = N(z)·N(w).

Сопряжение. Сопряжение z̄ = a − bi. Тогда z · z̄ = N(z).

Единицы Z[i]. Это элементы с нормой 1: {1, −1, i, −i}. Все остальные элементы не являются обратимыми.

Делимость. z делит w (z | w) в Z[i], если существует q ∈ Z[i] такое, что w = q·z.

Евклидово деление. Для z, w ∈ Z[i], w ≠ 0, существуют q, r такие, что z = qw + r и N(r) < N(w). Это доказывается геометрически: точки кратные w образуют квадратную решётку с шагом √N(w), и любая точка z находится на расстоянии < √N(w) от ближайшего узла. Отсюда следует, что Z[i] — евклидова область и, в частности, область с единственным разложением на простые (UFD).

Простые числа в Z[i]

Элемент π ∈ Z[i] называется гауссовым простым, если он не является единицей и не раскладывается в произведение двух неединичных элементов.

Теорема (классификация гауссовых простых). Каждое гауссово простое ассоциировано ровно с одним из следующих элементов:

1. (1 + i) — единственный гауссов простой с нормой 2. Действительно, 2 = −i(1+i)², а N(1+i) = 2.
2. a + bi, где a² + b² = p ≡ 1 (mod 4) — гауссово простое раскладывает обычное простое p ≡ 1 (mod 4) на два сопряжённых гауссовых простых: p = (a+bi)(a−bi).
3. p, где p ≡ 3 (mod 4) — обычные простые p ≡ 3 (mod 4) остаются простыми в Z[i] (называются инертными).

Доказательство. Рассмотрим обычное простое p.

Случай p = 2. 2 = −i(1+i)² — не просто в Z[i].

Случай p ≡ 1 (mod 4). По теореме о первообразном корне и теореме Вильсона, −1 — квадратный вычет mod p. Значит, существует a с a² ≡ −1 (mod p), т.е. p | a² + 1 = (a+i)(a−i). Если бы p было простым в Z[i], оно делило бы a+i или a−i. Но N(p) = p², а N(a+i) = a²+1. Из p | a+i следует a/p + i/p ∈ Z[i] — невозможно, т.к. i/p ∉ Z[i]. Значит, p не простое в Z[i], и поскольку Z[i] — UFD, p раскладывается на два гауссовых простых нормы p.

Случай p ≡ 3 (mod 4). Предположим, p = (a+bi)(c+di). Тогда p² = N(p) = (a²+b²)(c²+d²). Если a²+b² = p, то p является суммой двух квадратов. Но по лемме Ферма (которую мы докажем ниже), p ≡ 3 (mod 4) не является суммой двух квадратов — противоречие. Значит, разложение тривиально, и p простое в Z[i]. □

Теорема Ферма о суммах двух квадратов

Теорема (Ферма, 1640; доказательство Гаусса). Натуральное число n представимо в виде суммы двух квадратов n = a² + b² тогда и только тогда, когда в разложении n на простые каждый простой делитель вида p ≡ 3 (mod 4) входит в чётной степени.

Доказательство (через Z[i]).

Направление "если". Каждый простой p ≡ 1 (mod 4) раскладывается в Z[i] как p = (a+bi)(a−bi), и N(a+bi) = p — сумма двух квадратов. Простое 2 = 1² + 1². Простое p ≡ 3 (mod 4) в чётной степени: p^{2k} = (p^k)² + 0². Поскольку N мультипликативна, произведение сумм двух квадратов — снова сумма двух квадратов (формула Брахмагупты–Фибоначчи):

(a²+b²)(c²+d²) = (ac−bd)² + (ad+bc)²

Направление "только если". Предположим, n = a² + b² и p | n, p ≡ 3 (mod 4). Тогда в Z[i]: p | a²+b² = (a+bi)(a−bi). Поскольку p простое в Z[i], p | a+bi или p | a−bi. Оба случая дают p | a и p | b. Значит, p² | n. Индукция по степени p в n. □

Следствие. Число способов представить n = a² + b² (с учётом порядка и знаков) равно:

r₂(n) = 4 · (d₁(n) − d₃(n))

где d₁(n) — число делителей n вида 4k+1, d₃(n) — число делителей n вида 4k+3.

Альтернативная формула. Если n = 2^a₀ · ∏ p_i^{a_i} · ∏ q_j^{b_j}, где p_i ≡ 1 (mod 4) и q_j ≡ 3 (mod 4), то:

• n представимо как сумма двух квадратов ⟺ все b_j чётны
• Число представлений r₂(n) = 4 · ∏(a_i + 1), если все b_j чётны; иначе 0

Алгоритм нахождения разложения p = a² + b²

Для p ≡ 1 (mod 4) нужно явно найти a, b с a² + b² = p.

Алгоритм Корнакки–Смита (основан на алгоритме Гаусса):

1. Найти квадратный корень из −1 по модулю p: найти x с x² ≡ −1 (mod p).
   • Это можно сделать так: взять любой генератор g мультипликативной группы F_p^*, тогда x = g^{(p−1)/4} (mod p).
2. Применить алгоритм Евклида к паре (p, x), остановившись, когда остаток r₀ < √p.
   • Тогда: r₀² + r₁² = p (или r₁ = 1 при вырожденном случае).

Почему это работает? Из x² ≡ −1 (mod p) следует p | x²+1 = (x+i)(x−i) в Z[i]. Так как p не простое в Z[i], существует нетривиальный гауссов делитель. Алгоритм Евклида в Z[i], применённый к (p, x+i), находит этот делитель — и его норма равна p, т.е. это a+bi с a²+b² = p.

// Нахождение a² + b² = p для простого p ≡ 1 (mod 4)
pair<long long, long long> sum_of_squares(long long p) {
    // Шаг 1: найти x с x² ≡ -1 (mod p)
    long long x = 0;
    // Перебираем квадратные невычеты и берём корень
    for (long long g = 2; ; g++) {
        // Проверяем, что g — примитивный корень
        if (powmod(g, (p-1)/2, p) != p-1) continue;
        x = powmod(g, (p-1)/4, p);
        if (x * x % p == p - 1) break; // x² ≡ -1 (mod p)
    }
    // Шаг 2: алгоритм Евклида для (p, x), останавливаемся при r < √p
    long long r0 = p, r1 = x;
    while (r1 * r1 > p) {
        long long r2 = r0 % r1;
        r0 = r1;
        r1 = r2;
    }
    // Проверяем: r1² + r0_new² = p? Нет, нужен также второй остаток
    // Корректная реализация через шаг Евклида:
    long long a = r1;
    long long b_sq = p - a * a;
    long long b = (long long)round(sqrt((double)b_sq));
    // b может быть вычислен из цепочки: при остановке r1 < √p,
    // выполняется r1² + (r0 % r1)² = p не всегда. Используем:
    // После остановки последний r1 = a, предпоследний r0.
    // b = r0 % r1 правильно только в специальных случаях.
    // Правильный вариант:
    long long s0 = 1, s1 = 0;
    r0 = p; r1 = x; s0 = 0; s1 = 1;
    while ((long long)r1 * r1 > p) {
        long long q = r0 / r1;
        long long r2 = r0 - q * r1;
        long long s2 = s0 - q * s1; // не нужно, но для полноты
        r0 = r1; r1 = r2;
    }
    a = r1;
    // Вычисляем b из a² + b² = p
    b = (long long)sqrt((double)(p - a*a) + 0.5);
    while (b * b + a * a != p) b++;
    return {a, b};
}

Разложение составных чисел в Z[i]

Для составного n сначала факторизуем n в Z (Pollard ρ + Miller-Rabin), затем каждый простой делитель раскладываем в Z[i]:

• p = 2: 2 = −i(1+i)²
• p ≡ 1 (mod 4): p = (a+bi)(a−bi), где a²+b² = p (алгоритм выше)
• p ≡ 3 (mod 4): p остаётся простым в Z[i]

Затем перемножаем гауссовы делители с учётом правил кольца.

Подсчёт решётчатых точек

Число решётчатых точек (a, b) ∈ Z² с a²+b² ≤ R приближается к πR (площадь круга). Точная формула:

#{(a,b) : a²+b² ≤ R} = πR + O(R^{1/3+ε})

Для подсчёта решётчатых точек на окружности (a²+b² = n) используем r₂(n).

---

Ключевые формулы

---

Реализация на C++20

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

// ==================== Арифметика гауссовых целых ====================

struct GaussInt {
    ll re, im;
    GaussInt(ll r = 0, ll i = 0) : re(r), im(i) {}
    GaussInt operator+(const GaussInt& o) const { return {re+o.re, im+o.im}; }
    GaussInt operator-(const GaussInt& o) const { return {re-o.re, im-o.im}; }
    GaussInt operator*(const GaussInt& o) const {
        return {re*o.re - im*o.im, re*o.im + im*o.re};
    }
    GaussInt conj() const { return {re, -im}; }
    ll norm() const { return re*re + im*im; }
    bool operator==(const GaussInt& o) const { return re==o.re && im==o.im; }
    // Деление в Z[i] (округление к ближайшему)
    GaussInt operator/(const GaussInt& o) const {
        ll n = o.norm();
        ll q_re = (ll)round((double)(re*o.re + im*o.im) / n);
        ll q_im = (ll)round((double)(im*o.re - re*o.im) / n);
        return {q_re, q_im};
    }
    GaussInt operator%(const GaussInt& o) const {
        return *this - (*this / o) * o;
    }
};

// НОД в Z[i] (алгоритм Евклида)
GaussInt gcd_gaussian(GaussInt a, GaussInt b) {
    while (b.norm() != 0) {
        a = a % b;
        swap(a, b);
    }
    return a; // нормализовать при необходимости
}

// ==================== Нахождение a² + b² = p ====================

ll powmod(ll a, ll b, ll m) {
    ll res = 1;
    a %= m;
    while (b > 0) {
        if (b & 1) res = (__int128)res * a % m;
        a = (__int128)a * a % m;
        b >>= 1;
    }
    return res;
}

// Находим x: x² ≡ -1 (mod p), p ≡ 1 (mod 4)
ll sqrt_minus1(ll p) {
    // По теореме, g^{(p-1)/4} — корень из -1, если g — примитивный корень
    for (ll g = 2; ; g++) {
        ll x = powmod(g, (p-1)/4, p);
        if ((__int128)x * x % p == p - 1) return x;
    }
}

// Разложение простого p ≡ 1 (mod 4) в сумму двух квадратов
// Возвращает (a, b) с a² + b² = p, a > 0, b > 0, a ≤ b
pair<ll,ll> fermat_two_squares(ll p) {
    if (p == 2) return {1, 1};
    assert(p % 4 == 1);
    ll x = sqrt_minus1(p);
    // Алгоритм Корнакки-Смита: Евклид до r < √p
    ll r0 = p, r1 = x;
    while (r1 * r1 > p) {
        ll r2 = r0 % r1;
        r0 = r1;
        r1 = r2;
    }
    ll a = r1;
    ll b = (ll)round(sqrt((double)(p - a*a)));
    // Точная проверка
    while (a*a + b*b < p) b++;
    while (a*a + b*b > p) b--;
    assert(a*a + b*b == p);
    if (a > b) swap(a, b);
    return {a, b};
}

// ==================== Подсчёт представлений n = a² + b² ====================

// Факторизация (простая версия для демонстрации)
map<ll,int> factorize_small(ll n) {
    map<ll,int> f;
    for (ll p = 2; p*p <= n; p++) {
        while (n % p == 0) { f[p]++; n /= p; }
    }
    if (n > 1) f[n]++;
    return f;
}

// r₂(n) / 4 = d₁(n) - d₃(n), где d_k(n) = #{делители n ≡ k (mod 4)}
// Через мультипликативность: для p^e:
//   p=2:          вклад = 1 (нейтрально)
//   p≡1(4), p^e:  вклад d₁ - d₃ = e+1 (все делители ≡ 1 или i, норма нейтральна)
//   p≡3(4), p^e:  вклад = 1 если e чётно, 0 если нечётно (т.к. (-1)^e)
ll count_representations(ll n) {
    auto f = factorize_small(n);
    ll result = 1;
    for (auto [p, e] : f) {
        if (p == 2) {
            // нейтральный вклад
        } else if (p % 4 == 1) {
            result *= (e + 1);
        } else { // p ≡ 3 (mod 4)
            if (e % 2 != 0) return 0;
            // вклад = 1
        }
    }
    return 4 * result; // r₂(n) = 4 * (d₁ - d₃)
}

// Все представления n = a² + b², a ≥ 0, b ≥ 0
// Используем: если n = ∏ (a_i + b_i·i) · ∏ (a_i - b_i·i) · ∏ q_j^{2k_j}
// то все гауссовы делители нормы n дают все представления
vector<pair<ll,ll>> all_representations(ll n) {
    vector<pair<ll,ll>> res;
    for (ll a = 0; a*a <= n; a++) {
        ll b2 = n - a*a;
        ll b = (ll)sqrt((double)b2);
        while (b*b < b2) b++;
        while (b*b > b2) b--;
        if (b*b == b2 && a <= b) {
            res.push_back({a, b});
        }
    }
    return res;
}

// ==================== Арифметика в Z[i] mod (a + bi) ====================
// Для вычисления в факторкольце Z[i]/(π), где π — гауссово простое

struct GaussIntMod {
    ll re, im;
    ll p; // реальное простое (норма π)
    GaussInt pi; // гауссов простой
    // ... (сложная реализация, опускается для краткости)
};

// ==================== Пример: задача на решётки ====================

// Число целых точек (a, b) с 1 ≤ a ≤ n, gcd(a, b) = 1, a² + b² = p для простого p
// = r₂(p) / 4 (с учётом знаков и порядков)

int main() {
    // Тест: представления 50 = 2 · 5² = 2 · (2²+1²)²
    // 5 = 1² + 2², 2 = 1² + 1²
    // 50 = 1² + 7² = 5² + 5²
    cout << "r₂(50) = " << count_representations(50) << "\n"; // должно быть 12
    // (±1, ±7), (±7, ±1), (±5, ±5) = 4+4+4 = 12 ✓

    // Разложение 13 = 2² + 3²
    auto [a, b] = fermat_two_squares(13);
    cout << "13 = " << a << "² + " << b << "² = "
         << a*a << " + " << b*b << "\n"; // 4 + 9 = 13 ✓

    // Разложение 997 ≡ 1 (mod 4)
    auto [a2, b2] = fermat_two_squares(997);
    cout << "997 = " << a2 << "² + " << b2 << "²\n";

    return 0;
}

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

• fermat_two_squares(p): поиск g — O(log p) итераций в среднем (генераторов ≈ φ(p-1)/(p-1) ≈ 1/4), каждая — O(log p). Алгоритм Евклида — O(log p). Итого: O(log² p).
• count_representations(n): O(√n) с пробным делением, O(n^{1/4} log n) с Полларда ρ.
• all_representations(n): O(√n) — тривиальный перебор.

---

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

Задача. Дано n ≤ 10^9. Найти количество способов записать n = a² + b² (с учётом порядка и знаков).

Решение. Используем формулу r₂(n) = 4(d₁(n) − d₃(n)). Факторизуем n за O(√n). Для каждого простого p:

• p = 2: не влияет на разность
• p ≡ 1 (mod 4), степень e: умножаем результат на (e+1)
• p ≡ 3 (mod 4), степень e: если e нечётно, ответ 0; если чётно, не меняем

ll solve(ll n) {
    return count_representations(n); // из реализации выше
}

Пример. n = 325 = 5² · 13. 5 ≡ 1 (mod 4), e=2 → вклад 3. 13 ≡ 1 (mod 4), e=1 → вклад 2. r₂(325) = 4 · 3 · 2 = 24. Проверка: 325 = 1²+18² = 6²+17² = 10²+15² = 18²+1² = ... действительно 24 пары.

---

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

Задача. Дано простое p ≡ 1 (mod 4). Найти a, b > 0 с a² + b² = p. p ≤ 10^15.

Решение. Применяем алгоритм Корнакки-Смита. Детали реализации уже даны выше.

Тонкость: при поиске корня из −1 по mod p нужен примитивный корень g. Простой способ: перебираем g = 2, 3, 4, ... и проверяем powmod(g, (p−1)/2, p) == p−1 (тест примитивного корня), затем x = powmod(g, (p−1)/4, p).

Но перебор примитивных корней может быть долгим. Оптимизация: достаточно проверить powmod(g, (p−1)/4, p) ≠ 1 и powmod(g, (p−1)/4, p) ≠ p−1. Если ни то ни другое — тогда x = powmod(g, (p−1)/4, p) работает.

Время: O(log² p) ≈ O(50²) = 2500 операций для p ≤ 10^15. Мгновенно.

---

Разбор задачи 3 (ICPC WF / CF Div.1 E)

Задача (CF 424E, адаптация). Дано N точек на плоскости. Для каждой пары (i, j) расстояние d_{ij} = √(Δx² + Δy²). Подсчитать число пар с целым расстоянием.

Анализ. Расстояние целое ⟺ Δx² + Δy² — точный квадрат. Пусть S_k = {(a,b) : a²+b² = k²} для k = 1, 2, .... Нам нужно для каждой пары проверить, является ли Δx² + Δy² точным квадратом.

Ключевое наблюдение. Δx² + Δy² = k² ⟺ в Z[i]: (Δx + Δy·i)(Δx − Δy·i) = k². Это выполняется тогда и только тогда, когда в разложении (Δx + Δy·i) на гауссовы простые все гауссовы простые πα имеют степени, дающие в произведении чистый квадрат.

Алгоритм:

1. Для каждой пары вычислить d² = Δx² + Δy².
2. Факторизовать d² за O(d^{1/2}).
3. Проверить, является ли d² точным квадратом: √d² ∈ Z.

Более эффективно для N = 2000, N² = 4·10^6 пар: при max coord ≤ 10^6, d² ≤ 2·10^12 — Pollard ρ за O(d^{1/4}) ≈ O(10^3) итераций. Итого 4·10^9 — слишком много.

Оптимизация. Используем тот факт, что d² = a² + b² — надо проверить, является ли √(a²+b²) целым. Это проще: вычислить d = (ll)sqrt(a²+b²) и проверить d*d == a²+b². Это O(1) на пару! Итого O(N²).

int count_integer_distances(vector<pair<ll,ll>>& pts) {
    int cnt = 0;
    int n = pts.size();
    for (int i = 0; i < n; i++) {
        for (int j = i+1; j < n; j++) {
            ll dx = pts[i].first - pts[j].first;
            ll dy = pts[i].second - pts[j].second;
            ll d2 = dx*dx + dy*dy;
            ll d = (ll)sqrt((double)d2);
            // Корректируем округление
            while (d * d < d2) d++;
            while (d * d > d2) d--;
            if (d * d == d2) cnt++;
        }
    }
    return cnt;
}

---

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

---

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

1. Путаница с единицами. В Z[i] числа 1, −1, i, −i — единицы. Два элемента называются ассоциатами, если один получается из другого умножением на единицу. При подсчёте представлений n = a² + b² пары (a, b) и (−a, b) считаются разными! Формула r₂(n) учитывает ВСЕ 4 знака и оба порядка.
2. Ошибка в алгоритме Корнакки-Смита. Цикл Евклида нужно остановить при r₁² ≤ p (не строго < p). При p = a² + 0² (когда p = k², что невозможно для простых, но возможно для составных) нужно обрабатывать отдельно.
3. Неверная проверка p ≡ 1 (mod 4). Помните: p = 2 — особый случай (2 = 1² + 1²). Не включать его в цикл.
4. Переполнение. При проверке aa + bb == p для p ≤ 10^15 нужно использовать long long или __int128.
5. Игнорирование p ≡ 3 (mod 4) в чётной степени. p^{2k} = (p^k)² + 0², это законное представление нуля как одного из слагаемых! Но в формуле r₂(n) они учтены через нейтральный вклад ≡ 1.

---

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

Гауссовы целые позволяют мощно решать задачи, которые внешне кажутся не связанными с алгеброй. Ключевой паттерн: когда видите выражение вида a² + b², сразу думайте о Z[i]. Норма N(a + bi) = a² + b² мультипликативна — это нередко позволяет разложить сложную задачу на независимые части. Ещё один трюк: если задача формулируется как "найти число решёток точек на окружности a² + b² = n", ответ r₂(n) вычисляется через факторизацию n за O(n^{1/4} log n) без перебора точек. Пример: для n = 10^18 с Pollard ρ факторизация занимает ~1 мс, и мы мгновенно получаем r₂(n) — по сравнению с наивным O(√n) = O(10^9) перебором.

---

Итог

Кольцо гауссовых целых Z[i] — мощный инструмент для задач, связанных с суммами квадратов и решётками. Ключевые результаты: простые p ≡ 1 (mod 4) раскладываются в Z[i] как p = (a+bi)(a−bi); простые p ≡ 3 (mod 4) остаются простыми; 2 = −i(1+i)². Теорема Ферма описывает, когда n = a² + b², а формула r₂(n) = 4(d₁(n)−d₃(n)) считает число таких представлений. Алгоритм Корнакки-Смита находит конкретное разложение p = a²+b² за O(log² p). Понимание Z[i] открывает путь к более глубоким результатам теории чисел.