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;
    }
}

Answer 1

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.

Answer 2

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.