Seeding a user supplied random number generator in R

Question

I'm having some trouble with seeding a user defined RNG in R. It seems that

set.seed(123, kind='user', normal.kind='user')

Does not actually pass 123 to the user defined RNG initialization.

I went back to the documentation available at ?Random.user and tried the example code given there, with the minor modification that I print the seed passed to the user_unif_init function (full code below).

Steps to reproduce:

1. Paste code below in urand.c
2. Run R CMD SHLIB urand.c
3. Open R

4. Run the following commands:

   > dyn.load('urand.so')
   > set.seed(123, kind='user', normal.kind='user')
   Received seed: 720453763
   Received seed: 303482705 // any other numbers than 123
   

Here's the full code I used in urand.c:

// ##  Marsaglia's congruential PRNG

#include <stdio.h>
#include <R_ext/Random.h>

static Int32 seed;
static double res;
static int nseed = 1;

double * user_unif_rand()
{
    seed = 69069 * seed + 1;
    res = seed * 2.32830643653869e-10;
    return &res;
}

void  user_unif_init(Int32 seed_in) {
    printf("Received seed: %u\n", seed_in);
    seed = seed_in;
}
int * user_unif_nseed() { return &nseed; }
int * user_unif_seedloc() { return (int *) &seed; }

/*  ratio-of-uniforms for normal  */
#include <math.h>
static double x;

double * user_norm_rand()
{
    double u, v, z;
    do {
        u = unif_rand();
        v = 0.857764 * (2. * unif_rand() - 1);
        x = v/u; z = 0.25 * x * x;
        if (z < 1. - u) break;
        if (z > 0.259/u + 0.35) continue;
    } while (z > -log(u));
    return &x;
}

Any help would be greatly appreciated!

Answer 1

It appears that R scrambles the user supplied seed in RNG.c as follows:

for(j = 0; j < 50; j++)
    seed = (69069 * seed + 1)

(link to source)

Trying to unscramble this would be a way to get the original seed back.

UPDATE

Unscrambling can be done through the multiplicative inverse of 69069 as follows:

Int32 unscramble(Int32 scrambled)
{
        int j;
        Int32 u = scrambled;
        for (j=0; j<50; j++) {
                u = ((u - 1) * 2783094533);
        }
        return u;
}

Plugging this in my user_unif_init() function solves the problem.