Здесь я буду переводить на Сишный код алгоритмы из Кнута (том второй, глава 3 - случайные числа) и всех остальных самых популярных источников.

Итак, Кнут.
Алгоритм K мы пропускаем ввиду его полной бесполезности. Алгоритм можно реализовать через switch.

Код аддитивного генератора - Алгоритм А, взятый из NCBI TK:

Генератор основан на формуле вида:

X_n = (X_n-24 + X_n-55) mod m, n > = 55;

Можно заменить сложение на умножение как предложено Марсалья.

#define DIM(A) (sizeof(A)/sizeof((A)[0]))

/*
	Additive random number generator

	Modelled after "Algorithm A" in
	Knuth, D. E. (1981). The art of computer programming, volume 2, page 27.

	7/26/90 WRG
*/

static long	state[33] = {
	(long)0xd53f1852,  (long)0xdfc78b83,  (long)0x4f256096,  (long)0xe643df7,
	(long)0x82c359bf,  (long)0xc7794dfa,  (long)0xd5e9ffaa,  (long)0x2c8cb64a,
	(long)0x2f07b334,  (long)0xad5a7eb5,  (long)0x96dc0cde,  (long)0x6fc24589,
	(long)0xa5853646,  (long)0xe71576e2,  (long)0xdae30df,  (long)0xb09ce711,
	(long)0x5e56ef87,  (long)0x4b4b0082,  (long)0x6f4f340e,  (long)0xc5bb17e8,
	(long)0xd788d765,  (long)0x67498087,  (long)0x9d7aba26,  (long)0x261351d4,
	(long)0x411ee7ea,  (long)0x393a263,  (long)0x2c5a5835,  (long)0xc115fcd8,
	(long)0x25e9132c,  (long)0xd0c6e906,  (long)0xc2bc5b2d,  (long)0x6c065c98,
	(long)0x6e37bd55 };

#define r_off	12

static long	*rJ = &state[r_off],
			*rK = &state[DIM(state)-1];

NLM_EXTERN void LIBCALL Nlm_RandomSeed(long x)
{
	register size_t i;

	state[0] = x;
	/* linear congruential initializer */
	for (i=1; i<DIM(state); ++i) {
		state[i] = 1103515245*state[i-1] + 12345;
	}

	rJ = &state[r_off];
	rK = &state[DIM(state)-1];

	for (i=0; i<10*DIM(state); ++i)
		(void) Nlm_RandomNum();
}


/*
	Nlm_RandomNum --  return value in the range 0 <= x <= 2**31 - 1
*/
NLM_EXTERN long LIBCALL Nlm_RandomNum(void)
{
	register long	r;

	r = *rK;
	r += *rJ--;
	*rK-- = r;

	if (rK < state)
		rK = &state[DIM(state)-1];
	else
		if (rJ < state)
			rJ = &state[DIM(state)-1];
	return (r>>1)&0x7fffffff; /* discard the least-random bit */
}

Значит так, этот код я просто спер из TK и особо его не анализировал. Теперь займемся кодингом с самого начала, читая второе издание.
Алгоритм не предполагает цикла, предусматривается, что массив будет заполнен заранее, следовательно, нужна процедура инициализации и процедура возвращения нужного числа.
Заполнять можно динамически или заранее подготовить таблицу, как, например, это и сделано в NCBI TK. Однако в Кнуте говорится о "... который имеет длину периода 2⁵⁵-1 ... последовательность <X_n> должна иметь период по крайней мере такой же длины", а в TK таблица отведена только под 33 элемента, что уже является ошибкой как мне кажется!
Замечания по шагу A1.
"Запишем Y[k] <- (Y[k] + Y[j])mod 2^e"
mod 2^e тут наибольший возможный делитель - длина компьютерного слова. Обычно в качестве значения берется RAND_MAX + 1 из CRT (если последовательность заполняется функцией rand(), т.к. функция возвращает число в диапазоне от 0 до RAND_MAX), т.е. объяснения сводятся к коду вида:

int Y[55];	
int m = RAND_MAX + 1;
Y[k] = (Y[j] + Y[k]) % m;

volodya
А разве на асме этот алгоритм выглядит не лучше?

rnd_index dd 0
rnd_data rd 64
rand_maxp1 dd rand_max+1

rnd_init:
        mov ecx,64
    .l0:
        rdtsc
        mov [rnd_data+ecx*4-4],eax
        loop .l0
        ret

rnd:
        mov ecx,[rnd_index]
        add cl,40
        mov eax,[rnd_data+ecx]
        add cl,124
        add eax,[rnd_data+ecx]
        add cl,96
        xor edx,edx
        div [rand_maxp1]
        mov [rnd_data+ecx],edx
        mov [rnd_index],ecx
        ret 

Я не чувствую себя достаточно крутым для того, чтобы проанализировать твой вариант. Очевидно, ты инициализируешь последовательность динамически через rdtsc, но насколько это правильно? Очевидно, int64 в качестве seed для RNG - это очень неплохо и начальное значение eax:edx предсказать практически невозможно как на мой ламерский взгляд. Но вот последующие вызовы rdtsc будут вплотную ложиться друг к другу... Иначе говоря, последовательность будет строго инкрементированной и все ее члены будут отличаться друг от друга лишь на какое-то вполне постоянное число. Насколько это плохо - я сказать пока не могу, не хватает знаний. :(

Кстати, идея rdtsc в качестве начального значения (seed) используется, например в:
http://www.irisa.fr/caps/projects/hipsor/HAVEGE1.0.html

И Рэндалл Хайд в своем HLA пишет следующее:

"Because of the nature of the RDTSC instruction, you should not call rand.randomize frequently or you will compromise the quality of the random numbers (indeed, it's generally not a good idea to "randomize" a random number generator more than once per program invocation). Similarly, you should avoid calling this function from a fixed script after power-on since that may also degrade the quality of the randomization. (These two suggestions are only important to those who are extremely concerned about the quality of the randomness of the generated numbers)."

Вам, вероятно, доводилось слышать, что Intel встроила в PIII какой-то RNG. Так вот - технически эта фраза некорректна. RNG встроен в ЧИПСЕТ, а не в процессор.

Ссылка:
http://home.comcast.net/~andrex/hardware-RNG/index.html

Практическое использование можно найти в кодах линукса - файл:
linux-2.6.6\drivers\char\hw_random.c

Алгоритм М - рандомизация перемешиванием

Алгоритм М/Алгоритм В мне помог разобрать Аркадий, за что я ему жутко признателен. Он не поленился напечатать 4 страницы текста, поясняя мне те моменты, которые я не догнал в Кнуте.

#include <stdlib.h>
#include <stdio.h>

#define arraysize(a) (sizeof (a) / sizeof *(a))

/* пример реализации генераторов Xi и Yi */
unsigned rand_X () { return rand (); }
unsigned rand_Y () { return rand (); }

/* поскольку rand() генерит числа в диапазоне 0..RAND_MAX,
   соответственно, делителем будет RAND_MAX+1. */
/* ВАЖНО: нужно быть уверенным, что RAND_MAX меньше максимального
   целого числа, а значит, RAND_MAX+1 не "обмодулится" в 0. Для
   (16-битного) BC++ 3.1 RAND_MAX=0x7FFF, а это меньше 0xFFFF. */

#define Y_m_div (RAND_MAX+1)

/* инициализация вспомогательного массива */

unsigned mix_V [50]; /* k==50 */

void init_mix_rand () {
        for (unsigned i = 0; i < arraysize (mix_V); i++)
                mix_V [i] = rand_X ();
}

/* собственно генератор */

unsigned mix_rand () 
{
        unsigned k = 50;
		unsigned j = (unsigned)((unsigned long) rand_Y () * k / Y_m_div);
        unsigned ret = mix_V [j];
        mix_V [j] = rand_X ();
        return ret;
}

/* использование */

void main () {

int m = 0;
long n = 0;

        init_mix_rand ();
        printf ("random numbers: #1 = %u", mix_rand ());
        printf (", #2 = %u", mix_rand ());
        printf (", #3 = %u\n", mix_rand ());
}

Комментарии.
1) Массив V - вспомогательный статический массив, заполняющийся k значениями последовательности <X_n>
2) Наибольшие проблемы для меня представлял шаг M2, где было непонятно, что, собственно, Кнут хотел сказать.
"Присвоим j <- floor(kY/m), где m - модуль, используемый в последовательности <Y_n>".
floor - это операция взятия "пола" - наибольшего целого, которое меньше или равно данному числу. В CRT такая функция определена для работы с FPU-числами, здесь не используется, я взял ее просто потому, что математический значок взятия "пола" :) напечатать на клаве не могу.
Проблемы тут две:
1) Непонятно, что делать с модулем m
2) Непонятно, что делать при переполнении разрядной сетки и потере старших бит.

Разберем по порядку.
Т.к. мы решили использовать rand-функции CRT, а rand возвращает число в диапазоне от 0 до RAND_MAX (0x7fff), то модуль взят как RAND_MAX+1. Кнут ограничил максимальное число в последовательностях X_n и Y_n как m-1. Если бы мы выбрали заполнение последовательностей X и Y НЕ через rand(), то пришлось бы брать другой модуль.
Теперь вопрос о переполнении. Дело в том, что данный генератор ОЧЕНЬ архитектурно-зависим. Ассемблерщикам (т.е. всем нам тут) наверняка известно понятие "перенос" и "переполнение". На сайте была статья и народ долго спорил касательно многих положений. Так вот, вполне возможно возникновение "переполнения" и отброса старших значащих бит в операции:

unsigned j = (unsigned)((unsigned long) rand_Y () * k / Y_m_div);

что есть невероятно паскудно, т.к. Кнут утверждает, что именно перемешивание старших бит и является наиболее важным.
Отрывок из письма Аркадия ко мне:
"Другой вопрос, почему для индекса предложены именно старшие, а не
младшие разряды. Полгаю, ответ в ответе на упражнение 4. Также, выше, при
обсуждении формулы 7' (стр. 49 в новом издании), говорится "в то время как
старшие двоичные разряды более тщательно перемешаны", так что, видимо,
предполагается, что старшие биты Y "более случайны", чем младшие.
"
В нашем конкретном данном случае, а именно Win32 под x86, переполнения не произойдет и старшие значащие биты отброшены не будут, однако это надо тщательно проверять при переносе на другую систему, иначе генератор потеряет смысл.
Аркадий предлагает следующую проверку:

#include <limits.h>

#define V_size 64 /* степень двойки */

#if ((RAND_MAX + 1) & RAND_MAX) == 0 /* RAND_MAX+1 - степень двойки */

        /* практически любой компилятор эти деления превратит в сдвиг */
        unsigned j = rand_Y () / (Y_m_div / V_size);

#elif UINT_MAX / V_SIZE - 1 >= RAND_MAX

        unsigned j = rand_Y () * V_size / Y_m_div;

#elif ULONG_MAX / V_SIZE - 1 >= RAND_MAX

        unsigned j = (unsigned)((unsigned long) rand_Y () * V_size / Y_m_div);

#else

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

#endif

И еще одна цитата Аркадия:
"В случае, если k*(RAND_MAX+1) <= UINT_MAX (где тип результата rand() -
int или unsigned), то тут никакие кастинги не нужны (смотри пример выше), в
противном случае на некоторых числах, генеримых rand() (точнее, больших
UINT_MAX/k), умножение отбросит сташие биты.

Для иллюстрации:

unsigned r = INT_MAX;
unsigned j = 4 * r / (INT_MAX+1);

Если бы здесь не было переполнения, мы получилим бы j=3, однако переполнение
тут есть."

Алгоритм B - рандомизация перемешиванием

Как нам кажется, Кнут (или переводчики) допустил неточность при описании алгоритма - не ясно, что происходит с переменной Y. Однако MIX-код и комментарии к нему прояснили ситуацию.

unsigned rand_Y;
void init_mix_rand () 
{
        for ...
        rand_Y = rand_X ();
}
unsigned mix_rand () 
{
        unsigned j = rand_Y * k / m;
        unsigned ret = rand_Y = V[j];
        V[j] = rand_X ();
        return ret;
}

Неплохой линк:

http://choices.cs.uiuc.edu/~akapadia/project2/project2.html

Прошу прощения, что так часто стираю собственные темы в треде, т.к. посижу-подумаю, потом смотрю - господи, что я написал???

В книге Седжвика (ROBERT SEDGEWICK) приводится самопальный пример на прохождение PRNG теста хи-квадрат. Вот код на паскале:

function сhisquare(N, r, s: integer) : real;
var i, t: integer;
f: array [O..rmax] of integer;
begin
randnit(s);
for i:=O to rmax do f[i] :=O;
for i:=l to N do
begin
t:=randomint(r);
f[t]::=f[t]+l;
end ;
t:=0; for i:=0 to r-l do t:=t+f[i]*f[i];
chisquare:= ((r*t/N) - N);
end ;

Если внимательно посмотреть на строчку

t:=0; for i:=0 to r-l do t:=t+f[i]*f[i];

то видно, что t определена как integer, в нее постоянно добавляются другие integer из массива f, что приводит к переполнению (см. топ "Integer overflow"). Автор и сам осознал свою ошибку, поэтому в С-версии тихо переписал код на такой:

#include <math.h>
#include <stdlib.h>
typedef int numType;
#define R 1000
numType randNum()
  { return rand() % R; }
main(int argc, char *argv[])
  { int i, N = atoi(argv[1]);
    int *f = malloc(R*sizeof(int));
    float m1 = 0.0, m2 = 0.0, t = 0.0;
    numType x;
    for (i = 0; i < R; i++) f[i] = 0;
    for (i = 0; i < N; i++)
      {
        f[x = randNum()]++;
        m1 += (float) x/N; 
        m2 += (float) x*x/N;
      }
    for (i = 0; i < R; i++) t += f[i]*f[i];
    printf("       Average: %f\n", m1);
    printf("Std. deviation: %f\n", sqrt(m2-m1*m1));
    printf("    Chi-square: %f\n", (R*t/N)-N);
 }

Здесь тип t уже float, т.е. переполнение уже не грозит.

Ну хорошо, душа моя успокоилась.

Black_mirror

Я прогнал тест с последовательностью rdtsc. Формула в Седжвике по которой он построил свой тест chisquare - это формула №6 в Кнуте в разделе 3.3.1 (второе издание). Последовательность, генерируемая rdtsc, по твоему предложению (я отдаю себе отчет, что это только пример), теста хи-квадрат не проходит. Более того, даже на глаз видно (если создать гистограмму частот с наложенным нормальным распределением - любой нормальный пакет позволяет делать такое), что и тесту Колмогорова-Смирнова тут делать нечего. Словом, rdtsc можно использовать лишь раз - для затравки. В качестве последовательности применять его нельзя.

Теперь перехожу к спектральному тесту.

volodya
У меня rdtsc и используется только для инициализации массива. Хотя конечно инициализировать лучше так:

rnd_init:
        mov ecx,64
    .l0:
        rdtsc
        mov [rnd_data+ecx*4-4],eax
        loop .l0

        mov ecx,100000
    .l1:
        call rnd
        loop .l1 
        ret

Чтобы избавиться от наследия rdtsc.

А div [rand_maxp1] лучше наверно выкинуть нафиг.

У меня rdtsc и используется только для инициализации массива.

А я и клоню к тому, что этого делать НЕ стоит. Никакой линейный конгруэнтный генератор этого барахла НЕ вытянет. Последовательность должна быть сгенерирована как-то иначе. А rdtsc годится лишь для разового применения, например, вместо

srand(time(0));
rand();

лучше как раз использовать rdtsc. Вот так, например:

int get_rand()
{
	int seed;
	__asm
	{
		rdtsc
		mov seed, eax
	}
	return seed;
}

srand(get_rand());
rand();

А в остальных случаях о ней лучше забыть.

А div [rand_maxp1] лучше наверно выкинуть нафиг.

Ни в коем случае! Она спасает от переполнения. В топе о переполнении я давал линк на exploits о integer overflow. Там чел очень хорошо расписал, что такое "modulo".

volodya

X_n = (X_n-24 + X_n-55) mod m, n > = 55

У Кнута написано что длина периода данного генератора, при m равном 2^s, равна 2^s-1*(2⁵⁵-1). То есть для однобитного генератора у нас фактически получается две последовательности: одна соответствует случаю когда генератор генерирует только нули, так как все 55 чисел нулевые, а вот длина периода второй последовательности как раз и совпадает с количеством оставшихся состояний генератора. То есть в процессе её генерирования будут пройдены все 2⁵⁵-1 состояния генератора(кроме случая когда все числа нули), поэтому любое из них можно взять в качестве начального.

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


Ни в коем случае! Она спасает от переполнения
Отсутствие div эквивалентно div 2³² и Кнут пишет что m должно быть четным поэтому выкинуть div нафиг!

rand:
        mov edx,[rnd_index]
        add dl,40
        mov eax,[rnd_data+edx]
        add dl,124
        add eax,[rnd_data+edx]
        add dl,96
        mov [rnd_data+edx],eax
        mov [rnd_index],edx
        ret 

subscribe

но так как я не занимаюсь статистическим моделированием, мне его вполне хватает.

Да и я им не занимаюсь :)
И данного генератора действительно хватит для подавляющего числа прикладных программ (за исключением криптографических и каких-нибудь сложных мат. пакетов - например, для NASA этот PRNG не подойдет :)).
Но я к чему клоню.

Формула, что ты привел - она массив требует. Я клоню к тому, что rdtsc заполнять его не стоит. А для такой вот формулы:

Xn+1 = (aXn+ c)mod m

rdtsc в качестве затравки (X₀) вполне пойдет.

Отсутствие div эквивалентно div 2^32

Не понял!

volodya
Формула, что ты привел - она массив требует. Я клоню к тому, что rdtsc заполнять его не стоит.

Еще раз смотрим на однобитный генератор: длина его периода 2⁵⁵-1. Общее число его состояний 2⁵⁵. Состояние которое перехожит само в себя - это когда все числа нули. Следовательно какую бы лажу выдала rdtsc(щас нас интересуют только младшие биты, если конечно они не все нулевые), эта последовательность для инициализации массива вполне подойдет. Потому, что эта лажа на некотором шаге обязательно будет состоянием генератора, а если так, то почему бы его не взять в качестве самого первого состояния?

А что касается div, то ведь сложение у нас происходит по модулю 2^32, поэтому будем считать, что мой последний вариант генерирует числа в диапазоне 0 .. FFFFFFFFh 8)