Re: C code for pure float PRNG

Discussion in 'C Programming' started by Francois Grieu, Feb 5, 2004.

  1. /* this message intends to be a valid ANSI C program

    In article <>,
    Adam Ierymenko <> 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 */
     
    Francois Grieu, Feb 5, 2004
    #1
    1. Advertising

Want to reply to this thread or ask your own question?

It takes just 2 minutes to sign up (and it's free!). Just click the sign up button to choose a username and then you can ask your own questions on the forum.
Similar Threads
  1. Jean-Michel Besnard

    init and set state of PRNG with VC++

    Jean-Michel Besnard, Jul 27, 2004, in forum: C++
    Replies:
    0
    Views:
    633
    Jean-Michel Besnard
    Jul 27, 2004
  2. John Schutkeker

    PRNG Algorithm for RAN()

    John Schutkeker, Jul 4, 2003, in forum: C Programming
    Replies:
    4
    Views:
    2,039
    John Schutkeker
    Jul 6, 2003
  3. bd
    Replies:
    0
    Views:
    665
  4. Good binary PRNG

    , May 15, 2006, in forum: C Programming
    Replies:
    31
    Views:
    1,016
    CBFalconer
    May 19, 2006
  5. Carsten Fuchs
    Replies:
    45
    Views:
    1,651
    James Kanze
    Oct 8, 2009
Loading...

Share This Page