Applying Linear Algebra to Solve For Fixed Seeds in Pokémon Sword and Shield
...

Preamble
...

2x64-Bit State Implementation of the Xoroshiro PRNG
...

import numpy as np

class Xoroshiro128PlusGF:
	"""Xoroshiro128+ Generator used in GameFreak games"""

	XOROSHIRO_CONST = np.uint64(0x82A2B175229D6A5B)

	def __init__(
		self,
		seed_0: np.uint64,
		seed_1: np.uint64 = XOROSHIRO_CONST
	) -> None:
		self.state = np.empty(2, np.uint64)
		self.state[:] = seed_0, seed_1

	@staticmethod
	def rotl(integer: np.uint64, shift: np.uint64) -> np.uint64:
		"""Rotate the bits of a 64-bit unsigned integer left by
		a specified shift"""
		return (
			(np.uint64(integer) << np.uint64(shift))
			| (np.uint64(integer) >> np.uint64(64 - shift))
		)

	@staticmethod
	def bit_mask(integer: np.uint32) -> np.uint32:
		"""Generate a bitmask for a 32-bit integer"""
		return np.uint32((1 << (int(integer).bit_length())) - 1)

	def next_full(self) -> np.uint64:
		"""Generate the next 64-bit random number"""
		seed_0, seed_1 = self.state
		result = seed_0 + seed_1
		seed_1 ^= seed_0
		self.state[:] = (
			self.rotl(seed_0, 24)
			^ seed_1
			^ (seed_1 << np.uint64(16)),
			self.rotl(seed_1, 37)
		)

		return result

	def next(self) -> np.uint32:
		"""Generate the next 32-bit random number"""
		return np.uint32(self.next_full())

	def rand_max(self, maximum: np.uint32) -> np.uint32:
		"""Generate the next random number < maximum"""
		maximum = np.uint32(maximum)
		mask = self.bit_mask(maximum - 1)
		result = self.next() & mask

		while result >= maximum:
			result = self.next() & mask

		return result

-Matrix Implementation
...

Notation
...

Bit Vector
...

The bit vector representation of an -bit binary integer is defined as follows

Where each entry of the vector is its respective bit in the binary integer.

Shifted Matrices
...

The "-th shifted matrix" of A will be denoted as and is defined by shifting the elements of a matrix column-wise to the right times, discarding any elements that would be shifted past the bounds of the matrix, and shifting in zeros from the left, i.e. for any matrix , is defined as follows:

e.g.

Multiplying a bit vector by a shifted identity matrix is equivalent to the bit-wise left shift operation on the binary number where it is restricted to -bit storage
e.g. for an 5-bit integer

Rotated Matrices
...

The "-th rotated matrix" of A will be denoted as and very similar to the shifted matrix described above. It is defined by shifting the elements of a matrix column-wise to the right times, rotating in any elements that would be shifted past the bounds of the matrix to the left side of the matrix, i.e. for any matrix , is defined as follows:

e.g.

Multiplying this vector by a rotated identity matrix is equivalent to the bit-wise rotate-left operation on the binary number where it is restricted to -bit storage
e.g. for an 5-bit integer

Xoroshiro128+ Matrix Representation
...

The Xoroshiro128+ state transition function is -linear, this means it can be represented as multiplying the state's bit vector representation by some matrix . The bit-wise algorithm for the state transition function can be boiled down to the following operations:

# 1.
seed_1 ^= seed_0
# 2.
seed_0 = rotl(seed_0, 24)
# 3.
seed_0 ^= seed_1
# 4.
seed_0 ^= seed_1 << 16
# 5.
seed_1 = rotl(seed_1, 37)

where seed_0 and seed_1 are the two 64-bit halves of the internal state. Due to the relations between the two halves of the state, it is useful to refer to them each as 64-bit bit vectors . Due to the properties of addition and subtraction in , addition between bit vectors is equivalent to the exclusive-or operation on binary integers. Because of this and the properties of the shifted and rotated matrices defined earlier, we can represent the state transition function as follows:

By looking at each half of the state individually, we can build two matrices which when we multiply the entire state by result in and respectively. The matrices can be built by applying each step as follows:

  1. seed_0 = ...; seed_1 = ...
  2. seed_1 ^= seed_0
  3. seed_0 = rotl(seed_0, 24)
  4. seed_0 ^= seed_1
  5. seed_0 ^= seed_1 << 16
  6. seed_1 = rotl(seed_1, 37)
    Since maps the current state to the next seed_0, and maps the current state to the next seed_1, horizontally concatenating the matrices gives our final target matrix
    that describes the full state transition function.

Initial State Reconstruction From Partial States
...

Consider a modified form of the Xoroshiro128+ output function that when returning a random number in the interval would simply return the value of the first bit of the state. That function would look like this:

def rand2(self):
	result = self.seed0 & 1
	self.next_full()
	return result

This is evidently not how rand_max(2) would be handled in the actual implementation, but it is useful to consider for solving the problem of calculating the initial state of the PRNG from only its outputs. If we look at the first column for any given power of the state transition matrix , we can get a column vector that maps the initial state vector to the value of its first entry after multiplications of the state transition matrix, i.e.

What this means then, is if we construct an augmented square matrix made up of these column vectors for consecutive values of , i.e.

This matrix would map the initial state vector to some vector made up of the first bit of consecutive states

This vector , then is trivially able to be constructed by recording the outputs of the rand2 function described earlier:

prng = ... # prng with an unknown initial state
a = [prng.rand2() for _ in range(128)]

As we have now found a linear mapping from the initial state to this trivially observable vector , if this mapping can be shown to be invertible, the matrix inverse could be used to directly map from back to the initial state , and thereby expose the full initial state based only on the observable outputs. In the case of this implementation of Xoroshiro128+, this matrix can be shown to be invertible.

The Problem
...

In the Nintendo Switch console games Pokémon Sword and Shield, "overworld pokémon" have randomly generated attributes determined by the specific Xoroshiro128PlusGF implementation detailed earlier.

The specific attributes we will be focusing on are determined by the first few outputs of this PRNG when seeded by a 32-bit (for all intents and purposes) randomly generated number. In the most common case, the random number generator is tasked with generating upwards of 94 bits of information from this "fixed seed."

Due to the input size being significantly smaller than the potential output size, the possible results are a small subsection of the theoretical outcomes. This means that when arbitrarily choosing values for these random attributes, you are incredibly likely to pick a combination outside of the actually possible results.

The 94 bits of note are split into three categories:

  1. 32 bits for the "Personality ID" or "PID" for short
  2. 32 bits for the "Encryption Constant" or "EC" for short
  3. 5 bits for each of the 6 "Individual Values" or "IVs" for short

#3 is of particular note as the particular IVs a pokémon have determine its in-game stats, which then determine how well it can perform in battle. EC and PID are in most circumstances entirely hidden from the player and not of much use.

The goal then is to choose values for the IVs that are desirable, e.g. "perfect," having a 31 for each stat, and from there calculate which of the possible seeds can actually generate this result. Once we have a seed that can generate desirable IVs, it can be run through the same algorithm the game uses to generate the EC and PID to ensure all of the characteristics of the pokémon are valid and could actually be generated by the game.

Fixed Attribute Generation Algorithm
...

The game generates the fixed attributes via the following algorithm:

fixed_seed = ... # 32-bit value randomly decided
guaranteed_ivs = ... # a number [0,6] that determines how many IVs are guaranteed to be "perfect"
rng = Xoroshiro128PlusGF(fixed_seed)

ec = rng.rand_max(0xFFFFFFFF)
pid = rng.rand_max(0xFFFFFFFF)

# default values
ivs = [None, None, None, None, None, None]

for _ in range(guaranteed_ivs):
	# choose a random iv to set to "perfect" (31)
	random_iv_index = rng.rand_max(6)
	# if the iv has already been set to 31, choose a new iv until one is found that is not already perfect
	while ivs[random_iv_index] is not None:
		random_iv index = rng.rand_max(6)
	ivs[random_iv_index] = 31

# check each iv
for iv_index in range(6):
	# if the iv has not yet been set, set it to a random value [0,31]
	if ivs[iv_index] is None:
		ivs[iv_index] = rng.rand_max(32)

The thing that really complicates this algorithm and the progression of the PRNG is the guaranteed_ivs variable. This variable is only set to a non-0 value in the following specific cases:

  1. The pokémon is a "brilliant aura pokémon"
    • up to a 3% chance of occurring
    • sets the value to either 2 or 3 with a 50% chance of each
  2. The pokémon is a specific one-time "legendary pokémon encounter"
    • only happens with a specific set of pokémon
    • sets the value to 3

guaranteed_ivs being set can massively change the outcome of the generation algorithm, but due to it only happening in these two cases, for simplicity's sake we will only consider the case when it is set to 0. With guaranteed_ivs set to 0, the algorithm condenses to the following:

fixed_seed = ... # 32-bit value randomly decided
rng = Xoroshiro128PlusGF(fixed_seed)

ec = rng.rand_max(0xFFFFFFFF)
pid = rng.rand_max(0xFFFFFFFF)

ivs = [rng.rand_max(32) for _ in range(6)]

Obtaining Partial States from Individual Values
...

The input value for the rand_max that returns the EC and PID stands out as it is 1 less than , meaning the possible values for each includes every possible 32-bit value except . This is done by the game because it uses 0xFFFFFFFF which is -1 when interpreted as a signed 32-bit integer to denote when a value has not yet been set. If we replace rand_max with its definition we can more clearly see what this means for the relation between the seed and these results

fixed_seed = ... # 32-bit value randomly decided
rng = Xoroshiro128PlusGF(fixed_seed)

mask_0xFFFFFFFF = rng.bit_mask(0xFFFFFFFF - 1)
mask_32 = rng.bit_mask(32 - 1)

ec = rng.next() & mask_0xFFFFFFFF
while ec >= 0xFFFFFFFF:
	ec = rng.next() & mask_0xFFFFFFFF

pid = rng.next() & mask_0xFFFFFFFF
while pid >= 0xFFFFFFFF:
	pid = rng.next() & mask_0xFFFFFFFF

# default values
ivs = [None, None, None, None, None, None]
for iv_index in range(6):
	iv = rng.next() & mask_32
	while iv >= 32:
		iv = rng.next() & mask_32
	ivs[iv_index] = iv

mask_0xFFFFFFFF and mask_32 are constants which can be computed ahead of time

mask_0xFFFFFFFF = Xoroshiro128PlusGF.bit_mask(0xFFFFFFFF - 1)
mask_32 = Xoroshiro128PlusGF.bit_mask(32 - 1)
print(f"{hex(mask_0xFFFFFFFF)=} {mask_32=}")
fixed_seed = ... # 32-bit value randomly decided
rng = Xoroshiro128PlusGF(fixed_seed)

ec = rng.next() & 0xFFFFFFFF
while ec >= 0xFFFFFFFF:
	ec = rng.next() & 0xFFFFFFFF

pid = rng.next() & 0xFFFFFFFF
while pid >= 0xFFFFFFFF:
	pid = rng.next() & 0xFFFFFFFF

# default values
ivs = [None, None, None, None, None, None]
for iv_index in range(6):
	iv = rng.next() & 31
	while iv >= 32:
		iv = rng.next() & 31
	ivs[iv_index] = iv

You'll then notice that the loop that checks if iv >= 32 is redundant, as a bit-wise AND operation by 31 will never return a value greater than 31.

fixed_seed = ... # 32-bit value randomly decided
rng = Xoroshiro128PlusGF(fixed_seed)

ec = rng.next() & 0xFFFFFFFF
while ec >= 0xFFFFFFFF:
	ec = rng.next() & 0xFFFFFFFF

pid = rng.next() & 0xFFFFFFFF
while pid >= 0xFFFFFFFF:
	pid = rng.next() & 0xFFFFFFFF

ivs = [rng.next() & 31 for _ in range(6)]

The loops for EC and PID, however, are not full redundant as a bit-wise AND operation by 0xFFFFFFFF can produce a number 0xFFFFFFFF solely in the case where the value is already 0xFFFFFFFF. This means that if rng.next() returns 0xFFFFFFFF for either the EC or the PID it will cause an extra call to be needed to avoid it. This case, however, is 1 out of possible values that rng.next() can return, making it a very unlikely edge-case which can be accounted for easily down the line. For the sake of simplicity, we will be treating it as redundant.

fixed_seed = ... # 32-bit value randomly decided
rng = Xoroshiro128PlusGF(fixed_seed)

ec = rng.next() & 0xFFFFFFFF
pid = rng.next() & 0xFFFFFFFF
ivs = [rng.next() & 31 for _ in range(6)]

What we can now see is that under these conditions the IVs are generated using bits directly from the result of rng.next() consistently 2-8 "advances" away from the original fixed_seed. Let's now discard EC and PID completely and expand the rng.next() call used for generating the IVs

fixed_seed = ... # 32-bit value randomly decided
rng = Xoroshiro128PlusGF(fixed_seed)

rng.next_full()
rng.next_full()

# default values
ivs = [None, None, None, None, None, None]
for iv_index in range(6):
	seed_0, seed_1 = rng.state
	iv = (seed_0 + seed_1) & 31
	seed_1 ^= seed_0
	rng.state[:] = (
		rng.rotl(seed_0, 24)
		^ seed_1
		^ (seed_1 << np.uint64(16)),
		rng.rotl(seed_1, 37)
	)
	ivs[iv_index] = iv

As 31 in binary is 0b11111, when applying the bit-wise AND operation between it and another number, you simply mask the 5 least significant bits of that number. This means that each IV is simply the 5 least significant bits of the addition between the two halves of the internal state. Additionally, the operations done on the state afterwards are just the state transition algorithm so they will be abbreviated as rng.next_state()

fixed_seed = ... # 32-bit value randomly decided
rng = Xoroshiro128PlusGF(fixed_seed)

rng.next_full()
rng.next_full()

# default values
ivs = [None, None, None, None, None, None]
for iv_index in range(6):
	seed_0, seed_1 = rng.state
	iv = ((seed_0 & 31) + (seed_1 & 31)) & 31
	rng.next_state()
	ivs[iv_index] = iv

Thinking back to the example of initial state reconstruction, we now almost have an observable sequence of outputs that we can build a matrix to map the fixed_seed to. The one glaring issue here is that instead of directly having part of the state as output like we did in the example, we are adding two different portions of the state together as output. If we could somehow determine seed_0 & 31 and seed_1 & 31 from the IV, we would have our target observable sequence.

Unfortunately, addition with carry is a lossy non-linear operation and thus we cannot directly map from the initial state to the IVs. The solution to this problem is to introduce the element of brute-force. We know the value of iv, and if we could know the value of seed_1 & 31, we could simply subtract it from the value of iv to determine seed_0 & 31. We would then have our full observable sequence.

The idea then is to brute-force all 32 possible values of seed_1 & 31 for each IV, thereby giving us the value of both seed_0 & 31 and seed_1 & 31 for each IV. If we were to brute-force this for all 6 IVs, that would be a brute-force of or total values. This already is only 2 times more efficient than if we simply brute-forced all possible seeds, but the key here is that we are getting in return bits of information over the 6 states we are measuring. is way more information than we theoretically need, if everything worked out perfectly, we could build a 32-by-32 matrix to map the fixed_seed directly to 32 bits of observed information, and then invert that to map only 32 bits of information back to the fixed_seed.

If we then only focus on 3 of the IVs rather than all 6, we only need to brute-force possible values for theoretically up to 30 bits of the fixed_seed. We are then left with the task of finding the last bits, which we will later find is fairly easy.

Building an Invertible Matrix
...

The 5 least significant bits of seed_0 are represented by the 1st to 5th entries in the state bit vector, and the 5 least significant bits of seed_1 are represented by the 65th to 69th entries. We can map from the initial state to these bits with the following matrix:

where is an augmented vector containing the observable bits of interest after multiplications by the state transition matrix. We can then concatenate these matrices again to map from an initial state to the bits of interest from all 3 measured IVs

We start from here because of the two state transitions from generating the EC and then the PID prior to the generation of the first 3 IVs. now is 128-by-30 matrix mapping the full initial state to the 30 bits that make up the first 3 IVs. Looking at the initial state, the first 32 bits are made up of the fixed_seed, the next 32 bits are 0, and the last 64 bits are made up of the "Xoroshiro constant" 0x82A2B175229D6A5B. The Xoroshiro constant introduces a constant vector that is added to the result and the zero bits have no effect on the result so we can truncate from the 33rd row onward and handle the constant factor later.

As is a 32-by-30 matrix, it is not square and therefore is not directly invertible. In this case, we can use the Moore-Penrose inverse to compute a particular solution. Once we have a particular solution, we can combine it with any linear combination of the vectors that form the basis of the cokernel of to find all possible solutions. This works because the cokernel of consists of all vectors such that , this means that for any particular solution where , , and therefore also constitutes a solution.

If we find that the cokernel of has a nullity of , that is, the number of linearly independent vectors in the cokernel, we can determine that the entire cokernel is made up of vectors, as under there are only two possible scalars each vector could be multiplied by. This means there are unique solutions to the equation .

Solving the Problem
...

We have now detailed the process for finding all solutions to the equation , where is the bit vector representation of a "fixed seed", and is a augmented vector of "observable" bits based on the IVs of the pokémon we want to find the possible fixed seeds for:

  1. Build the matrix based on concatenating the columns of powers of the state transition matrix which correspond to the observable bits we are measuring, and truncating it to only include the rows that correspond to the bits of the fixed seed
  2. Calculate the Moore-Penrose inverse of , and multiply it our observable bits by it to get a particular solution
  3. Compute the cokernel of and combine each vector with our particular solution to find all unique solutions where is the nullity of the cokernel

What we will then find is that for the problem we have detailed, the nullity of is equal to 4, just greater than the minimum of 2 bits predicted earlier. The linearly independent vectors themselves are the bit vectors representing the integers 0x97a11084, 0x7addf188, 0x852536b0, and 0x2e4b6840 respectively.

We should now build a target vector based on the IVs we want to test our solution. What we need to recall is that half of the target vector is entirely free and each possible value needs to be brute-forced to get all possible solutions. As there are 15 free bits, there are more possible solutions, and therefore fixed seeds that will potentially produce the first three IVs we are using to calculate with.

These seeds can then be manually filtered for the ones that also have the target value for the last 3 IVs. For our example vector, we will use a trivial case where the first 10 free bits are unset, the last 5 bits are set, and our target first 3 IVs will be "perfect", each being 31. This example is chosen because it is known to produce a valid particular solution when multiplying it by , something easy to test when brute-forcing.

Working Through the Solution
...

We can build the vector as follows:

# only the first 3 ivs
target_ivs = [31, 31, 31]
# each are 5 bit integers
free_variables = [0, 0, 31]

# observable bits
l = []
for iv, seed_1_value in zip(target_ivs, free_variables):
	seed_0_value = (iv - seed_1_value) & 31
	seed_0_value_bit_vector = []
	seed_1_value_bit_vector = []
	for bit_index in range(5):
		seed_0_value_bit_vector.append(
			(seed_0_value >> bit_index) & 1
		)
		seed_1_value_bit_vector.append(
			(seed_1_value >> bit_index) & 1
		)
	# bits 1,2,3,4,5
	l.extend(seed_0_value_bit_vector)
	# bits 65,66,67,68,69
	l.extend(seed_1_value_bit_vector)

This gives the vector l = [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1]. To check that this vector produces a valid solution, we need to verify that , which in this case it does. Next we add the solution to every vector in the cokernel to get every unique solution:

v = l @ L_plus
solutions = []
for vector in co_kernel:
	solutions.append(v + vector)

After this we find 16 solutions, the bit vectors for the following integers:
159636180, 2653080144, 1935551324, 3841932248, 2359749732, 453508320, 4135263724, 1641680232, 667464340, 2959712784, 1561674524, 3400815512, 2733474852, 894211232, 3627054508, 1334928680
These bit vectors are all unique solutions to the equation , and as fixed seeds it seems as though they would generate pokémon with our target first three ivs, but we have not yet accounted for the Xoroshiro constant.

Accounting for The Xoroshiro Constant
...

To account for the Xoroshiro constant we need to take another look at the last 64 rows we truncated from our matrix that we will call . These rows map the bits of the Xoroshiro constant to its effect on the result vector . Due to the properties of matrix multiplication, we can show that where is the bit vector representation of the Xoroshiro constant. We had previously been assuming and thus using it interchangeably, to fix this, we simply need to replace with . is a constant scalar multiplied by a constant matrix, and we can calculate it to be the bit vector representing the integer 0x2c10b46c, and this constant is what will refer to moving forward.

Our solutions then become:

v = (l - c) @ L_plus
solutions = []
for vector in co_kernel:
	solutions.append(v + vector)

and we instead need to verify that , which does not end up being true for the originally chosen example, so lets choose a new example that works:

# only the first 3 ivs
target_ivs = [31, 31, 31]
# each are 5 bit integers
free_variables = [0, 0, 8]

# observable bits
l = []
for iv, seed_1_value in zip(target_ivs, free_variables):
	seed_0_value = (iv - seed_1_value) & 31
	seed_0_value_bit_vector = []
	seed_1_value_bit_vector = []
	for bit_index in range(5):
		seed_0_value_bit_vector.append(
			(seed_0_value >> bit_index) & 1
		)
		seed_1_value_bit_vector.append(
			(seed_1_value >> bit_index) & 1
		)
	# bits 1,2,3,4,5
	l.extend(seed_0_value_bit_vector)
	# bits 65,66,67,68,69
	l.extend(seed_1_value_bit_vector)

l = [1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0]

We then find = [1 1 1 1 0 1 1 0 1 1 1 1 0 1 1 0 1 0 1 0 0 0 0 0 0 1 0 1 1 0 0 0] as our particular solution. After adding this to all of the cokernel vectors, we get our 16 unique solutions as the bit vector representations of the following integers:
436563823, 2376368107, 1624809191, 4151938659, 2669697503, 142690651, 3858606167, 1918679251, 877528879, 2750355371, 1318319783, 3643991587, 2976592287, 650780955, 3417751575, 1545064595

Running these fixed seeds through the generation algorithm, we can see that we have in fact solved for 16 unique solutions that give the correct value for our first 3 IVs:

seeds = [436563823, 2376368107, 1624809191, 4151938659, 2669697503, 142690651, 3858606167, 1918679251, 877528879, 2750355371, 1318319783, 3643991587, 2976592287, 650780955, 3417751575, 1545064595]
for seed in seeds:
	rng = Xoroshiro128PlusGF(seed)
	rng.next_full()
	rng.next_full()
	ivs = [rng.next() & 31 for _ in range(6)]
	assert ivs[:3] == [31, 31, 31]
	print(seed, ivs)

Solving For All Solutions
...

The 16 solutions we have just found are only for the case where the free variables are set to 0, 0, and 8. In order to find every solution, we must repeat the same process on every possible value for a free variable. Doing this is as simple as iterating over all values, checking if they produce a valid particular solution, and computing the 16 unique solutions from that initial one. On modern hardware, checking values is fairly fast and so too is storing the potential valid solutions. Even so, looping over the values is an easy target for multi-processing as there is no reason they cannot be run in parallel, and the potential solutions is whittled down by invalid solutions as well as those that simply do not match the target value for the last 3 IVs.

Iterating over all possible values would look like this:

from itertools import product

target_ivs = [31, 31, 31]
fixed_seeds = []

for seed_1_values in product(range(32), range(32), range(32)):
	observable_bits = []
	for iv_index, seed_1_value in enumerate(seed_1_values):
		seed_0_value = (target_ivs[iv_index] - seed_1_value) & 31
		# bits 1,2,3,4,5
		observable_bits.extend(
			(seed_0_value >> bit_index) & 1
			for bit_index in range(5)
		)
		# bits 65,66,67,68,69
		observable_bits.extend(
			(seed_1_value >> bit_index) & 1
			for bit_index in range(5)
		)

	# account for the xoroshiro constant
	observable_bit_vector = (np.array(observable_bits, np.uint8) - c) % 2

	# if the vector does not produce a valid solution, skip it
	if any(
		((observable_bit_vector @ L_plus @ L) % 2)
		!= observable_bit_vector
	):
		continue

	principal_solution = (observable_bit_vector @ L_plus) % 2
	for co_kernel_vector in co_kernel:
		solution = (principal_solution + co_kernel_vector) % 2

		# convert solution from bit vector to integer
		solution_integer = 0
		for bit_index, bit in enumerate(solution):
			solution_integer |= int(bit) << bit_index

		fixed_seeds.append(solution_integer)

From this we find fixed seeds that generate pokémon with their first 3 IVs perfect. If we were to filter this down to only those that have all 6 IVs perfect, statistically, we would expect around or 16 solutions. Unfortunately, we find that it is actually impossible for all 6 IVs to generate perfect (in the guaranteed_ivs = 0 case) as none of the valid seeds found can actually produce this. Our search algorithm is not limited to only this case however, we can plug in any target IVs and find seeds that generate it. Let's try a fun case where all IVs generate as 0:

0x889603f7
0xfa153f53
0xf445cbac
0xc6c083fe
0x17033091
0x9f642e63
0x5249c75b
0x8a1e82e8

How about IVs that generate as 1, 2, 3, 4, 5, 6:

0xac8661b
0x1d281e69
0x39a54ef8

Or even all 9s:

0xaff58c97
0x62f3309e
0x9508d6a8
0xf03e8396

Conclusion
...

What we have shown is that using bit vectors to represent binary integers and matrices to represent binary operations on those integers, we can apply linear algebra concepts to solve problems involving the inputs and outputs of -linear pseudo random number generators. This general concept can be applied to any scenario where you are able to figure out portions of the internal state via its outputs and are looking to find the entire internal state at one point in time. This idea is the backbone behind programs I've written/aided in writing:

as well as the amazing work by Japanese researchers that laid the foundation for them:

Though I'm sure that throughout this document there is notation and formatting that may not align with how linear algebra would formally be presented, I hope that it is digestible enough for any curious mind to comprehend.

Full Functional Code
...

import sys
from itertools import product

# pyodide support
if sys.platform == "emscripten":
    import micropip
    await micropip.install("numpy")

import numpy as np

class Xoroshiro128PlusGF:
    """Xoroshiro128+ Generator used in GameFreak games"""

    XOROSHIRO_CONST = np.uint64(0x82A2B175229D6A5B)

    def __init__(
        self,
        seed_0: np.uint64,
        seed_1: np.uint64 = XOROSHIRO_CONST
    ) -> None:
        self.state = np.empty(2, np.uint64)
        self.state[:] = seed_0, seed_1

    @staticmethod
    def rotl(integer: np.uint64, shift: np.uint64) -> np.uint64:
        """Rotate the bits of a 64-bit unsigned integer left by
        a specified shift"""
        return (
            (np.uint64(integer) << np.uint64(shift))
            | (np.uint64(integer) >> np.uint64(64 - shift))
        )

    @staticmethod
    def bit_mask(integer: np.uint32) -> np.uint32:
        """Generate a bitmask for a 32-bit integer"""
        return np.uint32((1 << (int(integer).bit_length())) - 1)

    def next_full(self) -> np.uint64:
        """Generate the next 64-bit random number"""
        seed_0, seed_1 = self.state
        result = seed_0 + seed_1
        seed_1 ^= seed_0
        self.state[:] = (
            self.rotl(seed_0, 24)
            ^ seed_1
            ^ (seed_1 << np.uint64(16)),
            self.rotl(seed_1, 37)
        )

        return result

    def next(self) -> np.uint32:
        """Generate the next 32-bit random number"""
        return np.uint32(self.next_full())

    def rand_max(self, maximum: np.uint32) -> np.uint32:
        """Generate the next random number < maximum"""
        maximum = np.uint32(maximum)
        mask = self.bit_mask(maximum - 1)
        result = self.next() & mask

        while result >= maximum:
            result = self.next() & mask

        return result


def generate_ivs(fixed_seed: int) -> list[int]:
    """Generate ivs from a fixed seed"""

    rng = Xoroshiro128PlusGF(fixed_seed)
    rng.next_full()
    rng.next_full()

    return [rng.rand_max(32) for _ in range(6)]

def int_to_bit_vector(integer: int, bit_length: int) -> np.ndarray:
    """Convert an integer to its bit vector representation"""

    bit_vec = np.empty(bit_length, np.uint8)
    bit_vec[:] = list(((int(integer) >> i) & 1 for i in range(bit_length)))

    return bit_vec

def bit_vector_to_int(bit_vector: np.ndarray) -> int:
    """Convert a bit vector to its integer representation"""

    integer = 0
    for i, bit in enumerate(bit_vector):
        integer |= int(bit) << i
    return integer

def shift_matrix(matrix: np.ndarray, k: int) -> np.ndarray:
    """Build the k-th shifted matrix of matrix"""
    shifted_matrix = np.roll(matrix, -k, axis = 0)
    shifted_matrix[:, :k] = 0
    return shifted_matrix

def rotate_matrix(matrix: np.ndarray, k: int) -> np.ndarray:
    """Build the k-th rotated matrix of matrix"""
    return np.roll(matrix, -k, axis = 0)

def xoroshiro128plus_mat() -> np.ndarray:
    """Build the 128x128 Xoroshiro128+ state transition matrix"""

    s0_mat = np.zeros((128, 64), np.uint8)
    s1_mat = np.zeros((128, 64), np.uint8)

    s0_mat[0:64] = np.identity(64, np.uint8)
    s1_mat[64:128] = np.identity(64, np.uint8)

    s1_mat ^= s0_mat

    s0_mat = (s0_mat @ rotate_matrix(np.identity(64, np.uint8), 24)) % 2
    s0_mat ^= s1_mat
    s0_mat ^= (s1_mat @ shift_matrix(np.identity(64, np.uint8), 16)) % 2

    s1_mat = (s1_mat @ rotate_matrix(np.identity(64, np.uint8), 37)) % 2

    return np.hstack((s0_mat, s1_mat))

def build_map_mats() -> tuple[np.ndarray, np.ndarray]:
    """Build the 32x30 and 64x30 matrices that map
    from fixed_seed -> ivs and xoroshiro constant -> ivs
    respectively"""

    X = xoroshiro128plus_mat()

    L = np.empty((32, 0), np.uint8)
    C = np.empty((64, 0), np.uint8)

    for k in range(2, 5):
        X_k = np.identity(128, np.uint8)
        for _ in range(k):
            X_k = (X_k @ X) % 2
        for col in (1, 2, 3, 4, 5, 65, 66, 67, 68, 69):
            L = np.hstack((L, X_k[:32, col - 1].reshape(-1, 1)))
            C = np.hstack((C, X_k[64:, col - 1].reshape(-1, 1)))

    return L, C

def resize(matrix: np.ndarray, new_shape: tuple) -> np.ndarray:
    """Resize a matrix, truncating and filling with zeros"""

    mat_rows, mat_cols = matrix.shape
    new_rows, new_cols = new_shape
    new_mat = np.zeros(new_shape, np.uint8)
    new_mat[: min(mat_rows, new_rows), : min(mat_cols, new_cols)] = matrix[
        : min(mat_rows, new_rows), : min(mat_cols, new_cols)
    ]
    return new_mat

def reduced_row_echelon_form(
    matrix: np.ndarray
) -> tuple[np.ndarray, np.ndarray, int, list[int]]:
    """Convert a matrix to reduced row echelon form"""

    rows, columns = matrix.shape
    reduced_form = np.copy(matrix)
    inverse_form = np.identity(rows, np.uint8)
    rank = 0
    pivots = []

    for j in range(columns):
        for i in range(rank, rows):
            if reduced_form[i, j]:
                for k in range(rows):
                    if (k != i) and reduced_form[k, j]:
                        reduced_form[k] ^= reduced_form[i]
                        inverse_form[k] ^= inverse_form[i]
                reduced_form[[i, rank]] = reduced_form[[rank, i]]
                inverse_form[[i, rank]] = inverse_form[[rank, i]]
                pivots.append(j)
                rank += 1
                break
    return reduced_form, inverse_form, rank, pivots

def generalized_inverse(matrix: np.ndarray) -> np.ndarray:
    """Compute the generalized inverse of a matrix"""

    _, inverse_form, rank, pivots = reduced_row_echelon_form(matrix)
    inverse_form = resize(inverse_form, (matrix.shape[1], matrix.shape[0]))
    for i in range(rank - 1, -1, -1):
        column_index = pivots[i]
        inverse_form[[i, column_index]] = inverse_form[[column_index, i]]
    return inverse_form

def co_kernel_basis(matrix: np.ndarray) -> np.ndarray:
    """Compute the basis for the cokernel of a matrix"""

    matrix_inverse = generalized_inverse(matrix)
    basis = (matrix @ matrix_inverse) % 2
    basis = (basis + np.identity(basis.shape[0], np.uint8)) % 2
    basis, _, rank, _ = reduced_row_echelon_form(basis)
    return basis[:rank]

def co_kernel(matrix: np.ndarray) -> np.ndarray:
    """Compute the cokernel of a matrix"""

    basis = co_kernel_basis(matrix)
    space = np.zeros((2 ** basis.shape[0], basis.shape[1]), np.uint8)
    for k in range(space.shape[0]):
        vector = np.zeros(basis.shape[1], np.uint8)
        for i in range(basis.shape[0]):
            if (k >> i) & 1:
                vector ^= basis[i]
        space[k] = vector
    return space

L, C = build_map_mats()
L_plus = generalized_inverse(L)
c = (
    int_to_bit_vector(Xoroshiro128PlusGF.XOROSHIRO_CONST, 64) @ C
) % 2
L_co_kernel = co_kernel(L)

TARGET_IVS = [0, 0, 0, 0, 0, 0]
fixed_seeds = []

for seed_1_values in product(range(32), range(32), range(32)):
    observable_bits = []
    for iv_index, seed_1_value in enumerate(seed_1_values):
        seed_0_value = (TARGET_IVS[iv_index] - seed_1_value) & 31
        # bits 1,2,3,4,5
        observable_bits.extend(
            (seed_0_value >> bit_index) & 1
            for bit_index in range(5)
        )
        # bits 65,66,67,68,69
        observable_bits.extend(
            (seed_1_value >> bit_index) & 1
            for bit_index in range(5)
        )

    # account for the xoroshiro constant
    observable_bit_vector = (np.array(observable_bits, np.uint8) - c) % 2

    # if the vector does not produce a valid solution, skip it
    if any(
        ((observable_bit_vector @ L_plus @ L) % 2)
        != observable_bit_vector
    ):
        continue

    principal_solution = (observable_bit_vector @ L_plus) % 2
    for co_kernel_vector in L_co_kernel:
        solution = (principal_solution + co_kernel_vector) % 2
        solution = bit_vector_to_int(solution)

        solution_ivs = generate_ivs(solution)
        if solution_ivs == TARGET_IVS:
            fixed_seeds.append(solution)

for fixed_seed in fixed_seeds:
    print(f"Valid Fixed Seed: {fixed_seed:08X}")

Matrix-less Implementation
...

import sys
from itertools import product

# pyodide support
if sys.platform == "emscripten":
    import micropip
    await micropip.install("numpy")

import numpy as np

class Xoroshiro128PlusGF:
    """Xoroshiro128+ Generator used in GameFreak games"""

    XOROSHIRO_CONST = np.uint64(0x82A2B175229D6A5B)

    def __init__(
        self,
        seed_0: np.uint64,
        seed_1: np.uint64 = XOROSHIRO_CONST
    ) -> None:
        self.state = np.empty(2, np.uint64)
        self.state[:] = seed_0, seed_1

    @staticmethod
    def rotl(integer: np.uint64, shift: np.uint64) -> np.uint64:
        """Rotate the bits of a 64-bit unsigned integer left by
        a specified shift"""
        return (
            (np.uint64(integer) << np.uint64(shift))
            | (np.uint64(integer) >> np.uint64(64 - shift))
        )

    @staticmethod
    def bit_mask(integer: np.uint32) -> np.uint32:
        """Generate a bitmask for a 32-bit integer"""
        return np.uint32((1 << (int(integer).bit_length())) - 1)

    def next_full(self) -> np.uint64:
        """Generate the next 64-bit random number"""
        seed_0, seed_1 = self.state
        result = seed_0 + seed_1
        seed_1 ^= seed_0
        self.state[:] = (
            self.rotl(seed_0, 24)
            ^ seed_1
            ^ (seed_1 << np.uint64(16)),
            self.rotl(seed_1, 37)
        )

        return result

    def next(self) -> np.uint32:
        """Generate the next 32-bit random number"""
        return np.uint32(self.next_full())

    def rand_max(self, maximum: np.uint32) -> np.uint32:
        """Generate the next random number < maximum"""
        maximum = np.uint32(maximum)
        mask = self.bit_mask(maximum - 1)
        result = self.next() & mask

        while result >= maximum:
            result = self.next() & mask

        return result


def generate_ivs(fixed_seed: int) -> list[int]:
    """Generate ivs from a fixed seed"""

    rng = Xoroshiro128PlusGF(fixed_seed)
    rng.next_full()
    rng.next_full()

    return [rng.rand_max(32) for _ in range(6)]

def matmul(vector: int, matrix: list[int]) -> int:
    """Multiply vector by matrix"""
    result = 0
    for i, row_vector in enumerate(matrix):
        if (vector >> i) & 1:
            result ^= row_vector

    return result

# precomputed matrices
L = [
    1049601,
    35653634,
    71307268,
    176169000,
    352338000,
    671088768,
    301990144,
    603980288,
    134218752,
    268437504,
    536875008,
    34612256,
    69224512,
    138416256,
    311435520,
    622871040,
    138413057,
    311461890,
    622923780,
    172105736,
    344211472,
    688390144,
    268435456,
    536870912,
    1025,
    2050,
    4100,
    34644009,
    69288018,
    138543236,
    277119240,
    554238480,
]

L_plus = [
    75040,
    37749248,
    75498496,
    16779529,
    155325210,
    172190824,
    71349924,
    25240576,
    184680776,
    75506852,
    16852256,
    4194816,
    8389632,
    16793665,
    155353482,
    134217993,
    163980590,
    172705856,
    72398516,
    27337760,
    16777217,
    155320586,
    9248,
    776693824,
    310650420,
    188875016,
    67118116,
    16852000,
    4194304,
    8388608,
]

L_co_kernel = [
    0,
    2543915140,
    2061365640,
    3984384268,
    2233808560,
    310650420,
    4294494008,
    1750718396,
    776693824,
    3119151300,
    1419155912,
    3275196748,
    2876137200,
    1020218996,
    3518214008,
    1175633916
]

c = 0x2c10b46c

TARGET_IVS = [0, 0, 0, 0, 0, 0]
fixed_seeds = []

for seed_1_values in product(range(32), range(32), range(32)):
    observable_bits = 0
    for iv_index, seed_1_value in enumerate(seed_1_values):
        seed_0_value = (TARGET_IVS[iv_index] - seed_1_value) & 31
        observable_bits |= (seed_0_value | (seed_1_value << 5)) << (iv_index * 10)

    # account for the xoroshiro constant
    observable_bits ^= c

    # if the vector does not produce a valid solution, skip it
    if matmul(matmul(observable_bits, L_plus), L) != observable_bits:
        continue

    principal_solution = matmul(observable_bits, L_plus)
    for co_kernel_vector in L_co_kernel:
        solution = principal_solution ^ co_kernel_vector

        solution_ivs = generate_ivs(solution)
        if solution_ivs == TARGET_IVS:
            fixed_seeds.append(solution)

for fixed_seed in fixed_seeds:
    print(f"Valid Fixed Seed: {fixed_seed:08X}")