Xorshift

Initialization

There exists a function UnityEngine.Random::set_seed which seeds the rng from a 32-bit integer as follows:

uint32_t state[4];
void UnityEngine_Random_set_seed(uint32_t, seed) {
    state[0] = seed;
    state[1] = state[0] * 0x6c078965 + 1;
    state[2] = state[1] * 0x6c078965 + 1;
    state[3] = state[2] * 0x6c078965 + 1;
}

This function uses the LCG known as ARNG to fill out the state. This function is not actually utilized in BDSP and the state is simply filled by a call to nn::os::GenerateRandomBytes(state, 16) which is a per-process TinyMT.

State Transition Function

The state transition function can be implemented as:

void next_state(uint32_t state[4]) {
    uint32_t t = state[0];
    uint32_t s = state[3];
    s ^= s >> 19;
    t ^= t << 11;
    t ^= t >> 8;
    t ^= s;

    state[0] = state[1];
    state[1] = state[2];
    state[2] = state[3];
    state[3] = t;
}

This can be represented as a matrix by looking at each element of the state individually. Take the full state to be \(v =\) s[0] | (s[1] << 32) | (s[2] << 64) | (s[3] << 96). The state transition function \(f(S)\) is defined such that \({v_{0-31}}' = v_{32-63}\), \({v_{32-63}}' = v_{63-91}\), \({v_{63-91}}' = v_{91-128}\), and \({v_{91-128}}' = vT\). \(v_{0-91}\) is easily represented by a right shift of the entire state by 32 bits \(A_0 = R^{32}\). The matrix \(T\) can be constructed operation by operation as follows:

Isolate the initial values for s and t: \(S_0 = R^{96}\), \(T_0 = M\) where \(M\) is a masking matrix for the constant \(c =\) 0xFFFFFFFF.
Perform the xorshifts in order: \(S_1 = S_0 X_{19}\), \(T_1 = T_0 X_{-11} M\), \(T_2 = T_1 X_{8}\)
Perform the final XOR: \(T = T_2 + S_1\)

The full state transition matrix can be constructed by shifting \(T\) into place and adding it to \(A_0\): \(A = TL^{96} + A_0\). In its full form: \(A = (M X_{-11} M X_{8} + R^{96} X_{19})L^{96} + R^{32}\)

Alternatively, the transition matrix can be constructed directly.

Output Functions

The internal output function simply returns the t from next_state which is the same as truncating the state to just state[3].

uint32_t next_uint(uint32_t state[4]) {
    next_state(state);
    return state[3];
}

Modulo

When provided a range, n = max - min + 1, the result is simply next_uint reduced \(\bmod n\). When a range is not provided, default values min = -0x80000000, max = 0x7FFFFFFE are used:

uint32_t next_rand(uint32_t state[4], uint32_t min = -0x80000000, uint32_t max = 0x7FFFFFFE) {
    return (next_uint(state) % (max - min + 1)) + min;
}

The default range produces all outputs \([0, 2^{32} - 1]\) excluding \(2^{31} - 1\) as the required output from next_uint maps to \(2^{31}\) instead (making it twice as likely as all other outputs). Specified ranges work as expected but have slightly uneven distributions (see: LCG's modulo output).

Alternate Modulo

Many of the calls to this function do not provide a range and deal with the range reduction manually, resulting in an effective second output function:

uint32_t next_alt_rand(uint32_t state[4], uint32_t min, uint32_t max) {
    return (next_rand(state) % (max - min + 1)) + min;
}

Random Float

To generate a uniform float \(\in [0, 1]\), next_uint is truncated to the lower \(23\) bits and divided by \(2^{23} - 1\). To convert to a given range \([\text{min}, \text{max}]\), this value is linear interpolated in reverse:

float next_float(uint32_t state[4]) {
    float t = (next_uint(state) & 0x7FFFFF) / 8388607.0f;
    return (1.0f - t) * max + t * min
}

Stepping Backwards

Stepping backwards requires finding some \(f^{-1}(x)\) such that \(f^{-1}(xA) = x\). This is given by the matrix inverse \(A^{-1}\) as \(xAA^{-1} = x\) by definition. This can be found by matrix inversion algorithms or by using the fact that the RNG is cyclical. Since the period of the RNG is \(2^{128} - 1\), \(2^{128} - 1\) applications of the state transition function must be equivalent to the identity, or, \(A^{2^{128} - 1} = I\). This implies \(A^{2^{128} - 2} = A^{-1}\) and the inverse can be found via matrix exponentiation.

An easier to compute method can be found by realizing that the xorshift operation is easy to invert using only bitwise operations. With right xorshift of the form x^(x>>a), the upper a bits of x>>a must be \(0\) and thus the upper a bits of the input are fully preserved. In other words, (x>>(d-a)) == (f(x)>>(d-a)) where d is the bitlength of x. Since these bits are a part of x, they are also a part of x>>a and can be XORed out of the output: f(x)^((f(x) >> (d-a)) << (d-2*a)) or, if d-2*a < 0 then f(x)^((f(x) >> (d-a)) >> (2*a-d)). This operation restores a bits of x and can be repeated until all d are restored. A nearly identical argument can be made for a left xorshift x^(x<<a) which preserves the lower bits instead.

With this in mind, start by noticing that the initial value for s in next_state is simply the previous value of state[3] which is directly the current value of state[2]. This lets t ^= s be reversed by XORing state[3] by state[2] ^ (state[2] >> 19), or:

uint32_t t = state[3];
uint32_t s = state[2];
s ^= s >> 19;
t ^= s;

Then, to undo the right xorshift by \(8\):

t ^= (t >> (32 - 8)) << (32 - 16);
t ^= ((t >> (32 - 8 - 8)) << (32 - 16 - 8)) & 0xFFFF;
t ^= ((t >> (32 - 8 - 8 - 8)) << (32 - 16 - 8 - 8)) & 0xFF;

or

t ^= (t >> 24) << 16;
t ^= ((t >> 16) << 8) & 0xFFFF;
t ^= (t >> 8) & 0xFF;

which can be shown to be equivalent to

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

Finally, to undo the left xorshift by \(11\):

t ^= (t << (32 - 11)) >> (32 - 22);
t ^= ((t << (32 - 11 - 11)) >> (32 - 22 - 11)) & 0xFFC00000;

or

t ^= (t << 21) >> 10;
t ^= (t << 11) & 0xFFC00000;

which can be shown to be equivalent to

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

The full function can be implemented as:

void previous_state(uint32_t state[4]) {
    uint32_t t = state[3];
    uint32_t s = state[2];
    s ^= s >> 19;
    t ^= s;
    t ^= t >> 8;
    t ^= t >> 16;
    t ^= t << 11;
    t ^= t << 22;

    state[3] = state[2];
    state[2] = state[1];
    state[1] = state[0];
    state[0] = t;
}

Jumping

Efficient jumping of an \(\mathbb{F}_2\)-linear PRNG is described here.

Distance Function

There is no known efficient distance function other than those for the generic discrete logarithm over \(\mathbb{F}_2\). Because the period \(2^{128}-1\) is composite, the Silver-Pohlig-Hellman algorithm can be applied to reduce the size of the problem. An example implementation can be found here.

State Recovery

State recovery can be achieved as described here.