Джулиан Бакнелл - Фундаментальные алгоритмы и структуры данных в Delphi

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 и более будут учитываться в последней категории;