Velocity Reviews > Re: C code for pure float PRNG

Re: C code for pure float PRNG

Francois Grieu
Guest
Posts: n/a

 02-05-2004
/* this message intends to be a valid ANSI C program

In article <(E-Mail Removed) om>,

> I'm looking for a pure floating point random number
> generator in C.

My 2 cents worth:
*/

/* <--- put what follows in a header file ---- */

/* doubleFastRng version 2
A random number generator using floating point in standard C.
NOT crypto quality, but will still pass many tests.

2004-02-05 V2.1 Better rationale
2004-02-05 V2 Lower odds that gFastRngA enters a short cycle
2004-02-04 V1 Creation

Public domain, by Francois Grieu
*/

/* Meat of the generator
Use the following doubleFastRng() macro as you would use:
double doubleFastRng(void);

Result is random-like and near-uniform in [0 to 1]

Note: low-order bits of mantissa are NOT very good
see below for more serious limitations
Note: speed depends a lot on the speed of floor()
Note: results are architecture dependent; this shows
after like a hundred iterations.

Rationale:
- some chaotic generator A is built
- A is turned into a uniform generator B; consecutive
pairs of B are NOT independent.
- B is turned into a better uniform generator C;
consecutive pairs of B sampled at random points should
be well distributed, but triplets of C are NOT,
this may matter in some applications (e.g. selecting
random points in a cube); men this is a SIMPLE
generator, not a perfect one !
After the first pass, then from pass to pass,
- A is in range [0.07438 to 1.62009], non-uniform
- B is in range [0 to 1], near-uniform
- C is in range [0 to 1], random-like
Also
- ignoring a C*k term and precision considerations,
the feedback for A is
A <- ((A mod 1) + 3/11) ^ 2
- this transformation has no stationary point
- order 2 cycles are unstable, and it is conjectured
that higher order cycles are unstable, excluding
numerical precision issues
- but it is conceivable that a particlar seed will give
a short cycle on a given architecture, and worst that
such cycle is entered by accident
- in a (perhap silly) attempt to counter this, generator
C adds a small feedback on A; this probably makes the
distribution of the output slightly less uniform;
a binning test would reveal this, I guess.

*/
#define doubleFastRng() ( \
(gFastRngA += gFastRngC*(1./9973)+3./11-floor(gFastRngA)), \
(gFastRngB += (gFastRngA *= gFastRngA)), \
(gFastRngC += (gFastRngB -= floor(gFastRngB))), \
(gFastRngC -= floor(gFastRngC)) \
)

/* Optional: seed the RNG
Use the following doubleFastRng() macro as you would use:
double doubleFastRngSeed(double x);

Seed value x can be any non-negative, non-huge double
*/
#define doubleFastRngSeed(x) do { \
gFastRngA = (x); gFastRngB = gFastRngC = 0; \
} while(0)

/* Optional: define an additional state of the random generator
in a local context, which might be faster.
Seed value x can be any non-negative, non-huge double
*/
#define NEWFASTRNG(x) \
double gFastRngA = (x), gFastRngB = 0, gFastRngC = 0

/* Define a default global context using above macro
Note: will harmlessly define multiple equivalent contexts
if the header file is included many times
*/
NEWFASTRNG(0);

#include <math.h> /* needed for floor(), used in doubleFastRng */

/* --- put the above in a header file ----> */

/* Test code for the above; values shown on second and third
columns are expected to be of either sign, and often
(but not allways) in the range [-1 to +1] */
#include <stdio.h>
#include <stdlib.h>
int main(int argc, char *argv[])
{
NEWFASTRNG(0); /* define a local context for faster code */
double s=0,t=0,r;
unsigned long j;
if (argc>1 && argv[1]) /* seed the RNG if an argument is given */
doubleFastRngSeed(atof(argv[1]));
for(j=1; j<=1<<28; ++j)
{
s += r = doubleFastRng() - 1./2;
t += r*r - 1./12;
if ((j&(j-1))==0) /* j a power of two */
printf("%10lu %-+12g %-+12g\n",j,s*2./sqrt(j),t*12./sqrt(j));
}
return 0;
}

/* results of the above on a PowerPC G4 with MrC under MPW
1 -0.85124 +1.17383
2 -0.928248 +0.574721
4 -0.88665 +0.221971
8 -0.32914 +0.705442
16 -0.630065 +0.736995
32 -0.0536323 +1.27563
64 -0.291745 +0.51509
128 -0.37266 +0.745991
256 -0.758367 +0.8036
512 -0.23023 +0.40877
1024 +0.322027 +0.863811
2048 +0.686953 +0.403408
4096 +0.274433 +1.66491
8192 -0.109999 +0.515839
16384 -0.487148 +0.28429
32768 -0.368597 +0.917197
65536 -0.364967 +1.70257
131072 +0.0397136 +1.73647
262144 +0.21818 +1.96605
524288 -0.369138 +0.767671
1048576 -0.804561 -0.926822
2097152 -0.154766 -1.0003
4194304 +0.440308 -1.4173
8388608 +0.443935 -2.18429
16777216 -0.507174 -0.808045
33554432 -0.250194 +0.234557
67108864 -0.384728 -1.01272
134217728 -0.40035 -1.04248
268435456 +0.576571 -0.562746
*/

/* results of the above on an Athlon XP with GCC under Win32
1 -0.85124 +1.17383
2 -0.928248 +0.574721
4 -0.88665 +0.221971
8 -0.32914 +0.705442
16 -0.630065 +0.736995
32 -0.0536323 +1.27563
64 -0.291745 +0.51509
128 -0.514703 +1.52217
256 -0.212286 +1.24088
512 -0.422173 -0.272143
1024 -1.14005 +0.293329
2048 -0.952017 +0.819124
4096 -0.583247 +0.982798
8192 -0.589248 +1.13498
16384 -0.215161 +2.21373
32768 -0.308418 +0.281441
65536 -0.402067 +0.0375496
131072 +0.191811 +0.447173
262144 +0.25173 +0.156406
524288 +0.345341 +0.834386
1048576 +0.240871 +0.985605
2097152 +0.0712098 +0.489712
4194304 +0.224906 +0.160624
8388608 -0.661534 +0.23568
16777216 -0.73162 -0.421065
33554432 -0.581571 -0.299703
67108864 +0.677287 -0.984173
134217728 +0.150005 -0.23203
268435456 +0.352258 +0.0239879
*/

/* Francois Grieu */