Very fast floating point PRNG

Discussion in 'Scientific Statistics Math' started by Ribeiro Alvo, Aug 25, 2009.

  1. Ribeiro Alvo

    Phoenix Guest

    Luis A. Afonso

    Your code forgotten to put the variable (w) as double precision.

    DEFDBL S,X,w

    This makes much difference.

    Or, if you prefer use

    x=s + n * x-INT(s + n * x)

    without the temp variable (w)
     
    Phoenix, Aug 31, 2009
    #21
    1. Advertisements

  2. ************************************
    Phoemix said:

    Date: Aug 30, 2009 8:07 PM
    Author: Phoenix
    Subject: Re: Very fast floating point PRNG
    Your code forgotten to put the variable (w) as double precision.
    DEFDBL S,X,w
    This makes much difference. Or, if you prefer use
    x=s + n * x-INT(s + n * x)
    without the variable x
    ****************************************

    My response:

    Thank you for your concern. On contrary you state I DO NOT believe that DEFBL S,X,w will change the situation. It seems rather an insinuation in order to *save* in extremis the damned RNG.
    Surely I´ll check.

    Luis A. Afonso
     
    Luis A. Afonso, Aug 31, 2009
    #22
    1. Advertisements

  3. Ribeiro Alvo

    Vermute Guest

    Thomas Richter

    In your toy-example we have 16 iterations (i).

    x = 2^2
    n = 2^2
    i = x * n
    i = 2^2 * 2^2
    i = 2^4
    i = 16

    Ok.

    Let
    x = 2^64
    and
    n = 2^64

    with this we have

    i = 2^64 * 2^64

    i = 2^128

    No.

    But the algorithm have s.
    If we do s = 2^64 n=2^64 and x=2^64

    We have i = 2^64 + 2^128 or periods


    The values of s, n and x, can be bigger, dependending only on precission you can use in the machine.
    Pease read Periodicity and Precission on the page.
     
    Vermute, Aug 31, 2009
    #23
  4. **************************
    _________________________content/class
    ________chi-2___________min_______max______
    ___1____1478___________486_______1099______
    ___2____1474___________508_______1110______
    ___3____1461___________493_______1099______
    ___4____69.4E6__________0______105000______
    ___5____1498___________494_______1094______
    ___6____1450___________485_______1093______
    ___7____1450___________485_______1093______
    ___8____1456___________501_______1098______
    ___9____1478___________505_______1098______
    __10____1491___________499_______1106______
    __11____1553___________466_______1100______
    __12____1482___________477_______1090______
    __13____34.2E6__________0 _______62501______
    __14____1481___________485_______1097______
    __15____1537___________446_______1101______
    __16____1557___________490_______1116______
    __17____1507___________513_______1131______
    __18____1649___________488_______1127______
    __19____69.4E6__________0______125000______
    __20____1552___________469_______1104______

    On contrary I supposed the statement DEFDBL x modifies the outputs; however there still are a great possibility to provide very large Test Value as at 4, 13 and 19 trials above The general conclusion that the RNG is invalid in practical grounds is still true. The existence of empty cases is the most evident feature.

    Luis A. Afonso



    REM "puff"
    CLS
    DEFDBL S, W
    DEFDBL X
    DEFLNG N, Q
    DIM h(1000), chi(1000)
    q = 1000000: REM one million
    seed = TIMER - INT(TIMER)
    s = seed / 7
    mxxu = 0: ixu = 10 ^ 10
    FOR n = 1 TO q: w = s + n * x
    x = w - INT(w)
    xu = INT(1000 * x + .5)
    h(xu) = h(xu) + 1
    NEXT n
    xptd = q / 1001
    FOR j = 0 TO 1000
    chi(j) = (h(j) - xptd) ^ 2 / xptd
    PRINT USING "#### ######## "; j; h(j);
    IF h(j) > mxxu THEN mxxu = h(j)
    IF h(j) < ixu THEN ixu = h(j)
    chi = chi + chi(j)
    NEXT j
    PRINT
    PRINT chi
    PRINT ixu, mxxu
     
    Luis A. Afonso, Aug 31, 2009
    #24
  5. Ribeiro Alvo

    Vermute Guest

    Luis...

    That's funny. With your exactly code, I have this very normal results:

    _________________________content/class
    _____________chi-2_________min____max______

    ---1 - 1388.4794921875 --- 503 - 1096
    ---2 - 1496.721923828125 - 473 - 1111
    ---3 - 1442.65673828125 -- 515 - 1110
    ---4 - 1439.90380859375 -- 512 - 1127
    ---5 - 1504.160888671875 - 502 - 1128
    ---6 - 1468.486694335938 - 511 - 1117
    ---7 - 1501.528076171875 - 479 - 1102
    ---8 - 1486.96142578125 -- 501 - 1098
    ---9 - 1584.62158203125 -- 464 - 1090
    --10 - 1433.734008789062 - 454 - 1086
    --11 - 1550.087768554688 - 482 - 1134
    --12 - 1488.90625 -------- 482 - 1098
    --13 - 1515.394409179688 - 511 - 1091
    --14 - 1459.79248046875 -- 511 - 1091
    --15 - 1489.69775390625 -- 472 - 1101
    --16 - 1517.39697265625 -- 491 - 1100
    --17 - 1450.38134765625 -- 508 - 1111
    --18 - 1410.481201171875 - 537 - 1090
    --19 - 1571.742553710938 - 454 - 1101
    --20 - 1541.77783203125 -- 470 - 1102
    --21 - 1570.730712890625 - 491 - 1093
    --22 - 1539.788330078125 - 497 - 1104
    --23 - 1520.435668945312 - 483 - 1093
    --24 - 1498.64404296875 -- 484 - 1095
    --25 - 1416.042846679688 - 501 - 1107
    --26 - 1556.181640625 ---- 456 - 1102
    --27 - 1483.075927734375 - 487 - 1105
    --28 - 1391.3681640625 --- 518 - 1093
    --29 - 1503.323120117188 - 481 - 1092
    --30 - 1485.71337890625 -- 506 - 1092
    --31 - 1489.336547851562 - 499 - 1101
    --32 - 1483.072143554688 - 468 - 1116
    --33 - 1544.552734375 ---- 475 - 1120
    --34 - 1582.952514648438 - 471 - 1095
    --35 - 1526.141723632812 - 483 - 1094
    --36 - 1534.191528320312 - 504 - 1112
    --37 - 1409.592407226562 - 499 - 1088
    --38 - 1492.681640625 ---- 473 - 1115
    --39 - 1468.852294921875 - 490 - 1095
    --40 - 1494.030517578125 - 487 - 1090

    I do it more times, and don't find strange values.
    The RND function built in QBASIC, is an LCG with single precision and don't pass tests with ENT and Diehard.
    This algorithm pass ENT and Diehard battery of tests.

    If your values still different from my, I suggest a neutral person to verify that.
     
    Vermute, Aug 31, 2009
    #25
  6. As you had shown the proposed RNG is fast and trusty as well: so one is, possibly, before a bug in my software. A third opinion concerning the RNG quality is no more needed, I presume.
    Thank you by your concern.

    Luis A. Afonso
     
    Luis A. Afonso, Sep 1, 2009
    #26
  7. Ribeiro Alvo

    Vermute Guest

    Thank you (obrigado) Luis Afonso
     
    Vermute, Sep 1, 2009
    #27
  8. Oh, sure enough it can be bigger, depending on the precision of the machine, I certainly don't disagree with that.
    As far as I read it, s remains constant throughout the iteration, so you don't need to consider it as a state variable.

    Nevertheless, it is certainly not aperiodic.

    So long,
    Thomas
     
    Thomas Richter, Sep 1, 2009
    #28
  9. Ribeiro Alvo

    AverageJoe Guest

    /*
    rng_alvo.cpp - v1.00 - by AverageJoe

    Implementation of Alvo's RNG in conventional C++
    Added: Seed() via time(0) etc.
    The RNG returns floating point (float, double, long double) random numbers
    between 0 and 1 (excluding I think)

    Tested in g++
    g++ -Wall -O2 -o rng_alvo.exe rng_alvo.cpp -lm
    */

    #include <cstdio>
    #include <cstdlib>
    #include <ctime>
    #include <cstring>


    template <typename RealType = double> class Rng01_Alvo
    {
    public:
    Rng01_Alvo()
    : nbytes(sizeof(RealType))
    {
    uLoop = 0u;
    realX = 0.0;
    Seed(4711u);
    }

    void Seed(unsigned uSeed)
    { // uSeed=0 means random seed
    if (!uSeed)
    uSeed = time(0);
    char sz[2u + nbytes * 3u + 1u];
    strcpy(sz, "0.");
    unsigned ix = 2u;
    unsigned t = uSeed + !!uSeed;
    for (unsigned i = 0u; i < nbytes; ++i)
    {
    unsigned b = t & 0xffu;
    t >>= sizeof(unsigned char);

    char tmp[3 + 1];
    sprintf(tmp, "%u", b);
    for (int j = int(strlen(tmp)) - 1; j >= 0; --j)
    sz[ix++] = tmp[j];
    }
    sz[ix] = 0;
    realSeed = atof(sz) / 3.0;

    // printf("%u %s %lf\n", uSeed, sz, double(realSeed));

    // "calibrate" by eating the first two generated random numbers
    for (unsigned i = 0u; i < 2u; ++i)
    Rand();
    }

    RealType Rand()
    {
    realX = realSeed + realX * double(uLoop++);
    return realX -= unsigned(realX);
    }

    private:
    const unsigned nbytes;
    unsigned uLoop;
    RealType realSeed, realX;
    };


    int test_Rng01_Alvo()
    {
    Rng01_Alvo<> R; // using default type (ie. double)
    // Rng01_Alvo<float> R; // using float type
    // Rng01_Alvo<double> R; // using double
    // Rng01_Alvo<long double> R; // using long double

    R.Seed(0); // passing 0 or time(0) means random seed

    for (int i = 1; i <= 50; ++i)
    printf("%6d %.32lf\n", i, double(R.Rand()));

    return 0;
    }

    int main(int argc, char* argv[])
    {
    return test_Rng01_Alvo();
    }


    /*
    Sample output:

    1 0.98033557007968663565833367101732
    2 0.22875216131254050289101087400923
    3 0.69958727878935034194540776297799
    4 0.22418526243955971821719685976859
    5 0.12055641043058507033691739707137
    6 0.41208200478937295763159909256501
    7 0.44789435880783690091533344457275
    8 0.25391225780031323466801040922292
    9 0.37402140451786614061546742959763
    10 0.18935256964636715260041910369182
    11 0.97801028596324746100520997060812
    12 0.99722860036281923967749207804445
    13 0.28273191671921515943921576763387
    14 0.13351689514076503506601056869840
    15 0.73037842344787795578753275549388
    16 0.20426435642378848456246487330645
    17 0.41236143938987701584863998505170
    18 0.98666453810152465742078220500844
    19 0.03369414914764978785655102910823
    20 0.10203974131005244530001618841197
    21 0.77602810235536601179973104081000
    22 0.72777686979825761959261853917269
    23 0.24927292566880543400742453741259
    24 0.05744330397185404279980502906255
    25 0.93548682212602052299388333267416
    26 0.45619215256870249142195916647324
    27 0.01283888797518206725811751311994
    28 0.72507983814475140960098542564083
    29 0.52962913154615243183087613942917
    30 0.80483565193713757501825512008509
    31 0.69148647677678487522001660181559
    32 0.52910099810157795729992358246818
    33 0.37471024257091378739659148777719
    34 0.19135328603791745827322756667854
    35 0.59849918035497073898199005270726
    36 0.66854234475796614489695457450580
    37 0.06876810123275312047752549915458
    38 0.14656646145698526417788798426045
    39 0.48286569210748842007063785786158
    40 0.09029907153610938763677040697075
    41 0.30023345850292038061724042563583
    42 0.83757994354552445237516167253489
    43 0.85575171400823213385677945552743
    44 0.54740486930101761586797692871187
    45 0.60250803736295288626223509709234
    46 0.84996814169879875144886227644747
    47 0.82548139585404489704245634129620
    48 0.42662549947039707376461592502892
    49 0.51160028337475516035937062042649
    50 0.44188932977613082364598540152656
    */
     
    AverageJoe, Sep 2, 2009
    #29
  10. Ribeiro Alvo

    Vermute Guest

    AverageJoe

    Thank you for your great job.

    I don't no nothing about C/C++ programming language, but let me ask you.

    In your code, you avoid the value s=0.5?

    On the algorithm: x_n+1=frac(s+nx_n), the (s) value, never can be 0.5.

    See Seed mechanism at http://www.number.com.pt/index.html
     
    Vermute, Sep 2, 2009
    #30
  11. Ribeiro Alvo

    AverageJoe Guest

    Hi Vermute / Alvo (?),

    s = ...
    s = s / 3.0;

    s is a constant. It is computed only in the Seed() function.
    There, first a fractiun between 0 and 1 (excluding 0 and 1) is computed,
    and then that divided by 3.0. By this mechanism, s never can become 0.5.

    Seed(u) has same interface like the default C/C++ seed(u).
    Ie. accepting an unsigned integer provided by the user.
    But since in this algorithm the value 0 cannot be used, then
    passing 0 to Seed() is prevented by convention as follows:
    passing a 0 means taking the current time(0) value as input.
    Should even time(0) return 0 then it will made 1 in the line
    "unsigned t = uSeed + !!uSeed;"

    And the default internal Seed parameter is 4711 (see constructor;
    it is arbitrarily, one could have taken any value > 0).

    Ie. calling Seed() is optional. In a program one should call
    Seed() normally only once at program start.

    Hope this helps.
     
    AverageJoe, Sep 2, 2009
    #31
  12. Ribeiro Alvo

    AverageJoe Guest

    BTW, I myself haven't done any performance measurements,
    because you already did such tests and published the results
    at the web site http://www.number.com.pt/index.html
    Of cource, when making performance measurements then
    one should measure only the Rand() calls, not Seed().
     
    AverageJoe, Sep 2, 2009
    #32
  13. Ribeiro Alvo

    AverageJoe Guest

    <...>

    Bugfixes in v1.01:

    - "t >>= sizeof(unsigned char);" changed to
    "t >>= 8u;"

    - "realX = realSeed + realX * double(uLoop++);" changed to
    "realX = realSeed + realX * RealType(uLoop++);"
     
    AverageJoe, Sep 2, 2009
    #33
  14. Ribeiro Alvo

    AverageJoe Guest

    I did some performance tests against this method:

    double random_uniform_0_1b()
    { // 0 < r < 1
    // this uses the C library random number generator
    unsigned rnd = rand();
    while (!rnd || (rnd == RAND_MAX))
    rnd = rand();
    return double(rnd) / double(RAND_MAX);
    }

    Alvo's algorithm is about 25% faster than this.
     
    AverageJoe, Sep 2, 2009
    #34
  15. Ribeiro Alvo

    Vermute Guest

    AverageJoe:

    Vermute = Phoenix = Alvo :)

    On the page Seed mechanism, I only write on the problem of s=0.5.

    But the real problem is wend s=2**-n

    n=1 => s=0.5
    n=2 => s=0.25
    n=3 => s=0125
    ...
    and so on.

    A way to prevent this and s=0, is in BASIC code

    newseed:
    s = TIMER-INT(timer)
    if s=0 or -log(s)/log(2)=int(-log(s)/log(2)) then newseed

    Divide seed by 3.0 is not sufficent.

    See:
    if seed=0.75

    s=0.75/3.0

    s=0.25 and

    -log(0.25)/log(2)=2 = int(-log(s)/log(2))

    With this seed machanism, I think, we have no problem on small periods.
    I don't write this on the page yet, but I will, later

    About the test speed.

    I ask you if the time value, increase with the precision you use?
     
    Vermute, Sep 2, 2009
    #35
  16. Ribeiro Alvo

    AverageJoe Guest

    Ok, I understand.
    The timings for "float" and "double" are practically the same.
    "long double" takes about 14% more than double or float.

    (this is on i686 CPU (ie. AMD 32bit), Linux, and g++;
    generating 100 million Rand()'s, and using 2 clock() calls
    outside the loop).
     
    AverageJoe, Sep 2, 2009
    #36
  17. Ribeiro Alvo

    AverageJoe Guest

    Correction:
    the timings for float, double, and long double are practically all the same.
     
    AverageJoe, Sep 2, 2009
    #37
  18. Ribeiro Alvo

    Vermute Guest

    AverageJoe:
    Thats good news.
    This is not correct, becase can not avoid s=0.375, s=0.875 s=0.75


    The correction is:

    newseed:
    s = TIMER-INT(timer)
    if s/2^-30 = int(s/2^-30) then newseed

    With that we avoid 0, 1, 0.5, 0.25, 0.375, 0.75, 0.875 etc.
    I think the -30 value is a good number. But if we want, we can put a smaller value, deppending on the precision we use.

    Can you give the value for the pseudo-random numbers generated per second, on your system?

    If you agree, I put your C++ code on the page in plain text Made by (your name).ok.
     
    Vermute, Sep 2, 2009
    #38
  19. Ribeiro Alvo

    AverageJoe Guest

    Alvo, I'm sorry, I don't have time and have lost interest in this.
    I think this is getting messy.
    I would just recommend to use simply among the following 4 functions.
    They use the well tested standard rand(), here for generating uniform
    floating point random numbers between 0 and 1, with or without 0 and 1.
    And with them there are no such problems with seed.

    This is C++:


    #include <cstdlib>
    #include <ctime>

    // - Use srand(uSeed) once at pgmstart to set the seed, for example "srand(time(0));"

    // - Set this to a real type, ie. float, double, or long double
    typedef double REALTYPE;


    inline REALTYPE rand01_yy()
    { // 0 <= rnd <= 1
    return REALTYPE(rand()) / REALTYPE(RAND_MAX);
    }

    inline REALTYPE rand01_yn()
    { // 0 <= rnd < 1
    unsigned r = 0u;
    while (r == RAND_MAX)
    r = rand();
    return REALTYPE(r) / REALTYPE(RAND_MAX);
    }

    inline REALTYPE rand01_ny()
    { // 0 < rnd <= 1
    unsigned r = 0u;
    while (!r)
    r = rand();
    return REALTYPE(r) / REALTYPE(RAND_MAX);
    }

    inline REALTYPE rand01_nn()
    { // 0 < rnd < 1
    unsigned r = 0u;
    while (!r || (r == RAND_MAX))
    r = rand();
    return REALTYPE(r) / REALTYPE(RAND_MAX);
    }
     
    AverageJoe, Sep 2, 2009
    #39
    1. Advertisements

Ask a Question

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

You'll need to choose a username for the site, which only take a couple of moments (here). After that, you can post your question and our members will help you out.