Commit d77a7f24 authored by Benoit Bayol's avatar Benoit Bayol

initial commit from version 3.05 of http://maths-people.anu.edu.au/~brent/random.html

parents
xorgens version 3.05 (20 September 2008)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
Contents:
test-xorgens.c Simple test program.
xorgens.h typedef statements etc.
xorgens.c xor4096i (for integer random numbers) for
xor4096r (for real random numbers).
xorgensyz.out Output from the test program in four cases:
yz = 44: 32-bit integer, 32-bit real (float);
yz = 48: 32-bit integer, 64-bit real (double);
yz = 84: 64-bit integer, 32-bit real;
yz = 88: 64-bit integer, 64-bit real.
rpb224.pdf A preprint describing the theory.
R. P. Brent <xorgens@rpbrent.com>
File added
/* test-xorgens.c
A simple test program for xorgens version 3.04.
Thanks to K. M. Briggs for assistance in making this portable.
Needs xorgens.h, xorgens.c
R. P. Brent, 20060628
===================================================================
Compilation hints:
On some machines (e.g. Sun Ultrasparc) can compile
in 32-bit integer mode with:
gcc -O4 -o test-xorgens -lm test-xorgens.c
or 64-bit integer mode with:
cc -Ofast -xtarget=ultra -xarch=v9 -lm -o test-xorgens test-xorgens.c
On other machines (e.g. AMD Opteron) 64-bit mode is the default:
gcc -O4 -o test-xorgens -lm test-xorgens.c
but if you really want 32-bit integers you can get them by changing
the definition of UINT in xorgens.h
With Fedora core 4 on AMD64, and Fedora core 5 on Pentiums:
gcc -O3 -Wall test-xorgens.c -o test-xorgens -lm
*/
#include <stdio.h>
#include <math.h>
#include "xorgens.h"
#include "xorgens.c"
const int bins=64;
const int trials=10000000;
const int seed0=1;
int main() {
UINT i, n=trials, seed=seed0, h[bins];
double x, sum, e, xerr;
printf("sizeof(UINT) = %d (should be 4 or 8)\n", (int)sizeof(UINT));
printf("sizeof(UREAL) = %d (should be 4 or 8)\n", (int)sizeof(UREAL));
printf("%d trials and %d bins with seed %ld\n",(int)n,(int)bins,(long)seed);
for (i=0; i<bins; i++) h[i]=0;
x = (double)xor4096r(seed);
printf("First random number %8g\n",x);
for (i=0; i<n; i++) h[(unsigned int)(bins*(double)xor4096r(0))]++;
for (i=0; i<bins; i++) printf("%3d %d\n",(int)i,(int)h[i]);
e = (double)n/(double)bins; /* Expected number per bin */
for (i=0, sum=0; i<bins; i++) /* Compute chi-squared statistic */
sum += ((double)h[i] - e)*((double)h[i] - e);
sum = sum/e;
printf ("Chi-squared_%d %g\n", bins-1, sum);
xerr = sqrt(2.0*sum) - sqrt(2.0*(bins-1)-1.0);
/* If bins is not too small then xerr should have normal distribution */
printf ("Following should be roughly N(0,1): %g\n", xerr);
return 0;
}
/* xorgens.c version 3.05, R. P. Brent, 20080920.
==========================================================================
| |
| Copyright (C) 2004, 2006, 2008 R. P. Brent. |
| |
| This program is free software; you can redistribute it and/or |
| modify it under the terms of the GNU General Public License, |
| version 2, June 1991, as published by the Free Software Foundation. |
| For details see http://www.gnu.org/copyleft/gpl.html . |
| |
| If you would like to use this software but the GNU GPL creates legal |
| problems, then please contact the author to negotiate a special |
| agreement. |
| |
==========================================================================
For type definitions see xorgens.h */
UINT xor4096i(UINT seed)
{
/* 32-bit or 64-bit integer random number generator
with period at least 2**4096-1.
It is assumed that "UINT" is a 32-bit or 64-bit integer
(see typedef statements in xorgens.h).
xor4096i should be called exactly once with nonzero seed, and
thereafter with zero seed.
One random number uniformly distributed in [0..2**wlen) is returned,
where wlen = 8*sizeof(UINT) = 32 or 64.
R. P. Brent, 20060628.
*/
/* UINT64 is TRUE if 64-bit UINT,
UINT32 is TRUE otherwise (assumed to be 32-bit UINT). */
#define UINT64 (sizeof(UINT)>>3)
#define UINT32 (1 - UINT64)
#define wlen (64*UINT64 + 32*UINT32)
#define r (64*UINT64 + 128*UINT32)
#define s (53*UINT64 + 95*UINT32)
#define a (33*UINT64 + 17*UINT32)
#define b (26*UINT64 + 12*UINT32)
#define c (27*UINT64 + 13*UINT32)
#define d (29*UINT64 + 15*UINT32)
#define ws (27*UINT64 + 16*UINT32)
static UINT w, weyl, zero = 0, x[r];
UINT t, v;
static int i = -1 ; /* i < 0 indicates first call */
int k;
if ((i < 0) || (seed != zero)) { /* Initialisation necessary */
/* weyl = odd approximation to 2**wlen*(3-sqrt(5))/2. */
if (UINT32)
weyl = 0x61c88647;
else
weyl = ((((UINT)0x61c88646)<<16)<<16) + (UINT)0x80b583eb;
v = (seed!=zero)? seed:~seed; /* v must be nonzero */
for (k = wlen; k > 0; k--) { /* Avoid correlations for close seeds */
v ^= v<<10; v ^= v>>15; /* Recurrence has period 2**wlen-1 */
v ^= v<<4; v ^= v>>13; /* for wlen = 32 or 64 */
}
for (w = v, k = 0; k < r; k++) { /* Initialise circular array */
v ^= v<<10; v ^= v>>15;
v ^= v<<4; v ^= v>>13;
x[k] = v + (w+=weyl);
}
for (i = r-1, k = 4*r; k > 0; k--) { /* Discard first 4*r results */
t = x[i = (i+1)&(r-1)]; t ^= t<<a; t ^= t>>b;
v = x[(i+(r-s))&(r-1)]; v ^= v<<c; v ^= v>>d;
x[i] = t^v;
}
}
/* Apart from initialisation (above), this is the generator */
t = x[i = (i+1)&(r-1)]; /* Assumes that r is a power of two */
v = x[(i+(r-s))&(r-1)]; /* Index is (i-s) mod r */
t ^= t<<a; t ^= t>>b; /* (I + L^a)(I + R^b) */
v ^= v<<c; v ^= v>>d; /* (I + L^c)(I + R^d) */
x[i] = (v ^= t); /* Update circular array */
w += weyl; /* Update Weyl generator */
return (v + (w^(w>>ws))); /* Return combination */
#undef UINT64
#undef UINT32
#undef wlen
#undef r
#undef s
#undef a
#undef b
#undef c
#undef d
#undef ws
}
UREAL xor4096r(UINT seed)
{
/* 64-bit or 32-bit real random number generator
with period at least 2**4096-1.
It is assumed that "UINT" is a 32-bit or 64-bit integer and "UREAL"
is "double" or "float". If "double" this is an IEEE standard
floating-point number with 53 bits in the fraction; if "single" it
has 24 bits in the fraction (including 1 implicit bit in each case).
In the 64-bit integer case, the method used is to call xor4096i to get
64 random bits, then the high 53 (for double) or 24 bits (for float)
are scaled to the open interval (0.0, 1.0), except that they are
discarded if all zero.
In the 32-bit integer case, one or two calls to xor4096i are made to
get 32 or 64 random bits, some are discarded, and the remaining bits
(if nonzero) are scaled to the open interval (0.0, 1.0).
xor4096r should be called exactly once with nonzero seed, and
thereafter with zero seed.
One random number of type UREAL is returned per call.
The results be should uniformly distributed in (0.0, 1.0) to the
resolution of the floating-point system (0.5**53 or 0.5**24).
The results are never exactly 0.0 or 1.0.
R. P. Brent, 20060628.
*/
#define UINT64 (sizeof(UINT)>>3)
#define UINT32 (1 - UINT64)
#define UREAL64 (sizeof(UREAL)>>3)
#define UREAL32 (1 - UREAL64)
/* sr = number of bits discarded = 11 for double, 40 or 8 for float */
#define sr (11*UREAL64 +(40*UINT64 + 8*UINT32)*UREAL32)
/* ss (used for scaling) is 53 or 21 for double, 24 for float */
#define ss ((53*UINT64 + 21*UINT32)*UREAL64 + 24*UREAL32)
/* SCALE is 0.5**ss, SC32 is 0.5**32 */
#define SCALE ((UREAL)1/(UREAL)((UINT)1<<ss))
#define SC32 ((UREAL)1/((UREAL)65536*(UREAL)65536))
UREAL res;
res = (UREAL)0;
while (res == (UREAL)0) /* Loop until nonzero result. */
{ /* Usually only one iteration . */
res = (UREAL)(xor4096i(seed)>>sr); /* Discard sr random bits. */
seed = (UINT)0; /* Zero seed for next time. */
if (UINT32 && UREAL64) /* Need another call to xor4096i. */
res += SC32*(UREAL)xor4096i(seed);/* Add low-order 32 bits. */
}
return (SCALE*res); /* Return result in (0.0, 1.0). */
#undef UINT64
#undef UINT32
#undef UREAL64
#undef UREAL32
#undef SCALE
#undef SC32
#undef sr
#undef ss
}
/* xorgens.h for xorgens version 3.04, R. P. Brent 20060628 */
typedef unsigned long UINT; /* Type for random 32 or 64-bit integer,
e.g. unsigned long, unsigned long long,
uint64_t, unsigned int or uint32_t */
typedef double UREAL; /* Type for random 32 or 64-bit real,
e.g. double or float */
UINT xor4096i(UINT seed); /* integer random number generator */
UREAL xor4096r(UINT seed); /* real random number generator */
sizeof(UINT) = 4 (should be 4 or 8)
sizeof(UREAL) = 4 (should be 4 or 8)
10000000 trials and 64 bins with seed 1
First random number 0.152044
0 156548
1 156110
2 156504
3 156602
4 156970
5 156436
6 155900
7 156727
8 156658
9 155868
10 156124
11 155515
12 156798
13 156250
14 155770
15 155819
16 156391
17 155706
18 155823
19 156022
20 155710
21 156219
22 156269
23 156436
24 156284
25 156579
26 156233
27 155779
28 156319
29 156681
30 155089
31 156178
32 156740
33 156482
34 156222
35 156279
36 155760
37 156849
38 157158
39 156095
40 155684
41 156569
42 156182
43 155784
44 155734
45 156097
46 156152
47 156329
48 156087
49 156177
50 156652
51 156581
52 156435
53 155538
54 156563
55 157009
56 155717
57 155879
58 156104
59 156722
60 156651
61 156806
62 155842
63 156804
Chi-squared_63 72.6872
Following should be roughly N(0,1): 0.876793
sizeof(UINT) = 4 (should be 4 or 8)
sizeof(UREAL) = 8 (should be 4 or 8)
10000000 trials and 64 bins with seed 1
First random number 0.152044
0 156300
1 155777
2 156604
3 156302
4 156564
5 156262
6 156060
7 156663
8 156463
9 155459
10 156573
11 156275
12 156558
13 156263
14 156077
15 156422
16 156366
17 155690
18 155632
19 156391
20 156316
21 156617
22 155585
23 156520
24 156454
25 156410
26 155977
27 156105
28 156552
29 156439
30 155311
31 155913
32 156558
33 155972
34 156403
35 156153
36 155538
37 157160
38 156358
39 156325
40 155498
41 156168
42 156208
43 155855
44 156149
45 155627
46 156627
47 156479
48 156111
49 155987
50 156341
51 156607
52 156809
53 155525
54 156450
55 156703
56 155779
57 156227
58 156236
59 157060
60 156279
61 156632
62 156216
63 157060
Chi-squared_63 64.8234
Following should be roughly N(0,1): 0.205918
sizeof(UINT) = 8 (should be 4 or 8)
sizeof(UREAL) = 4 (should be 4 or 8)
10000000 trials and 64 bins with seed 1
First random number 0.720041
0 156665
1 156070
2 155227
3 155928
4 156837
5 157179
6 156354
7 156131
8 156709
9 156443
10 155702
11 156149
12 156443
13 156537
14 156246
15 155864
16 156568
17 155851
18 155921
19 155961
20 155975
21 156580
22 155896
23 155848
24 156131
25 157103
26 156436
27 157050
28 156251
29 156364
30 156434
31 156329
32 156818
33 156383
34 155872
35 156568
36 156346
37 156018
38 155516
39 156711
40 156285
41 156853
42 155904
43 155817
44 156055
45 156345
46 156614
47 156402
48 156355
49 155692
50 156099
51 155948
52 155846
53 156373
54 155943
55 156121
56 156049
57 156556
58 155836
59 155906
60 156959
61 156521
62 155672
63 156435
Chi-squared_63 65.1177
Following should be roughly N(0,1): 0.231731
sizeof(UINT) = 8 (should be 4 or 8)
sizeof(UREAL) = 8 (should be 4 or 8)
10000000 trials and 64 bins with seed 1
First random number 0.720041
0 156665
1 156070
2 155227
3 155928
4 156837
5 157179
6 156354
7 156131
8 156709
9 156443
10 155702
11 156149
12 156443
13 156537
14 156246
15 155864
16 156568
17 155851
18 155921
19 155961
20 155975
21 156580
22 155896
23 155848
24 156131
25 157103
26 156436
27 157050
28 156251
29 156364
30 156434
31 156329
32 156818
33 156383
34 155872
35 156568
36 156346
37 156018
38 155516
39 156711
40 156285
41 156853
42 155904
43 155817
44 156055
45 156345
46 156614
47 156402
48 156355
49 155692
50 156099
51 155948
52 155846
53 156373
54 155943
55 156121
56 156049
57 156556
58 155836
59 155906
60 156959
61 156521
62 155672
63 156435
Chi-squared_63 65.1177
Following should be roughly N(0,1): 0.231731
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment