Path: g2news1.google.com!postnews.google.com!k11g2000vbg.googlegroups.com!not-for-mail
From: Jonathan Lee <cho...@shaw.ca>
Newsgroups: sci.math,comp.lang.c,comp.lang.c++
Subject: Re: RNGs: A double KISS
Date: Wed, 14 Apr 2010 12:55:42 -0700 (PDT)
Organization: http://groups.google.com
Lines: 120
Message-ID: <8ee39635-b420-4592-a84d-e2c38788f35b@k11g2000vbg.googlegroups.com>
References: <a054eef2-5d89-47ff-bda0-ef5ce2325a86@e7g2000yqf.googlegroups.com> 
	<MPG.262fa3d3db369ed09896e8@news.eternal-september.org> <b36dec89-1f40-4849-8fdd-c3b619c60fcf@e21g2000vbb.googlegroups.com>
NNTP-Posting-Host: 64.231.111.105
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: posting.google.com 1271274943 5590 127.0.0.1 (14 Apr 2010 19:55:43 GMT)
X-Complaints-To: groups-abuse@google.com
NNTP-Posting-Date: Wed, 14 Apr 2010 19:55:43 +0000 (UTC)
Complaints-To: groups-abuse@google.com
Injection-Info: k11g2000vbg.googlegroups.com; posting-host=64.231.111.105; 
	posting-account=-ysU1woAAABnrwtj-u1ZwHgijoTNBPAW
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/5.0 (X11; U; Linux x86_64; en-US) AppleWebKit/533.2 
	(KHTML, like Gecko) Chrome/5.0.342.9 Safari/533.2,gzip(gfe)

On Apr 14, 2:53=A0pm, Jonathan Lee <cho...@shaw.ca> wrote:
> On Apr 14, 2:06=A0pm, Dann Corbit <dcor...@connx.com> wrote:
>
> Thanks for sharing.
>
> Isn't this dangerous for y:
>

Just confirmed: if y isn't masked to 32-bits the output
is wrong.

Anyway, here's my take on the code:

// kiss2.hpp
-------------------------------------------------------------------
#include <cstddef>

class kiss2 {
    double Q[1220];         // meh.. I guess this could be a
std::vector
    const std::size_t qlen; // length of Q array
    std::size_t indx;
    double c;               // current CSWB
    double zc;              // current SWB `borrow`
    double zx;              // SWB seed1
    double zy;              // SWB seed2
public:
    kiss2(unsigned long =3D 123456789UL, unsigned long =3D 362436069UL);
    double operator()();    // For getting a random number
};

// kiss2.cpp
-------------------------------------------------------------------
//#include "kiss2.hpp"
#include <cfloat>
using std::size_t;
/**
\todo? Replace constant initialization of indx and qlen w/ sizeof Q/
sizeof Q[0]
\todo Magic numbers make me uneasy (re: zx, zy)
*/
kiss2::kiss2(
    unsigned long x,  // seed1 value
    unsigned long y   // seed2 value
):  qlen(1220), indx(1220), c(0.0), zc(0.0),
    zx(5212886298506819.0 / 9007199254740992.0),
    zy(2020898595989513.0 / 9007199254740992.0)
{
    #if (FLT_RADIX !=3D 2 || DBL_MIN_EXP > -53)
    #error "Machine double not supported"
    #endif

    // fill in Q[i] "using 9th bit from Cong+Xorshift"
    for (size_t i =3D 0; i < qlen; i++) {
        double s =3D 0.0, t =3D 1.0;

        // Build Q[i] one bit at a time
        for (int j =3D 0; j < 52; j++) {
            t *=3D 0.5; // make t=3D.5/2^j
            x =3D (69069 * x + 123);
            y =3D (y ^ (y << 13)) & 0xFFFFFFFFUL;
            y =3D (y ^ (y >> 17)) & 0xFFFFFFFFUL;
            y =3D (y ^ (y << 5)) & 0xFFFFFFFFUL;

            if (((x + y) >> 23) & 1UL) // conditionally set bit of s
                s +=3D t;
        }
        Q[i] =3D s;
    }
}

/**
\todo Separate refill function?
\todo I guess we sorta need a guarantee that qlen >=3D 30
*/
double kiss2::operator()() { // Takes 14 nanosecs, Intel Q6600,2.40GHz
    static const double cc =3D 1.0 / 9007199254740992.0; // 2^-53

    double t; // t: first temp, then next CSWB value
    // First get zy as next lag-2 SWB
    t =3D zx - zy - zc;
    zx =3D zy;
    if (t < 0) {
        zy =3D t + 1.0, zc =3D cc;
    } else zy =3D t, zc =3D 0.0;

    // Then get t as the next lag-1220 CSWB value
    if (indx >=3D qlen) { // refill Q[n] via Q[n-1220]-Q[n-1190]-c
        for (size_t i =3D 0; i < qlen; i++) {
            size_t j =3D (i < 30) ? i + (qlen - 30) : i - 30;
            t =3D Q[j] - Q[i] + c; // Get next CSWB element
            if (t > 0) {
                t =3D t - cc, c =3D cc;
            } else t =3D t - cc + 1.0, c =3D 0.0;
            Q[i] =3D t;
        }
        indx =3D 0; // reset indx
    }

    // return Q[indx] - zy (mod 1)
    t =3D Q[indx++] - zy;
    return (t < 0.0 ? 1.0 + t : t);
}

// main.cpp
--------------------------------------------------------------------
#include <cstdio>
#include <cstdlib>
// #include "kiss2.hpp"
int main(void)
{ /* You must provide at least two 32-bit seeds */
    kiss2 krng; // Use the default seeds
    double t;
    int i;
    for (i =3D 0; i < 1000000000; i++)
        t =3D krng();
    printf ("%.16f\n", krng());
}

--Jonathan