Краткий материал по алгоритмам генерации псевдослучайных последовательностей чисел

Последовательности случайных чисел (СЧ) являются неотъемлемым инструментом решения многих математических задач, и в том числе задач имитационного моделирования. Чаще всего необходимы последовательности случайных чисел, которые равномерно и равновероятно распределены на некотором отрезке. Большинство природных процессов являются случайными, и, согласно теории математической статистики, они подчиняются различным законам распределения случайных величин. Для порождения случайных чисел, равномерно распределенных на отрезке [0..M], достаточно для некоторого случайного процесса, подчиняющегося закону равномерного распределения, ввести меру случайной величины. А затем последовательно проводить эксперименты, измерять значения случайной величины и после нормирования получать необходимую равномерную случайную последовательность чисел. Например, можно положить в корзину десять шаров, помеченных десятью десятичными цифрами. Затем, каждый раз перемешивать тщательно шары, случайно доставать один из них, записывать его цифру и класть шар обратно в корзину.

Подобные методы дадут очень хорошие статистические результаты, но потребуют колоссального времени для получения сколько-нибудь длинной последовательности. К счастью, математики разработали рекуррентные формулы получения псевдослучайных числовых последовательностей. Название «псевдослучайный» обусловлено хотя бы тем фактом, что если известно некоторое i-е число, то по формуле мы однозначно вычислим (i + 1)-ый элемент последовательности. Более того, из рекуррентной природы найденной формулы получается, что при одном и том же x0 мы получим при повторной генерации ту же самую последовательность. К тому же, все арифметические алгоритмы генерации псевдослучайных последовательностей периодичны, то есть существует некоторое pN, для которого выполняется xi = xi + pn при любом натуральном n. Поэтому любой генератор псевдослучайных последовательностей должен быть инициализирован случайной величиной, неким внешним источником случайных значений. Такими «аппаратными» генераторами в ПК могут выступать, например, электрические шумы в полупроводниковых устройствах, а также текущее количество выполненных процессором тактов или текущее значение времени.

Недостатки генераторов псевдослучайных чисел

Инициализация генератора псевдослучайной последовательности

Чтение текущего значения миллисекунд

Как уже было сказано выше, рекуррентную формулу вычисления i-го члена псевдослучайной последовательности при реализации на ПК можно инициализировать значением миллисекунд текущего времени в момент запуска программы. Для этого можно вызвать функцию «2Ch» прерывания 21h операционной системы DOS и получить «почти» случайное число на отрезке [0..99].

MOV   AH, 2Ch        ; номер функции получения времени
INT   21h            ; получаем время
MOV   HUNDREDS, DL   ; получаем сотые доли текущей секунды из DL
 

Чтение значения счетчика тактов процессора

Начиная с линейки процессоров Pentium в архитектуре x86 появилась инструкция, позволяющая прочитать счетчик тактов процессора с момента последнего сброса. Для иструкции rdtsc (Read Time Stamp Counter) задан машинный код 0F 31. 64-битное значение счетчика возвращается в паре регистров <EDX:EAX>. Если для инициализации генератора достаточно 32-битного значения, то необходимо использовать наиболее «чувствительные» младшие 32 бита счетчика тактов. При использовании старого компилятора языка ассемблера для обращения к 32-битному регистру EAX понадобится вручную проставить префикс с кодом 66h.

; Чтение значения младших 32 битов счетчика тактов процессора в пару <DX:AX>
DB    0Fh, 31h       ; машинный код инструкции RDTSC, результат в <EDX:EAX>
MOV   CL, 16         ; 32-битное значение EAX помещаем в пару <DX:AX>
DB    66h            ; = ROR EAX, CL - циклический сдвиг старшей части EAX в младшую
ROR   AX, CL
MOV   DX, AX         ; старшие 16 бит копируются в DX
DB    66h            ; = ROR EAX, CL - возвращаем младшие биты на место
ROL   AX, CL
 

Линейный конгруэнтный метод

Данный алгоритм был предложен Д. Х. Лемером в 1948 году. Рекуррентная формула выглядит следующим образом:

xn = (axn − 1 + cmod m.

Период порождаемой последовательности не превышает m. Естественно, что от выбора параметров acm и значения первого члена x0 существенно зависят основные свойства порождаемой последовательности. Длина периода будет максимальной (равной m) только в том случае, когда:

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

Алгоритм Блюма-Блюма-Шуба

В 1986 году трое авторов Ленор Блюм, Мануэль Блюм и Майкл Шуб предложили алгоритм генерации псевдослучайной последовательности, стойкий к обратным преобразованиям. Основная рекуррентная формула алгоритма:

xn = (xn − 1)2 mod pq,
zn = parity (xn),

где p и q — два больших простых числа. Для повышения качества получаемой последовательности на очередном шаге выбираются не все биты xn, а только младшие, или даже только бит четности. Из полученных «случайных битов» формируются двоичные псевдослучайные числа произвольной разрядности. Одной из особенностей вычислительной формулы является сквозная возможность вычислить xn без генерации всех предыдущих членов последовательности.

xn = (x0)2n mod (p − 1)(q − 1) mod pq,

Данный алгоритм более требователен к вычислительным ресурсам, но, с другой стороны, обладает хорошими статистическими характеристиками.

Алгоритм xor-shift

Профессором университета Флориды Джорджем Марсаглией был разработан очень быстрый генератор, который был назван «xor-shift». «Псевдокод» на Турбо Паскале выглядит следующим образом:

var
  x : LongWord = 3456789;
  y : LongWord = 62436069,
  z : LongWord = 1288629,
  w : LongWord = 675123;
 
function xor128: LongWord;
var
  t : LongWord;
begin
    t := (x xor (x shl 11));
    x := y;
    y := z;
    z := w;
    xor128 := w = (w xor (w shr 19)) xor (t xor (t shr 8)));
end;
 

Комбинация данного алгоритма с линейным конгруэнтным алгоритмом и алгоритмом Фибоначчи с запаздываниями представлена ниже.

program Sample;
var
    x : LongWord = 123456789;
    y : LongWord = 362436000;
    z : LongWord = 521288629;
    c : LongWord = 7654321;
 
function Random() : LongWord;
var
    t : int64;
begin
    x := int64(69069)*x + 12345;
    y := y xor (y shl 13);
    y := y xor (y shr 17);
    y := y xor (y shl 5);
    t := int64(698769069)*z + c;
    c := t shr 32;
    z := t;
    Random := x + y + z;
end;
 
begin
    {...}
    writeln ( Random() );
end.
 

Литература