Генератор случайных чисел с большим периодом для параллельных программ

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

Сер. 10. 2010. Вып. 1

УДК 519.6:533.6.011

Ю. П. Галюк, В. П. Мемнонов

ГЕНЕРАТОР СЛУЧАЙНЫХ ЧИСЕЛ С БОЛЬШИМ ПЕРИОДОМ ДЛЯ ПАРАЛЛЕЛЬНЫХ ПРОГРАММ

1. Введение. Для численного моделирования различных задач в науке и технике широко применяются статистические методы Монте-Карло. Их важной составной частью является интенсивное использование последовательностей случайных чисел для моделирования различных статистических распределений и исходов событий. Такие последовательности вырабатываются специальными компьютерными программами, генераторами случайных чисел (ГСЧ). Они являются поэтому только некоторыми псевдослучайными приближениями к последовательностям настоящих случайных чисел. Степень данного приближения проверяется с помощью разных локальных и глобальных тестов [1], и хороший генератор должен выдерживать их с необходимостью. В то же время как достаточные в некотором практическом смысле можно рассматривать специальные проверки путем решения задач, допускающих независимую оценку результатов с помощью аналитических или численных методов и относящихся к области дальнейшего применения ГСЧ [2, 3]. Поэтому, например, сравнение результатов численного моделирования задачи с точным решением той же задачи дает оценку качества применявшегося генератора, которую можно распространить для практического использования и в некотором круге других подобных задач, так или иначе связанных с тестовой задачей.

C этой целью в работе методом прямого статистического моделирования (ПСМ) численно моделируется на параллельном кластере задача о релаксации вращательной энергии молекул из первоначального неравновесного состояния посредством межмоле-кулярных столкновений, для расчета которых применяются псевдослучайные подпоследовательности испытуемого генератора. В конечном, равновесном состоянии должно установиться распределение Больцмана с известной и одинаковой для всех степеней свободы температурой Tt. Именно по величине отклонений от этого теоретического значения с учетом ошибок самого метода ПСМ, в наибольшей степени связанных со статистическим рассеянием из-за неизбежной ограниченности величины используемых статистических выборок, и оценивалось качество тестировавшегося генератора. Для оценки степени независимости случайных величин, моделируемых с помощью

Галюк Юрий Парфенич — заведующий отделом сетевых серверных систем и информационновычислительных ресурсов Петродворцового телекоммуникационного центра Санкт-Петербургского государственного университета. Количество опубликованных работ: 19. Научные направления: информационные технологии, цифровая обработка сигналов. E-mail: galyuck@ptc.spbu.ru.

Мемнонов Владимир Павлович — кандидат физико-математических наук, старший научный сотрудник аэродинамической лаборатории математико-механического факультета Санкт-Петербургского государственного университета. Количество опубликованных работ: 44. Научные направления: численные методы Монте-Карло, параллельное программирование распределенных систем, асимптотические методы. E-mail: pokusa@star.math.spbu.ru.

© Ю. П. Галюк, В. П. Мемнонов, 2010

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

2. Генератор случайных чисел с большим периодом для применения в параллельных программах. В настоящей работе тестируется линейный мультипликативный генератор вида

Zi+l = A  Zi

причем для моделирования с помощью методов Монте-Карло последовательность Zi нормируется и превращается в последовательность U:

U =[Ui = Zi/M; i = 1, 2,...],

элементы которой представляют собой как бы результат моделирования из равномерного распределения на интервале [0,1). В работе [4] уже выполнена проверка подобного генератора Фишмана (F) с модулем M = 231 — 1 и 207 множителями Aj, прошедших всесторонние теоретические тесты в статье [1]. В результате были получены оценки различных множителей Aj и способы их применения в параллельных программах. Однако для очень больших и продолжительных расчетов с использованием многих сотен процессоров нужны генераторы с гораздо большими и модулем, и периодом. Одним из них является генератор Дядькина-Хэмильтона (DH) с модулем M = 2128 и множителем A = 5lool°9. Его период P очень большой, P = 2126. Этот генератор предложен ими в статье [5]. В ней также проведены над ним некоторые тесты, в которых он проявил себя с наилучшей стороны.

2.1. Реализация генератора DH. Для реализации операции умножения по модулю M с очень большими числами в генераторе DH удобно представлять каждый из соответствующих множителей U и V в виде двух частей, выделяя отдельно биты левее n = 128, которые не повлияют на результат, Ul, V1, и остальные справа Uo, Vo, которые как раз очень существенны. Тогда при перемножении U = 2nUl + Uo на V = 2nVl + Vo получаем выражение

UV = 22nUlVl + 2n

в котором первое слагаемое просто отбрасывается при этой операции, а в остальных после умножений в зависимости от разрядности машины, в свою очередь, могут потребоваться дополнительное разбиение на части и умножение их столбиками и перекрестно. На языках высокого уровня можно описать операнды длиной В байт (integer * 8 Фортран, long long С). Но все же арифметика по модулю 2128 может быть осуществлена только программным путем. В частности, если на 32-разрядной машине аппаратно реализуются только арифметические операции (умножение и сложение целых чисел) над 4-байтными операндами (integer * 4, long), то для реализации генератора DH на такой машине нужно из 4 целых чисел типа integer * 4 образовать целое число integer * 16, а при перемножении двух подобных больших целых чисел достаточно провести 10 операций умножения из 16 и сложить результаты (6 умножений можно не делать, так как они не войдут в сумму из-за взятия результата по модулю 2128). Но лучше, конечно, использовать целые числа типа integer * 8, если они имеются. Тогда достаточно двух

таких чисел для получения целого числа integer * 16, а при перемножении подобной пары уже нужно выполнить только 3 операции умножения (см. (2)). Все частичные произведения должны делаться как беззнаковые (старший бит является значимым, а не знаковым, на языке С аналогом служат операнды unsigned long). При суммировании длинных операндов обязательно надо учитывать возможность появления бита переполнения, который надо добавить к старшей части операнда. Проще всего обойти эти проблемы, используя программирование на Ассемблере, где существует беззнаковый вариант операции умножения (команда MUL). Более того, одна из форм такой команды позволяет получать неусеченный результат (двойная размерность от размерности операндов), чего нет в языках высокого уровня. Также в Ассемблере имеется форма операции сложения с учетом бита переноса (команда ADC). Такая программа на Ассемблере была создана Ю. П. Галюком и приводится в п. 4. Получаемый после компиляции файл rand128.o эквивалентен фортрановской подпрограмме с именем rand128, для чего к ее имени внутри ассемблеровского текста программы rand128.s добавлен символ подчеркивания, и после этого она может использоваться при вызове из программ на Фортране или С. Существенно также то, что ассемблеровский генератор rand128.s оказался в 2 раза более быстрым, чем генератор DH на Фортране, реализованный с помощью специальных приемов для преодоления вышеотмеченных трудностей в работе [6]. Отметим также, что совпадение с последовательностью случайных чисел, вырабатываемых реализацией генератора DH этой работы, послужила дополнительной проверкой ассемблеровского умножения в rand128.s. Теперь сравним среднее время rrn для выработки одного случайного числа с помощью рассматриваемой реализации генератора DH, а также генераторов (F) [1] и Mersenne twister (MT) [7], используя процессор Xeon 5335 с тактовой частотой 2 ГГц при прогонки для каждого из них последовательности из 500 млн случайных чисел:

Генераторы ... F DH MT

тгп  108, с ... 1.3637 0.8802 1.5875

Как видно, генератор DH оказался заметно более быстрым.

2.2. Разбиение на подпоследовательности и назначение их на процессоры. Для получения независимых реализаций на отдельных процессорах в параллельной программе необходимо иметь возможность распределять на них большие независимые подпоследовательности генератора. Для этого была создана специальная программа, в которой заранее с помощью умножения посредством генератора rand128.s осуществляются большие скачки по его значениям для перехода на такие подпоследовательности. Начальные значения данных подпоследовательностей в виде двух частей z(1), z(2) с целыми типа integer * 8 для представления целого числа типа integer * 16 запоминаются в некотором двумерном массиве nsmf() и по позиции In в нем с помощью своего значения myrank, которое предварительно должно быть передано в рабочие подпрограммы, каждый процессор может уникально их использовать. Действительно, для каждого процесса myrank получаем начальное значение его подпоследовательности по формуле

z(1) = nsmf (2 * myrank + 1 + 2 * In * numprocs),

z

где numprocs обозначает число запущенных процессов в параллельной MPI-задаче.

А с помощью изменения целого числа In достигается и более широкое применение

формулы (3), например при переходе на новые подпоследовательности, при распределенных вычислениях с одновременным запуском на разных кластерах отдельных MPI-задач. Кроме того, также с помощью умножения посредством rand128.s в этой программе получается множитель A = 5100109 целого типа integer * 16 в виде двух частей двумерного массива ml(1), ml(2) целых типа integer * 8, которые применяются уже во всех процессорах по формуле (1).

3. Проверка генератора DH с помощью численного моделирования точно решаемой задачи. Как уже отмечалось выше, сравнение результатов численного моделирования задачи с ее точным решением дает оценку качества применявшегося генератора, которую можно распространить для практического использования и в некотором круге других близких задач. В нашей точно решаемой задаче рассматриваются молекулы, равномерно распределенные в объеме с зеркально отражающими стенками и имеющие только поступательные и вращательные степени свободы. Причем в начальный момент создается неравновесное состояние системы, когда вся энергия сосредоточена только на поступательных степенях свободы. В конечном состоянии молекулы должны иметь распределение Больцмана для всех степеней свободы. Отклонения поступательной Ttr и вращательной Tr температур от известного теоретического значения в равновесии Tt могут существовать, но при идеальных генераторах это возможно только вследствие ограниченного размера статистической выборки N. Связанную с этим ошибку г^ метода Монте-Карло [8], например, для вращательной температуры Tr можно записать в виде

rN = x

где D(Tr) - дисперсия вращательной температуры; N - величина выборки. Параметр x принимался равным 3, соответствуя коэффициенту доверия 0.997. Поэтому, если в результате релаксации со статистической выборкой N получаются отклонения от Tt,

STr =| Tr - Tt \,

больше rN, то это означает, что испытуемый генератор с вероятностью 0.997 нехорош. Относительно взаимодействия частиц предполагалось, что молекулы сталкиваются как твердые сферы, а в неупругих вероятностях вращательных переходов используется значение 5 для параметра числа столкновений, что соответствует вероятности неупругого перехода 0.2. Расчетное поле разбивается на одномерные ячейки, в которых производятся процедуры столкновений метода ПСМ [2]. Временной шаг т был равным половине

среднего времени между столкновениями для одной молекулы, а вся расчетная область

содержала Nm молекул, причем во всех расчетах Nm = 100 001. После достижения равновесного состояния через промежутки времени At, равные

At = Ntf х т,

выполнялись еще Kt дополнительных выборок параметров системы. Причем целое число временных шагов между выборками было большим, чаще всего Ntf = 24. В работе [4] для линейного мультипликативного генератора с множителями из статьи [1] при использовании той же самой точно решаемой задачи было установлено, что этот выбор для данной задачи оптимален, обеспечивая отсутствие корреляций для хороших множителей. В итоге суммарный объем статистической выборки N получается равным

N = Nm х Kt х numprocs.

В начальном состоянии, как уже говорилось выше, вся энергия была сосредоточена только на поступательных степенях свободы молекул. Поступательная температура Т полагалась равной 500 К, вращательная Тг была равна нулю. Предполагалось, что молекулы обладают только двумя вращательными степенями свободы. В использованной модели вращательной релаксации Ларсена-Боргнаке (ЛБ) [9] строго выполняются законы сохранения энергии и импульса и условие детального баланса. Поэтому в процессе вращательной релаксации должно устанавливаться конечное равновесное состояние с одинаковыми значениями температуры Т4, равными 300 К по закону сохранения полной энергии для данного случая. Несмотря на то, что модель ЛБ только приближенно описывает во времени ход вращательной релаксации, конечное равновесное состояние получается точным. И сравнение результатов численного моделирования проводится именно по нему.

3.1. Результаты численного моделирования. Численное моделирование вращательной релаксации выполнялось на 40-процессорном кластере Санкт-Петербургского государственного университета. Однако практически можно было использовать одновременно не более 16 процессоров. Поэтому в приводимых ниже данных сначала накапливалась необходимая выборка шагами по 16, а затем она обрабатывалась специальной программой. В частности, именно таким путем были получены представленные на рис. 1 результаты для 208 процессоров, применяя только одну выборку К = 1 на каждом из них. Начальные значения подпоследовательностей для них выбирались вышеописанным способом (см. п. 2.2). Причем для оценки влияния выбора величины промежутка между соседними подпоследовательностями были выполнены расчеты, когда они следовали друг за другом через интервалы в 2КМ чисел при КМ = 29, 32, 64. Эти значения КМ отложены по оси абсцисс на рис. 1. При численном моделировании задачи с К = 1 процессоры использовали от 4, 95 • 108 до 5, 06 • 108 случайных чисел. Поскольку 229 = 536 870 912, то значение КМ = 29 выбрано как минимально возможный интервал между соседними подпоследовательностями для данной задачи. Тем не менее максимальное отклонение \5Т1Г,Г|тах (см. (5)) конечных поступательной Ти вращательной Тг температур от известного теоретического значения в равновесии Т = 300 К, представленное на рис. 1 кривой 2, для всех значений КМ гораздо меньше теоретической ошибки Монте-Карло гN (см. (4)), равной для рассматриваемого варианта ГN = 0.197 и изображенной прямой 1. Так что генератор DH может быть рекомендован для численных расчетов методами ПСМ в подобных задачах при Кг = 1.

В работе [4] было показано, что после суммирования на каждом процессоре конечных температур по всем Мт молекулам и деления результата на Мт получаются асимптотически нормально рапределенные случайные величины, например, для вращательных температур обозначим их Тг. Дело в том, что именно в случае двумерного нормального распределения обращение в нуль коэффициента корреляции эквивалентно независимости рассматриваемой паре случайных величин. Причем Мт = 105 по оценкам в [4] уже обеспечивает хорошее асимптотическое приближение к нормальности, и на рис. 1 в виде кривых 3 и 5 представлены выборочные коэффициенты корреляции между этими конечными температурами Тг - соответственно гааеу и гъе, полученные на соседних процессорах с нечетными и четными значениями тугапк и, наоборот, из разных половин их, последовательно начиная с максимально удаленных. В частности, входящие в определение этих выборочных коэффициентов корреляции

СОУ

го<1еу — ^ ^ ,

е

о.з г

о.о-

25 30 35 40 45 50 55 60 65

КМ

Рис. 1. Значения Q кривых: для ошибки гN = 0,197

выборочные ковариации СОУ (Х,У) записываются для г0^РЛ как

П2

СОУ(X, У) = (£)(ТГ(23 - 1) - Хау) • (Тг(2Я - Уау))/п2,

3=1

а для гье как

П2

СОУье (Х,У) = ^(Тг 3) - Хау ) • (Тг (п - 3 +1) - Уау ))/ш,

3 = 1

где п - число процессоров, п2 = п/2, индекс суммирования 3 для удобства полагался равным 3 = тугапк +1, а выборочные дисперсии и Ьу и выборочные средние Хау и Уау записываются в виде

ЬХ = ( Е(^ (23 - 1) - Хау )2 /(п2 - 1), = £(ТГ (23) - Уау )2 | /(п2 - 1),

31 ) \3=1

ХаУ = Тг (23 - 1)^ /п2, Уау = ^ Тг (23)^ /п2.

Истинный коэффициент корреляции рху между нормально распределенными случайными величинами Х и У можно оценить [10] с помощью г0^РЛ или гье. Для больших выборок среднее квадратичное отклонение аг коэффициента, например г0аеу (см. (8)), от рху, согласно [10], может быть представлено следующим образом:

о-г =

при этом г0аеУ приближенно следует нормальному закону с параметрами (рху,&г). И для процентного уровня значимости ц имеем неравенство

гойеу tq • &т < рХУ < гойеу + tq • &т,

где для ц = 1% и п2 > 100 по [10] получаем tq = 2.58. Таким образом, при выполнении неравенства

| гойеу > tq •

коэффициент корреляции рху отличен от нуля и между рассматриваемыми случайными величинами с вероятностью 0.99 есть корреляционная связь. На рис. 1 кривая 4 как раз и представляет слившиеся в одну линию граничные значения из (10), (11) для г0^РЛ и гье. Как видно, кривые 3 и 5 проходят намного ниже ее, поэтому неравенство (11) нигде не выполняется, коэффициент корреляции между Тг не отличим от нуля и, значит, эти случайные величины независимы. Порождающие их псевдослучайные подпоследовательности, видимо, также независимы.

2

Рис. 2. Значения Q кривых: для ошибки гN (1), максимального отклонения \5Т1г,г\тв.х (2), коэффициента корреляции \гое \ (3), его границы гЬое (4), коэффициента корреляции \гЬе\ (5) и его границы гЬЬе (6) при К = 14

В практических расчетах методом ПСМ часто для увеличения суммарной выборки (см. (7)) делают дополнительные выборки через временные интервалы Дt (см. (6)) в конечном состоянии. Результаты численного моделирования для такого случая с К = 14 иллюстрирует рис. 2. Здесь приведены результаты решений 13 МР1-задач по 16 процессов в каждой с К = 14, начальные значения независимых подпоследовательностей генератора DH для них выбирались с помощью переменных тугапк и 1п

в формуле (4). Монте-карловская ошибка гN была постоянной и равнялась 0.190. Она изображена прямой 1. Кривой 2 на рис. 2 представлено максимальное отклонение для конечных температур бТг = | Тг — Т | или бТг = | Т^ — Т |. Как видно, оно часто оказывается выше прямой 1 ошибки ГN, что указывает, возможно, на неадекватную работу генератора в таких условиях.

Поэтому такой способ увеличения суммарной выборки нужно применять с осторожностью, используя, например, предварительную проверку с помощью многократного расчета одного типичного варианта задачи при К = 1 для накопления статистической выборки N, равной выборкам будущих экономных прогонов с К >> 1. Собственно нужно сделать К независимых расчетов, образовать файл суммарных результатов и сравнить ошибки с соответствующим экономным прогоном, затем, считая расчет с К = 1 наиболее близким к точному, оценить возможность применения К >> 1.

О

Рис. 3. Значения Q кривых: для ошибки гN = 0,411 (1), максимального отклонения \5Т1г,г\тв.х (2), коэффициента корреляции \'Гог1£ъ\ (3), его границы гЪогЫъ (4 ) при 1п = 1, и при 1п = 8 для тех же величин: \5Т1г,г\-та.х (5), \г

odev \ (6), гЪ

odev

Коэффициенты | и для применявшегося четного числа процессоров

и при четном К оценивают соответственно попарные корреляции между последовательными по времени соседними промоделированными значениями Тг для каждого из 16 процессоров и, наоборот, последовательно удаленными парами на разных процессорах. Они изображаются кривыми 3 и 5. И, наконец, граничные значения неравенства (11) гЪ04еь и гЪъе изображаются кривыми 4 и 6. Эти величины определяют границу, при преодолении которой выборочными коэффициентами корреляции |г0dev | и \гъе | появляются уже отличные от нуля истинные коэффициенты корреляции. И теперь из рис. 2 следует, что соседние пары коррелированы, поскольку кривая 3 везде выше своей граничной кривой 4. Возможно, эти корреляции каким-то образом связаны и с вышеупомянутыми выбросами кривой 2 выше прямой 1 монте-карловской ошибки ГN. Но данный вопрос требует дополнительных исследований. В то же время кривая 5,

представляющая коэффициент корреляции \тье|, почти всюду меньше своей граничной кривой 6. Так что пары, содержащие случайные величины Tr с разных процессоров, не имеют корреляционной связи.

Как видно из рис. 2, для Kt = 14 при некоторых значениях параметра In и, следовательно, при соответствующих подпоследовательностях псевдослучайных чисел генератора DH численное моделирование рассматриваемой задачи сопровождается ошибкой, превосходящей теоретическую ошибку метода Монте-Карло tn . Интересно, а как будет меняться результат, если постепенно уменьшать величину Kt?

Некоторые результаты такого исследования представлены на рис. 3. Прямая 1 на рис. 3 соответствует поддерживаемой постоянной посредством изменения числа процессоров теоретической ошибке Монте-Карло tn = 0.411. Кривая 2 представляет максимальный модуль отклонения \STtr,r\max конечных поступательной и вращательных температур от теоретического значения Tt = 300 К для весьма плохой для Kt = 14, согласно рис. 2, подпоследовательности при In =1. Но при Kt = 2, 4 она проходит все-таки еще под контрольной прямой 1 с tn = 0.411, хотя при Kt = 8, 16 отклонения \STtr,r \max уже находятся выше этой прямой и соответственно больше теоретической ошибки. Кривая 3 представляет попарный коэффициент корреляции \Todev \, полученный по статистической выборке, состоящей из пар соседних значений температур Tr в пределах Kt на каждом процессоре, а кривая 4 - границу, при преодолении которой становится отличным от 0 истинный коэффициент корреляции, и между ними появляется корреляционная связь. Заметим, что это происходит в стабильном режиме приблизительно одновременно с выходом \STtr,r\max за теоретический предел tn = 0.411. На рис. 3 также показаны полученные при тех же условиях результаты для умеренно плохой для Kt = 14 подпоследовательности генератора с In = 8. В данной ситуации даже при Kt = 8, 16 отрезок 5 для отклонения \5Ttrr\max гораздо ниже прямой 1, правда, коэффициент корреляции \Todev \ (см. (8)) оказывается выше границы 7 в точке Kt = 8. Таким образом, при некотором снижении требований к точности расчета генератор может работать и при Kt >> 1. Но, конечно, необходима предварительная проверка с прогонками типичного варианта при Kt = 1, о которой говорилось выше.

4. Приложение. Приведем программу DH на Ассемблере для процессоров с архитектурой x86-64, созданную Ю. П. Галюком. Программа написана в нотации ATT для компилятора с языка Ассемблер: as (Linux/Unix), но может быть переработана для Windows (компиляторы FASM, ML64). Функция rand128 не имеет параметров; все необходимые для работы генератора переменные (множитель A и начальное случайное число Z из (1)) находятся в соттоп-блоке тт16 и доступны из вызывающей программы, что позволяет использовать программу, в частности, для реализации больших скачков при образовании независимых подпоследовательностей генератора.

program r128.s

.textliteral:.long 0x00000000,0x3c000000 .comm rr16_,32,8 .comm varrnd_,8,8 .global rand128_ rand128_:

// common /rr/z,m// integer*4 z(4),m(4) push % rcx

incq varrnd_

// 0x1

mov rr16_,%rcx push % rbx

mov rr16_+24,%rax mul %rcx mov %rax,%rbx // 1x0

mov rr16_+8,%rax mull rr16_+16

add %rax,%rbx // 0x0

mov rr16_+16,%rax mul %rcx // z

add %rbx,%rdx mov %rax,rr16_ mov %rdx,rr16_+8 // ppr16 shr %rdx

cvtsi2sdq %rdx,%xmm0

5. Заключение. Создана новая быстродействующая реализация генератора DH Дядькина и Хэмилтона на Ассемблере. Получен удобный способ разбиения на независимые подпоследовательности и назначения их на процессоры. Проведено численное моделирование точно решаемой задачи с помощью генератора DH. Когда на каждом процессоре производилась только одна выборка, ошибки были всегда меньше теоретического предела ошибки метода Монте-Карло. Были также оценены истинные коэффициенты корреляции между конечными температурами на разных процессорах, которые оказались не отличимыми от нуля и подтвердили независимость используемых подпоследовательностей. В то же время, если на каждом процессоре производятся дополнительные выборки, как это часто делают на практике для увеличения суммарной статистической выборки, то появляются ошибки, превосходящие возможные теоретические ошибки метода Монте-Карло. Численная величина их продемонстрирована в 13 разных параллельных расчетах для 16 процессоров с 14 подобными дополнительными выборками на каждом. В этом случае оценки истинных коэффициентов корреляции оказались уже отличными от нуля и выявили наличие корреляционной связи между конечными температурами, полученными в дополнительных выборках. Для оценки возможной полезности выполнения дополнительных выборок предложена специальная предварительная проверка с расчетом типичного варианта задачи двумя способами.

Литература

1. Fishman G. S., Moore L. R. An exhaustive analysis of multiplicative congruential random number generators with modulus 231 — 1 // SIAM J. Sci. Stat. Comput. 1986. Vol. 7, N 1. P. 24-45.

2. Bird G. A. Molecular gas dynamics and the direct simulation of gas flows. Oxford: Clarendon Press, 1994. 458 p.

3. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. М.: Наука, 1982. 296 с.

4. Мемнонов В. П. Тестирование генераторов случайных чисел с помощью численного моделирования точно решаемой задачи // Вестн. С.-Петерб. ун-та. Сер. 1: Математика, механика, астрономия. 2009. Вып. 1. С. 25-32.

5. Dyadkin I. G., Hamilton K. G. Study of 128-bit multiplyers for congruent pseudorandom number generators // Comput. Phys. Commun. 2000. Vol. 125. P. 239-258.

6. Marchenko M. A., Mikhailov G. A. Parallel realization of statistical simulation and random number generators // Russ. J. Numer. Anal. Math. Modelling. 2002. Vol. 17, N 1. P. 113-124.

7. Matsumoto M., Nishimura T. Mersenne Twister: A 623-dimensional equidistributed uniform pseudorandom number generator // ACM Trans. on Modeling and Computer Simulation. 1998. Vol. 8, N 1. P. 3-30.

8. Соболь И. М. Численные методы Монте-Карло. М.: Наука, 1973. 311 с.

9. Borgnakke C., Larsen P. S. Statistical collision model for Monte Carlo simulation of polyatomic gas mixture // J. Comput. Phys. 1975. Vol. 18. P. 405-420.

10. Смирнов Н. В., Дунин-Барковский И. В. Курс теории вероятностей и математической статистики. М.: Наука, 1969. 511 с.

Статья рекомендована к печати проф. Л. А. Петросяном.

Статья принята к печати 24 сентября 2009 г.