KISS4691, a potentially top-ranked RNG.
in [Fortran]
|
Prev: solutions book
Next: real kind declaration
From: Uno on 2 Aug 2010 21:26 sturlamolden wrote: > I did a speed comparison on my laptop of KISS4691 against Mersenne > Twister 19937. KISS4691 produced about 110 million random numbers per > second, whereas MT19937 produced 118 million. If the compiler was > allowed to inline KISS4691, the performance increased to 148 millon > per second. The function call overhead is imporant, and without taking > it into consideration, MT19937 will actually be faster than KISS4691. > > Also note that a SIMD oriented version of MT19937 is twice as fast as > the one used in my tests. There is also a version of Mersenne Twister > suitable for parallel processing and GPUs. So a parallel or SIMD > version of Mersenne Twister is likely to be a faster PRNG than > KISS4691. > > Then for the question of numerical quality: Is KISS4691 better than > MT19937? I don't feel qualified to answer this. Marsaglia is the > world's foremost authority on random number generators, so I trust his > advice, but MT19937 is claimed to give impeckable quality except for > crypto. > > For those who need details: > > The speed test was basically running the PRNGs for five seconds or > more, querying Windows' performance counter on every ten million call, > and finally subtracting an estimate of the timing overhead. The C > compiler was Microsoft C/C++ version 15.00.30729.01 (the only > optimization flag I used was /O2, i.e. maximize for speed). The > processor is AMD Phenom N930 @ 2.0 Gz. I'd be curious to see the source you used for this purpose. This is very-slightly adapted from Dann Corbitt in ch 13 of unleashed. $ gcc -Wall -Wextra cokus2.c -o out $ ./out 3510405877 4290933890 2191955339 564929546 152112058 4262624192 2687398418 268830360 1763988213 578848526 4212814465 3596577449 4146913070 950422373 1908844540 1452005258 3029421110 142578355 1583761762 1816660702 2530498888 1339965000 3874409922 3044234909 1962617717 2324289180 310281170 981016607 908202274 3371937721 2244849493 675678546 3196822098 1040470160 3059612017 3055400130 2826830282 2884538137 3090587696 2262235068 3506294894 2080537739 1636797501 4292933080 2037904983 2465694618 1249751105 30084166 112252926 1333718913 880414402 334691897 3337628481 17084333 1070118630 .... 3009858040 3815089086 2493949982 3668001592 1185949870 2768980234 3004703555 1411869256 2625868727 3108166073 3689645521 4191339889 1933496174 1218198213 3716194408 1148391246 1345939134 3517135224 3320201329 4292973312 3428972922 1172742736 275920387 617064233 3754308093 842677508 529120787 1121641339 $ cat cokus2.c #include <stdio.h> #include <stdlib.h> #include "mtrand.h" /* uint32 must be an unsigned integer type capable of holding at least 32 bits; exactly 32 should be fastest, but 64 is better on an Alpha with GCC at -O3 optimization so try your options and see what's best for you */ typedef unsigned long uint32; /* length of state vector */ #define N (624) /* a period parameter */ #define M (397) /* a magic constant */ #define K (0x9908B0DFU) /* mask all but highest bit of u */ #define hiBit(u) ((u) & 0x80000000U) /* mask all but lowest bit of u */ #define loBit(u) ((u) & 0x00000001U) /* mask the highest bit of u */ #define loBits(u) ((u) & 0x7FFFFFFFU) /* move hi bit of u to hi bit of v */ #define mixBits(u, v) (hiBit(u)|loBits(v)) /* state vector + 1 extra to not violate ANSI C */ static uint32 state[N + 1]; /* next random value is computed from here */ static uint32 *next; /* can *next++ this many times before reloading */ static int left = -1; /* ** ** We initialize state[0..(N-1)] via the generator ** ** x_new = (69069 * x_old) mod 2^32 ** ** from Line 15 of Table 1, p. 106, Sec. 3.3.4 of Knuth's ** _The Art of Computer Programming_, Volume 2, 3rd ed. ** ** Notes (SJC): I do not know what the initial state requirements ** of the Mersenne Twister are, but it seems this seeding generator ** could be better. It achieves the maximum period for its modulus ** (2^30) iff x_initial is odd (p. 20-21, Sec. 3.2.1.2, Knuth); if ** x_initial can be even, you have sequences like 0, 0, 0, ...; ** 2^31, 2^31, 2^31, ...; 2^30, 2^30, 2^30, ...; 2^29, 2^29 + 2^31, ** 2^29, 2^29 + 2^31, ..., etc. so I force seed to be odd below. ** ** Even if x_initial is odd, if x_initial is 1 mod 4 then ** ** the lowest bit of x is always 1, ** the next-to-lowest bit of x is always 0, ** the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... , ** the 3rd-from-lowest bit of x 4-cycles ... 0 1 1 0 0 1 1 0 ... , ** the 4th-from-lowest bit of x has the 8-cycle ... 0 0 0 1 1 1 1 0 ... , ** ... ** ** and if x_initial is 3 mod 4 then ** ** the lowest bit of x is always 1, ** the next-to-lowest bit of x is always 1, ** the 2nd-from-lowest bit of x alternates ... 0 1 0 1 0 1 0 1 ... , ** the 3rd-from-lowest bit of x 4-cycles ... 0 0 1 1 0 0 1 1 ... , ** the 4th-from-lowest bit of x has the 8-cycle ... 0 0 1 1 1 1 0 0 ... , ** ... ** ** The generator's potency (min. s>=0 with (69069-1)^s = 0 mod 2^32) is ** 16, which seems to be alright by p. 25, Sec. 3.2.1.3 of Knuth. It ** also does well in the dimension 2..5 spectral tests, but it could be ** better in dimension 6 (Line 15, Table 1, p. 106, Sec. 3.3.4, Knuth). ** ** Note that the random number user does not see the values generated ** here directly since reloadMT() will always munge them first, so maybe ** none of all of this matters. In fact, the seed values made here could ** even be extra-special desirable if the Mersenne Twister theory says ** so-- that's why the only change I made is to restrict to odd seeds. */ void mtsrand(uint32 seed) { register uint32 x = (seed | 1U) & 0xFFFFFFFFU, *s = state; register int j; for (left = 0, *s++ = x, j = N; --j; *s++ = (x *= 69069U) & 0xFFFFFFFFU); } uint32 reloadMT(void) { register uint32 *p0 = state, *p2 = state + 2, *pM = state + M, s0, s1; register int j; if (left < -1) mtsrand(4357U); left = N - 1, next = state + 1; for (s0 = state[0], s1 = state[1], j = N - M + 1; --j; s0 = s1, s1 = *p2++) *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); for (pM = state, j = M; --j; s0 = s1, s1 = *p2++) *p0++ = *pM++ ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); s1 = state[0], *p0 = *pM ^ (mixBits(s0, s1) >> 1) ^ (loBit(s1) ? K : 0U); s1 ^= (s1 >> 11); s1 ^= (s1 << 7) & 0x9D2C5680U; s1 ^= (s1 << 15) & 0xEFC60000U; return (s1 ^ (s1 >> 18)); } uint32 mtrand(void) { uint32 y; if (--left < 0) return (reloadMT()); y = *next++; y ^= (y >> 11); y ^= (y << 7) & 0x9D2C5680U; y ^= (y << 15) & 0xEFC60000U; return (y ^ (y >> 18)); } #define UNIT_TEST #ifdef UNIT_TEST int main(void) { int j; /* you can seed with any uint32, but the best are odds in 0..(2^32 - 1) */ mtsrand(4357U); /* print the first 2,002 random numbers seven to a line as an example */ for (j = 0; j < 2002; j++) printf(" %10lu%s", (unsigned long) mtrand(), (j % 7) == 6 ? "\n" : ""); return (EXIT_SUCCESS); } #endif // gcc -Wall -Wextra cokus2.c -o out $ cat mtrand.h /* ** The proper usage and copyright information for ** this software is covered in DSCRLic.TXT ** This code is Copyright 1999 by Dann Corbit */ /* ** Header files for Mersenne Twister pseudo-random number generator. */ extern void mtsrand(unsigned long seed); extern unsigned long reloadMT(void); extern unsigned long mtrand(void); $ -- Uno
From: sturlamolden on 2 Aug 2010 23:00 On 3 Aug, 03:26, Uno <merrilljen...(a)q.com> wrote: > I'd be curious to see the source you used for this purpose. This is > very-slightly adapted from Dann Corbitt in ch 13 of unleashed. http://folk.uio.no/sturlamo/prngtest.zip
From: Dann Corbit on 3 Aug 2010 00:37 In article <4d49afc5-15d7-4e5e-bd65- 6f0708d614ee(a)q35g2000yqn.googlegroups.com>, sturlamolden(a)yahoo.no says... > http://folk.uio.no/sturlamo/prngtest.zip > I got this: Mersenne Twister 19937 speed, single (Hz): 1.63278e+008 Mersenne Twister 19937 speed, array (Hz): 1.3697e+008 KISS 4691 speed, single (Hz): 1.86338e+008 KISS 4691 speed, array (Hz): 1.87675e+008
From: robin on 3 Aug 2010 06:41 "Uno" <merrilljensen(a)q.com> wrote in message news:8bfos5FakcU1(a)mid.individual.net... | robin wrote: | > I have already posted a PL/I version using unsigned arithmetic. | > | > Here is another version, this time using signed arithmetic :-- | > | > (NOSIZE, NOFOFL): | > RNG: PROCEDURE OPTIONS (MAIN, REORDER); | > | > declare (xs initial (521288629), xcng initial (362436069), | > Q(0:4690) ) static fixed binary (31); | > | > MWC: procedure () returns (fixed binary (31)); | > declare (t,x,i) fixed binary (31); | > declare (c initial (0), j initial (4691) ) fixed binary (31) static; | > declare (t1, t2, t3) fixed binary (31); | > | > if j < hbound(Q,1) then j = j + 1; else j = 0; | > x = Q(j); | > t = isll(x,13)+c+x; | > t1 = iand(x, 3) - iand(t, 3); | > t2 = isrl(x, 2) - isrl(t, 2); | > if t2 = 0 then t2 = t1; | > if t2 > 0 then t3 = 1; else t3 = 0; | > c = t3 + isrl(x, 19); | > Q(j)=t; | > return (t); | > end MWC; | > | > CNG: procedure returns (fixed binary (31)); | > xcng=bin(69069)*xcng+bin(123); | > return (xcng); | > end CNG; | > | > XXS: procedure returns (fixed binary (31)); | > xs = ieor (xs, isll(xs, 13) ); | > xs = ieor (xs, isrl(xs, 17) ); | > xs = ieor (xs, isll(xs, 5) ); | > return (xs); | > end XXS; | > | > KISS: procedure returns (fixed binary (31)); | > return ( MWC()+CNG+XXS ); | > end KISS; | > | > declare (i,x) fixed binary (31); | > declare y fixed decimal (11); | > | > Q = CNG+XXS; /* Initialize. */ | > do i = 1 to 1000000000; x=MWC(); end; | > put skip edit (" Expected MWC result = 3740121002", 'computed =', x) | > (a, skip, x(12), a, f(11)); | > y = iand(x, 2147483647); | > if x < 0 then y = y + 2147483648; | > put skip edit (y) (x(11), f(22)); put skip; | > do i = 1 to 1000000000; x=KISS; end; | > put skip edit ("Expected KISS result = 2224631993", 'computed =', x) | > (a, skip, x(12), a, f(11)); | > y = iand(x, 2147483647); | > if x < 0 then y = y + 2147483648; | > put skip edit (y) (x(11), f(22)); | > | > end RNG; | If you were to comment out the PL/I command line that compiled this, | what would it be? ???
From: James Waldby on 3 Aug 2010 13:15
On Tue, 03 Aug 2010 20:41:15 +1000, robin wrote: > "Uno" <merrilljensen> wrote: [snip code] >> If you were to comment out the PL/I command line that compiled this, >> what would it be? > > ??? Does that mean you don't understand Uno's question, or don't know the answer? In case you don't understand the question, it appears to be: "What command is used to compile the code?" -- jiw |