Генератор псевдослучайных чисел (a * seed + b) % m: можно ли распараллелить?

Данный очерк затрагивает проблему генерации последовательности псевдослучайных чисел вида x[n+1] = (a * x[n] + b) % m. "a % b" используется для обозначения остатка от целочисленного деления a на b. Этот способ генерирования предложен в конкурсе параллельного программирования Acceler8 2011 для заполнения исходной матрицы. На первый взгляд, такая генерация не может быть распараллелена, что вызывает опасения при масштабировании решения на случай нескольких ядер. Не затрагивая вопрос о целесообразности распараллеливания генерации матрицы в контексте задачи конкурса, здесь будет описан способ ускорения и распараллеливания выполнения этой задачи.

Описание проблемы


Итак, необходимо сгенерировать последовательность

x[n+1] = (a * x[n] + b) % m, где x[0] = (a * seed0 + b) % m.


Тут a, b, seed0 и m - параметры генерации. Рассмотрим здесь лишь случай НОД(a-1, m) = 1. Случай, когда НОД(a-1, m) делит b, может быть сведён к предыдущему, а все остальные дают вырожденную длину цикла. Последний случай тоже интересен, но его анализ будет предложен к самостоятельному анализу читателя, дабы не загромождать пост. Более подробно о различных параметрах генерации и аспектах практической реализации можно прочесть в работе [1]. Далее обсуждаемый алгоритм генерации будем называть алгоритмом Лемера (согласно [1] впервые он был предложен Д. Х. Лемером).

Цель данного поста - показать возможность параллелизации генерации, а также указать практические аспекты реализации.

Параллельная генерация


Дальнейшие соображения верны в случае НОД(a-1, m) = 1. Более того, все выкладки приведены для простого m. В случае составного m, m-1 необходимо заменить значением функции Эйлера от m - это также сделано для простоты.

Модификация задачи
Здесь будут использованы результаты из прикладной алгебры без формулировок и доказательств, однако источника [1] вполне достаточно, чтобы разобраться. Впрочем, автор надеется, что читатель сможет ограничиться нижеприведенными короткими рассуждениями.

Рассмотрим следующее элементарное алгебраическое преобразование:

y(n) = (x(n) + c) % m


При этом выберем такое c, чтобы выполнялось

y[n+1] = (a * y[n]) % m


Из алгоритма Лемера получим

(x[n+1] + c) % m = (a * x[n] + a * c) % m


откуда сразу получим выражение для с:

c = (inv(a - 1) * b) % m


Здесь inv(a - 1) - означает обратный элемент к a - 1 по модулю m: ((a - 1) * inv(a - 1)) % m = 1. Несложно показать, что

inv(a - 1) = ((a - 1)^(m - 1)) % m


Таким образом, c может быть найдено за log(m), поскольку inv(a - 1) может быть найден за log(m). Эти очевидные размышления приведены для того, чтобы показать, что подготовительные вычисления не сводят на нет предложенную вычислительную схему.

Наконец, окончательная формула для вычислений x[n] следующая:

y[n+1] = (y[n] * a) % m = (y[1] * a^n) % m


и

x[n] = (y[n] + (-c)) % m = (y[n] + m - c) % m


Тут (-c) = m - c. А c = ((a - 1)^(m - 1) * b) % m, как показано выше.

Параллелизация
Пусть мы хотим распараллелить вычисления на k процессоров. Вычислим m - c, а также y[1], y[2], ..., y[k] и ak = (a^k) % m.

Разделим всю задачу на k задач: i-тая задача заключается в вычислении x[i+k], x[i+2*k], ... . Для подсчёта каждого следующего элемента необходимо единственное умножение по модулю m: y[i+k] = (y[i] * ak) % m, и сложение x[i+k] = (y[i+k] + (m - c)) % m.

Сложность
Сложность каждой задачи порядка m/k. Заметьте, что длина цикла не более m, поэтому достаточно сгенерировать m первых элементов, а затем скопировать повторения до конца массива. При этом, среднее значение элементов массива может быть посчитано для всего массива через сумму элементов цикла (и тех, разумеется, что, возможно, идут перед его первым вхождением в массиве), и отнято от элементов цикла перед копированием (см. детали генерации в условиях конкурса). Такой подход позволяет ускорить процесс генерации даже без распараллеливания.

Практические аспекты


В предложенной параллельной вычислительной схеме тонким местом является обращение процессов к общей памяти для записи. Для решения этой проблемы можно писать каждым процессом в отдельную строку исходной матрицы. Такой порядок записи позволяет при использовании определённого количества задач не только быстро генерировать данные, но и расположить их в удобном для дальнейшей обработки порядке. Однако, эти размышления выходят за рамки обсуждения параллельного алгоритма генерации.

Целесообразность реализации описанной выше параллельной вычислительной схемы в контексте задачи конкурса оставляю на суд читателя.

Используемая литература


1. Кнут Д. Искусство программирования для ЭВМ: пер. с англ.– М.: Мир, 1977.– Т.2.
For more complete information about compiler optimizations, see our Optimization Notice.