Прямые генератора
Линейный конгруэнтный генератор выдаёт числа xn+1 = (a·xn + c) mod m. Если отметить на квадрате все пары подряд (xn/m, xn+1/m), они лягут на параллельные прямые. Целые s1, s2, не оба нулевые, задают семейство прямых s1·u + s2·v = k + s2·c/m, k — любое целое, и это семейство накрывает все пары, если s1 + a·s2 делится на m. Расстояние между соседними прямыми семейства — 1/√(s12 + s22).
Спектральный тест ищет семейство с самыми редкими прямыми, то есть наименьшее s12 + s22: чем оно больше, тем гуще прямые и тем равномернее пары заполняют квадрат. Найдите его для каждого генератора. От c ответ не зависит.
Вход
В первой строке — число генераторов q (1 ≤ q ≤ 2·104). Дальше q строк «a m»: 2 ≤ m ≤ 231, 1 ≤ a < m.
Выход
Для каждого генератора — наименьшее значение s12 + s22 по всем целым s1, s2, не обоим нулевым, у которых s1 + a·s2 делится на m.
В примере — генераторы из статьи: RANDU, MINSTD с множителями 16807 и 48271 и три генератора из сцены с модулем 4096. При a = 2049 ответ 8: вектор (−2, 2), все пары лежат на четырёх прямых.
Примеры
Пример 1
6 65539 2147483648 16807 2147483647 48271 2147483647 5 4096 2049 4096 1697 4096
2147221514 282475250 1990735345 26 8 4562
Ограничения
- Время · Python
- 4000 мс
- Время · Go
- 1000 мс
- Память
- 256 МБ
- Проверка
- 6 тестов
Разбор
Открывайте по ступеням: каждая следующая подсказывает больше. Можно сразу к решению — но задача, решённая с одной подсказкой, запоминается лучше.
1Подсказка 1
Все пары (s1, s2) с s1 + a·s2 ≡ 0 (mod m) образуют решётку на плоскости. Какие два вектора её порождают? И как в решётке найти кратчайший ненулевой вектор, похоже на алгоритм Евклида?
2Идея
Любое решение записывается как s2·(−a, 1) + j·(m, 0): второе слагаемое добавляет к s1 кратное m. Значит, решётку порождают векторы (m, 0) и (−a, 1), и нужен её кратчайший ненулевой вектор.
На плоскости его находит редукция Лагранжа — Гаусса — двумерный алгоритм Евклида. Держим два вектора, порождающих решётку: u длиннее, v короче. Вычитаем из u кратное v, ближайшее к проекции: u ← u − k·v, где k — округлённое (u·v)/(v·v). Решётка от этого не меняется. Если новый u стал короче v, меняем их ролями и продолжаем; если нет — v и есть кратчайший. Почему: после остановки |u| ≥ |v| и |u·v| ≤ v·v/2, и тогда любая комбинация αu + βv с целыми α ≠ 0 не короче v, а при α = 0 это кратное v.
3Разбор
Начинаем с u = (m, 0), v = (−a, 1), переставив их, если u короче. Округление k = ⌊(2·u·v + v·v) / (2·v·v)⌋ — к ближайшему целому. Ответ — v·v в момент остановки.
Длины в квадрате и скалярные произведения не больше m2 ≤ 262 и помещаются в 64-битное знаковое целое. А вот 2·u·v уже может не поместиться, поэтому в Go округляют через остаток: k = ⌊d/vv⌋, и ещё единица, если остаток не меньше половины vv.
3+Решение
4Сложность
Как у алгоритма Евклида: длины падают в геометрической прогрессии, шагов O(log m) на генератор, всего O(q log m). Больше всего шагов — при a, близком к m·(√5 − 1)/2, как у Евклида на соседних числах Фибоначчи.
5Типичные ошибки
- Перебирать s2 до небольшой границы. Кратчайший вектор бывает длиной до (4/3)1/4·√m, и при m около 231 s2 доходит до десятков тысяч.
- Остановиться, когда u короче, но не переставить. Шаг должен всегда вычитать короткий вектор из длинного, иначе редукция застрянет или выдаст не кратчайший.
- Округлять отбрасыванием дробной части. Отбрасывание к нулю вместо ближайшего целого оставляет u длиннее, чем можно, и цикл кончается раньше времени.