jeudi 1 juillet 2021

Creating Numerical Recipes random number generator in Python

I would like to generate random numbers from the U[0,1] distribution. I am trying to recreate the random number generator from Numerical Recipes Chapter 7.1. Unfortunately, my C++ is not particularly strong and so I am having some issues following the code! I was wondering if someone could help clarify my understanding of really what is happening here. In C++ the code is provided as:

struct Ran {
    Ullong u, v, w;
    Ran(Ullong j) : v(4101842887655102017LL), w(1) {
        u = j ^ v; int64();
        v = u; int64();
        w = v; int64();
    }
    inline Ullong int64() {
        u = u * 2862933555777941757LL + 7046029254386353087LL;
        v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
        w = 4294957665U*(w & 0xffffffff) + (w >> 32);
        Ullong x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
        return (x + v) ^ w;
    }
    inline Doub doub() { return 5.42101086242752217E-20 * int64(); }
    inline Uint int32() { return (Uint)int64(); }
};

My understanding is this:

Ran(Ullong j) : v(4101842887655102017LL), w(1) {
        u = j ^ v; int64();
        v = u; int64();
        w = v; int64();
    }

This is our initialisation. j here is our seed. v is set to be a really big prime number. ^ is a bit-wise XOR operator and so this provides our first level of randomness as we shuffle u, v and w.

Beyond this I start getting confused.

inline Ullong int64() {
        u = u * 2862933555777941757LL + 7046029254386353087LL;
        v ^= v >> 17; v ^= v << 31; v ^= v >> 8;
        w = 4294957665U*(w & 0xffffffff) + (w >> 32);
        Ullong x = u ^ (u << 21); x ^= x >> 35; x ^= x << 4;
        return (x + v) ^ w;
    }

In particular, I don't understand what

    v ^= v >> 17;

does. I think it is bit-shifting v by 17 digits to the right and then performing XOR of the original v against this and then saving this output as v.

I also don't understand

w = 4294957665U*(w & 0xffffffff) + (w >> 32);

at all. A quick search here says 0xffffffff is either -1 or a big number but I still don't really understand what this line is doing (other than some bit-wise operations).

Finally, I believe that this line:

inline Doub doub() { return 5.42101086242752217E-20 * int64(); }

is returning the random number between 0 and 1.

If anyone could correct, clarify or confirm what I have here I would greatly appreciate it.




Aucun commentaire:

Enregistrer un commentaire