1. Идея

Нужно решить сравнение a * x ≡ b (mod n) и вывести все решения x из диапазона 0..n-1.

Ключевой факт из теории чисел:

• пусть g = gcd(a, n);
• если b не делится на g, то решений нет;
• если b делится на g, то решений ровно g.

После этого сравнение можно упростить, разделив всё на g:

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

Тогда получаем сравнение

a1 * x ≡ b1 (mod n1),

причём gcd(a1, n1) = 1. Значит, у числа a1 существует обратный элемент по модулю n1, и можно найти одно базовое решение. Из него уже строятся все остальные.

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

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

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

Известный факт: такое сравнение разрешимо тогда и только тогда, когда gcd(a, n) делит b.

Поэтому:

• если b % g != 0, ответ — 0 решений.

Наблюдение 2. После деления на НОД коэффициенты становятся взаимно простыми

Если g = gcd(a, n), то после деления:

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

получаем gcd(a1, n1) = 1.

Это важно, потому что тогда у a1 есть обратный элемент по модулю n1.

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

Если gcd(a1, n1) = 1, то существует число inv, такое что

a1 * inv ≡ 1 (mod n1).

Тогда одно решение:

x0 ≡ b1 * inv (mod n1).

Обратный элемент находится через расширенный алгоритм Евклида.

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

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

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

Именно таких решений будет ровно g, и все они различны в диапазоне 0..n-1.

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;
   • с помощью расширенного алгоритма Евклида найти коэффициент x, такой что
     a1 * x + n1 * y = 1;
   • тогда inv = x mod n1 — обратный элемент к a1 по модулю n1;
   • найти базовое решение: x0 = (b1 * inv) % n1;
   • для всех t от 0 до g - 1 добавить в ответ число x0 + t * n1.
5. Вывести количество решений и сами решения.

Так как x0 < n1, а далее мы прибавляем n1, ответы сразу идут по возрастанию.

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

Пусть дано сравнение a * x ≡ b (mod n).

Это означает, что существует целое k, для которого

a * x - b = k * n.

То есть нужно решить линейное диофантово уравнение. Для таких уравнений известно:

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

Обозначим g = gcd(a, n).

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

Тогда уравнение неразрешимо, а значит и сравнение не имеет решений.

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

Разделим всё на g. Получим:

a1 * x ≡ b1 (mod n1),

где a1 = a/g, b1 = b/g, n1 = n/g.

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

a1 * inv ≡ 1 (mod n1).

Умножим сравнение на inv:

x ≡ b1 * inv (mod n1).

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

x = x0 + t * n1, где x0 = (b1 * inv) % n1.

Почему нужно брать t именно от 0 до g-1?

Потому что диапазон всех ответов — 0..n-1, а n = g * n1.
Числа

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

все лежат в диапазоне 0..n-1 и различны. Следующее число, x0 + g*n1, уже отличается на n, то есть совпадает с первым по модулю n.

Следовательно, мы получаем все и только все решения, причём их ровно g.

5. Сложность

• Нахождение gcd занимает O(log n).
• Расширенный алгоритм Евклида тоже работает за O(log n).
• Построение всех решений занимает O(g), где g = gcd(a, n).

Итого:

• время: O(log n + g);
• память: O(g) для хранения ответов.

С учётом ограничения n <= 10^5 это работает быстро.

6. Код на C++17

#include <iostream>
#include <vector>
#include <numeric>

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 g = gcd(a, n);

    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 x, y;
    extended_gcd(a1, n1, x, y);

    long long 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(a, b):
    if b == 0:
        return a, 1, 0
    g, x1, y1 = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return g, x, y

g = __import__("math").gcd(a, n)

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

    _, x, y = extended_gcd(a1, n1)
    inv = x % n1
    x0 = (b1 * inv) % n1

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

    print(g)
    print(" ".join(ans))