Prev: Square Roots and Sequences
Next: Cycles of 1D Dynamical System
From: Jonathan Lee on 14 Apr 2010 15:55
On Apr 14, 2:53 pm, Jonathan Lee <cho...(a)shaw.ca> wrote:
> On Apr 14, 2:06 pm, Dann Corbit <dcor...(a)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 = 123456789UL, unsigned long = 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 != 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 = 0; i < qlen; i++) {
double s = 0.0, t = 1.0;

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

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

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

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

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

// return Q[indx] - zy (mod 1)
t = 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 = 0; i < 1000000000; i++)
t = krng();
printf ("%.16f\n", krng());
}

--Jonathan
From: rossum on 15 Apr 2010 10:04
On Wed, 14 Apr 2010 03:41:05 -0700 (PDT), geo <gmarsaglia(a)gmail.com>
wrote:

>1. Simple and very fast, some 70 million double-precision
> uniform(0,1] random variables/second, in IEEE 754 standard
> form
This gives an issue with the Java conversion. Java expects the return
from Random.nextDouble() to be in the range [0.0, 1.0), while this
code produces numbers in the range (0.0, 1.0].

I know I can make the switch by using:

return (1.0 - DoubleKISS());

Does anyone have any alternative suggestions that are faster and do
not break George's excellent RNG.

rossum

From: steve on 15 Apr 2010 10:35
On Apr 15, 7:04 am, rossum <rossu...(a)coldmail.com> wrote:
> On Wed, 14 Apr 2010 03:41:05 -0700 (PDT), geo <gmarsag...(a)gmail.com>
> wrote:
>
> >1. Simple and very fast, some 70 million double-precision
> >   uniform(0,1] random variables/second, in IEEE 754 standard
> >   form
>
> This gives an issue with the Java conversion.  Java expects the return
> from Random.nextDouble() to be in the range [0.0, 1.0), while this
> code produces numbers in the range (0.0, 1.0].
>
> I know I can make the switch by using:
>
>   return (1.0 - DoubleKISS());
>
> Does anyone have any alternative suggestions that are faster and do
> not break George's excellent RNG.

Not sure if it would be faster, but one could do

return(nextafter(DoubleKISS(), -1.));

Two items that George did not comment on are subnormal
numbers and how he determined the distribution is normal
over (0,1]. Does his generator create subnormals? If
I understand Figure D-1 in Goldberg, "What every computer
scientist should know about floating-point arithmetic,"
the distribution cannot be normal.
From: rossum on 16 Apr 2010 10:27
On Wed, 14 Apr 2010 03:41:05 -0700 (PDT), geo <gmarsaglia(a)gmail.com>
wrote:

>for(i=0;i<1000000000;i++) t=dUNI();
>printf("%.16f\n",dUNI());
>}
>
>/*
>Output, after 10^9 random uniforms:
> 0.6203646342357479
Small mistake here George. You get 0.6203646342357479 after (10^9) +
1 random uniforms because there is an extra call to dUNI() in the
printf statement.

Took me a while to find, as my code had the equivalent of:
printf("%.16f\n",t);
and I couldn't find my mistake anywhere in the dUNI() code.

I get:
0.7976863625406643 after 1,000,000,000 calls
0.6203646342357480 after 1,000,000,001 calls

rossum

From: geo on 16 Apr 2010 12:41
On Apr 15, 10:35 am, steve <kar...(a)comcast.net> wrote:
> On Apr 15, 7:04 am, rossum <rossu...(a)coldmail.com> wrote:
>
>
>
>
>
> > On Wed, 14 Apr 2010 03:41:05 -0700 (PDT), geo <gmarsag...(a)gmail.com>
> > wrote:
>
> > >1. Simple and very fast, some 70 million double-precision
> > >   uniform(0,1] random variables/second, in IEEE 754 standard
> > >   form
>
> > This gives an issue with the Java conversion.  Java expects the return
> > from Random.nextDouble() to be in the range [0.0, 1.0), while this
> > code produces numbers in the range (0.0, 1.0].
>
> > I know I can make the switch by using:
>
> >   return (1.0 - DoubleKISS());
>
> > Does anyone have any alternative suggestions that are faster and do
> > not break George's excellent RNG.
>
> Not sure if it would be faster, but one could do
>
>    return(nextafter(DoubleKISS(), -1.));
>
> Two items that George did not comment on are subnormal
> numbers and how he determined the distribution is normal
> over (0,1].  Does his generator create subnormals?  If
> I understand Figure D-1 in Goldberg, "What every computer
> scientist should know about floating-point arithmetic,"
> the distribution cannot be normal.- Hide quoted text -
>
> - Show quoted text -

Users who insist on getting every possible float
in the unit interval could form dUNI()*.5^j;
with j taking values
0, 1, 2, 3,...
with probabilities
q, qp, qp^2, qp^3, ...
and p = 2^{-53) and q = 1-p..
Even at the rate of 70 million per second,
you wouild have to wait an average of over
four years before you had to multiply by .5,
8 years to multiply by .5^2, etc..

George Marsaglia
First  |  Prev  |  Next  |  Last
Pages: 1 2 3
Prev: Square Roots and Sequences
Next: Cycles of 1D Dynamical System