rs.ts

1/**2 * Reed-Solomon codec — GF(2^8) with primitive polynomial 0x11D3 * (x^8 + x^4 + x^3 + x^2 + 1, the standard CCSDS / NASA polynomial).4 *5 * Two API surfaces:6 *7 * 1. Per-block primitives (Mayan v2): rsEncodeBlock / rsDecodeBlock with8 *    arbitrary parity count. Used for RS(255,204) — 204 data + 51 parity —9 *    where the Mayan layer does its own blocking and interleaving.10 *11 * 2. Legacy block-mode RS(255,178): reedSolomonEncode / reedSolomonDecode12 *    auto-block a payload into 178/77 blocks. Used by Mayan v1 codes only.13 *    New v2 codes do not go through this path.14 */1516const FIELD_SIZE = 256;17const PRIMITIVE_POLY = 0x11d;1819const BLOCK_DATA = 178;20const BLOCK_PARITY = 77;21const BLOCK_TOTAL = BLOCK_DATA + BLOCK_PARITY; // 2552223// ── GF(2^8) tables ──────────────────────────────────────────────────────────2425const EXP_TABLE = new Uint8Array(FIELD_SIZE * 2);26const LOG_TABLE = new Uint8Array(FIELD_SIZE);2728(function initTables() {29  let x = 1;30  for (let i = 0; i < 255; i++) {31    EXP_TABLE[i] = x;32    LOG_TABLE[x] = i;33    x <<= 1;34    if (x & 0x100) x ^= PRIMITIVE_POLY;35  }36  // Double-length exp table avoids modulo on the hot path.37  for (let i = 255; i < 512; i++) EXP_TABLE[i] = EXP_TABLE[i - 255];38})();3940function gfMul(a: number, b: number): number {41  if (a === 0 || b === 0) return 0;42  return EXP_TABLE[LOG_TABLE[a] + LOG_TABLE[b]];43}4445function gfDiv(a: number, b: number): number {46  if (b === 0) throw new Error("rs: divide by zero");47  if (a === 0) return 0;48  return EXP_TABLE[(LOG_TABLE[a] + 255 - LOG_TABLE[b]) % 255];49}5051function gfPow(x: number, power: number): number {52  if (x === 0) return 0;53  return EXP_TABLE[(LOG_TABLE[x] * power) % 255];54}5556function gfInv(x: number): number {57  if (x === 0) throw new Error("rs: inverse of zero");58  return EXP_TABLE[(255 - LOG_TABLE[x]) % 255];59}6061// ── Polynomial helpers (coefficients high→low; out[0] = x^(n-1)) ────────────6263function polyScale(p: Uint8Array, x: number): Uint8Array {64  const r = new Uint8Array(p.length);65  for (let i = 0; i < p.length; i++) r[i] = gfMul(p[i], x);66  return r;67}6869function polyAdd(p: Uint8Array, q: Uint8Array): Uint8Array {70  const len = Math.max(p.length, q.length);71  const r = new Uint8Array(len);72  for (let i = 0; i < p.length; i++) r[i + len - p.length] = p[i];73  for (let i = 0; i < q.length; i++) r[i + len - q.length] ^= q[i];74  return r;75}7677function polyMul(p: Uint8Array, q: Uint8Array): Uint8Array {78  const r = new Uint8Array(p.length + q.length - 1);79  for (let j = 0; j < q.length; j++) {80    if (q[j] === 0) continue;81    for (let i = 0; i < p.length; i++) {82      if (p[i] === 0) continue;83      r[i + j] ^= gfMul(p[i], q[j]);84    }85  }86  return r;87}8889function polyEval(p: Uint8Array, x: number): number {90  let y = p[0];91  for (let i = 1; i < p.length; i++) y = gfMul(y, x) ^ p[i];92  return y;93}9495function generatorPoly(parity: number): Uint8Array {96  let g: Uint8Array = new Uint8Array([1]);97  for (let i = 0; i < parity; i++) {98    g = polyMul(g, new Uint8Array([1, EXP_TABLE[i]]));99  }100  return g;101}102103// Cache generators per parity length — only one is used in practice.104const generatorCache = new Map<number, Uint8Array>();105function getGenerator(parity: number): Uint8Array {106  let g = generatorCache.get(parity);107  if (!g) {108    g = generatorPoly(parity);109    generatorCache.set(parity, g);110  }111  return g;112}113114// ── Block encode/decode ─────────────────────────────────────────────────────115116function encodeBlock(data: Uint8Array, parity: number): Uint8Array {117  const gen = getGenerator(parity);118  const out = new Uint8Array(data.length + parity);119  out.set(data);120  for (let i = 0; i < data.length; i++) {121    const coef = out[i];122    if (coef !== 0) {123      for (let j = 1; j < gen.length; j++) {124        out[i + j] ^= gfMul(gen[j], coef);125      }126    }127  }128  // First `data.length` bytes still hold the original message; replace with129  // it (it was used as scratch above) and copy the parity tail unchanged.130  out.set(data);131  return out;132}133134function calcSyndromes(received: Uint8Array, parity: number): Uint8Array {135  const synd = new Uint8Array(parity + 1);136  synd[0] = 0;137  for (let i = 0; i < parity; i++) {138    synd[i + 1] = polyEval(received, EXP_TABLE[i]);139  }140  return synd;141}142143function findErrorLocator(synd: Uint8Array, parity: number): Uint8Array | null {144  let errLoc: Uint8Array = new Uint8Array([1]);145  let oldLoc: Uint8Array = new Uint8Array([1]);146147  for (let i = 0; i < parity; i++) {148    oldLoc = new Uint8Array([...oldLoc, 0]);149    let delta = synd[i + 1];150    for (let j = 1; j < errLoc.length; j++) {151      delta ^= gfMul(errLoc[errLoc.length - 1 - j], synd[i + 1 - j]);152    }153    if (delta !== 0) {154      if (oldLoc.length > errLoc.length) {155        const newLoc = polyScale(oldLoc, delta);156        oldLoc = polyScale(errLoc, gfInv(delta));157        errLoc = newLoc;158      }159      errLoc = polyAdd(errLoc, polyScale(oldLoc, delta));160    }161  }162163  // Strip leading zeros164  let lead = 0;165  while (lead < errLoc.length && errLoc[lead] === 0) lead++;166  errLoc = errLoc.slice(lead);167168  const errs = errLoc.length - 1;169  if (errs * 2 > parity) return null; // too many errors170  return errLoc;171}172173function findErrors(errLoc: Uint8Array, length: number): number[] | null {174  // Λ has roots at α^(-poly_degree). Chien finds Λ(α^j)=0 at j = (b+1) mod 255175  // for byte index b in a length-255 codeword. So b = (j - 1) mod 255.176  // For our use this is always full RS(255,178) blocks; the modulus is 255.177  const errs = errLoc.length - 1;178  const positions: number[] = [];179  const N = 255;180  for (let i = 0; i < N; i++) {181    if (polyEval(errLoc, gfPow(2, i)) === 0) {182      const byteIdx = (i - 1 + N) % N;183      if (byteIdx < length) positions.push(byteIdx);184    }185  }186  if (positions.length !== errs) return null;187  return positions;188}189190function correctErrors(191  received: Uint8Array,192  synd: Uint8Array,193  errPositions: number[],194): Uint8Array {195  const coefPos = errPositions.map((p) => received.length - 1 - p);196  // Forney syndrome197  const errLocPoly = (() => {198    let p: Uint8Array = new Uint8Array([1]);199    for (const pos of coefPos) {200      p = polyMul(p, polyAdd(new Uint8Array([1]), new Uint8Array([gfPow(2, pos), 0])));201    }202    return p;203  })();204205  const syndRev = new Uint8Array(synd.length - 1);206  for (let i = 0; i < syndRev.length; i++) syndRev[i] = synd[synd.length - 1 - i];207  let omega: Uint8Array = polyMul(syndRev, errLocPoly);208  omega = omega.slice(omega.length - errLocPoly.length + 1);209210  // Formal derivative of errLocPoly (skip even-degree terms)211  const errLocDeriv = new Uint8Array(errLocPoly.length - 1);212  for (let i = 0; i < errLocDeriv.length; i++) {213    const deg = errLocPoly.length - 1 - i;214    if (deg % 2 === 1) errLocDeriv[i] = errLocPoly[i];215  }216217  const out = new Uint8Array(received);218  for (let i = 0; i < coefPos.length; i++) {219    const pos = coefPos[i];220    const xi = gfPow(2, pos);221    const xiInv = gfInv(xi);222    const num = polyEval(omega, xiInv);223    const den = polyEval(errLocDeriv, xiInv);224    if (den === 0) throw new Error("rs: denominator zero in Forney");225    const magnitude = gfMul(xi, gfDiv(num, den));226    out[received.length - 1 - pos] ^= magnitude;227  }228  return out;229}230231function decodeBlock(received: Uint8Array, parity: number): Uint8Array | null {232  const synd = calcSyndromes(received, parity);233  let hasError = false;234  for (let i = 1; i < synd.length; i++) {235    if (synd[i] !== 0) { hasError = true; break; }236  }237  if (!hasError) return received.slice(0, received.length - parity);238239  const errLoc = findErrorLocator(synd, parity);240  if (!errLoc) return null;241242  const positions = findErrors(errLoc, received.length);243  if (!positions) return null;244245  let corrected: Uint8Array;246  try {247    corrected = correctErrors(received, synd, positions);248  } catch {249    return null;250  }251252  // Verify by recomputing syndromes253  const synd2 = calcSyndromes(corrected, parity);254  for (let i = 1; i < synd2.length; i++) {255    if (synd2[i] !== 0) return null;256  }257  return corrected.slice(0, received.length - parity);258}259260// ── Public block-mode API ───────────────────────────────────────────────────261262/**263 * Encode `data` with Reed-Solomon. Splits into BLOCK_DATA-byte data blocks264 * (zero-padded), appends BLOCK_PARITY parity bytes per block, returns the265 * concatenation. `_advisoryParityBytes` is accepted for spec compatibility266 * and ignored — the actual parity volume is fixed by the block geometry.267 */268export function reedSolomonEncode(269  data: Uint8Array,270  _advisoryParityBytes?: number,271): Uint8Array {272  if (data.length === 0) return new Uint8Array(0);273  const blocks = Math.ceil(data.length / BLOCK_DATA);274  const out = new Uint8Array(blocks * BLOCK_TOTAL);275  for (let b = 0; b < blocks; b++) {276    const start = b * BLOCK_DATA;277    const slice = data.slice(start, start + BLOCK_DATA);278    const padded = slice.length === BLOCK_DATA279      ? slice280      : (() => { const p = new Uint8Array(BLOCK_DATA); p.set(slice); return p; })();281    const enc = encodeBlock(padded, BLOCK_PARITY);282    out.set(enc, b * BLOCK_TOTAL);283  }284  return out;285}286287/**288 * Decode an RS-encoded stream. Splits into BLOCK_TOTAL-byte received blocks289 * and decodes each. `dataBytes` is the original payload length; the result290 * is trimmed accordingly. Throws if any block fails (too many errors).291 *292 * `_parityBytes` is accepted for spec compatibility; the per-block parity293 * count is fixed by the block geometry.294 */295export function reedSolomonDecode(296  received: Uint8Array,297  dataBytes: number,298  _parityBytes?: number,299): Uint8Array {300  if (dataBytes === 0) return new Uint8Array(0);301  if (received.length % BLOCK_TOTAL !== 0) {302    throw new Error(303      `rs: received length ${received.length} not a multiple of ${BLOCK_TOTAL}`,304    );305  }306  const blocks = received.length / BLOCK_TOTAL;307  const out = new Uint8Array(blocks * BLOCK_DATA);308  for (let b = 0; b < blocks; b++) {309    const slice = received.slice(b * BLOCK_TOTAL, (b + 1) * BLOCK_TOTAL);310    const decoded = decodeBlock(slice, BLOCK_PARITY);311    if (!decoded) throw new Error(`rs: block ${b} uncorrectable`);312    out.set(decoded, b * BLOCK_DATA);313  }314  return out.slice(0, dataBytes);315}316317/** Constants the Mayan layer needs for capacity math. */318export const RS_BLOCK_DATA = BLOCK_DATA;319export const RS_BLOCK_PARITY = BLOCK_PARITY;320export const RS_BLOCK_TOTAL = BLOCK_TOTAL;321322// ── Per-block API (Mayan v2 uses these directly with parity=51) ─────────────323324/**325 * Encode a single shortened RS codeword: append `parity` parity bytes to326 * `data`, returning a `data.length + parity`-byte block. Caller controls327 * the block geometry; for RS(255,204) pass 204-byte data + parity=51.328 */329export function rsEncodeBlock(data: Uint8Array, parity: number): Uint8Array {330  return encodeBlock(data, parity);331}332333/**334 * Decode a single shortened RS codeword. Returns the data portion (without335 * parity) on success, or null if uncorrectable. Block length and parity must336 * match what was used at encode time.337 */338export function rsDecodeBlock(received: Uint8Array, parity: number): Uint8Array | null {339  return decodeBlock(received, parity);340}