Оптимизация PRNG

Код:

Код (Text):

1.  
2. #define m 100000000L
3. #define m1 10000L
4. #define b 31415821L
5.  
6. unsigned long mult(long p, long q) {
7.         unsigned long p1=p/m1, p0=p%m1, q1=q/m1, q0=q%m1;
8.         return (((p0*q1+p1*q0) % m1) * m1+p0*q0) % m;
9. };
10.  
11. unsigned long NextRandomRange(long range, long* seed) {
12.         *seed=(mult(*seed, b)+1) % m;
13.         return (((*seed/m1)*range)/m1);
14. };
15.  
16. unsigned long NextRandomNumber(long* seed) {
17.         long n1=NextRandomRange(256, seed);
18.         long n2=NextRandomRange(256, seed);
19.         long n3=NextRandomRange(256, seed);
20.         long n4=NextRandomRange(256, seed);
21.         return (n1<<24)|(n2<<16)|(n3<<8)|n4;
22. };
23.  

Цель - максимальная оптимизация по скорости выполнения. Ну просто максимальная. На асме. Всяких там SSE2 и прочих я не знаю, но если их можно применить - то пусть будет Если такая оптимизация отожрет гиг памяти, но алго будет выполнятся быстрее - мя это устраивает.



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

 

/m1 встречается часто, можно медленную операцию деления заменить на обработку фиксированного предвычисленного значения, известного как RECIPROCAL_DIVISOR (A.Fog).



А что говорит cl /Ox ? Очень медленно?

 

Как это ни странно - да. Результаты в "попугаях" (меньше-лучше)



1. CPU = Intel Celeron 2GHz

1.1 cl 13.10.3077

Флаги: /G7 /Ox

Попугаи: 25

1.2 gcc 3.4.4

Флаги: -O3 -march=pentium4 -ftracer

Попугаи: 19



2. CPU = AMD Mobile Athlon XP 2500+

2.1 cl 13.10.3077

Флаги: /G7 /Ox

Попугаи: 14

2.2 gcc 3.4.4

Флаги: -O3 -march=athlon -ftracer

Попугаи: 8



/me был несколько удивлен. Пошел ставить Intel C Compiler

 

volodya

У тебя seed не превышает m, то есть начав с некоторого значения seed, не более чем через m вызовов функций NextRandomRange значение seed станет равно одному из уже встретившихся. То есть в данном случае не более чем через 100 млн. итераций значения начнут повторяться. Если длина последовательности генерируемой функцией NextRandomRange кратна 2 или 4, то длина последовательности, генерируемой функцией NextRandomNumber, будет соответственно в 2 или в 4 раза меньше. Так что в пол гига памяти поместится вся последовательность (можно сразу записывать результаты вызова NextRandomNumber).

 

самое медленное тут /m1, с ним можно поступить так, избавившись от div:

Код (Text):

1. #define RD -776530087
2.  
3.     unsigned long p1=p/m1
4.  
5.     превратить в
6.  
7. __asm {
8.     mov     eax,p             // eax - dividend
9.     mul     RD
10.     shr     edx,13            // edx - result (p/10000)
11.     }

Точно так же и с q1=q/m1 ; *seed/m1 ; (((*seed/m1)*range)/m1)

 

Это все верно. Но функция mult написана весьма хитро. И для того, чтобы действительно все влезло в полгига надо уметь из 32-битного seed получить соответствующий ему seed, меньший m.

 

cresta



Это делает оптимизирующий компилятор.



Black_mirror



Вот это уже очень светлая идея. Мы думали о таблице предпосчитаных значений, сорри, я просто не говорю всех деталей. В общем и в частности, там 15 гигов надо для прямого lookup. Это слишком дофига получается

 

volodya

Набросал тут замену unsigned long mult без использования div на 78 тактов против 185 у vc 7.1 с ключами оптимизации

cl /Ox /O2 /Ot

cl.exe вставил 4 div'а с этими ключами. Возможно есть другие ключи, с которыми будет без деления.

Если надо, приаттачу.

 

volodya

А функция mult это просто запись средствами С выражения (seed*b+1) mod m или в ней есть скрытый смысл? И почему её параметры не unsigned, получение отрицательного seed возможно?

Для функции NextRandomNumber вычисления можно распараллелить, при условии что между её вызовами не вызывается NextRandomRange.



До вызова функции считаем:



s1=(seed*b+1) mod m

s2=(s1*b+1) mod m

s3=(s2*b+1) mod m

s4=(s3*b+1) mod m



При вызове считаем:

n1=s1/m1*256/m1

n2=s2/m1*256/m1

n3=s3/m1*256/m1

n4=s4/m1*256/m1

и

s1=(s1*(b^4 mod m)+((b^3+b^2+b+1) mod m)) mod m

s2=(s2*(b^4 mod m)+((b^3+b^2+b+1) mod m)) mod m

s3=(s3*(b^4 mod m)+((b^3+b^2+b+1) mod m)) mod m

s4=(s4*(b^4 mod m)+((b^3+b^2+b+1) mod m)) mod m

ну а return оставляем как есть.



А вообще задача стоит ускорить именно эти функции, или работу вызывающего их кода?

 

volodya

Можно уменьшить число умножений/делений если хранить seed ввиде SH*m1+SL, где SH и SL<m1:

Код (Text):

1. #define m 100000000L
2. #define m1 10000L
3. #define bh 3141L
4. #define bl 5821L
5.  
6. struct seed_t
7. {
8.     long h,l;
9. }
10.  
11. unsigned long NextRandomRange(long range, seed_t* seed) {
12.         long t=seed->l*bl+1;
13.         seed->h=(seed->h*bl+seed->l*bh+t/m1)%m1;
14.         seed->l=t%m1;
15.         return (seed->h*range/m1);
16. };

 

А вообще задача стоит ускорить именно эти функции, или работу вызывающего их кода?



Хороший вопрос!

Я бы с удовольствием избавился бы от этих функций вообще. Смысл прост - с помощью данных функций наполняется некий буфер - char buff[]. Это узкое горлышко программы.



Если рассматривать NextRandomNumber как некую f(), на вход которой подается некоторое число, а на выходе имеем 1 < result < m, то я бы с удовольствием заменил бы эту f() на ту, которая дает аналогичный результат, но жрет меньше тактов.

 

cresta



Приаттачь, пожалуйста.

 

Black_mirror



Спасибо, очень интересные предложения! Буду разбираться



оффтоп: у тебя аси, случаем, нет?

 

Вообще говоря для положительных seed это действительно просто (seed*b+1) mod m. Надо лишь аккуратно разобраться, что происходит при отрицательных....



Если решить этот вопрос, то тогда действительно можно будет создать буфер из m 8 байтовых кусков (выход NextRandomNumber и следующий сид) и использовать его для lookup. Можно создавать буфер пятибайтовых кусков (выход NextRandomRange(256) и следующий сид).

 

volodya

Буфер заполняется последовательно при помощи NextRandomNumber?



NextRandomRange используется только в функции NextRandomNumber?



Кто еще и как часто кроме этих двух функций изменяет значение seed?



Для NextRandomRange можно возвращать seed*range/m или то, что используются только старшие цифры очень важно? И что нужно записать в буфер: просто случайные числа или числа сгенерированные именно этим генератором?



offtop: 341367720

 

RD и RD_1 рассчитаны для фиксированных m и m1 (100000000 и 10000).

Процедура stdcall, возможно переход к fastcall ещё что-то даст.



827258944__unsigned_long_mult.asm

 

cresta



спасибо

 

А вы верной дорогой идете, товариcчи ?

Насколько я понимаю если нач.значение seed неотрицательное и range < m1 все сводится к стандартным вычислениям:

Код (Text):

1. seed = (seed*b+1) % m
2. number = seed*range/m

Если вычислять в лоб через div, то суммарное время составит порядка 90-100 тиков (на P3,P4 и атлоне div кушает ~40 тактов + умножение и мелкие расходы). Прикинем насколько получится быстрее, если заменить div'ы умножениями. Проблема в том, что m и seed достаточно большие 32-битные числа и => сделать подряд два умножения не получится - нужно дробить числа на части как предлагает Black_mirror. При таком разбиении требуется минимум 4 умножения и 3 взятия остатка (по 2 умножения + сдвиг) - итого 10 умножений. P4 Northwood (~14 тиков на mul и 4 на shr) может отдыхать даже с поправкой на распараллеливание, Prescott (~10 на mul) - фифти-фифти. Из реальных претендентов, которые могут дать выигрыш остаются атлон (6 тиков/mul) и P3\PM (4 тика/mul).

Вот такие пироги - следует определиться с целевым процессором, иначе можно не выиграть, а проиграть



Хотя если оптимизировать NextRandomNumber в целом или считать NextRandomRange не по новому, а по предыдущему значению seed (предвычисление первого значения), то те же 50-60 тиков можно легко и просто получить за счет распараллеливания вычислений seed и number. Запускаем вычисление seed через mul и div и пока осуществляется деление вычисляем number по пред. значению seed на fpu (fild + fmul + fistp)



PS: Можно попробовать вариант с разбиением и заменой деления на умножения на SSE реализовать, но на вскидку вроде как существенного выигрыша я пока не вижу

 

Именно. Небольшие проблемы с отрицательными сидами ну и еще с конкретными деталями применения этого PRNG

Как показывает практика, на mobile Athlon XP 2500+ метод Black_mirror с разбитием дает примерно двукратный выигрыш по времени.



Как пока мне представляется - единственным существенным продвижением в этом деле может быть только полное табулирование значений генератора (~500Mb) и потом уже прямое вытаскивание всего из памяти... Хотя вопрос в скорости тут тоже остается открытым....

 

aSL

> "Небольшие проблемы с отрицательными сидами"

Да вроде как и при отрицательных сидах получается одинаковый результат. В данном случае q=b положительное. Поэтому при p=seed < 0 все частные и остатки в mult будут отрицательными (p1 и p0 отрицательные), как и при простом idiv.



> "дает примерно двукратный выигрыш по времени"

По сравнению с чем ? С простым div или с mult ?