Random Numbers with OpenCL using Random123

Question

I have been looking at this lib Random123 and associated quote:

One mysterious man came to my booth and asked what I knew about generating random numbers with OpenCL. I told him about implementations of the Mersenne Twister, but he wasn't impressed. He told me about a new technical paper that explains how to generate random numbers on GPUs by combining integer counters and block ciphers. In reverential tones, he said that counter-based random number generators (CBRNGs) produce numbers with greater statistical randomness than the MT and with much greater speed.

I was able to get a demo running using this kernel:

__kernel void counthits(unsigned n, __global uint2 *hitsp) {
    unsigned tid = get_global_id(0);
    unsigned hits = 0, tries = 0;
    threefry4x32_key_t k = {{tid, 0xdecafbad, 0xfacebead, 0x12345678}};
    threefry4x32_ctr_t c = {{0, 0xf00dcafe, 0xdeadbeef, 0xbeeff00d}};
    while (tries < n) {
        union {
            threefry4x32_ctr_t c;
            int4 i;
        } u;
        c.v[0]++;
        u.c = threefry4x32(c, k);
        long x1 = u.i.x, y1 = u.i.y;
        long x2 = u.i.z, y2 = u.i.w;
        if ((x1*x1 + y1*y1) < (1L<<62)) {
            hits++;
        }
        tries++;
        if ((x2*x2 + y2*y2) < (1L<<62)) {
            hits++;
        }
        tries++;
    }
    hitsp[tid].x = hits;
    hitsp[tid].y = tries;
}

My questions are now, will this not generate the same random numbers every time its run, a random number is based on the global id ? How can I generate new random numbers each time. Possible to provide a seed as a parameter for the kernel and then use that somehow?

Anyone who have been using this lib and can give me some more insight in the use of it?

Answer 1

Yes. The example code generates the same sequences of random numbers every time it is called.

To get different streams of random numbers, just initialize any of the values k[1..3] and/or c[1..3] differently. You can initialize them from command line arguments, environment variables, time-of-day, saved state, /dev/urandom, or any other source. Just be aware that:

a) if you initialize all of them exactly the same way in two different runs, then those two runs will get the same stream of random numbers

b) if you initialize them differently in two different runs, then those two runs will get different streams of random numbers.

Sometimes you want property a). Sometimes you want property b). Take a moment to think about which you want and be sure that you're doing what you intend.

More generally, the functions in the library, e.g., threefry4x32, have no state. If you change any bit in the input (i.e., any bit in any of the elements of c or k), you'll get a completely different random, statistically independent, uniformly distributed output.

P.S. I'm one of the authors of the library and the paper "Parallel Numbers: As Easy as 1, 2, 3": http://dl.acm.org/citation.cfm?id=2063405

If you're not a subscriber to the ACM digital library, the link above may hit a pay-wall. Alternatively, you can obtain the paper free of charge by following the link on this page:

http://www.thesalmons.org/john/random123/index.html

Answer 2

I can't help you with the library per se, but I can tell you that the most common way to generate random numbers in OpenCL is to save some state between calls to the kernel.

Random number generators usually use a state, from which a new state and a random number are generated. In practice, this isn't complicated at all: you just pass an extra array that holds state. In my codes, I implement random numbers as follows:

uint rand_uint(uint2* rvec) {  //Adapted from http://cas.ee.ic.ac.uk/people/dt10/research/rngs-gpu-mwc64x.html
    #define A 4294883355U
    uint x=rvec->x, c=rvec->y; //Unpack the state
    uint res = x ^ c;          //Calculate the result
    uint hi = mul_hi(x,A);     //Step the RNG
    x = x*A + c;
    c = hi + (x<c);
    *rvec = (uint2)(x,c);      //Pack the state back up
    return res;                //Return the next result
    #undef A
}
inline float rand_float(uint2* rvec) {
    return (float)(rand_uint(rvec)) / (float)(0xFFFFFFFF);
}
__kernel void my_kernel(/*more arguments*/ __global uint2* randoms) {
    int index = get_global_id(0);
    uint2 rvec = randoms[index];

    //Call rand_uint or rand_float a number of times with "rvec" as argument.
    //These calls update "rvec" with new state, and return a random number

    randoms[index] = rvec;
}

. . . then, all you do is pass an extra array that holds the RNG's state into random. In practice, you'll want to seed this array differently for each work item.