Velocity Reviews - Computer Hardware Reviews

Velocity Reviews > Newsgroups > Programming > C Programming > Re: C code for pure float PRNG

Reply
Thread Tools

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>,
Adam Ierymenko <(E-Mail Removed)> wrote:

> 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
Added doubleFastRng
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 */
 
Reply With Quote
 
 
 
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are Off


Similar Threads
Thread Thread Starter Forum Replies Last Post
float to string to float, with first float == second float Carsten Fuchs C++ 45 10-08-2009 09:47 AM
Good binary PRNG pinkfloydhomer@gmail.com C Programming 31 05-19-2006 09:08 PM
init and set state of PRNG with VC++ Jean-Michel Besnard C++ 0 07-27-2004 05:12 PM
Re: float->byte->float is same with original float image. why float->ubyte->float is different??? bd C Programming 0 07-07-2003 12:09 AM
PRNG Algorithm for RAN() John Schutkeker C Programming 4 07-06-2003 11:43 AM



Advertisments