diff --git a/CMakeLists.txt b/CMakeLists.txt
index d611add..b46c5e4 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -303,7 +303,12 @@ add_library(aaruformat
src/ngcw/ngcw_junk.h
src/ngcw/wii_crypto.c
src/ngcw/wii_crypto.h
- src/compression/zstd.c)
+ src/compression/zstd.c
+ src/lib/gf256.c
+ src/lib/gf256.h
+ src/lib/reed_solomon.c
+ src/lib/reed_solomon.h
+ include/aaruformat/structs/erasure.h)
# Set up include directories for the target
target_include_directories(aaruformat
diff --git a/src/lib/gf256.c b/src/lib/gf256.c
new file mode 100644
index 0000000..549c579
--- /dev/null
+++ b/src/lib/gf256.c
@@ -0,0 +1,357 @@
+/*
+ * This file is part of the Aaru Data Preservation Suite.
+ * Copyright (c) 2019-2026 Natalia Portillo.
+ *
+ * This library is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation; either version 2.1 of the
+ * License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, see .
+ */
+
+/**
+ * @file gf256.c
+ * @brief GF(2^8) arithmetic with SIMD-accelerated region multiply.
+ *
+ * Galois Field GF(2^8) with irreducible polynomial x^8 + x^4 + x^3 + x^2 + 1
+ * (0x11D). Uses log/antilog tables for scalar operations and 4-bit nibble split
+ * tables with SIMD shuffle for vectorized region multiply-accumulate.
+ */
+
+#include
+#include
+#include
+
+#include "gf256.h"
+#include "aaruformat/simd.h"
+
+/* -------------------------------------------------------------------------
+ * Log / anti-log tables for GF(2^8) with polynomial 0x11D
+ * Generator element: 2
+ * ------------------------------------------------------------------------- */
+
+/** Log table: gf256_log[x] = discrete log base 2 of x in GF(2^8). gf256_log[0] is undefined. */
+static uint8_t gf256_log_table[256];
+/** Anti-log (exp) table: gf256_exp[i] = 2^i mod P. Extended to 512 entries to avoid modular reduction. */
+static uint8_t gf256_exp_table[512];
+
+/** Flag to ensure tables are initialized exactly once. */
+static int gf256_tables_initialized = 0;
+
+/**
+ * @brief Initialize log/antilog tables for GF(2^8) with polynomial 0x11D.
+ */
+static void gf256_init_tables(void)
+{
+ if(gf256_tables_initialized) return;
+
+ unsigned x = 1;
+ for(int i = 0; i < 255; i++)
+ {
+ gf256_exp_table[i] = (uint8_t)x;
+ gf256_exp_table[i+255] = (uint8_t)x; /* wrap-around for easy mod 255 */
+ gf256_log_table[x] = (uint8_t)i;
+
+ /* Multiply by generator 2 in GF(2^8) */
+ x <<= 1;
+ if(x & 0x100) x ^= 0x11D;
+ }
+ gf256_log_table[0] = 0; /* Convention: log(0) = 0, unused in mul since we short-circuit */
+ gf256_exp_table[510] = gf256_exp_table[0]; /* Complete wrap */
+ gf256_exp_table[511] = gf256_exp_table[1];
+
+ gf256_tables_initialized = 1;
+}
+
+/* -------------------------------------------------------------------------
+ * Scalar operations
+ * ------------------------------------------------------------------------- */
+
+uint8_t gf256_mul(uint8_t a, uint8_t b)
+{
+ if(a == 0 || b == 0) return 0;
+ gf256_init_tables();
+ return gf256_exp_table[gf256_log_table[a] + gf256_log_table[b]];
+}
+
+uint8_t gf256_div(uint8_t a, uint8_t b)
+{
+ if(a == 0) return 0;
+ /* b must be non-zero */
+ gf256_init_tables();
+ return gf256_exp_table[gf256_log_table[a] + 255 - gf256_log_table[b]];
+}
+
+uint8_t gf256_inv(uint8_t a)
+{
+ /* a must be non-zero */
+ gf256_init_tables();
+ return gf256_exp_table[255 - gf256_log_table[a]];
+}
+
+/* -------------------------------------------------------------------------
+ * SIMD region multiply-accumulate: dst[i] ^= GF_mul(src[i], coeff)
+ *
+ * Technique: 4-bit nibble decomposition.
+ * For a given coeff, precompute:
+ * low_tbl[i] = GF_mul(i, coeff) for i = 0..15
+ * hi_tbl[i] = GF_mul(i<<4, coeff) for i = 0..15
+ * Then for each byte b:
+ * GF_mul(b, coeff) = low_tbl[b & 0x0F] ^ hi_tbl[b >> 4]
+ * This maps to SIMD shuffle (pshufb / vpshufb / vqtbl1q_u8).
+ * ------------------------------------------------------------------------- */
+
+/**
+ * @brief Build the two 16-byte nibble lookup tables for a given coefficient.
+ */
+static void gf256_build_mul_tables(uint8_t coeff, uint8_t low_tbl[16], uint8_t hi_tbl[16])
+{
+ gf256_init_tables();
+ for(int i = 0; i < 16; i++)
+ {
+ low_tbl[i] = gf256_mul((uint8_t)i, coeff);
+ hi_tbl[i] = gf256_mul((uint8_t)(i << 4), coeff);
+ }
+}
+
+/* ---------- Scalar fallback ---------- */
+
+static void gf256_mul_region_scalar(uint8_t *dst, const uint8_t *src, uint8_t coeff, size_t len)
+{
+ uint8_t low_tbl[16], hi_tbl[16];
+ gf256_build_mul_tables(coeff, low_tbl, hi_tbl);
+
+ for(size_t i = 0; i < len; i++)
+ dst[i] ^= low_tbl[src[i] & 0x0F] ^ hi_tbl[src[i] >> 4];
+}
+
+static void gf256_xor_region_scalar(uint8_t *dst, const uint8_t *src, size_t len)
+{
+ size_t i = 0;
+
+ /* Process 8 bytes at a time */
+ for(; i + 8 <= len; i += 8)
+ {
+ uint64_t d, s;
+ memcpy(&d, dst + i, 8);
+ memcpy(&s, src + i, 8);
+ d ^= s;
+ memcpy(dst + i, &d, 8);
+ }
+
+ for(; i < len; i++)
+ dst[i] ^= src[i];
+}
+
+/* ---------- x86 SSSE3 ---------- */
+
+#if defined(__x86_64__) || defined(__amd64) || defined(_M_AMD64) || defined(_M_X64) || \
+ defined(__I386__) || defined(__i386__) || defined(__THW_INTEL) || defined(_M_IX86)
+
+#include /* SSSE3: _mm_shuffle_epi8 */
+
+SSSE3 static void gf256_mul_region_ssse3(uint8_t *dst, const uint8_t *src, uint8_t coeff, size_t len)
+{
+ uint8_t low_tbl[16], hi_tbl[16];
+ gf256_build_mul_tables(coeff, low_tbl, hi_tbl);
+
+ const __m128i low_v = _mm_loadu_si128((const __m128i *)low_tbl);
+ const __m128i hi_v = _mm_loadu_si128((const __m128i *)hi_tbl);
+ const __m128i mask = _mm_set1_epi8(0x0F);
+
+ size_t i = 0;
+ for(; i + 16 <= len; i += 16)
+ {
+ __m128i s = _mm_loadu_si128((const __m128i *)(src + i));
+ __m128i d = _mm_loadu_si128((const __m128i *)(dst + i));
+ __m128i s_lo = _mm_and_si128(s, mask);
+ __m128i s_hi = _mm_and_si128(_mm_srli_epi16(s, 4), mask);
+ __m128i lo = _mm_shuffle_epi8(low_v, s_lo);
+ __m128i hi = _mm_shuffle_epi8(hi_v, s_hi);
+ __m128i r = _mm_xor_si128(_mm_xor_si128(lo, hi), d);
+ _mm_storeu_si128((__m128i *)(dst + i), r);
+ }
+
+ /* Tail */
+ for(; i < len; i++)
+ dst[i] ^= low_tbl[src[i] & 0x0F] ^ hi_tbl[src[i] >> 4];
+}
+
+SSSE3 static void gf256_xor_region_ssse3(uint8_t *dst, const uint8_t *src, size_t len)
+{
+ size_t i = 0;
+ for(; i + 16 <= len; i += 16)
+ {
+ __m128i d = _mm_loadu_si128((const __m128i *)(dst + i));
+ __m128i s = _mm_loadu_si128((const __m128i *)(src + i));
+ _mm_storeu_si128((__m128i *)(dst + i), _mm_xor_si128(d, s));
+ }
+ for(; i < len; i++) dst[i] ^= src[i];
+}
+
+/* ---------- x86 AVX2 ---------- */
+
+#include /* AVX2: _mm256_shuffle_epi8 */
+
+AVX2 static void gf256_mul_region_avx2(uint8_t *dst, const uint8_t *src, uint8_t coeff, size_t len)
+{
+ uint8_t low_tbl[16], hi_tbl[16];
+ gf256_build_mul_tables(coeff, low_tbl, hi_tbl);
+
+ /* Broadcast 16-byte tables to both 128-bit lanes of 256-bit register */
+ const __m128i low_128 = _mm_loadu_si128((const __m128i *)low_tbl);
+ const __m128i hi_128 = _mm_loadu_si128((const __m128i *)hi_tbl);
+ const __m256i low_v = _mm256_broadcastsi128_si256(low_128);
+ const __m256i hi_v = _mm256_broadcastsi128_si256(hi_128);
+ const __m256i mask = _mm256_set1_epi8(0x0F);
+
+ size_t i = 0;
+ for(; i + 32 <= len; i += 32)
+ {
+ __m256i s = _mm256_loadu_si256((const __m256i *)(src + i));
+ __m256i d = _mm256_loadu_si256((const __m256i *)(dst + i));
+ __m256i s_lo = _mm256_and_si256(s, mask);
+ __m256i s_hi = _mm256_and_si256(_mm256_srli_epi16(s, 4), mask);
+ __m256i lo = _mm256_shuffle_epi8(low_v, s_lo);
+ __m256i hi = _mm256_shuffle_epi8(hi_v, s_hi);
+ __m256i r = _mm256_xor_si256(_mm256_xor_si256(lo, hi), d);
+ _mm256_storeu_si256((__m256i *)(dst + i), r);
+ }
+
+ /* Tail: SSSE3 for remaining 16-byte chunks, then scalar */
+ const __m128i low_v2 = low_128;
+ const __m128i hi_v2 = hi_128;
+ const __m128i mask2 = _mm_set1_epi8(0x0F);
+ for(; i + 16 <= len; i += 16)
+ {
+ __m128i s = _mm_loadu_si128((const __m128i *)(src + i));
+ __m128i d = _mm_loadu_si128((const __m128i *)(dst + i));
+ __m128i s_lo = _mm_and_si128(s, mask2);
+ __m128i s_hi = _mm_and_si128(_mm_srli_epi16(s, 4), mask2);
+ __m128i lo = _mm_shuffle_epi8(low_v2, s_lo);
+ __m128i hi = _mm_shuffle_epi8(hi_v2, s_hi);
+ __m128i r = _mm_xor_si128(_mm_xor_si128(lo, hi), d);
+ _mm_storeu_si128((__m128i *)(dst + i), r);
+ }
+
+ for(; i < len; i++)
+ dst[i] ^= low_tbl[src[i] & 0x0F] ^ hi_tbl[src[i] >> 4];
+}
+
+AVX2 static void gf256_xor_region_avx2(uint8_t *dst, const uint8_t *src, size_t len)
+{
+ size_t i = 0;
+ for(; i + 32 <= len; i += 32)
+ {
+ __m256i d = _mm256_loadu_si256((const __m256i *)(dst + i));
+ __m256i s = _mm256_loadu_si256((const __m256i *)(src + i));
+ _mm256_storeu_si256((__m256i *)(dst + i), _mm256_xor_si256(d, s));
+ }
+ for(; i + 16 <= len; i += 16)
+ {
+ __m128i d = _mm_loadu_si128((const __m128i *)(dst + i));
+ __m128i s = _mm_loadu_si128((const __m128i *)(src + i));
+ _mm_storeu_si128((__m128i *)(dst + i), _mm_xor_si128(d, s));
+ }
+ for(; i < len; i++) dst[i] ^= src[i];
+}
+
+/* Forward declarations for CPUID-based detection (from simd.c) */
+int have_ssse3(void);
+int have_avx2(void);
+
+#endif /* x86 */
+
+/* ---------- ARM NEON ---------- */
+
+#if defined(__aarch64__) || defined(_M_ARM64) || defined(__arm__) || defined(_M_ARM)
+
+#include
+
+TARGET_WITH_SIMD static void gf256_mul_region_neon(uint8_t *dst, const uint8_t *src, uint8_t coeff, size_t len)
+{
+ uint8_t low_tbl[16], hi_tbl[16];
+ gf256_build_mul_tables(coeff, low_tbl, hi_tbl);
+
+ const uint8x16_t low_v = vld1q_u8(low_tbl);
+ const uint8x16_t hi_v = vld1q_u8(hi_tbl);
+ const uint8x16_t mask = vdupq_n_u8(0x0F);
+
+ size_t i = 0;
+ for(; i + 16 <= len; i += 16)
+ {
+ uint8x16_t s = vld1q_u8(src + i);
+ uint8x16_t d = vld1q_u8(dst + i);
+ uint8x16_t s_lo = vandq_u8(s, mask);
+ uint8x16_t s_hi = vandq_u8(vshrq_n_u8(s, 4), mask);
+ uint8x16_t lo = vqtbl1q_u8(low_v, s_lo);
+ uint8x16_t hi = vqtbl1q_u8(hi_v, s_hi);
+ uint8x16_t r = veorq_u8(veorq_u8(lo, hi), d);
+ vst1q_u8(dst + i, r);
+ }
+
+ for(; i < len; i++)
+ dst[i] ^= low_tbl[src[i] & 0x0F] ^ hi_tbl[src[i] >> 4];
+}
+
+TARGET_WITH_SIMD static void gf256_xor_region_neon(uint8_t *dst, const uint8_t *src, size_t len)
+{
+ size_t i = 0;
+ for(; i + 16 <= len; i += 16)
+ {
+ uint8x16_t d = vld1q_u8(dst + i);
+ uint8x16_t s = vld1q_u8(src + i);
+ vst1q_u8(dst + i, veorq_u8(d, s));
+ }
+ for(; i < len; i++) dst[i] ^= src[i];
+}
+
+int have_neon(void);
+
+#endif /* ARM */
+
+/* -------------------------------------------------------------------------
+ * Dispatch functions
+ * ------------------------------------------------------------------------- */
+
+void gf256_mul_region(uint8_t *dst, const uint8_t *src, uint8_t coeff, size_t len)
+{
+ if(coeff == 0) return;
+ if(coeff == 1) { gf256_xor_region(dst, src, len); return; }
+
+ gf256_init_tables();
+
+#if defined(__x86_64__) || defined(__amd64) || defined(_M_AMD64) || defined(_M_X64) || \
+ defined(__I386__) || defined(__i386__) || defined(__THW_INTEL) || defined(_M_IX86)
+ if(have_avx2()) { gf256_mul_region_avx2(dst, src, coeff, len); return; }
+ if(have_ssse3()) { gf256_mul_region_ssse3(dst, src, coeff, len); return; }
+#endif
+
+#if defined(__aarch64__) || defined(_M_ARM64) || defined(__arm__) || defined(_M_ARM)
+ if(have_neon()) { gf256_mul_region_neon(dst, src, coeff, len); return; }
+#endif
+
+ gf256_mul_region_scalar(dst, src, coeff, len);
+}
+
+void gf256_xor_region(uint8_t *dst, const uint8_t *src, size_t len)
+{
+#if defined(__x86_64__) || defined(__amd64) || defined(_M_AMD64) || defined(_M_X64) || \
+ defined(__I386__) || defined(__i386__) || defined(__THW_INTEL) || defined(_M_IX86)
+ if(have_avx2()) { gf256_xor_region_avx2(dst, src, len); return; }
+ if(have_ssse3()) { gf256_xor_region_ssse3(dst, src, len); return; }
+#endif
+
+#if defined(__aarch64__) || defined(_M_ARM64) || defined(__arm__) || defined(_M_ARM)
+ if(have_neon()) { gf256_xor_region_neon(dst, src, len); return; }
+#endif
+
+ gf256_xor_region_scalar(dst, src, len);
+}
diff --git a/src/lib/gf256.h b/src/lib/gf256.h
new file mode 100644
index 0000000..2d6f487
--- /dev/null
+++ b/src/lib/gf256.h
@@ -0,0 +1,72 @@
+/*
+ * This file is part of the Aaru Data Preservation Suite.
+ * Copyright (c) 2019-2026 Natalia Portillo.
+ *
+ * This library is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation; either version 2.1 of the
+ * License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, see .
+ */
+
+#ifndef LIBAARUFORMAT_GF256_H
+#define LIBAARUFORMAT_GF256_H
+
+#include
+#include
+
+/**
+ * @brief Multiply two elements in GF(2^8) with polynomial 0x11D.
+ * @param a First operand.
+ * @param b Second operand.
+ * @return Product a*b in GF(2^8).
+ */
+uint8_t gf256_mul(uint8_t a, uint8_t b);
+
+/**
+ * @brief Divide two elements in GF(2^8).
+ * @param a Dividend.
+ * @param b Divisor (must be non-zero).
+ * @return Quotient a/b in GF(2^8).
+ */
+uint8_t gf256_div(uint8_t a, uint8_t b);
+
+/**
+ * @brief Compute multiplicative inverse in GF(2^8).
+ * @param a Element (must be non-zero).
+ * @return Inverse a^(-1) in GF(2^8).
+ */
+uint8_t gf256_inv(uint8_t a);
+
+/**
+ * @brief Multiply-accumulate a region: dst[i] ^= GF_mul(src[i], coeff) for all i.
+ *
+ * Uses SIMD acceleration when available (AVX2 > SSSE3 > NEON > scalar).
+ * If coeff is 0, this is a no-op. If coeff is 1, this is XOR.
+ *
+ * @param dst Destination buffer (read-modify-write).
+ * @param src Source buffer (read-only).
+ * @param coeff GF(2^8) coefficient.
+ * @param len Number of bytes to process.
+ */
+void gf256_mul_region(uint8_t *dst, const uint8_t *src, uint8_t coeff, size_t len);
+
+/**
+ * @brief XOR a region: dst[i] ^= src[i] for all i.
+ *
+ * Uses SIMD acceleration when available.
+ *
+ * @param dst Destination buffer (read-modify-write).
+ * @param src Source buffer (read-only).
+ * @param len Number of bytes to process.
+ */
+void gf256_xor_region(uint8_t *dst, const uint8_t *src, size_t len);
+
+#endif /* LIBAARUFORMAT_GF256_H */
diff --git a/src/lib/reed_solomon.c b/src/lib/reed_solomon.c
new file mode 100644
index 0000000..a22d6a6
--- /dev/null
+++ b/src/lib/reed_solomon.c
@@ -0,0 +1,300 @@
+/*
+ * This file is part of the Aaru Data Preservation Suite.
+ * Copyright (c) 2019-2026 Natalia Portillo.
+ *
+ * This library is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation; either version 2.1 of the
+ * License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, see .
+ */
+
+/**
+ * @file reed_solomon.c
+ * @brief Reed-Solomon erasure codec over GF(2^8).
+ *
+ * Implements RS(K, M) encoding and decoding using a Vandermonde-derived
+ * generator matrix. Encoding is incremental (one data shard at a time).
+ * Decoding uses Gaussian elimination to reconstruct erased shards.
+ *
+ * The coding matrix is (K+M) x K where:
+ * - Top K rows form an identity matrix (data shards pass through unchanged)
+ * - Bottom M rows are the generator matrix (parity = G * data)
+ *
+ * Generator matrix construction:
+ * Start with a (K+M) x K Vandermonde matrix V where V[i][j] = i^j in GF(2^8).
+ * Invert the top K x K submatrix and multiply the entire matrix by the inverse
+ * so that the top K rows become identity. The bottom M rows are the generator.
+ */
+
+#include
+#include
+
+#include "reed_solomon.h"
+#include "gf256.h"
+
+struct rs_context
+{
+ uint16_t K; /**< Number of data shards. */
+ uint16_t M; /**< Number of parity shards. */
+ uint8_t *gen; /**< Generator matrix: M rows x K columns (row-major). */
+ uint8_t *coding; /**< Full coding matrix: (K+M) rows x K columns. */
+};
+
+/**
+ * @brief Build a Vandermonde matrix (K+M) x K in GF(2^8).
+ *
+ * V[i][j] = i^j in GF(2^8), where i is the row index and j is the column index.
+ * Row 0 is all zeros except column 0 (since 0^0 = 1 by convention, 0^j = 0 for j>0).
+ * We use row indices 0..K+M-1.
+ */
+static uint8_t *build_vandermonde(uint16_t K, uint16_t M)
+{
+ const uint16_t N = K + M;
+ uint8_t *V = calloc((size_t)N * K, sizeof(uint8_t));
+ if(!V) return NULL;
+
+ for(uint16_t i = 0; i < N; i++)
+ {
+ uint8_t val = 1; /* i^0 = 1 */
+ for(uint16_t j = 0; j < K; j++)
+ {
+ V[(size_t)i * K + j] = val;
+ val = gf256_mul(val, (uint8_t)i);
+ }
+ }
+ return V;
+}
+
+/**
+ * @brief Invert a K x K matrix in-place using Gaussian elimination over GF(2^8).
+ *
+ * @param mat The matrix to invert, stored row-major in K*K bytes.
+ * @param inv Output inverse matrix (must be pre-initialized to identity).
+ * @param K Matrix dimension.
+ * @return 0 on success, -1 if singular.
+ */
+static int invert_matrix(const uint8_t *mat, uint8_t *inv, uint16_t K)
+{
+ /* Work on a copy to avoid modifying input */
+ uint8_t *work = malloc((size_t)K * K);
+ if(!work) return -1;
+ memcpy(work, mat, (size_t)K * K);
+
+ /* Initialize inv to identity */
+ memset(inv, 0, (size_t)K * K);
+ for(uint16_t i = 0; i < K; i++)
+ inv[(size_t)i * K + i] = 1;
+
+ /* Forward elimination */
+ for(uint16_t col = 0; col < K; col++)
+ {
+ /* Find pivot */
+ uint16_t pivot = col;
+ while(pivot < K && work[(size_t)pivot * K + col] == 0)
+ pivot++;
+ if(pivot == K) { free(work); return -1; } /* Singular */
+
+ /* Swap rows if needed */
+ if(pivot != col)
+ {
+ for(uint16_t j = 0; j < K; j++)
+ {
+ uint8_t tmp = work[(size_t)col * K + j];
+ work[(size_t)col * K + j] = work[(size_t)pivot * K + j];
+ work[(size_t)pivot * K + j] = tmp;
+
+ tmp = inv[(size_t)col * K + j];
+ inv[(size_t)col * K + j] = inv[(size_t)pivot * K + j];
+ inv[(size_t)pivot * K + j] = tmp;
+ }
+ }
+
+ /* Scale pivot row to make diagonal element 1 */
+ uint8_t diag = work[(size_t)col * K + col];
+ if(diag != 1)
+ {
+ uint8_t inv_diag = gf256_inv(diag);
+ for(uint16_t j = 0; j < K; j++)
+ {
+ work[(size_t)col * K + j] = gf256_mul(work[(size_t)col * K + j], inv_diag);
+ inv[(size_t)col * K + j] = gf256_mul(inv[(size_t)col * K + j], inv_diag);
+ }
+ }
+
+ /* Eliminate column in all other rows */
+ for(uint16_t row = 0; row < K; row++)
+ {
+ if(row == col) continue;
+ uint8_t factor = work[(size_t)row * K + col];
+ if(factor == 0) continue;
+ for(uint16_t j = 0; j < K; j++)
+ {
+ work[(size_t)row * K + j] ^= gf256_mul(factor, work[(size_t)col * K + j]);
+ inv[(size_t)row * K + j] ^= gf256_mul(factor, inv[(size_t)col * K + j]);
+ }
+ }
+ }
+
+ free(work);
+ return 0;
+}
+
+rs_context *rs_create(uint16_t K, uint16_t M)
+{
+ if(K == 0 || M == 0 || (uint32_t)K + M > 255) return NULL;
+
+ rs_context *ctx = calloc(1, sizeof(rs_context));
+ if(!ctx) return NULL;
+ ctx->K = K;
+ ctx->M = M;
+
+ const uint16_t N = K + M;
+
+ /* Build Vandermonde matrix */
+ uint8_t *V = build_vandermonde(K, M);
+ if(!V) { free(ctx); return NULL; }
+
+ /* Invert top K x K submatrix */
+ uint8_t *top_inv = malloc((size_t)K * K);
+ if(!top_inv) { free(V); free(ctx); return NULL; }
+
+ if(invert_matrix(V, top_inv, K) != 0)
+ {
+ free(top_inv);
+ free(V);
+ free(ctx);
+ return NULL;
+ }
+
+ /* Compute coding matrix = V * top_inv^(-1) so top K rows become identity */
+ ctx->coding = calloc((size_t)N * K, sizeof(uint8_t));
+ if(!ctx->coding) { free(top_inv); free(V); free(ctx); return NULL; }
+
+ for(uint16_t i = 0; i < N; i++)
+ {
+ for(uint16_t j = 0; j < K; j++)
+ {
+ uint8_t val = 0;
+ for(uint16_t m = 0; m < K; m++)
+ val ^= gf256_mul(V[(size_t)i * K + m], top_inv[(size_t)m * K + j]);
+ ctx->coding[(size_t)i * K + j] = val;
+ }
+ }
+
+ free(top_inv);
+ free(V);
+
+ /* Generator matrix = bottom M rows of the coding matrix */
+ ctx->gen = ctx->coding + (size_t)K * K;
+
+ return ctx;
+}
+
+void rs_free(rs_context *ctx)
+{
+ if(!ctx) return;
+ free(ctx->coding); /* gen points inside coding, don't free separately */
+ free(ctx);
+}
+
+uint8_t rs_get_coefficient(const rs_context *ctx, uint16_t m, uint16_t k)
+{
+ return ctx->gen[(size_t)m * ctx->K + k];
+}
+
+void rs_encode_incremental(uint8_t coeff, const uint8_t *data, uint8_t *parity, size_t shard_size)
+{
+ gf256_mul_region(parity, data, coeff, shard_size);
+}
+
+int rs_decode(const rs_context *ctx, uint8_t **shards, const uint8_t *present, size_t shard_size)
+{
+ const uint16_t K = ctx->K;
+ const uint16_t M = ctx->M;
+ const uint16_t N = K + M;
+
+ /* Count erasures */
+ uint16_t num_erased = 0;
+ for(uint16_t i = 0; i < N; i++)
+ if(!present[i]) num_erased++;
+
+ if(num_erased == 0) return 0; /* Nothing to do */
+ if(num_erased > M) return -1; /* Too many erasures */
+
+ /* Build the submatrix from rows of the coding matrix corresponding to
+ * the K surviving shards. We need exactly K surviving shards to form
+ * a K x K system. */
+
+ /* Collect indices of surviving shards (pick first K) */
+ uint16_t *surviving = malloc((size_t)K * sizeof(uint16_t));
+ if(!surviving) return -2;
+
+ uint16_t s = 0;
+ for(uint16_t i = 0; i < N && s < K; i++)
+ {
+ if(present[i]) surviving[s++] = i;
+ }
+
+ if(s < K) { free(surviving); return -1; } /* Not enough surviving shards */
+
+ /* Build K x K submatrix from surviving rows of the coding matrix */
+ uint8_t *submat = malloc((size_t)K * K);
+ if(!submat) { free(surviving); return -2; }
+
+ for(uint16_t i = 0; i < K; i++)
+ memcpy(submat + (size_t)i * K, ctx->coding + (size_t)surviving[i] * K, K);
+
+ /* Invert the submatrix */
+ uint8_t *submat_inv = malloc((size_t)K * K);
+ if(!submat_inv) { free(submat); free(surviving); return -2; }
+
+ if(invert_matrix(submat, submat_inv, K) != 0)
+ {
+ free(submat_inv);
+ free(submat);
+ free(surviving);
+ return -1; /* Should not happen if coding matrix is MDS */
+ }
+
+ /* Reconstruct erased shards:
+ * For each erased shard e, compute:
+ * shard[e] = sum over j=0..K-1 of (coding[e][j] * decoded_data[j])
+ *
+ * But decoded_data = submat_inv * surviving_shards
+ * So: shard[e] = sum_j coding[e][j] * (sum_k submat_inv[j][k] * surviving_shards[k])
+ *
+ * Reorder: shard[e] = sum_k (sum_j coding[e][j] * submat_inv[j][k]) * surviving_shards[k]
+ * Let repair_row[e][k] = sum_j coding[e][j] * submat_inv[j][k]
+ */
+ for(uint16_t e = 0; e < N; e++)
+ {
+ if(present[e]) continue;
+
+ /* Compute repair coefficients for this erased shard */
+ memset(shards[e], 0, shard_size);
+
+ for(uint16_t k = 0; k < K; k++)
+ {
+ /* Compute combined coefficient: sum_j coding[e][j] * submat_inv[j][k] */
+ uint8_t coeff = 0;
+ for(uint16_t j = 0; j < K; j++)
+ coeff ^= gf256_mul(ctx->coding[(size_t)e * K + j], submat_inv[(size_t)j * K + k]);
+
+ if(coeff != 0)
+ gf256_mul_region(shards[e], shards[surviving[k]], coeff, shard_size);
+ }
+ }
+
+ free(submat_inv);
+ free(submat);
+ free(surviving);
+ return 0;
+}
diff --git a/src/lib/reed_solomon.h b/src/lib/reed_solomon.h
new file mode 100644
index 0000000..ae5d01c
--- /dev/null
+++ b/src/lib/reed_solomon.h
@@ -0,0 +1,96 @@
+/*
+ * This file is part of the Aaru Data Preservation Suite.
+ * Copyright (c) 2019-2026 Natalia Portillo.
+ *
+ * This library is free software; you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as
+ * published by the Free Software Foundation; either version 2.1 of the
+ * License, or (at your option) any later version.
+ *
+ * This library is distributed in the hope that it will be useful, but
+ * WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with this library; if not, see .
+ */
+
+#ifndef LIBAARUFORMAT_REED_SOLOMON_H
+#define LIBAARUFORMAT_REED_SOLOMON_H
+
+#include
+#include
+
+/**
+ * @brief Opaque Reed-Solomon codec context.
+ */
+typedef struct rs_context rs_context;
+
+/**
+ * @brief Create a Reed-Solomon codec for RS(K, M) over GF(2^8).
+ *
+ * K = number of data shards, M = number of parity shards.
+ * K + M must be <= 255 (GF(2^8) field size minus 1).
+ *
+ * The codec precomputes the Vandermonde-derived generator matrix
+ * (M rows x K columns) used for encoding.
+ *
+ * @param K Number of data shards (>= 1).
+ * @param M Number of parity shards (>= 1).
+ * @return Codec context, or NULL on error.
+ */
+rs_context *rs_create(uint16_t K, uint16_t M);
+
+/**
+ * @brief Free a Reed-Solomon codec context.
+ * @param ctx Context returned by rs_create().
+ */
+void rs_free(rs_context *ctx);
+
+/**
+ * @brief Get the generator matrix coefficient for parity shard m, data shard k.
+ *
+ * During incremental encoding, call this to get the coefficient, then call
+ * rs_encode_incremental() with it.
+ *
+ * @param ctx Codec context.
+ * @param m Parity shard index (0 .. M-1).
+ * @param k Data shard index (0 .. K-1).
+ * @return GF(2^8) coefficient.
+ */
+uint8_t rs_get_coefficient(const rs_context *ctx, uint16_t m, uint16_t k);
+
+/**
+ * @brief Incrementally accumulate one data shard's contribution to one parity shard.
+ *
+ * Computes: parity[i] ^= GF_mul(data[i], coeff) for all i in [0, shard_size).
+ *
+ * This is the core primitive for streaming write: for each data block written,
+ * call this M times (once per parity shard) with the appropriate coefficient.
+ *
+ * For M=1 (XOR-only), coeff is always 1, and this reduces to XOR.
+ *
+ * @param coeff GF(2^8) coefficient from rs_get_coefficient().
+ * @param data Data shard bytes (read-only, shard_size bytes).
+ * @param parity Parity shard accumulator (read-write, shard_size bytes, must be zeroed before first call).
+ * @param shard_size Number of bytes per shard.
+ */
+void rs_encode_incremental(uint8_t coeff, const uint8_t *data, uint8_t *parity, size_t shard_size);
+
+/**
+ * @brief Decode (reconstruct) erased shards.
+ *
+ * Given K+M shards where some are erased, reconstruct the erased ones using
+ * Gaussian elimination over GF(2^8).
+ *
+ * @param ctx Codec context.
+ * @param shards Array of K+M shard pointers (each shard_size bytes). Erased shards must
+ * point to allocated buffers of shard_size bytes (content will be overwritten).
+ * @param present Boolean array of K+M entries: 1 = shard is valid, 0 = shard is erased.
+ * @param shard_size Number of bytes per shard.
+ * @return 0 on success, -1 if too many erasures (> M), -2 on allocation failure.
+ */
+int rs_decode(const rs_context *ctx, uint8_t **shards, const uint8_t *present, size_t shard_size);
+
+#endif /* LIBAARUFORMAT_REED_SOLOMON_H */
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index f26e638..3a26bff 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -102,6 +102,7 @@ add_executable(tests_run
large_file_io.cpp
mode2_nocrc.cpp
mode2_errored.cpp
+ reed_solomon.cpp
${CMAKE_CURRENT_SOURCE_DIR}/../src/lib/aes128.c
${CMAKE_CURRENT_SOURCE_DIR}/../src/ps3/ps3_crypto.c
${CMAKE_CURRENT_SOURCE_DIR}/../src/ps3/ps3_encryption_map.c
@@ -111,6 +112,8 @@ add_executable(tests_run
../tool/ps3/ird.c
../tool/ps3/sfo.c
../tool/ps3/iso9660_mini.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/../src/lib/gf256.c
+ ${CMAKE_CURRENT_SOURCE_DIR}/../src/lib/reed_solomon.c
)
aaru_enable_large_file_support(tests_run)