Finding inverse operation to George Marsaglia's XorShift RNG

Question

Abstract

Hi, suppose you have 128 bit automata (represented by four 32 bit words X, Y, Z, W) that changes it's state according to a following rule:

X = ...
Y = ...
Z = ...
W = ...

void next()
{
    var t = X ^ (X << 11);

    X = Y;
    Y = Z;
    Z = W;

    W = W ^ (W >> 19) ^ (t ^ (t >> 8));
}

^ - denotes binary XOR operation

<< - denotes binary shift left operation

>> - denotys binary shift right operation

It is guaranteed that the above automata generates no collisions i.e. each state is a result of one (and only one) previous state. It is also guaranteed that the above state machine produces 2^128 unique states.

Question

For any given state (X,Y,Z,W) produce inverse to next, (i.e. prev) operation that would revert the state to previous one.

In other words, if you have the following state (X=1, Y=2, Z=3, W=4) and will call next, the state will change to (X=2, Y=3, Z=4, W=2061), it is supposed that after calling prev the state should be equal again to (X=1, Y=2, Z=3, W=4).

P.S.

The next operation is one of the implementations to XorShift pseudorandom number generators that was discovered by George Marsaglia

https://en.wikipedia.org/wiki/Xorshift

The inverse to this operation would be very useful in general, consider the implications of Guid.Next(...), Guid.Prev(...) availability

Edit

I have somewhat improved Niklas B.'s original answer and ported result to C#, so here's the final piece of code, hope someone will benefit from Random.Next() and Random.Prev() operations:

public class Xor128
{
    public UInt32 X { get; set; }
    public UInt32 Y { get; set; }
    public UInt32 Z { get; set; }
    public UInt32 W { get; set; }

    public Xor128()
    {

    }

    public Xor128(UInt32 x, UInt32 y, UInt32 z, UInt32 w)
    {
        X = x;
        Y = y;
        Z = z;
        W = w;
    }

    //private UInt32 UnXorShl(UInt32 x, Int32 shift)
    //{
    //    for (var i = shift; i < 32; i <<= 1) {
    //        x ^= x << i;
    //    }

    //    return x;
    //}

    //private UInt32 UnXorShr(UInt32 x, Int32 shift)
    //{
    //    for (var i = shift; i < 32; i <<= 1) {
    //        x ^= x >> i;
    //    }

    //    return x;
    //}

    //public UInt32 Prev()
    //{
    //    var t = UnXorShr(W ^ Z ^ (Z >> 19), 8);

    //    W = Z;
    //    Z = Y;
    //    Y = X;
    //    X = UnXorShl(t, 11);

    //    return W;
    //}

    public UInt32 Prev()
    {
        var t = W ^ Z ^ (Z >> 19);

        t ^= t >> 8;
        t ^= t >> 16;

        W = Z;
        Z = Y;
        Y = X;

        t ^= t << 11;
        t ^= t << 22;

        X = t;

        return W;
    }


    public UInt32 Curr()
    {
        return W;
    }

    public UInt32 Next()
    {
        UInt32 t = X ^ (X << 11);

        X = Y;
        Y = Z;
        Z = W;

        return W = W ^ (W >> 19) ^ (t ^ (t >> 8));
    }
}

btw. Here's a swift version:

public class Xor128 {
    public var X: UInt32
    public var Y: UInt32
    public var Z: UInt32
    public var W: UInt32
    
    public convenience init(uuid: uuid_t) {
        let xa = (UInt32(uuid.0 ) << 24)
        let xb = (UInt32(uuid.1 ) << 16)
        let xc = (UInt32(uuid.2 ) << 8 )
        let xd = (UInt32(uuid.3 ) << 0 )

        let ya = (UInt32(uuid.4 ) << 24)
        let yb = (UInt32(uuid.5 ) << 16)
        let yc = (UInt32(uuid.6 ) << 8 )
        let yd = (UInt32(uuid.7 ) << 0 )
        
        let za = (UInt32(uuid.8 ) << 24)
        let zb = (UInt32(uuid.9 ) << 16)
        let zc = (UInt32(uuid.10) << 8 )
        let zd = (UInt32(uuid.11) << 0 )

        let wa = (UInt32(uuid.12) << 24)
        let wb = (UInt32(uuid.13) << 16)
        let wc = (UInt32(uuid.14) << 8 )
        let wd = (UInt32(uuid.15) << 0)
        
        self.init(
            x: xa + xb + xc + xd,
            y: ya + yb + yc + yd,
            z: za + zb + zc + zd,
            w: wa + wb + wc + wd
        )
    }
    
    public convenience init(uuid: UUID) {
        self.init(uuid: uuid.uuid)
    }
    
    public init(x: UInt32, y: UInt32, z: uint32, w: UInt32) {
        X = x
        Y = y
        Z = z
        W = w
    }
    
    @discardableResult
    public func next() -> UInt32 {
        let t = X ^ (X << 11);
        
        X = Y;
        Y = Z;
        Z = W;
        
        W = W ^ (W >> 19) ^ (t ^ (t >> 8))
        
        return W;
    }
    
    public var curr: UInt32 {
        return W
    }
    
    @discardableResult
    public func prev() -> UInt32 {
        var t = W ^ Z ^ (Z >> 19);
        
        t ^= t >> 8;
        t ^= t >> 16;
        
        W = Z;
        Z = Y;
        Y = X;
        
        t ^= t << 11;
        t ^= t << 22;
        
        X = t;
        
        return W;
    }
}

Interesting problem. What have you tried? How close to answering it did you get? You are supposed to try to solve it yourself before asking for answers! — tucuxi, Jul 20 '15 at 09:43
I don't think this will be possible, since bit-shifting operations will delete the bits shifted out of the variable range => irreversable — xXliolauXx, Jul 20 '15 at 09:47
On purely theoretic grounds, you can do exhaustive search, find all 2^128 states, and therefore map it all out. Takes some (imposibly long & large) time and space, but in theory it works. — tucuxi, Jul 20 '15 at 09:49
@tucuxi - I think you only need `2^32` storage - we know the previous value of `W` (and `Y` and `Z`). We can use that to compute `(t ^ (t >> 8))` with the current `W` value and formula. `t` only depends on `X`, so if we have a mapping of `X` values to `(t ^ (t >>8 ))` values we should be able to reverse it. — Damien_The_Unbeliever, Jul 20 '15 at 09:56
It cannot have 2^128 reachable states since X=Y=Z=W=0 will loop forever — joop, Jul 20 '15 at 10:07
You are right, (0,0,0,0) is the special case, but the most interesting thing is that state (1,0,0,0) is not reachable from any other state, the above algorithm gives that (1,0,0,0) is produced from (0,0,0,0) which is obviously an ill-defined point — Lu4, Jul 20 '15 at 22:01

Niklas B. · Accepted Answer · 2015-07-20T11:44:49.857

The basic building block you need is an algorithm to reverse the XOR with left shift operation f(x) = x ^ (x << s) for some s > 0. Given f(x), you already know the lower s bits of x directly.

You can reconstruct the rest of the bits iteratively from low to high, because you already know at each point the two bits that have been XORed to get the bit of f(x). Here's an example in Python:

def reverse_xor_lshift(y, shift, w=32):
    x = y & ((1<<shift) - 1)
    for i in range(w - shift):
        x |= (1 if bool(x & (1<<i)) ^ bool(y & (1<<(shift+i))) else 0)<<(shift+i)
    return x

Now the rest becomes rather easy. Note that I'm reusing the left shift reversal for the right shift analogue:

def reverse_bin(x, w=32):
    return int(bin(x)[2:].rjust(w, '0')[::-1], 2)

def reverse_xor_rshift(y, shift, w=32):
    # for simplicity, we just reuse reverse_xor_lshift here
    return reverse_bin(reverse_xor_lshift(reverse_bin(y), shift))

def forward(X, Y, Z, W):
    t = (X ^ (X << 11)) & 0xffffffff
    X = Y
    Y = Z
    Z = W
    W = W ^ (W >> 19) ^ (t ^ (t >> 8))
    return (X, Y, Z, W)

def backward(X, Y, Z, W):
    t = reverse_xor_rshift(W ^ Z ^ (Z >> 19), 8)
    return (reverse_xor_lshift(t, 11), X, Y, Z)

backward is the function that reverses the state transition. Some random testing:

import random
for _ in range(1000):
    X, Y, Z, W = [random.randint(0,2**32-1) for _ in range(4)]
    assert backward(*forward(X,Y,Z,W)) == (X, Y, Z, W)

Seems to work.

Btw. Hope you'll find this following question http://stackoverflow.com/questions/31521910/simplify-the-inverse-of-z-x-x-y-function/31522122#31522122 interesting — Lu4, Jul 20 '15 at 19:31
@Lu4 actually that's exaxtly the same algorithm, with a much nicer formulation. Thanks :) — Niklas B., Jul 20 '15 at 20:39
Oh to be honest, it's asymptotically more efficient, because it's O(log w) rather than O(w), so that's great — Niklas B., Jul 20 '15 at 23:01
@Lu4 I was referring to harold's answer to that other question — Niklas B., Jul 21 '15 at 10:48

score 3 · Answer 2 · edited May 23 '17 at 11:54

3

For Y, Z and W, we can easily reverse it. For X, we need to make some observations:

W' = W ^ (W >> 19) ^ (t ^ (t >> 8)), -> t ^ (t >> 8) = W' ^ (W ^ (W >> 19))

So, now, we have t ^ (t >> 8) = W' ^ (W ^ (W >> 19)) = a

t = X ^ (X << 11) 

-> t ^ (t >> 8) = X ^ (X << 11) ^ ((X ^ (X <<11)) >> 8) 
                = X ^ (X << 11) ^ (X >> 8) ^ (X << 3)

Denoting each bit of X as x0, x1, x2, ... x31, and each bit of a as a0, a1, ... we can form following equation system:

x0 ^ x8 = a0
x1 ^ x9 = a1
.....

Or, equivalent to:

(x0 + x8) % 2 = a0
(x1 + x9) % 2 = a1
....

Which we can easily solve by applying Gaussian elimination.

edited May 23 '17 at 11:54

Community

1
1

answered Jul 20 '15 at 10:07

Pham Trung

11,204
2
24
43

@Lu4 Look at the equation `X ^ (X << 11) ^ (X >> 8) ^ (X << 3)`, for the first bit of a, we have `a0 = (First bit of X) ^ (First bit of X << 11) ^ (First bit of X >> 8) ^ (First bit of X << 3) = x0 ^ 0 ^ x8 ^ 0 = x0 ^ x8`. – Pham Trung Jul 20 '15 at 11:02
ok, now I get what you mend, I was initially thinking in big-endian format (first bit is the most significant bit) :))) – Lu4 Jul 20 '15 at 11:26
let me write-out all the equations, would it be a problem if I'll use your answer for that purpose? – Lu4 Jul 20 '15 at 11:27
@Lu4 No problem at all :), be careful, there are 32 equations, so it will be long, you can write some script to do that. – Pham Trung Jul 20 '15 at 11:30
And suppose we get the values of `A`, still it will be required for it to get decomposed into `X ^ (X << 11)`... – Lu4 Jul 20 '15 at 11:31
@Lu4 We have smt like this : V[32][32] * X[32] = A[32], with V is fixed, and A is varied depending on value of `A`, so you only need to create matrix V once. – Pham Trung Jul 20 '15 at 11:34

Finding inverse operation to George Marsaglia's XorShift RNG

Abstract

Question

P.S.

Edit

2 Answers2