From c16120ae5be03ae3b465b3bfdeafdd54830d1e2b Mon Sep 17 00:00:00 2001 From: Natalia Portillo Date: Sat, 11 Apr 2026 20:14:43 +0100 Subject: [PATCH] Add GF and Reed Solomon. --- CMakeLists.txt | 7 +- src/lib/gf256.c | 357 +++++++++++++++++++++++++++++++++++++++++ src/lib/gf256.h | 72 +++++++++ src/lib/reed_solomon.c | 300 ++++++++++++++++++++++++++++++++++ src/lib/reed_solomon.h | 96 +++++++++++ tests/CMakeLists.txt | 3 + 6 files changed, 834 insertions(+), 1 deletion(-) create mode 100644 src/lib/gf256.c create mode 100644 src/lib/gf256.h create mode 100644 src/lib/reed_solomon.c create mode 100644 src/lib/reed_solomon.h 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)