Примитивные и грязные генераторы случайных чисел
ОГЛАВЛЕНИЕ
В стандарте ANSI-C имеется функция rand(), описанная в файле <stdlib.h> и выдающая равномерно распределенное число от 0 до RAND_MAX. С ней связана функция srand(), производящая начальную установку счетчика. Почти все подобные генераторы используют рекуррентную последовательность
I(n+1)=(a*I(n)+c)(mod m).
Число a называется мультипликатором, число c инкрементом, а число m - модулем.
Такой выбор может послужить серьезным ограничением на статистические свойства последовательности случайных чисел. Кроме того, в большинстве случаев RAND_MAX=35767, что значительно меньше, чем диапазон изменения целых чисел. В некоторых испытаниях теория рекомендует 106 - 109 случайных проб, но, пользуясь подобным счетчиком, можно получить не более RAND_MAX одинаковых случайных чисел, т.е. в 30 тыс.- 30 млн. раз меньше рекомендуемого.
Генератор ANSI-C был опубликован комиссией как "пример". Мы его тоже приводим, но как "не рекомендованный" для серьезных приложений.
/* (в модуле stdlib.h) */
#define RAND_MAX 32767
/* "пример" от комитета ANSI-C */
unsigned long next=1;
int rand(void) {
next=next*1103515245+12345;
return((unsigned int)(next/65536)%32768);
}
void srand(unsigned int seed) {
next=seed;
}
Мультипликатор и инкремент этого примера (который, скорее всего и поставляется со стандартной библиотекой C) не являются оптимальными. Об оптимальных константах для такого рода последовательностей ниже.
Быстрые и грязные генераторы случайных чисел, отличные от стандартных
Иногда требуется произвести не слишком изысканную последовательных случайных действительных или целых чисел, при этом код генерации случайного числа желательно держать "в строке" (для скорости) либо производить выбор для различного возможного числа значащих битов беззнакового целого числа. Ниже приводятся фрагменты программ, осуществляющие подобную генерацию для нормального распределения действительных чисел между 0 и 1 и для целых чисел произвольного диапазона.
Генерация случайного действительного числа, равномерно распределенного
от 0 до 1 и целого, равномерно распределенного от jlow до jhigh.
static unsigned long iran;
unsigned long rand_a,rand_c,rand_m;
float fran;
int jran;
...................
/* число с плавающей точкой fran от 0 до 1 */
iran=(iran*rand_a+rand_c)%rand_m;
fran=(float)iran/(float)rand_m;
..................
/* целое jran на промежутке от jlow до jhigh */
iran=(iran*rand_a+rand_c)%rand_m;
jran=jlow+((jhigh-jlow+1)*iran)/im;
Ниже приводятся оптимальные значения коэффициентов rand_a, rand_c, rand_m для различного значения числа значащих бит в беззнаковом целом.
bits | rand_m | rand_a | rand_c |
20 | 6075 | 106 | 1283 |
21 | 7875 | 211 | 1663 |
22 | 7875 | 421 | 1663 |
23 | 6075 | 1366 | 1283 |
6635 | 936 | 1399 | |
11979 | 430 | 2531 | |
24 | 14406 | 967 | 3041 |
29282 | 419 | 6173 | |
53125 | 171 | 11213 | |
25 | 12960 | 1741 | 2731 |
14000 | 1541 | 2957 | |
21870 | 1291 | 4621 | |
31104 | 625 | 6571 | |
139968 | 205 | 29573 | |
26 | 29282 | 1255 | 6173 |
81000 | 421 | 17117 | |
134456 | 281 | 28411 | |
27 | 86436 | 1093 | 18257 |
121500 | 1021 | 25673 | |
259200 | 421 | 54773 | |
28 | 117128 | 1277 | 24749 |
121500 | 2041 | 25673 | |
312500 | 741 | 66037 | |
29 | 145800 | 3661 | 30809 |
175000 | 2661 | 36979 | |
233280 | 1861 | 49297 | |
244944 | 1597 | 51749 | |
30 | 139968 | 3877 | 29573 |
214326 | 3613 | 45289 | |
714025 | 1366 | 150889 | |
31 | 134456 | 8121 | 28411 |
259200 | 7141 | 54773 | |
32 | 233280 | 9301 | 49297 |
714025 | 4096 | 150889 |
Самый быстрый генератор для 32-битового представления целых и действительных чисел
В большинстве случаев, число типа unsigned long имеет 32 бита. В этом случае для генерации числа в диапазоне 0 - 232-1 достаточно простого умножения на мультипликатор и сложения с инкрементом. Деление по модулю будет произведено автоматически при переполнении. Значения мультипликатора и инкремента для этого случая получены в исследованиях Кнута и Льюиса.
/* генерация целого числа от 0 до 0xFFFFFFFF. */
static unsigned long iran;
...........
// целое случайное число
iran=1664525L*iran+1013904223L;
Для реализации на этой основе очень быстрого генератора равномерного распределения действительных чисел от 0 до 1 важно, что числа с плавающей точкой одинарной точности (тип float) в большинстве случаев представлены также 32 битами. Кроме того, во многих случаях (включая SUN, ALPHA, Silicon Graphics и IBM PC) представление этого числа отвечает одному стандарту IEEE. Для машины VAX это не так. Грязным трюком, позволяющим избегнуть деления на действительное число, является маскировка экспоненты и дальнейшее вычитание из числа единицы. Необходимые коэффициенты в программе приведены для IEEE (по умолчанию) и для VAX.
static unsigned long iran;
unsigned long temp;
float fran;
..........
#if !defined(VAX)
static unsigned long jflone=0x3f800000;
static unsigned long jflmsk=0x007fffff;
#else
static unsigned long jflone=0x00004080;
static unsigned long jflmsk=0xffff007f;
#endif
..........
iran=1664525L*iran+1013904223L;
temp=jflone|(jflmsk&iran);
// случайное число с плавающей точкой
fran=(*(float *)&temp)-1.F;
С точки зрения программиста, занимающегося численными методами, а не железом, этот генератор является "особо грязным", поскольку реализация зависит от представления чисел в компьютере. Впрочем, этот генератор в подавляющем большинстве практических случаев работает верно, а главное - очень быстро.