Числа Кармайкла в окне
Старая библиотека проверяет простоту тестом Ферма: берёт основание a и считает an−1 mod n. Её авторы собирают список чисел, которые этот тест обманывает при любом основании, взаимно простом с n, — чисел Кармайкла, — чтобы добавить их в чёрный список. Составное n — число Кармайкла, если an−1 ≡ 1 (mod n) для всех a, взаимно простых с n. Найдите все числа Кармайкла в нескольких окнах [L, R].
Проверять все основания не нужно: по признаку Корсельта составное n — число Кармайкла ровно тогда, когда оно не делится на квадрат простого числа и для каждого его простого делителя p число p − 1 делит n − 1. Например, 561 = 3 · 11 · 17, а 560 делится на 2, 10 и 16.
Вход
В первой строке число окон q (1 ≤ q ≤ 5). В каждой из следующих q строк — два целых L и R (1 ≤ L ≤ R ≤ 10¹², R − L < 50 000). Оба конца входят в окно.
Выход
Для каждого окна — одна строка: сначала количество чисел Кармайкла в окне, затем сами числа по возрастанию через пробел.
Примеры
Пример 1
1 1 3000
5 561 1105 1729 2465 2821
Ограничения
- Время · Python
- 6000 мс
- Время · Go
- 1500 мс
- Память
- 256 МБ
- Проверка
- 7 тестов
Разбор
Открывайте по ступеням: каждая следующая подсказывает больше. Можно сразу к решению — но задача, решённая с одной подсказкой, запоминается лучше.
1Подсказка 1
Раскладывать на множители каждое число окна перебором дорого: окно лежит у 10¹², корень — миллион. Каким быстрым тестом, который все числа Кармайкла заведомо проходят, можно отсеять почти всё окно сразу?
2Подсказка 2
Числа Кармайкла нечётные и проходят тест Ферма по основанию 2. После такого фильтра в окне из 50 000 чисел у 10¹² остаются около 1800 чисел — и почти все они простые. Как быстро и без ошибок отбросить простые?
3Подсказка 3
Следите за признаком целиком. Квадрат простого 1093 проходит тест Ферма по основанию 2, и 1092 делит 1093² − 1. Но это не число Кармайкла.
4Идея
Число Кармайкла нечётно (у чётного n основание n − 1 даёт (−1)n−1 = −1) и взаимно просто с двойкой, поэтому проходит тест Ферма по основанию 2. Значит, сначала фильтр: 2n−1 mod n = 1 — одно быстрое возведение в степень на число. Из выживших выбрасываем простые детерминированным тестом Миллера — Рабина: для n < 1,12·10¹² достаточно оснований 2, 13, 23 и 1662803. Остаются составные псевдопростые по основанию 2 — их единицы на всё окно, а чаще ни одного, — и только их раскладываем перебором до корня и проверяем признак Корсельта.
5Разбор
Перебираем нечётные n окна. Возводим 2 в степень d, где n − 1 = 2sd, и s раз в квадрат: последнее число цепочки — это 2n−1, и по нему работает фильтр Ферма, а вся цепочка сразу даёт и сильный тест по основанию 2. Эталон пользуется этим, чтобы для простого не возводить в степень лишний раз.
Если n проходит фильтр и сильный тест по основаниям 2, 13, 23, 1662803, оно простое. Иначе раскладываем его перебором: делители ищем до √n ≤ 10⁶, а делать это приходится редко. Число подходит, если в разложении нет повторяющихся простых, их не меньше трёх и каждое p − 1 делит n − 1.
Почему простых делителей хотя бы три? Если n = pq и p < q, то n − 1 = p(q − 1) + (p − 1), и из q − 1 | n − 1 следовало бы q − 1 | p − 1, а p − 1 меньше. Квадраты простых признак исключает отдельным условием.
5+Решение
6Сложность
На каждое число окна — одно возведение в степень, O(log n) умножений; на каждое простое — ещё три раунда. Разложение перебором, до 10⁶ делений, достаётся только псевдопростым по основанию 2, а их до 10¹² всего около ста тысяч на триллион чисел. Итого O((R − L) log R) на окно.
7Типичные ошибки
- Признак без «не делится на квадрат». 1093² = 1 194 649 и 3511² = 12 327 121 — псевдопростые по основанию 2, и p − 1 у них делит n − 1. Простые 1093 и 3511 — единственные известные простые Вифериха, у которых 2p−1 ≡ 1 даже по модулю p².
- Псевдопростые вместо чисел Кармайкла. 341 = 11 · 31 проходит тест Ферма по основанию 2, но 3340 ≢ 1 (mod 341).
- Сильный тест вместо теста Ферма в фильтре. Сильный тест по основанию 2 проходят лишь немногие числа Кармайкла — 43 из 646 до миллиарда, — и фильтр выбросит почти все нужные числа.
- Переполнение. В Go произведение двух остатков по модулю около 10¹² доходит до 10²⁴ и в 64 бита не помещается. Нужно 128-битное произведение —
bits.Mul64иbits.Div64, как в эталоне.