Neda/call/lib/core/fec/reed_solomon.dart

295 lines
11 KiB
Dart
Raw Blame History

This file contains ambiguous Unicode characters

This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.

import 'dart:typed_data';
/// Reed-Solomon codec RS(15,11) over GF(16).
///
/// GF(16) = GF(2^4) with primitive polynomial p(x) = x^4 + x + 1 (0x13).
/// Primitive element α has order 15 (α^15 = 1).
///
/// Parameters:
/// n = 15 (codeword length in nibbles)
/// k = 11 (data symbols per codeword)
/// t = 2 (corrects up to 2 nibble-errors per codeword)
///
/// Wire mapping:
/// 11 data nibbles + 4 parity nibbles = 15 nibbles = 7.5 bytes.
/// To keep byte alignment the codec works on full-byte streams:
/// encode(Uint8List data) where data.length must be a multiple of 5.5 B
/// Callers pack/unpack nibbles via [packNibbles] / [unpackNibbles].
///
/// Usage (encode 11-byte block → 15-byte codeword):
/// final cw = ReedSolomon.encodeBlock(data11Bytes); // Uint8List(15)
/// final ok = ReedSolomon.decodeBlock(cw); // corrects in-place
abstract final class ReedSolomon {
// ── GF(16) tables ─────────────────────────────────────────────────
//
// Computed offline for p(x) = x^4 + x + 1 = 0x13.
// α^0 = 1 α^1 = 2 α^2 = 4 α^3 = 8
// α^4 = 3 α^5 = 6 α^6 = 12 α^7 = 11
// α^8 = 5 α^9 = 10 α^10 = 7 α^11 = 14
// α^12 = 15 α^13 = 13 α^14 = 9 α^15 = 1 (wrap)
static const List<int> _exp = [
1, 2, 4, 8, 3, 6, 12, 11, 5, 10, 7, 14, 15, 13, 9,
// duplicate for modular index without % 15
1, 2, 4, 8, 3, 6, 12, 11, 5, 10, 7, 14, 15, 13, 9,
];
static const List<int> _log = [
// index 0 → undefined (−∞); we store 0 but never use it for non-zero
0, // log[0] = -∞ (sentinel)
0, // log[1] = 0 (α^0)
1, // log[2] = 1
4, // log[3] = 4
2, // log[4] = 2
8, // log[5] = 8
5, // log[6] = 5
10, // log[7] = 10
3, // log[8] = 3
14, // log[9] = 14
9, // log[10] = 9
7, // log[11] = 7
6, // log[12] = 6
13, // log[13] = 13
11, // log[14] = 11
12, // log[15] = 12
];
// Generator polynomial g(x) = (x-α)(x-α^2)(x-α^3)(x-α^4)
// = x^4 + 13x^3 + 12x^2 + 8x + 7 (coefficients in GF(16))
// g[0] is the constant term, g[4] is the leading coefficient (=1).
static const List<int> _gen = [7, 8, 12, 13, 1];
static const int _n = 15;
static const int _k = 11;
static const int _t = 2;
// ── GF(16) arithmetic ─────────────────────────────────────────────
// Addition = XOR in GF(2^m).
static int _add(int a, int b) => a ^ b;
// Multiplication using log/exp lookup tables.
static int _mul(int a, int b) {
if (a == 0 || b == 0) return 0;
return _exp[(_log[a] + _log[b]) % 15];
}
// Division a / b.
static int _div(int a, int b) {
if (a == 0) return 0;
assert(b != 0, 'GF division by zero');
return _exp[(_log[a] - _log[b] + 15) % 15];
}
// Power α^e.
static int _pow(int e) => _exp[e % 15];
// ── Nibble pack / unpack ──────────────────────────────────────────
/// Pack n pairs of nibbles into bytes. nibbles.length must be even.
static Uint8List packNibbles(List<int> nibbles) {
assert(nibbles.length.isEven, 'nibbles must have even length');
final out = Uint8List(nibbles.length >> 1);
for (int i = 0; i < out.length; i++) {
out[i] = ((nibbles[i * 2] & 0xF) << 4) | (nibbles[i * 2 + 1] & 0xF);
}
return out;
}
/// Unpack bytes into a List<int> of nibbles (high nibble first).
static List<int> unpackNibbles(List<int> bytes) {
final out = <int>[];
for (final b in bytes) {
out.add((b >> 4) & 0xF);
out.add(b & 0xF);
}
return out;
}
// ── Encoder ───────────────────────────────────────────────────────
/// Encode [data] (11-nibble list) into a 15-nibble RS codeword.
///
/// The encoding is systematic: the first 11 nibbles of the codeword
/// are the data, the last 4 are parity computed via polynomial division.
static List<int> _encodeNibbles(List<int> data) {
assert(data.length == _k, 'RS encode: input must be $_k nibbles');
// Work buffer: data shifted left by t positions.
final r = List<int>.filled(_k + _t, 0);
for (int i = 0; i < _k; i++) {
r[i] = data[i];
}
// Polynomial long division: r(x) = x^(n-k)·m(x) mod g(x)
for (int i = 0; i < _k; i++) {
final coeff = r[i];
if (coeff != 0) {
for (int j = 1; j <= _t * 2; j++) {
r[i + j] = _add(r[i + j], _mul(_gen[_t * 2 - j], coeff));
}
}
}
// Build codeword: data || parity (last 4 elements of r)
final cw = List<int>.from(data);
for (int i = _k; i < _n; i++) {
cw.add(r[i]);
}
return cw;
}
// ── Public API: byte-level encode (11 bytes → 15 bytes) ───────────
/// Encode an 11-byte data block into a 15-byte RS codeword.
///
/// The 11 bytes are treated as 22 GF(16) nibbles split into two
/// 11-nibble RS codewords. Output is 30 nibbles = 15 bytes.
///
/// [data] must have exactly 11 bytes. Returns Uint8List(15).
static Uint8List encodeBlock(Uint8List data) {
assert(data.length == 11, 'encodeBlock: need 11 bytes, got ${data.length}');
final nibs = unpackNibbles(data); // 22 nibbles
// Codeword 0: nibbles [0..10]
final cw0 = _encodeNibbles(nibs.sublist(0, 11));
// Codeword 1: nibbles [11..21]
final cw1 = _encodeNibbles(nibs.sublist(11, 22));
// Merge 30 nibbles back to 15 bytes
return packNibbles([...cw0, ...cw1]);
}
// ── Decoder (Berlekamp-Massey + Chien + Forney) ───────────────────
/// Decode and correct a 15-nibble codeword (in-place list<int>).
/// Returns true if no uncorrectable error was detected.
static bool _decodeNibbles(List<int> cw) {
assert(cw.length == _n, 'RS decode: codeword must be $_n nibbles');
// ── Step 1: compute syndromes S_i = c(α^i) for i=1..2t ─────────
final s = List<int>.filled(2 * _t, 0);
for (int i = 0; i < 2 * _t; i++) {
int ev = 0;
for (int j = 0; j < _n; j++) {
ev = _add(_mul(ev, _pow(i + 1)), cw[j]);
}
s[i] = ev;
}
// All-zero syndromes → no errors.
if (s.every((v) => v == 0)) return true;
// ── Step 2: Berlekamp-Massey error-locator polynomial Λ(x) ──────
var sigma = List<int>.filled(_t + 1, 0);
sigma[0] = 1;
var prev = List<int>.filled(_t + 1, 0);
prev[0] = 1;
int lfsr = 0, d = 1;
for (int n = 0; n < 2 * _t; n++) {
// Discrepancy Δ = S[n] + Λ_1·S[n-1] + … + Λ_L·S[n-L]
int delta = s[n];
for (int i = 1; i <= lfsr; i++) {
if (n - i >= 0) delta = _add(delta, _mul(sigma[i], s[n - i]));
}
if (delta == 0) {
d++;
} else {
final temp = List<int>.from(sigma);
final scale = _div(delta, _exp[0]); // delta / 1 = delta
// sigma(x) -= (delta / Δ_prev) · x^d · prev(x)
for (int i = d; i <= _t + d && i <= _t; i++) {
if (i - d < prev.length) {
sigma[i] = _add(sigma[i], _mul(scale, prev[i - d]));
}
}
if (2 * lfsr <= n) {
lfsr = n + 1 - lfsr;
prev = temp;
d = 1;
} else {
d++;
}
}
}
// ── Step 3: Chien search — find roots of Λ(x) ──────────────────
final errorPositions = <int>[];
for (int i = 0; i < _n; i++) {
// Evaluate Λ(α^(-i)) = Λ(α^(15-i))
int ev = 0;
for (int j = sigma.length - 1; j >= 0; j--) {
ev = _add(_mul(ev, _pow((15 - i) % 15)), sigma[j]);
}
if (ev == 0) errorPositions.add(i);
}
if (errorPositions.length > _t) return false; // uncorrectable
// ── Step 4: Forney algorithm — compute error magnitudes ─────────
for (final pos in errorPositions) {
// Evaluate error magnitude using Forney formula.
// For t≤2 we use the direct closed-form when |errors|≤2.
// Omega(x) = S(x)·Λ(x) mod x^(2t)
final omega = List<int>.filled(2 * _t, 0);
for (int i = 0; i < 2 * _t; i++) {
for (int j = 0; j <= i && j < sigma.length; j++) {
if (i - j < s.length) {
omega[i] = _add(omega[i], _mul(sigma[j], s[i - j]));
}
}
}
// Λ'(x) formal derivative (odd powers survive in GF(2))
// For binary GF, the derivative has only odd-indexed terms.
int sigmaDerivAtXi = 0;
for (int j = 1; j < sigma.length; j += 2) {
sigmaDerivAtXi = _add(sigmaDerivAtXi,
_mul(sigma[j], _pow(j * (15 - pos) % 15)));
}
if (sigmaDerivAtXi == 0) return false;
// Evaluate Omega at α^(-pos)
int omegaVal = 0;
for (int i = omega.length - 1; i >= 0; i--) {
omegaVal = _add(_mul(omegaVal, _pow((15 - pos) % 15)), omega[i]);
}
// Error value = -X^(-1) · Omega(X^(-1)) / Λ'(X^(-1))
// In GF(2^m): -a = a, so e = _mul(inv(α^pos), omegaVal) / sigmaDerivAtXi
final xi = _pow(pos); // α^pos
final e = _div(_mul(_div(1, xi), omegaVal), sigmaDerivAtXi);
cw[pos] = _add(cw[pos], e);
}
return true;
}
/// Decode and error-correct a 15-byte RS block in-place.
///
/// Returns true if the block is valid (or was corrected successfully).
/// Returns false if errors exceeded the correction capacity.
static bool decodeBlock(Uint8List block) {
assert(block.length == 15, 'decodeBlock: need 15 bytes, got ${block.length}');
final nibs = unpackNibbles(block); // 30 nibbles
final cw0 = nibs.sublist(0, 15);
final cw1 = nibs.sublist(15, 30);
final ok0 = _decodeNibbles(cw0);
final ok1 = _decodeNibbles(cw1);
if (!ok0 || !ok1) return false;
// Write corrected nibbles back to byte array.
final merged = [...cw0, ...cw1];
final packed = packNibbles(merged);
for (int i = 0; i < 15; i++) { block[i] = packed[i]; }
return true;
}
/// Extract just the 11 data bytes from a (possibly corrected) 15-byte block.
static Uint8List extractData(Uint8List block) {
assert(block.length == 15);
final nibs = unpackNibbles(block);
// Data nibbles: [0..10] from cw0 and [15..25] from cw1 → merge
final dataNibs = [...nibs.sublist(0, 11), ...nibs.sublist(15, 26)];
return packNibbles(dataNibs); // 22 nibbles → 11 bytes
}
}