diff --git a/SerialPrograms/Source/Pokemon/Pokemon_Gf2Matrix.cpp b/SerialPrograms/Source/Pokemon/Pokemon_Gf2Matrix.cpp
new file mode 100644
index 0000000000..800b212c9b
--- /dev/null
+++ b/SerialPrograms/Source/Pokemon/Pokemon_Gf2Matrix.cpp
@@ -0,0 +1,183 @@
+/*  GF(2) Linear Algebra
+ *
+ *  From: https://github.com/PokemonAutomation/
+ *
+ */
+
+#include <bit>
+#include "Common/Cpp/Exceptions.h"
+#include "Pokemon_Gf2Matrix.h"
+
+namespace PokemonAutomation{
+namespace Pokemon{
+
+
+
+bool Gf2Vec128::get(size_t index) const{
+    return index < 64
+        ? ((high >> (63 - index)) & 1) != 0
+        : ((low >> (127 - index)) & 1) != 0;
+}
+void Gf2Vec128::set(size_t index, bool value){
+    uint64_t& word = index < 64 ? high : low;
+    uint64_t mask = (uint64_t)1 << (index < 64 ? 63 - index : 127 - index);
+    if (value){
+        word |= mask;
+    }else{
+        word &= ~mask;
+    }
+}
+bool Gf2Vec128::dot(const Gf2Vec128& x) const{
+    int bits = std::popcount(high & x.high) + std::popcount(low & x.low);
+    return (bits & 1) != 0;
+}
+
+
+
+Gf2Matrix128 Gf2Matrix128::identity(){
+    Gf2Matrix128 ret;
+    for (size_t c = 0; c < 128; c++){
+        ret.m_rows[c].set(c, true);
+    }
+    return ret;
+}
+
+//  row i of the product is the XOR of every row k in x that has x[k][i] == 1.
+Gf2Matrix128 Gf2Matrix128::operator*(const Gf2Matrix128& x) const{
+    Gf2Matrix128 ret;
+    for (size_t i = 0; i < 128; i++){
+        Gf2Vec128 accumulator;
+        uint64_t bits = m_rows[i].high;
+        while (bits != 0){
+            //  index = 63 - position.
+            accumulator ^= x.m_rows[63 - (size_t)std::countr_zero(bits)];
+            bits &= bits - 1;
+        }
+        bits = m_rows[i].low;
+        while (bits != 0){
+            //  index = 127 - position.
+            accumulator ^= x.m_rows[127 - (size_t)std::countr_zero(bits)];
+            bits &= bits - 1;
+        }
+        ret.m_rows[i] = accumulator;
+    }
+    return ret;
+}
+Gf2Vec128 Gf2Matrix128::operator*(const Gf2Vec128& column) const{
+    Gf2Vec128 ret;
+    for (size_t c = 0; c < 128; c++){
+        ret.set(c, m_rows[c].dot(column));
+    }
+    return ret;
+}
+
+Gf2Matrix128 Gf2Matrix128::pow(uint64_t exponent) const{
+    Gf2Matrix128 ret = identity();
+    Gf2Matrix128 base = *this;
+    while (exponent != 0){
+        if ((exponent & 1) != 0){
+            ret = ret * base;
+        }
+        exponent >>= 1;
+        if (exponent != 0){
+            base = base * base;
+        }
+    }
+    return ret;
+}
+
+Gf2SolveResult gf2_solve_128(
+    const std::vector<Gf2Vec128>& equations,
+    const std::vector<bool>& rhs
+){
+    if (equations.size() != rhs.size()){
+        throw InternalProgramError(
+            nullptr, PA_CURRENT_FUNCTION,
+            "gf2_solve_128(): Coefficient and constant counts do not match."
+        );
+    }
+
+    //  Augmented system. We reduce to row echelon form, then back-substitute.
+    std::vector<Gf2Vec128> rows = equations;
+    std::vector<bool> constants = rhs;
+    const size_t height = rows.size();
+
+    //  pivot_row[c] is the row that owns column "c" as its pivot, or NO_PIVOT.
+    const size_t NO_PIVOT = (size_t)0 - 1;
+    std::array<size_t, 128> pivot_row;
+    pivot_row.fill(NO_PIVOT);
+
+    size_t next_row = 0;
+    for (size_t column = 0; column < 128 && next_row < height; column++){
+        size_t pivot = NO_PIVOT;
+        for (size_t row = next_row; row < height; row++){
+            if (rows[row].get(column)){
+                pivot = row;
+                break;
+            }
+        }
+        if (pivot == NO_PIVOT){
+            continue;
+        }
+
+        std::swap(rows[next_row], rows[pivot]);
+        {
+            //  std::vector<bool> has no swappable references.
+            bool tmp = constants[next_row];
+            constants[next_row] = constants[pivot];
+            constants[pivot] = tmp;
+        }
+
+        for (size_t row = 0; row < height; row++){
+            if (row != next_row && rows[row].get(column)){
+                rows[row] ^= rows[next_row];
+                constants[row] = constants[row] != constants[next_row];
+            }
+        }
+
+        pivot_row[column] = next_row;
+        next_row++;
+    }
+
+    Gf2SolveResult result;
+
+    //  Any all-zero row with a nonzero constant makes the system unsolvable.
+    for (size_t row = 0; row < height; row++){
+        if (rows[row].is_zero() && constants[row]){
+            return result;
+        }
+    }
+    result.consistent = true;
+
+    //  Free variables are set to zero, so each pivot variable is just its constant.
+    for (size_t column = 0; column < 128; column++){
+        if (pivot_row[column] != NO_PIVOT){
+            result.solution.set(column, constants[pivot_row[column]]);
+        }
+    }
+
+    //  One null space basis vector per free column.
+    for (size_t column = 0; column < 128; column++){
+        if (pivot_row[column] != NO_PIVOT){
+            continue;
+        }
+        Gf2Vec128 basis;
+        basis.set(column, true);
+        for (size_t pivot_column = 0; pivot_column < 128; pivot_column++){
+            size_t row = pivot_row[pivot_column];
+            if (row != NO_PIVOT && rows[row].get(column)){
+                basis.set(pivot_column, true);
+            }
+        }
+        result.null_space_basis.emplace_back(basis);
+    }
+    result.null_space_dimension = result.null_space_basis.size();
+
+    return result;
+}
+
+
+
+
+}
+}
diff --git a/SerialPrograms/Source/Pokemon/Pokemon_Gf2Matrix.h b/SerialPrograms/Source/Pokemon/Pokemon_Gf2Matrix.h
new file mode 100644
index 0000000000..10bedbdf8e
--- /dev/null
+++ b/SerialPrograms/Source/Pokemon/Pokemon_Gf2Matrix.h
@@ -0,0 +1,95 @@
+/*  GF(2) Linear Algebra
+ *
+ *  From: https://github.com/PokemonAutomation/
+ *
+ */
+
+#ifndef PokemonAutomation_Pokemon_Gf2Matrix_H
+#define PokemonAutomation_Pokemon_Gf2Matrix_H
+
+#include <stddef.h>
+#include <stdint.h>
+#include <array>
+#include <string>
+#include <vector>
+
+namespace PokemonAutomation{
+namespace Pokemon{
+
+
+//  A 128-bit vector over GF(2).
+//  Index 0 is the MSB of "high". Index 127 is the LSB of "low".
+struct Gf2Vec128{
+    uint64_t high = 0;
+    uint64_t low = 0;
+
+    Gf2Vec128() = default;
+    Gf2Vec128(uint64_t p_high, uint64_t p_low) : high(p_high), low(p_low) {}
+
+    bool operator==(const Gf2Vec128& x) const{ return high == x.high && low == x.low; }
+    bool operator!=(const Gf2Vec128& x) const{ return !(*this == x); }
+
+    bool get(size_t index) const;
+    void set(size_t index, bool value);
+
+    bool is_zero() const{ return (high | low) == 0; }
+
+    //  Parity of the bitwise AND. This is the GF(2) dot product.
+    bool dot(const Gf2Vec128& x) const;
+
+    Gf2Vec128 operator^(const Gf2Vec128& x) const{ return Gf2Vec128(high ^ x.high, low ^ x.low); }
+    Gf2Vec128& operator^=(const Gf2Vec128& x){ high ^= x.high; low ^= x.low; return *this; }
+};
+
+
+//  A 128x128 matrix over GF(2), stored as 128 row vectors.
+//  Vectors are treated as columns, so "matrix * vector" is the usual product.
+class Gf2Matrix128{
+public:
+    static Gf2Matrix128 identity();
+
+    const Gf2Vec128& operator[](size_t row) const{ return m_rows[row]; }
+    Gf2Vec128& operator[](size_t row){ return m_rows[row]; }
+
+    bool operator==(const Gf2Matrix128& x) const{ return m_rows == x.m_rows; }
+    bool operator!=(const Gf2Matrix128& x) const{ return !(*this == x); }
+
+    Gf2Matrix128 operator*(const Gf2Matrix128& x) const;
+    Gf2Vec128 operator*(const Gf2Vec128& column) const;
+
+    Gf2Matrix128 pow(uint64_t exponent) const;
+
+private:
+    std::array<Gf2Vec128, 128> m_rows;
+};
+
+
+struct Gf2SolveResult{
+    //  False means the system has no solution at all. This normally indicates
+    //  corrupt observations rather than a bug in the caller.
+    bool consistent = false;
+
+    //  One particular solution. Only meaningful if "consistent".
+    Gf2Vec128 solution;
+
+    //  Zero means the solution is unique. Anything larger means the
+    //  observations under-determine the state and more are needed.
+    size_t null_space_dimension = 0;
+
+    //  Basis of the null space. Adding any XOR-combination of these to
+    //  "solution" gives another valid solution.
+    std::vector<Gf2Vec128> null_space_basis;
+};
+
+//  Solve "equations * x = rhs" over GF(2), where each entry of "equations" is one
+//  row of coefficients and the matching entry of "rhs" is that row's constant.
+//  The system may have any number of rows; more than 128 is normal and desirable.
+Gf2SolveResult gf2_solve_128(
+    const std::vector<Gf2Vec128>& equations,
+    const std::vector<bool>& rhs
+);
+
+
+}
+}
+#endif
diff --git a/SerialPrograms/Source/Pokemon/Pokemon_Xorshift128.cpp b/SerialPrograms/Source/Pokemon/Pokemon_Xorshift128.cpp
new file mode 100644
index 0000000000..846766d311
--- /dev/null
+++ b/SerialPrograms/Source/Pokemon/Pokemon_Xorshift128.cpp
@@ -0,0 +1,156 @@
+/*  Xorshift128
+ *
+ *  From: https://github.com/PokemonAutomation/
+ *
+ */
+
+#include "Common/Cpp/PrettyPrint.h"
+#include "Pokemon_Xorshift128.h"
+
+namespace PokemonAutomation{
+namespace Pokemon{
+
+
+const uint64_t JUMP_THRESHOLD = 512;
+
+
+std::string Xorshift128State::to_string() const{
+    return "[0x" + tostr_hex_padded(8, s0)
+        + ", 0x" + tostr_hex_padded(8, s1)
+        + ", 0x" + tostr_hex_padded(8, s2)
+        + ", 0x" + tostr_hex_padded(8, s3) + "]";
+}
+
+
+Xorshift128State xorshift128_state_from_seed_pair(uint64_t seed0, uint64_t seed1){
+    return Xorshift128State(
+        (uint32_t)(seed0 >> 32), (uint32_t)seed0,
+        (uint32_t)(seed1 >> 32), (uint32_t)seed1
+    );
+}
+void xorshift128_state_to_seed_pair(const Xorshift128State& state, uint64_t& seed0, uint64_t& seed1){
+    seed0 = ((uint64_t)state.s0 << 32) | state.s1;
+    seed1 = ((uint64_t)state.s2 << 32) | state.s3;
+}
+
+Gf2Vec128 xorshift128_state_to_vector(const Xorshift128State& state){
+    Gf2Vec128 ret;
+    xorshift128_state_to_seed_pair(state, ret.high, ret.low);
+    return ret;
+}
+Xorshift128State xorshift128_state_from_vector(const Gf2Vec128& vector){
+    return xorshift128_state_from_seed_pair(vector.high, vector.low);
+}
+
+
+
+uint32_t Xorshift128::next(){
+    uint32_t temp = m_state.s0 ^ (m_state.s0 << 11);
+    uint32_t old_s3 = m_state.s3;
+
+    m_state.s0 = m_state.s1;
+    m_state.s1 = m_state.s2;
+    m_state.s2 = m_state.s3;
+    m_state.s3 = temp ^ (temp >> 8) ^ old_s3 ^ (old_s3 >> 19);
+
+    return m_state.s3;
+}
+
+void Xorshift128::prev(){
+    uint32_t temp = (m_state.s2 >> 19) ^ m_state.s2 ^ m_state.s3;
+
+    temp ^= temp >> 8;
+    temp ^= temp >> 16;
+
+    temp ^= temp << 11;
+    temp ^= temp << 22;
+
+    m_state.s3 = m_state.s2;
+    m_state.s2 = m_state.s1;
+    m_state.s1 = m_state.s0;
+    m_state.s0 = temp;
+}
+
+void Xorshift128::advance(uint64_t count){
+    if (count < JUMP_THRESHOLD){
+        for (uint64_t c = 0; c < count; c++){
+            next();
+        }
+        return;
+    }
+    m_state = xorshift128_state_from_vector(
+        xorshift128_transition_power(count) * xorshift128_state_to_vector(m_state)
+    );
+}
+void Xorshift128::rewind(uint64_t count){
+    if (count < JUMP_THRESHOLD){
+        for (uint64_t c = 0; c < count; c++){
+            prev();
+        }
+        return;
+    }
+    m_state = xorshift128_state_from_vector(
+        xorshift128_inverse_transition_matrix().pow(count) * xorshift128_state_to_vector(m_state)
+    );
+}
+
+
+
+template <typename StepFunction>
+static Gf2Matrix128 build_step_matrix(StepFunction&& step){
+    Gf2Matrix128 matrix;
+    for (size_t column = 0; column < 128; column++){
+        Gf2Vec128 basis;
+        basis.set(column, true);
+
+        Xorshift128 rng(xorshift128_state_from_vector(basis));
+        step(rng);
+        Gf2Vec128 image = xorshift128_state_to_vector(rng.state());
+
+        for (size_t row = 0; row < 128; row++){
+            if (image.get(row)){
+                matrix[row].set(column, true);
+            }
+        }
+    }
+    return matrix;
+}
+
+const Gf2Matrix128& xorshift128_transition_matrix(){
+    static Gf2Matrix128 matrix = build_step_matrix([](Xorshift128& rng){ rng.next(); });
+    return matrix;
+}
+const Gf2Matrix128& xorshift128_inverse_transition_matrix(){
+    static Gf2Matrix128 matrix = build_step_matrix([](Xorshift128& rng){ rng.prev(); });
+    return matrix;
+}
+
+//  T^(2^k) for every k that fits in a 64-bit count.
+static const std::array<Gf2Matrix128, 64>& xorshift128_transition_powers_of_two(){
+    static std::array<Gf2Matrix128, 64> table = [](){
+        std::array<Gf2Matrix128, 64> ret;
+        ret[0] = xorshift128_transition_matrix();
+        for (size_t c = 1; c < ret.size(); c++){
+            ret[c] = ret[c - 1] * ret[c - 1];
+        }
+        return ret;
+    }();
+    return table;
+}
+
+Gf2Matrix128 xorshift128_transition_power(uint64_t count){
+    const std::array<Gf2Matrix128, 64>& powers = xorshift128_transition_powers_of_two();
+    Gf2Matrix128 ret = Gf2Matrix128::identity();
+    for (size_t bit = 0; count != 0; count >>= 1, bit++){
+        if ((count & 1) != 0){
+            ret = powers[bit] * ret;
+        }
+    }
+    return ret;
+}
+
+
+
+
+}
+}
diff --git a/SerialPrograms/Source/Pokemon/Pokemon_Xorshift128.h b/SerialPrograms/Source/Pokemon/Pokemon_Xorshift128.h
new file mode 100644
index 0000000000..4c5a3ff0a7
--- /dev/null
+++ b/SerialPrograms/Source/Pokemon/Pokemon_Xorshift128.h
@@ -0,0 +1,123 @@
+/*  Xorshift128
+ *
+ *  From: https://github.com/PokemonAutomation/
+ *
+ */
+
+#ifndef PokemonAutomation_Pokemon_Xorshift128_H
+#define PokemonAutomation_Pokemon_Xorshift128_H
+
+#include <stddef.h>
+#include <stdint.h>
+#include <array>
+#include <string>
+#include "Pokemon_Gf2Matrix.h"
+
+namespace PokemonAutomation{
+namespace Pokemon{
+
+
+struct Xorshift128State{
+    uint32_t s0 = 0;
+    uint32_t s1 = 0;
+    uint32_t s2 = 0;
+    uint32_t s3 = 0;
+
+    Xorshift128State() = default;
+    Xorshift128State(uint32_t p_s0, uint32_t p_s1, uint32_t p_s2, uint32_t p_s3)
+        : s0(p_s0), s1(p_s1), s2(p_s2), s3(p_s3)
+    {}
+
+    bool operator==(const Xorshift128State& x) const{
+        return s0 == x.s0 && s1 == x.s1 && s2 == x.s2 && s3 == x.s3;
+    }
+    bool operator!=(const Xorshift128State& x) const{ return !(*this == x); }
+
+    //  "[0x........, 0x........, 0x........, 0x........]"
+    std::string to_string() const;
+};
+
+
+Xorshift128State xorshift128_state_from_seed_pair(uint64_t seed0, uint64_t seed1);
+void xorshift128_state_to_seed_pair(const Xorshift128State& state, uint64_t& seed0, uint64_t& seed1);
+
+
+Gf2Vec128 xorshift128_state_to_vector(const Xorshift128State& state);
+Xorshift128State xorshift128_state_from_vector(const Gf2Vec128& vector);
+
+
+class Xorshift128{
+public:
+    Xorshift128() = default;
+    explicit Xorshift128(const Xorshift128State& state) : m_state(state) {}
+    Xorshift128(uint32_t s0, uint32_t s1, uint32_t s2, uint32_t s3) : m_state(s0, s1, s2, s3) {}
+
+    const Xorshift128State& state() const{ return m_state; }
+    void set_state(const Xorshift128State& state){ m_state = state; }
+
+    uint32_t next();
+
+    //  Step backwards. Undoes exactly one next().
+    void prev();
+
+    void advance(uint64_t count);
+    void rewind(uint64_t count);
+
+private:
+    Xorshift128State m_state;
+};
+
+
+const Gf2Matrix128& xorshift128_transition_matrix();
+const Gf2Matrix128& xorshift128_inverse_transition_matrix();
+
+Gf2Matrix128 xorshift128_transition_power(uint64_t count);
+
+
+inline uint32_t bdsp_gen_transform(uint32_t raw){
+    return (uint32_t)(raw % 0xFFFFFFFF) + 0x80000000;
+}
+
+
+//  A sliding window of pre-generated outputs. Size must be a power of two.
+template <size_t Size>
+class Xorshift128List{
+    static_assert(Size != 0 && (Size & (Size - 1)) == 0, "Size must be a power of two.");
+
+public:
+    explicit Xorshift128List(const Xorshift128& rng)
+        : m_rng(rng)
+    {
+        for (uint32_t& value : m_buffer){
+            value = m_rng.next();
+        }
+    }
+
+    uint32_t next_raw(){ return m_buffer[m_index++ & (Size - 1)]; }
+    uint32_t next_gen(){ return bdsp_gen_transform(next_raw()); }
+
+    uint32_t next_raw_modulo(uint32_t modulo){ return next_raw() % modulo; }
+    uint32_t next_gen_modulo(uint32_t modulo){ return next_gen() % modulo; }
+
+    //  Skip values without reading them.
+    void advance(size_t count){ m_index += count; }
+
+    void advance_state(){
+        m_buffer[m_head++ & (Size - 1)] = m_rng.next();
+        m_index = m_head;
+    }
+
+    //  Rewind the read position to the start of the current window.
+    void reset_index(){ m_index = m_head; }
+
+private:
+    Xorshift128 m_rng;
+    std::array<uint32_t, Size> m_buffer;
+    size_t m_head = 0;
+    size_t m_index = 0;
+};
+
+
+}
+}
+#endif
diff --git a/SerialPrograms/cmake/SourceFiles.cmake b/SerialPrograms/cmake/SourceFiles.cmake
index 5ee82caae2..b7a288cefc 100644
--- a/SerialPrograms/cmake/SourceFiles.cmake
+++ b/SerialPrograms/cmake/SourceFiles.cmake
@@ -1347,6 +1347,8 @@ file(GLOB LIBRARY_SOURCES
     Source/Pokemon/Pokemon_DataTypes.h
     Source/Pokemon/Pokemon_EncounterStats.cpp
     Source/Pokemon/Pokemon_EncounterStats.h
+    Source/Pokemon/Pokemon_Gf2Matrix.cpp
+    Source/Pokemon/Pokemon_Gf2Matrix.h
     Source/Pokemon/Pokemon_IvJudge.cpp
     Source/Pokemon/Pokemon_IvJudge.h
     Source/Pokemon/Pokemon_NatureChecker.cpp
@@ -1365,6 +1367,8 @@ file(GLOB LIBRARY_SOURCES
     Source/Pokemon/Pokemon_Types.h
     Source/Pokemon/Pokemon_Xoroshiro128Plus.cpp
     Source/Pokemon/Pokemon_Xoroshiro128Plus.h
+    Source/Pokemon/Pokemon_Xorshift128.cpp
+    Source/Pokemon/Pokemon_Xorshift128.h
     Source/Pokemon/Pokemon_AdvRng.cpp
     Source/Pokemon/Pokemon_AdvRng.h
     Source/Pokemon/Resources/Pokemon_BerryNames.cpp