1. Идея

Нужно найти все числа x из диапазона [0, n - 1], для которых выполняется сравнение:

a * x ≡ b (mod n).

Это стандартное линейное сравнение. Его удобно решать через расширенный алгоритм Евклида.

Основная идея такая:

• сначала найти g = gcd(a, n);
• если b не делится на g, то решений нет;
• если делится, то сравнение можно сократить на g и получить более простое: a1 * x ≡ b1 (mod n1), где a1 = a / g, b1 = b / g, n1 = n / g;
• после сокращения числа a1 и n1 взаимно просты, значит у a1 существует обратный элемент по модулю n1;
• находим одно решение x0, а затем выписываем все решения в диапазоне [0, n - 1].

2. Наблюдения

Наблюдение 1. Когда решений нет

Сравнение

a * x ≡ b (mod n)

эквивалентно тому, что число a * x - b делится на n.

Из теории линейных сравнений известно:

• решение существует тогда и только тогда, когда gcd(a, n) делит b.

Если это не так, ответ сразу равен 0.

Наблюдение 2. После деления на НОД задача упрощается

Пусть g = gcd(a, n), и b делится на g.

Тогда можно разделить всё на g:

• a1 = a / g
• b1 = b / g
• n1 = n / g

Получим сравнение:

a1 * x ≡ b1 (mod n1).

Теперь gcd(a1, n1) = 1, а значит, у a1 есть обратный элемент по модулю n1.

Наблюдение 3. Как получить одно решение

Если inv — обратное к a1 по модулю n1, то

x0 ≡ b1 * inv (mod n1).

Это даёт одно решение x0 в диапазоне [0, n1 - 1].

Наблюдение 4. Как получить все решения

Если одно решение по модулю n1 найдено, то все решения по модулю n имеют вид:

x = x0 + t * n1, где t = 0, 1, ..., g - 1.

Почему именно столько?

• шаг между решениями равен n1;
• диапазон [0, n - 1] имеет длину n = g * n1;
• значит, туда помещается ровно g различных решений.

3. Алгоритм

1. Считать a, b, n.
2. С помощью расширенного алгоритма Евклида найти g = gcd(a, n).
3. Если b % g != 0, то:
   • вывести 0;
   • вывести пустую строку.
4. Иначе:
   • вычислить a1 = a / g, b1 = b / g, n1 = n / g;
   • снова применить расширенный алгоритм Евклида к a1 и n1, чтобы получить коэффициент, из которого можно восстановить обратный элемент inv для a1 по модулю n1;
   • вычислить одно решение: x0 = (b1 * inv) % n1;
   • для всех t от 0 до g - 1 добавить в ответ: x0 + t * n1.
5. Вывести количество решений и сами решения.

4. Почему это работает

Пусть g = gcd(a, n).

Случай 1. b не делится на g

Тогда линейное сравнение

a * x ≡ b (mod n)

неразрешимо. Это классический факт: сравнение a * x ≡ b (mod n) имеет решение тогда и только тогда, когда gcd(a, n) делит b.

Значит, в этом случае правильный ответ — отсутствие решений.

Случай 2. b делится на g

Тогда после деления на g получаем:

a1 * x ≡ b1 (mod n1),

где gcd(a1, n1) = 1.

Раз a1 и n1 взаимно просты, существует число inv, такое что

a1 * inv ≡ 1 (mod n1).

Умножим обе части сравнения на inv:

x ≡ b1 * inv (mod n1).

Значит, все решения по модулю n1 имеют вид:

x ≡ x0 (mod n1),

где x0 = (b1 * inv) % n1.

То есть все целые решения:

x = x0 + t * n1.

Теперь нужно оставить только те, которые лежат в диапазоне [0, n - 1].

Так как n = g * n1, то при t = 0, 1, ..., g - 1 получаем ровно g различных чисел:

• x0
• x0 + n1
• x0 + 2*n1
• ...
• x0 + (g-1)*n1

Все они меньше n, и все различны. Других решений в диапазоне [0, n - 1] быть не может, потому что все решения сравнимы с x0 по модулю n1.

Значит, алгоритм находит все и только все подходящие значения x.

5. Сложность

Расширенный алгоритм Евклида работает за O(log n).

Дальше мы выводим все решения. Их ровно g, где g <= n.

Итоговая сложность:

• по времени: O(log n + g)
• по памяти: O(g) из-за хранения ответов

В худшем случае это O(n) по времени и памяти, что подходит при n <= 10^5.

6. Код на C++17

#include <iostream>
#include <vector>
#include <algorithm>

using namespace std;

long long extended_gcd(long long a, long long b, long long &x, long long &y) {
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    long long x1, y1;
    long long g = extended_gcd(b, a % b, x1, y1);
    x = y1;
    y = x1 - (a / b) * y1;
    return g;
}

int main() {
    long long a, b, n;
    cin >> a >> b >> n;

    long long x, y;
    long long g = extended_gcd(a, n, x, y);
    if (g < 0) g = -g;

    if (b % g != 0) {
        cout << 0 << "\n\n";
        return 0;
    }

    long long a1 = a / g;
    long long b1 = b / g;
    long long n1 = n / g;

    long long inv_x, inv_y;
    extended_gcd(a1, n1, inv_x, inv_y);
    long long inv = ((inv_x % n1) + n1) % n1;

    long long x0 = (b1 * inv) % n1;

    vector<long long> ans;
    for (long long t = 0; t < g; t++) {
        ans.push_back(x0 + t * n1);
    }

    cout << ans.size() << "\n";
    for (int i = 0; i < (int)ans.size(); i++) {
        if (i > 0) cout << ' ';
        cout << ans[i];
    }
    cout << "\n";

    return 0;
}

7. Код на Python 3

a, b, n = map(int, input().split())

def extended_gcd(x, y):
    if y == 0:
        return x, 1, 0
    g, x1, y1 = extended_gcd(y, x % y)
    return g, y1, x1 - (x // y) * y1

g, x, y = extended_gcd(a, n)
g = abs(g)

if b % g != 0:
    print(0)
    print()
else:
    a1 = a // g
    b1 = b // g
    n1 = n // g

    g2, inv_x, inv_y = extended_gcd(a1, n1)
    inv = inv_x % n1

    x0 = (b1 * inv) % n1

    ans = []
    for t in range(g):
        ans.append(str(x0 + t * n1))

    print(len(ans))
    print(" ".join(ans))