aaru-dps/libaaruformat
1
0
Fork 0
mirror of https://github.com/aaru-dps/libaaruformat.git synced 2026-07-08 18:06:18 +00:00
Code Issues 3 Packages Projects Releases Wiki Activity

Add GF and Reed Solomon.

Browse Source
This commit is contained in:
Natalia Portillo
2026-04-11 20:14:43 +01:00
parent acfa2f0d2f
commit c16120ae5b
6 changed files with 834 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

7
CMakeLists.txt
View File

@@ -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

357
src/lib/gf256.c Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
* @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 <stddef.h>
#include <stdint.h>
#include <string.h>
#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 <tmmintrin.h> /* 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 <immintrin.h> /* 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 <arm_neon.h>
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);
}

72
src/lib/gf256.h Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef LIBAARUFORMAT_GF256_H
#define LIBAARUFORMAT_GF256_H
#include <stddef.h>
#include <stdint.h>
/**
* @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 */

300
src/lib/reed_solomon.c Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
/**
* @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 <stdlib.h>
#include <string.h>
#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;
}

96
src/lib/reed_solomon.h Normal file
View File

@@ -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 <http://www.gnu.org/licenses/>.
*/
#ifndef LIBAARUFORMAT_REED_SOLOMON_H
#define LIBAARUFORMAT_REED_SOLOMON_H
#include <stddef.h>
#include <stdint.h>
/**
* @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 */

3
tests/CMakeLists.txt
View File

@@ -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)