Текст книги "Фундаментальные алгоритмы и структуры данных в Delphi"

Автор книги: Джулиан Бакнелл

сообщить о нарушении

Текущая страница: 15 (всего у книги 36 страниц)

В случае если значение с равно 0, генератор называется мультипликативным линейным конгруэнтным генератором случайных чисел (multiplicative linear congruential generator). Чтобы гарантировать, что цикл повторения последовательности максимален, необходимо в качестве значения параметра m выбирать простое число. Самым известным генератором подобного рода является так называемый минимальный стандартный генератор случайных чисел (minimal standard random number generator), предложенный Стивеном Парком (Stephen Park) и Кейтом Миллером (Keith Miller) в 1988 году. Для него а = 16807, а m = 2147483647 (или 2(^31^) – 1). После разработки этого генератора было проведено большое количество статистических тестов, и генератор прошел большинство из них (несмотря на то что предложенный генератор обладает определенными нежелательными свойствами, которые мы рассмотрим чуть ниже).

Мультипликативные линейные конгруэнтные генераторы случайных чисел имеют одну аномалию: они никогда не дают числа 0. (Это объясняется тем, что, во-первых, m представляет собой простое число, во-вторых, a mod m не равно нулю, и, в-третьих, если начальное число не равно нулю, Х(_0_) mod m тоже не равно нулю.) Следовательно, если генераторы никогда не дают числа 0, их нельзя назвать случайными. На практике невозможность генерации нуля, как правило, игнорируется, – в конце концов, в 32-разрядной операционной системе это всего лишь отсутствие всего одного числа из примерно 2 миллиардов.

При реализации минимального стандартного генератора случайных чисел (как, в общем-то, и любого другого) особое внимание необходимо уделить исключению возможности возникновения переполнения, поскольку значение текущего начального числа, умноженное на а, может легко превысить максимально допустимое значение для 32-битного целого числа. Если не позаботиться об исключении переполнения, возможно возникновение ошибок, которые негативно скажутся на достаточно хорошем генераторе случайных чисел. Для обработки случаев переполнения используется метод Шрейга (Schrage) (его описание в этой книге не приводится, но его можно найти в статье Парка и Миллера [16]).

Для сравнения и тестирования различных генераторов случайных чисел будет создана иерархия классов, базовый класс которой будет содержать виртуальный метод, инкапсулирующий основные функциональные возможности генератора, в частности, генерация случайного числа с плавающей запятой в диапазоне от 0 до 1 (мы будем пользоваться переменными типа double). Этот виртуальный метод будет перекрываться в дочерних классах, что позволит генерировать случайное число в соответствии с алгоритмами дочерних классов. В базовом классе метод будет применяться для создания других типов случайных чисел, например, случайных чисел целого типа не больше определенного значения или случайного числа из определенного диапазона.

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

Листинг 6.2. Базовый класс генератора случайных чисел

type

TtdBasePRNG = class private

FName : TtdNameString;

protected procedure bError(aErrorCode : integer;

const aMethodName : TtdNameString);

public

function AsDouble : double; virtual;

abstract;

{вернуть случайное число из диапазона от 0 включительно до 1 исключительно}

function AsLimitedDouble(aLower, aUpper : double): double;

{-вернуть случайное число из диапазона от aLower включительно до aUpper исключительно}

function AsInteger(aUpper : integer): integer;

{-вернуть случайное число из диапазона от 0 включительно до aUpper исключительно}

property Name : TtdNameString read FName write FName;

end;

function TtdBasePRNG.AsLimitedDouble(aLower, aUpper : double): double;

begin

if (aLower < 0.0) or (aUpper < 0.0) or (aLower >= aUpper) then

bError(tdeRandRangeError, 'AsLimitedDouble');

Result := (AsDouble * (aUpper – aLower)) + aLower;

end;

function TtdBasePRNG.AsInteger(aUpper : integer): integer;

begin

if (aUpper <= 0) then

bError(tdeRandRangeError, 'AsInteger');

Result := Trunc(AsDouble * aUpper);

end;

procedure TtdBasePRNG.bError(aErrorCode : integer;

const aMethodName : TtdNameString);

begin

raise EtdRandGenException.Create(

FmtLoadStr(aErrorCode,

[UnitName, ClassName, aMethodName, Name]));

end;

В листинге 6.2 приведен код базового класса генератора случайных чисел. В нем определен виртуальный метод AsDouble, который возвращает случайное число X в диапазоне 0< х< 1. Кроме того, в классе объявлены два простых метода, один из которых возвращает случайное число с плавающей запятой из заданного диапазона значений, а второй – из диапазона значений от 0 до некоторой заданной верхней границы (аналогично тому, как функция Random (Limit) использует целое значение Limit). Теперь, когда базовый класс определен, для реализации алгоритма Парка и Миллера можно объявить дочерний класс.

Листинг 6.3. Минимальный стандартный генератор псевдослучайных чисел

type

TtdMinStandardPRNG = class(TtdBasePRNG) private

FSeed : longint;

protected

procedure msSetSeed(aValue : longint);

public

constructor Create(aSeed : longint);

function AsDouble : double; override;

property Seed : longint read FSeed write msSetSeed;

end;

constructor TtdMinStandardPRNG.Create(aSeed : longint);

begin

inherited Create;

Seed := aSeed;

end;

function TtdMinStandardPRNG.AsDouble : double;

const

a = 16807;

m = 2147483647;

q = 127773; {равно m diva}

r = 2836; {равно m mod a}

OneOverM : double = 1.0 / 2147483647.0;

var

k : longint;

begin

k := FSeed div q;

FSeed := (a * (FSeed – (k * q))) – (k * r);

if (FSeed <= 0) then

inc( FSeed, m);

Result := FSeed * OneOverM;

end;

function GetTimeAsLong : longint;

{$IFDEF Delphi1}

assembler;

asm

mov ah, $2С

call DOS3Call

mov ax, cx end;

{$ENDIF}

{$IFDEF Delph2Plus}

begin

Result := longint(GetTickCount);

end;

{$ENDIF}

{$IFDEF KylixlPlus}

var

T : TTime_t;

begin

_time(@T);

Result := longint(T);

end;

{$ENDIF}

procedure TtdMinStandardPRNG.msSetSeed(aValue : longint);

const

m = 2147483647;

begin

if (aValue > 0) then

FSeed := aValue

else

FSeed := GetTimeAsLong;

{убедиться, что значение начального числа находится в переделах от 0 до m-1 включительно}

if (FSeed >=m-1) then

FSeed := FSeed – (m – 1) + 1;

end;

Как несложно заметить в коде метода AsDouble, метод Шрейга выглядит гораздо сложнее, нежели простая формула X(_n+1_) = aX(_n_) mod m со значениями а = 16807 и m = 2(^31^) – 1. Тем не менее, используя достаточно сложные математические выкладки, можно доказать его равенство приведенной формуле.

Кроме того, как уже упоминалось, в генераторе случайных чисел подобного типа использование нуля в качестве начального числа нежелательно, поскольку тогда бы все генерируемые значения были бы нулевыми. Поэтому метод msSetSeed использует значение 0 в качестве флага при необходимости установки начального числа по значению системных часов. К сожалению, для выполнения этой операции в 16– и 32-разрядных системах Windows используется разный код.

Создадим класс случайных чисел, который будет использовать системный генератор случайных чисел – функцию Random. В листинге 6.4 показан код метода AsDouble для такого класса.

Листинг 6.4. Использование в классе системной функции Random

function TtdSystemPRNG.AsDouble : double;

var

OldSeed : longint;

begin

OldSeed := System.RandSeed;

System.RandSeed := Seed;

Result := System.Random;

Seed := System.RandSeed;

System.RandSeed := OldSeed;

end;

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

Тестирование

В основе всех тестов будут лежать одни и те же принципы. Мы будем генерировать большое количество случайных чисел из диапазона от 0.0 (включительно) до 1.0 (исключительно). Получаемые в результате работы генераторов значения будут разбиваться на несколько категорий, будет подсчитываться количество значений в каждой категории, а затем вероятность попадания значения в каждую категорию. На основе результатов вычислений будет определяться значение функции хи-квадрат, на основе которого будет прогоняться тест по критерию хи-квадрат. При этом количество степеней свободы будет на единицу меньше, чем количество категорий значений. Это было всего лишь краткое введение, но через несколько минут мы приступим к собственно тестированию.

Тест на однородность

Первый тест самый простой – проверка на однородность. О нем мы уже говорили. Фактически случайные числа будут проверяться на равномерность распределения по диапазону от 0.0 до 1.0. Разобьем весь диапазон на 100 поддиапазонов, сформируем набор из 1000000 случайных чисел и вычислим количество значений, попавших в каждый поддиапазон. В поддиапазоне 0 будут находиться значения от 0.00 до 0.01, в поддиапазоне 1 – значения от 0.01 до 0.02 и т.д. Вероятность попадания случайного числа в любой поддиапазон составляет 0.01. Для полученного распределения вычислим значение параметра хи-квадрат и сравним его с данными для стандартного распределения хи-квадрат, находящимися в строке, для 99 степеней свободы.

Листинг 6.5. Тест на однородность

procedure UnifomityTest(RandGen : TtdBasePRNG;

var ChiSquare : double; var DegsFreedo : integer);

var

BucketNumber, i : integer;

Expected, ChiSqVal : double;

Bucket : array [0..pred(Uniformitylntervals) ] of integer;

begin

{вычислить количество чисел в каждом поддиапазоне}

FillChar(Bucket, sizeof(Bucket), 0);

for i := 0 to pred(UniformityCount) do

begin

BucketNumber := trunc(RandGen.AsDouble * Uniformitylntervals);

inc (Bucket [BucketNumber]);

end;

{вычислить значение параметра xu-квадрат}

Expected := UniformityCount / Uniformitylntervals;

ChiSqVal := 0.0;

for i := 0 to pred(Uniformitylntervals) do

ChiSqVal := ChiSqVal + (Sqr (Expected – Bucket [i]) / Expected);

{вернуть значения}

ChiSquare := ChiSqVal;

DegsFreedom := pred(Uniformitylntervals);

end;

Тест на пропуски

Второй тест, который мы проведем, – тест на пропуски – несколько сложнее первого. Тест на пропуски гарантирует, что последовательность случайных чисел не будет попадать сначала в один поддиапазон, а затем в другой, третий и т.д., несмотря на то, что в целом значения будут распределены равномерно по всему диапазону. Определим в диапазоне поддиапазон, скажем, первую половину – от 0.0 до 0.5. Сформируем набор случайных чисел. Для каждого генерируемого числа будем проверять, попадает ли оно в выбранный поддиапазон (попадание) или нет (промах). В результате проверок будет получена последовательность попаданий и промахов. Найдите последовательности из одного и большего количества промахов (такие последовательности называются пропусками, отсюда и название теста – тест на пропуски). Вы получите последовательности из одного, двух и даже большего количества промахов. Разбейте длины пропусков на категории. Если известно, что вероятность попадания равна p (в нашем случае она будет равна длине выбранного поддиапазона), то вероятность промаха будет (1 -p). На основе этих данных можно определить вероятность возникновения пропуска из одного промаха – (1 -p)p, двух промахов – (1 -p)(^2^)p, n промахов – (1 -p)(^n^)p, а, следовательно, вычислить ожидаемое количество пропусков любой длины. После этого применим тест по критерию хи-квадрат. Будем использовать 10 категорий пропусков (поскольку вероятность возникновения пропусков длиной 11 и более промахов очень мала, все пропуски длиной 10 и более будут учитываться в последней категории;

при этом, конечно, следует учитывать реальную вероятность попадания длины пропуска в эту последнюю категорию), следовательно, мы получим девять степеней свободы. Как правило, тест на пропуски проводится пять раз: для первой и второй половины диапазона, а также для первой, второй и третьей третей диапазона.

Листинг 6.6. Тест на пропуски

procedure GapTest(RandGen : TtdBasePRNG;

Lower, Upper : double;

var ChiSquare : double;

var DegsFreedom : integer);

var

NumGaps : integer;

GapLen : integer;

i : integer;

p : double;

Expected : double;

ChiSqVal : double;

R : double;

Bucket : array [0..pred(GapBucketCount) ] of integer;

begin

{вычислить длины пропусков и определить количество пропусков в каждой категории}

FillChar(Bucket, sizeof(Bucket), 0);

GapLen := 0;

NumGaps := 0;

while (NumGaps < GapsCount) do

begin

R := RandGen.AsDouble;

if (Lower <= R) and (R < Upper) then begin

if (GapLen >= GapBucketCount) then

GapLen := pred(GapBucketCount);

inc(Bucket[GapLen]);

inc(NumGaps);

GapLen := 0;

end else

if (GapLen < GapBucketCount) then

inc(GapLen);

end;

p := Upper – Lower;

ChiSqVal := 0.0;

{обработать все категории, кроме последней}

for i := 0 to GapBucketCount-2 do

begin

Expected := p * IntPower(1-p, i) * NumGaps;

ChiSqVal := ChiSqVal + (Sqr (Expected – Bucket [i]) / Expected);

end;

{обработать последнюю категорию}

i := pred(GapBucketCount);

Expected IntPower (1-p, i) * NumGaps;

ChiSqVal := ChiSqVal + (Sqr (Expected – Bucket [i]) / Expected);

{вернуть значения}

ChiSquare := ChiSqVal;

DegsFreedom := pred(GapBucketCount);

end;

Тест «покер»

Третий тест известен под названием "покер" (poker test). Случайные числа группируются в наборы по пять, а затем преобразуются в "карты", которые представляют собой цифры от 0 до 9. После этого определяется количество разных карт в каждом наборе (оно будет равно от одного до пяти) и полученные результаты разбиваются на категории. Поскольку вероятность пятикратного повторения одной и той же карты достаточно низка, случай выпадения только одной карты, как правило, включается в категорию "две разные цифры". К полученным четырем категориям применятся тест по критерию хи-квадрат (три степени свободы). Вероятность возникновения события для каждой категории вычислить не так уж легко (к тому же математические выкладки основаны на использовании комбинаторных значений, называемых числами Стерлинга), поэтому вычисления в этой книге не приводятся. Если вам интересно, то подробное описание можно найти в [11].

Листинг 6.7. Тест "покер"

procedure PokerTest(RandGen : TtdBasePRNG;

var ChiSquare : double;

var DegsFreedom : integer);

var

i, j, jlBucketNumber, NumFives : integer;

Accum, Divisor, Expected, ChiSqVal : double;

Bucket : array [0..4] of integer;

Flag : array [0..9] of boolean;

p : array [0..4] of double;

begin

{подготовительные операции}

FillChar(Bucket, sizeof(Bucket), 0);

NumFives PokerCount div 5;

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

Accum := 1.0;

Divisor := IntPower(10.0, 5);

for i := 0 to 4 do

begin

Accum := Accum * (10.0 – i);

p[i] := Accum * Stirling(5, succ(i)) / Divisor;

end;

{для каждой группы из пяти случайных чисел преобразовать все значения и числа от 1 до 10, определить количество разных цифр}

for i := 1 to NumFives do

begin

FillChar(Flag, sizeof(Flag), 0);

for j := 1 to 5 do begin

Flag [trunc(RandGen.AsDouble * 10.0)] :=true;

end;

BucketNumber := -1;

for j := 0 to 9 do

if Flag[j] then

inc(BucketNumber);

inc(Bucket[BucketNumber]);

end;

{объединить две первые категории – это будет сумма категорий "все цифры одинаковы" и "две разные цифры"}

inc(Bucket[1], Bucket[0]);

Expected := (p[0]+p[1]) * NumFives;

ChiSqVal := Sqr(Expected – Bucket[1]) / Expected;

{обработать другие категории}

for i := 2 to 4 do

begin

Expected :=p[i] * NumFives;

ChiSqVal := ChiSqVal + (Sqr (Expected – Bucket [i]) / Expected);

end;

{вернуть значения}

ChiSquare := ChiSqVal;

DegsFreedom := 3;

end;

Тест «сбор купонов»

Четвертый тест, который мы будем проводить, называется "сбор купонов" (coupon collector's test). Случайные числа считываются по одному и преобразуются в "купоны" – числа от 0 до 4. Фиксируется длина последовательности до получения полного комплекта купонов (т.е. цифр от 0 до 4). Очевидно, что получаемые длины будут от пяти и выше. После набора полного комплекта сбор купонов начинается снова. Длины последовательностей разбиваются на категории, к которым затем применяется тест по критерию хи-квадрат. Как правило, используются категории для длин от 5 до 19 и еще одна дополнительная категория для больших длин. Таким образом, мы получаем 16 категорий, а, следовательно, 15 степеней свободы. Как и в тесте "покер", вычисление вероятностей для каждой из категорий включает сложные математические выкладки, которые в этой книге не приводятся. Соответствующие подробности можно найти в [11].

Листинг 6.8. Тест "сбор купонов"

procedure CouponCollectorsTest(RandGen : TtdBasePRNG;

var ChiSquare : double;

var DegsFreedom : integer);

var

NumSeqs, LenSeq, NumVals, NewVal, i : integer;

Expected, ChiSqVal : double;

Bucket : array [5..20] of integer;

Occurs : array [0..4] of boolean;

p : array [5..20] of double;

begin

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

p[20] := 1.0;

for i := 5 to 19 do

begin

p[i] := (120.0 * Stirling(i-1, 4)) / IntPower(5.0, i);

p[20] := p[20] – p[i];

end;

NumSeqs := 0;

FillChar(Bucket, sizeof(Bucket), 0);

while (NumSeqs < CouponCount) do

begin

{продолжать сбор купонов (т.е. случайных чисел) до получения полного набора из пяти купонов}

LenSeq := 0;

NumVals := 0;

FillChar (Occurs, sizeof(Occurs), 0);

repeat

inc(LenSeq);

NewVal := trune(RandGen.AsDouble * 5);

if not Occurs [NewVal] then begin

Occurs[NewVal] := true;

inc(NumVals);

end;

until (NumVals = 5);

{обновить значение для соответствующей категории в зависимости от количества собранных купонов}

if (LenSeq > 20) then

LenSeq := 20;

inc(Bucket[LenSeq]);

inc (NumSeqs);

end;

{вычислить значение xu-квадрат}

ChiSqVal := 0.0;

for i := 5 to 20 do

begin

Expected := p [ i ] * NumSeqs;

ChiSqVal := ChiSqVal + (Sqr(Expected – Bucket [i]) / Expected);

end;

{вернуть значения}

ChiSquare := ChiSqVal;

DegsFreedom := 15;

end;

Результаты выполнения тестов

В разделе сопровождающих эту книгу материалов, который расположен на Web-сайте издательства, можно найти тестовую программу, которая применяет все рассмотренные нами тесты к стандартному генератору случайных чисел Delphi и минимальному стандартному генератору случайных чисел. На рис. 6.1 приведены результаты проведения одного из тестов для генератора случайных чисел Delphi.

Рисунок. 6.1. Тестирование стандартного генератора Delphi

Как видите, данный конкретный тест свидетельствует о том, что стандартный генератор Delphi успешно прошел проверку. (под успешным прохождением теста понимается, что предложенная к тестированию последовательность случайных чисел не дает результатов, которые являются значимыми на уровне 5%.)

Окно в правой части окна программы представляет собой снимок случайных чисел, полученных на выходе генератора. Координаты точек вычисляются на основе двух случайных чисел: первого для оси х, а второго – для оси Y. После этого точки выводятся в окне, если они находятся в прямоугольнике (0.0, 0.0, 0.001, 1.0) или, другими словами, в прямоугольнике, нижний левый угол которого расположен в точке (0.0, 0.0), а верхний правый – в точке (0.001, 1.0). Чтобы распределение точек было удобнее изучать, прямоугольник растянут по оси х. Как видите, на рисунке точки случайным образом разбросаны по всему прямоугольнику. Никакой системы при этом не наблюдается.

У вас может возникнуть удивление, зачем мы так много говорим об этом окне. Посмотрите на рис. 6.2, на котором показана та же программа, но для минимального стандартного генератора случайных чисел. Как видите, и этот генератор успешно прошел все тесты, но посмотрите на распределение случайных точек. Очевидно, что генератор дает последовательность случайных чисел, которые при переносе их на график формируют определенный регулярный рисунок.

Регулярность минимального стандартного генератора не позволяет использовать его для некоторых приложений, особенно тех, которые требуют пар случайных чисел. Даже незначительной регулярности бывает достаточно для того, чтобы приложение давало неверные результаты. Кроме того, отсутствие регулярности в результатах стандартного генератора случайных чисел Delphi в двухмерной плоскости не означает, что регулярности не будет в гиперплоскостях более высокой размерности. Существуют тесты, которые проверяют случайные числа на наличие регулярности в А> мерном пространстве, но давайте не будем погружаться в изучение слишком сложных тестов, а рассмотрим методы использования двух уже известных нам генераторов для дальнейшей рандомизации их выходных данных.

Рисунок 6.2. Тестирование минимального стандартного генератора

Мы рассмотрим три метода: первый известен как комбинаторный, второй – аддитивный и третий – метод тасования.

Комбинирование генераторов

Комбинирование генераторов заключается в параллельном использовании двух (или большего количества) мультипликативных линейных конгруэнтных генераторов с различными длинами циклов. Случайные числа генерируются обоими генераторами, а затем вычисляется их разность. Если получено отрицательное число, необходимо сделать его положительным, сложив его с длиной цикла первого генератора.

Листинг 6.9. Комбинирование генераторов type

TtdCombinedPRNG = class (TtdBasePRNG) private

FSeed1 : longint;

FSeed2 : longint;

protected

procedure cpSetSeed1(aValue : longint);

procedure cpSetSeed2(aValue : longint);

public

constructor Create(aSeed1, aSeed2 : longint);

function AsDouble : double; override;

property Seed1 : longint read FSeed1 write cpSetSeed1;

property Seed2 : longint read FSeed2 write cpSetSeed2;

end;

constructor TtdCombinedPRNG.Create(aSeed1, aSeed2 begin

inherited Create;

Seed1 := aSeed1;

Seed2 := aSeed2;

end;

longint);

function TtdCombinedPRNG.AsDouble : double;

const

al = 40014;

m1 = 2147483563;

ql = 53668;

{равно m1 div al}

rl = 12211;

{равно m1 mod al}

a2 = 40692;

m2 = 2147483399;

q2 = 52774;

{равно m2 div a2}

r2 = 3791;

{равно m2 mod a2}

OneOverMl : double = 1.0 / 2147483563.0;

var k : longint;

Z : longint;

begin

{получить случайное число с помощью первого генератора}

k := FSeed1 div ql;

FSeed1 := (al * (FSeed1 – (k * ql))) – (k * rl);

if (FSeed1 <= 0) then

inc(FSeed1, m1);

{получить случайное число с помощью второго генератора}

k := FSeed2 divq2;

FSeed2 := (a2 * (FSeed2 – (k * q2))) – (k * r2);

if (FSeed2 <= 0) then

inc(FSeed2, m2);

{объединить два случайных числа}

Z := FSeed1 – FSeed2;

if (Z <= 0) then

Z := Z + m1 – 1;

Result := Z * OneOverMl;

end;

procedure TtdCombinedPRNG.cpSetSeed1(aValue : longint);

const

m1 = 2147483563;

begin

if (aValue > 0) then

FSeed1 := aValue

else

FSeed1 := GetTimeAsLong;

{убедиться, что случайное число находится в диапазоне от 1 до m-1 включительно}

if (FSeed1 > – m1-1) then

FSeed1 := FSeed1 – (m1-1) + 1;

end;

procedure TtdCombinedPRNG.cpSetSeed2(aValue : longint);

const

m2 = 2147483399;

begin

if (aValue > 0) then

FSeed2 := aValue else

FSeed2 := GetTimeAsLong;

{убедиться, что случайное число находится в диапазоне от 1 до m-1 включительно}

if (FSeed2 >=m2-1) then

FSeed2 := FSeed2 – (m2 – 1) + 1;

end;

Как видите, код метода AsDouble в листинге 6.9 содержит два мультипликативных линейных конгруэнтных генератора: первый с параметрами {а, m} = {40014,2147483563}

и второй с параметрами {а, m} = {40692, 2147483399}.

Циклы обоих генераторов отличаются, но, тем не менее, близки к 2(^31^). Для преобразования промежуточного значения типа longint в значение типа double используется генератор с более длинным циклом.

Приведенный в листинге 6.9 генератор исключает двухмерную регулярность простого мультипликативного линейного конгруэнтного генератора, в чем можно убедиться с помощью программы тестирования. Можно показать, что длина цикла полученного комбинированного генератора составляет примерно 2 * 10(^18^). (Для сравнения, длина цикла стандартного генератора Delphi примерно равна 4 * 10(^9^).) Последовательность, вычисляемая с помощью комбинированного генератора полностью, определяется двумя начальными числами – по одному для каждого внутреннего генератора, в то время как для простого мультипликативного генератора было достаточно одного числа.

Аддитивные генераторы

Второй стандартный метод получения "более случайных" чисел от простого генератора называется аддитивным.

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

Листинг 6.10. Аддитивный генератор

type

TtdAdditiveGenerator = class (TtdBasePRNG) private

FInx1 : integer;

FInx2 : integer;

FPRNG : TtdMinStandardPRNG;

FTable : array [0..54] of double;

protected

procedure agSetSeed(aValue : longint);

procedure agInitTable;

public

constructor Create(aSeed : longint);

destructor Destroy; override

function AsDouble : double; override

property Seed : longint write agSetSeed;

end;

constructor TtdAdditiveGenerator.Create(aSeed : longint);

begin

inherited Create;

FPRNG := TtdMinStandardPRNG.Create(aSeed);

agInitTable;

FInx1 := 54;

FInx2 := 23;

end;

destructor TtdAdditiveGenerator.Destroy;

begin

FPRNG.Free

inherited Destroy;

end;

procedure TtdAdditiveGenerator.agSetSeed(aValue : longint);

begin

FPRNG.Seed := aValue;

agInitTable;

end;

procedure TtdAdditiveGenerator.agInitTable;

var

i : integer;

begin

for i := 54 downto 0 do

FTable[i] := FPRNG.AsDouble;

end;

function TtdAdditiveGenerator.AsDouble : double;

begin

Result := FTable[FInx1] + FTable[FInx2];

if (Result >= 1.0) then

Result := Result – 1.0;

FTable[FInx1] := Result;

inc(FInx1);

if (FInx1 >= 55) then

FInx1 := 0;

inc(FInx2);

if (FInx2 >= 55) then

FInx2 := 0;

end;

Если внимательно изучить код, показанный в листинге 6.10, можно обратить внимание, что для формирования массива, используемого при работе аддитивного генератора, применяется минимальный стандартный генератор случайных чисел. Несмотря на то что мы не можем определить «начальное число» для аддитивного генератора (фактически по истечении некоторого времени начальное число эквивалентно всему массиву;

внутренний генератор псевдослучайных чисел вызывается только 55 раз), мы можем его установить. При установке начального значения вызывается внутренний генератор, который заполняет массив, предназначенный для инициализации аддитивного генератора.

Длина массива, 55, и значения индексов, 54 и 23, – это не просто взятые наугад значения. Было показано, что они дают хорошие последовательности случайных чисел при генерации целых значений. (В книге [11] можно найти таблицы других удачных значений длины массива и индексов.)

Самым хорошим свойством данного генератора является длина цикла. Она просто огромна (при реализации на основе значений типа longint длина цикла будет составлять 230(255– 1), или приблизительно 3 * 1025). Даже если бы вы генерировали каждую секунду триллион случайных чисел, то для того, чтобы пройти весь цикл, потребовались бы годы.

Тасующие генераторы

И последний тип рассматриваемых нами генераторов, позволяющих получать "более случайные" числа, принадлежит к алгоритмам тасования. Здесь мы опишем генератор, реализованный на основе одного внутреннего генератора, хотя существуют и другие генераторы, аналогичным образом использующие два внутренних генератора.

Как и для аддитивного генератора, на первом этапе создается массив случайных чисел с плавающей запятой. Количество элементов в массиве не имеет особого значения. Кнут (Knuth) предложил использовать длины порядка 100. В нашем примере будет использоваться массив из 97 элементов – простое число, близкое к 100 [11]. (Кстати, применение простого числа не обязательно, оно просто выбрано в качестве примера.) Заполним массив случайными числами, полученными с помощью минимального стандартного генератора случайных чисел. Введем новую вспомогательную переменную и установим ее значение равным следующему случайному числу в последовательности.

При необходимости генерации следующего случайного числа с помощью тасующего генератора, вспомогательная переменная используется для вычисления случайного числа из диапазона от 0 до 96. Устанавливаем значение вспомогательной переменной равным значению элемента с вычисленным индексом и заменяем элемент новым случайным числом, полученным от внутреннего генератора случайных чисел. В качестве результата тасующего генератора используется значение вспомогательной переменной.

Листинг 6.11. Тасующий генератор

type

TtdShuffleGenerator = class(TtdBasePRNG) private

FAux : double;

FPRNG : TtdMinStandardPRNG;

FTable : array [0..96] of double;

protected

procedure sgSetSeed(aValue : longint);

procedure sgInitTable;

public

constructor Create(aSeed : longint);

destructor Destroy; override;

function AsDouble : double; override;

property Seed : longint write sgSetSeed;

end;

constructor TtdShuffleGenerator.Create(aSeed : longint);

begin

inherited Create;

FPRNG := TtdMinStandardPRNG.Create(aSeed);

sgInitTable;

end;

destructor TtdShuffleGenerator.Destroy;

begin

FPRNG.Free;

inherited Destroy;

end;

function TtdShuffleGenerator.AsDouble : double;

var

Inx : integer;

begin

Inx := Trunc(FAux * 97.0);

Result := FTable[Inx];

FAux := Result;

FTable[Inx] := FPRNG.AsDouble;

end;

procedure TtdShuffleGenerator.sgSetSeed(aValue : longint);

begin

FPRNG.Seed := aValue;

sgInitTable;

end;

procedure TtdShuffleGenerator.sgInitTable;

var

i : integer;

begin

for i := 96 downto 0 do

FTable[i] := FPRNG.AsDouble;

FAux := FPRNG.AsDouble;

end;

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