v) Значение с не существенно, когда а - хороший множитель, за исключением того, что с не должно иметь общего множителя с т, когда т - размер компьютерного слова. Таким образом, можно выбрать с = 1 или с = а. Многие используют с = О вместе с m = 2, но они жертвуют двумя двоичными разрядами точности и половиной начальных значений, чтобы сэкономить всего несколько наносекунд счета (см. упр. 3.2.1.2-9).

vi) Младшие значащие цифры (справа) X не очень случайны, так что решения, основанные на числе X, всегда должны опираться, главным образом, на старшие значащие цифры. Обычно лучше считать X случайной дробью Х/т между О и 1, т. е. представлять себе X с десятичной точкой слева, а не относиться к X как к случайному целому числу между О и m - 1. Чтобы подсчитать случайное целое число между О и fc - 1, нужно умножить его на fc (но не делить на fc; см. упр. 3.4.1-3) и округлить результат.

vii) Важное ограничение случайности последовательности (1) обсуждалось в разделе 3.3.4, в котором показано, что "точность" при размерности t будет только порядка /m. Применяя метод Монте-Карло, необходимо использовать случайные последовательности высокой надежности. Их можно получить с помощью технических приемов, описанных в разделе 3.2.2.

viii) Можно генерировать не больше т/1000 чисел, иначе последующие будут вести себя подобно предыдущим. Если т = 2, значит, новая схема (например, новый множитель а) должна использоваться после генерирования нескольких миллионов случайных чисел.

Комментарии, приведенные выше, относятся, главным образом, к машинному языку программирования. Некоторые понятия прекрасно работают также на языках высокого уровня, например (1) превращается в Х=а*Х+с на языке С, если X имеет тип беззнаковое длинное и если т равно модулю беззнаково длинной арифметики (обычно или 2). Но С не дает возможности рассматривать X в качестве дроби, как того требует (vi), если не обращаться к числам с плавающей точкой двойной точности.

Поэтому для языка программирования, подобного С, часто используют другой вариант (1): выбирается простое число т, близкое к наибольшему легко вычисляемому целому числу, а полагается равным первообразному корню m и приращение с в этом случае равняется нулю. Тогда (1) может полностью выполняться простой арифметикой над числами, остающимися между -т и +т, так, как в упр. 3.2.1.1-9. Например, когда а = 48271 и m = 2 - 1 (см. строку 20 табл. 3.3.4-1), можно вычиашть А - аХ mod т на языке С.

#define ММ 2147483647 /♦ простое число Мерсена ♦/

#define АА 48271 /♦ здесь хорош спектральный критерий ♦/ #define QQ 44488 /♦ (long)(ММ/АА) ♦/ #define RR 3399 /♦ MM /. АА; важно, что RR<QQ ♦/

Х=АА* (X/.QQ) -RR* (long) (X/QQ) ; if (X<0) X+=MM;

Здесь X имеет тип long и X можно инициализировать, как ненулевое начальное значение, меньшее ММ. Поскольку ММ - простое число, наименьший значащий дво-

ичный разряд X так же случаен, как и самый старший, поэтому в предостережении (vi) нет необходимости.

Чтобы получить миллионы и миллионы случайных чисел, можно скомбинировать эту программу с другой, как в 3.3.4-(38), дописав дополнительно несколько операций.

#define МММ 2147483399 /♦ простое число, но не Мерсена ♦/

#define AAA 40692 /♦ другой удачный спектральный критерий */

#define QQQ 52774 /♦ (long)(МММ/AAA) ♦/

#define RRR 3791 /* МММ I, AAA; снова меньше, чем QQQ ♦/

Y=AAA*(Y/.QQQ)-RRR*(long) (Y/QQQ) ; if (Y<0) Y+=MMM; Z=X-Y; if (Z<=0) Z+=MM;

Подобно X, случайная величина Y должна быть установлена не равной нулю. Эти операции несколько отличаются от операций, привс 1,енных в 3.3.4-(38): выходное Z никогда не равно нулю и Z всегда лежит точно между О и 2. Длина периода последовательности Z приблизительно равна 74 квадриллионам, и ее числа сейчас имеют точность, приблизительно вдвое большую, че.м числа X.

Этот метод мобильный и совершенно простой, но не очень устойчивый. Альтернативная схема, основанная на последовательности Фибоначчи с запаздьтанием и вычитанием (упр. 3.2.2-23), даже более привлекательна, так как не только легко преобразуется для использования на разных компьютерах, но значительно быстрее работает и поставляет случайные числа лучшего качества, поскольку -мерная точность, возможно, хороша для t < 100. Ниже приведена программа на языке С гап.мггау(long аа[], int п), которая генерирует п новых случайных чисел и засылает их в заданный массив аа, используя рекуррентное соотношение

Xj = (Х, 1оо - Х,-з7) mod 2=*°. (2)

Это соотношение, в частности, хорошо подходит для современных компьютеров. Значение п должно быть равно по крайней мере 100; рекомендуются большие значения, например 1 ООО.

#define КК 100 /♦ длинное запаздывание ♦/

#define LL 37 /♦ короткое запаздывание ♦/

#define ММ (1L«30) /♦ модуль ♦/

#define mod.diff (х,у) ((.(х)-(у))&(ММ-1)) /♦ (х-у) mod ММ ♦/

long ran x[KK]; /♦ состояние генератора */

void ran array(long aa[],int n) { register int i,j;

for (j=0;j<KK;j++) aa[j]=ran x[j];

for (;j<n;j++) aa[j]=mod diff(aa[j-KK],aa[j-LL]);

for (i=0;i<LL;i++,j++) ran x[i]=mod diff(aa[j-KK],aa[j-LL]);

for (;i<KK;i++,j++) ran x[i]=mod diff(aa[j-KK],ran x[i-LL]);

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

пировать во время вычислений, чтобы позднее повторить запуск программы с такими же значениями, не проходя весь путь назад, к началу последовательности. Особенностью использования рекуррентных соотношений, подобных (2), является, конечно, необходимость получения сразу всех начальных значений путем присвоения подходящих значений Xq, Х99. Следующая программа ranstart {long seed) хорошо запускает генератор, когда задано любое начальное число, находящееся между О и 2° - 3 = 1,073,741,821 включительно.

#define ТТ 70 /♦ гарантирует разделение потоков */

#define is odd(x) ((х)&1) /♦ блок двоичных разрядов х */

#define evenize(x) ((х)&(ММ-2)) /* делаем х четным ♦/

void ran start(long seed) { /* используется для размещения гап array */ register int t,j;

long x[KK+KK-l]; /* подготовка буфера ♦/

register long ss=evenize(seed+2);

for (j=0;j<KK;j++) {

x[j]=ss; /* самозагрузка буфера */

ss«=l; if (ss>=MM) ss-=MM-2; /* Щ1клический сдвиг 29 двоичных

разрядов */

for (;j<KK+KK-l;j++) x[j]=0;

x[l]++; /♦ делаем x[l] (и только x[l]) нечетным */

ss=seed&(MM-l); t=TT-l; while (t) {

for (j=KK-l;j>0;j-) x[j+j]=xtj]; /* "квадрат" ♦/

for (j=KK+KK-2;j>KK-LL;j-=2) x[КК+КК-1-j]=evenize(x[j]); for (j=KK+KK-2;j>=KK;j-) if (is odd(x [j] )) { x[j-(KK-LL)]=mod diff(x[j-(KK-LL)],x[j]); x[j-KK]=mod diff (x[j-KK] ,x[j]);

if (is odd(ss)) { /* "умножаем на z" */

for (j=KK;j>0;j--) x[j]=x[j-l];

x[0]=x[KK]; /♦ сдвигаем буфер щжлично */

if (is odd(x[KK])) x[LL]=mod diff(x[LL],x[KK]);

if (ss) ss»=l; else t-;

for (j=0;j<LL;j++) ran x[j+KK-LL]=x[j] ; for (;j<KK;j++) ran x[j-LL]=x[j];

Довольно любопытное маневрирование ran.start приведено в упр. 9, в котором доказывается, что последовательности чисел, генерируемых с различными начальными значениями, независимы: каждый блок 100 последовательных значений Х„, Хп+1, Хп+99 В последующем выходном массиве ran.array будет отличаться от блоков, появляющихся с другим начальным значением. (Строго говоря, известно, что это справедливо только тогда, когда п < 2™, но меньше 2 не в год.) Некоторые процессы можно, следовательно, запускать параллельно с различными начальными