rs.ts: Reed-Solomon GF(2^8) RS(255,178) block codec with Berlekamp-Massey + Forney; corrects up to 38 byte errors per block (30.2% ECC); 8 round-trip tests pass

@@ -0,0 +1,84 @@1+import { describe, it, expect } from "vitest";2+import {3+  reedSolomonEncode,4+  reedSolomonDecode,5+  RS_BLOCK_DATA,6+  RS_BLOCK_PARITY,7+  RS_BLOCK_TOTAL,8+} from "./rs";9+10+function randomBytes(n: number): Uint8Array {11+  const out = new Uint8Array(n);12+  for (let i = 0; i < n; i++) out[i] = (Math.random() * 256) | 0;13+  return out;14+}15+16+function corruptBytes(buf: Uint8Array, n: number): Uint8Array {17+  const out = new Uint8Array(buf);18+  const positions = new Set<number>();19+  while (positions.size < n) {20+    positions.add((Math.random() * buf.length) | 0);21+  }22+  for (const p of positions) out[p] ^= ((Math.random() * 255) | 0) + 1;23+  return out;24+}25+26+describe("Reed-Solomon RS(255,178)", () => {27+  it("constants are coherent", () => {28+    expect(RS_BLOCK_DATA + RS_BLOCK_PARITY).toBe(RS_BLOCK_TOTAL);29+    expect(RS_BLOCK_TOTAL).toBe(255);30+  });31+32+  it("round-trips a small payload (no errors)", () => {33+    const data = new Uint8Array([0x4d, 0x01, 0x05, 0x00, 1, 2, 3, 4, 5]);34+    const encoded = reedSolomonEncode(data);35+    expect(encoded.length).toBe(RS_BLOCK_TOTAL); // single block36+    const decoded = reedSolomonDecode(encoded, data.length);37+    expect(Array.from(decoded)).toEqual(Array.from(data));38+  });39+40+  it("round-trips a multi-block payload (no errors)", () => {41+    const data = randomBytes(RS_BLOCK_DATA * 3 + 17);42+    const encoded = reedSolomonEncode(data);43+    expect(encoded.length).toBe(RS_BLOCK_TOTAL * 4); // 4 blocks44+    const decoded = reedSolomonDecode(encoded, data.length);45+    expect(Array.from(decoded)).toEqual(Array.from(data));46+  });47+48+  it("corrects 1 byte error per block", () => {49+    const data = randomBytes(RS_BLOCK_DATA);50+    const encoded = reedSolomonEncode(data);51+    const corrupted = corruptBytes(encoded, 1);52+    const decoded = reedSolomonDecode(corrupted, data.length);53+    expect(Array.from(decoded)).toEqual(Array.from(data));54+  });55+56+  it("corrects up to floor(parity/2) byte errors per block", () => {57+    const data = randomBytes(RS_BLOCK_DATA);58+    const encoded = reedSolomonEncode(data);59+    // RS_BLOCK_PARITY=77, can correct up to 38 errors per block60+    const corrupted = corruptBytes(encoded, 38);61+    const decoded = reedSolomonDecode(corrupted, data.length);62+    expect(Array.from(decoded)).toEqual(Array.from(data));63+  });64+65+  it("fails when more than floor(parity/2) errors per block", () => {66+    const data = randomBytes(RS_BLOCK_DATA);67+    const encoded = reedSolomonEncode(data);68+    const corrupted = corruptBytes(encoded, 50); // > 3869+    expect(() => reedSolomonDecode(corrupted, data.length)).toThrow(/uncorrectable/);70+  });71+72+  it("handles empty payload", () => {73+    const encoded = reedSolomonEncode(new Uint8Array(0));74+    expect(encoded.length).toBe(0);75+    const decoded = reedSolomonDecode(encoded, 0);76+    expect(decoded.length).toBe(0);77+  });78+79+  it("rejects malformed received length", () => {80+    expect(() =>81+      reedSolomonDecode(new Uint8Array(RS_BLOCK_TOTAL + 1), 10),82+    ).toThrow(/not a multiple/);83+  });84+});

@@ -0,0 +1,328 @@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+ * Block-based RS(255, 178): 178 data bytes + 77 parity bytes per block.6+ * Yields 30.2% ECC, matching the Mayan spec's ECC_RATE = 0.30.7+ *8+ * The Mayan spec (peach_mayan_spec.md §4) cites "30% ECC" and9+ * "RS(255,223)" in the same sentence, but those numbers contradict —10+ * RS(255,223) is 12.5% ECC. The encoder code in the spec uses11+ * ECC_RATE = 0.30 consistently, so we honor the rate (30%) and pick12+ * standard 255-byte blocks at the matching ratio.13+ *14+ * The exposed API matches the call sites in the Mayan spec:15+ *16+ *   reedSolomonEncode(data, totalParityBytes)17+ *   reedSolomonDecode(received, dataBytes, parityBytes)18+ *19+ * `totalParityBytes` is the spec's `Math.ceil(data.length * ECC_RATE)`;20+ * we treat it as advisory and instead generate 77 parity bytes per21+ * 178-byte block so the symbol-size math matches what the spec expected.22+ */23+24+const FIELD_SIZE = 256;25+const PRIMITIVE_POLY = 0x11d;26+27+const BLOCK_DATA = 178;28+const BLOCK_PARITY = 77;29+const BLOCK_TOTAL = BLOCK_DATA + BLOCK_PARITY; // 25530+31+// ── GF(2^8) tables ──────────────────────────────────────────────────────────32+33+const EXP_TABLE = new Uint8Array(FIELD_SIZE * 2);34+const LOG_TABLE = new Uint8Array(FIELD_SIZE);35+36+(function initTables() {37+  let x = 1;38+  for (let i = 0; i < 255; i++) {39+    EXP_TABLE[i] = x;40+    LOG_TABLE[x] = i;41+    x <<= 1;42+    if (x & 0x100) x ^= PRIMITIVE_POLY;43+  }44+  // Double-length exp table avoids modulo on the hot path.45+  for (let i = 255; i < 512; i++) EXP_TABLE[i] = EXP_TABLE[i - 255];46+})();47+48+function gfMul(a: number, b: number): number {49+  if (a === 0 || b === 0) return 0;50+  return EXP_TABLE[LOG_TABLE[a] + LOG_TABLE[b]];51+}52+53+function gfDiv(a: number, b: number): number {54+  if (b === 0) throw new Error("rs: divide by zero");55+  if (a === 0) return 0;56+  return EXP_TABLE[(LOG_TABLE[a] + 255 - LOG_TABLE[b]) % 255];57+}58+59+function gfPow(x: number, power: number): number {60+  if (x === 0) return 0;61+  return EXP_TABLE[(LOG_TABLE[x] * power) % 255];62+}63+64+function gfInv(x: number): number {65+  if (x === 0) throw new Error("rs: inverse of zero");66+  return EXP_TABLE[(255 - LOG_TABLE[x]) % 255];67+}68+69+// ── Polynomial helpers (coefficients high→low; out[0] = x^(n-1)) ────────────70+71+function polyScale(p: Uint8Array, x: number): Uint8Array {72+  const r = new Uint8Array(p.length);73+  for (let i = 0; i < p.length; i++) r[i] = gfMul(p[i], x);74+  return r;75+}76+77+function polyAdd(p: Uint8Array, q: Uint8Array): Uint8Array {78+  const len = Math.max(p.length, q.length);79+  const r = new Uint8Array(len);80+  for (let i = 0; i < p.length; i++) r[i + len - p.length] = p[i];81+  for (let i = 0; i < q.length; i++) r[i + len - q.length] ^= q[i];82+  return r;83+}84+85+function polyMul(p: Uint8Array, q: Uint8Array): Uint8Array {86+  const r = new Uint8Array(p.length + q.length - 1);87+  for (let j = 0; j < q.length; j++) {88+    if (q[j] === 0) continue;89+    for (let i = 0; i < p.length; i++) {90+      if (p[i] === 0) continue;91+      r[i + j] ^= gfMul(p[i], q[j]);92+    }93+  }94+  return r;95+}96+97+function polyEval(p: Uint8Array, x: number): number {98+  let y = p[0];99+  for (let i = 1; i < p.length; i++) y = gfMul(y, x) ^ p[i];100+  return y;101+}102+103+function generatorPoly(parity: number): Uint8Array {104+  let g: Uint8Array = new Uint8Array([1]);105+  for (let i = 0; i < parity; i++) {106+    g = polyMul(g, new Uint8Array([1, EXP_TABLE[i]]));107+  }108+  return g;109+}110+111+// Cache generators per parity length — only one is used in practice.112+const generatorCache = new Map<number, Uint8Array>();113+function getGenerator(parity: number): Uint8Array {114+  let g = generatorCache.get(parity);115+  if (!g) {116+    g = generatorPoly(parity);117+    generatorCache.set(parity, g);118+  }119+  return g;120+}121+122+// ── Block encode/decode ─────────────────────────────────────────────────────123+124+function encodeBlock(data: Uint8Array, parity: number): Uint8Array {125+  const gen = getGenerator(parity);126+  const out = new Uint8Array(data.length + parity);127+  out.set(data);128+  for (let i = 0; i < data.length; i++) {129+    const coef = out[i];130+    if (coef !== 0) {131+      for (let j = 1; j < gen.length; j++) {132+        out[i + j] ^= gfMul(gen[j], coef);133+      }134+    }135+  }136+  // First `data.length` bytes still hold the original message; replace with137+  // it (it was used as scratch above) and copy the parity tail unchanged.138+  out.set(data);139+  return out;140+}141+142+function calcSyndromes(received: Uint8Array, parity: number): Uint8Array {143+  const synd = new Uint8Array(parity + 1);144+  synd[0] = 0;145+  for (let i = 0; i < parity; i++) {146+    synd[i + 1] = polyEval(received, EXP_TABLE[i]);147+  }148+  return synd;149+}150+151+function findErrorLocator(synd: Uint8Array, parity: number): Uint8Array | null {152+  let errLoc: Uint8Array = new Uint8Array([1]);153+  let oldLoc: Uint8Array = new Uint8Array([1]);154+155+  for (let i = 0; i < parity; i++) {156+    oldLoc = new Uint8Array([...oldLoc, 0]);157+    let delta = synd[i + 1];158+    for (let j = 1; j < errLoc.length; j++) {159+      delta ^= gfMul(errLoc[errLoc.length - 1 - j], synd[i + 1 - j]);160+    }161+    if (delta !== 0) {162+      if (oldLoc.length > errLoc.length) {163+        const newLoc = polyScale(oldLoc, delta);164+        oldLoc = polyScale(errLoc, gfInv(delta));165+        errLoc = newLoc;166+      }167+      errLoc = polyAdd(errLoc, polyScale(oldLoc, delta));168+    }169+  }170+171+  // Strip leading zeros172+  let lead = 0;173+  while (lead < errLoc.length && errLoc[lead] === 0) lead++;174+  errLoc = errLoc.slice(lead);175+176+  const errs = errLoc.length - 1;177+  if (errs * 2 > parity) return null; // too many errors178+  return errLoc;179+}180+181+function findErrors(errLoc: Uint8Array, length: number): number[] | null {182+  // Λ has roots at α^(-poly_degree). Chien finds Λ(α^j)=0 at j = (b+1) mod 255183+  // for byte index b in a length-255 codeword. So b = (j - 1) mod 255.184+  // For our use this is always full RS(255,178) blocks; the modulus is 255.185+  const errs = errLoc.length - 1;186+  const positions: number[] = [];187+  const N = 255;188+  for (let i = 0; i < N; i++) {189+    if (polyEval(errLoc, gfPow(2, i)) === 0) {190+      const byteIdx = (i - 1 + N) % N;191+      if (byteIdx < length) positions.push(byteIdx);192+    }193+  }194+  if (positions.length !== errs) return null;195+  return positions;196+}197+198+function correctErrors(199+  received: Uint8Array,200+  synd: Uint8Array,201+  errPositions: number[],202+): Uint8Array {203+  const coefPos = errPositions.map((p) => received.length - 1 - p);204+  // Forney syndrome205+  const errLocPoly = (() => {206+    let p: Uint8Array = new Uint8Array([1]);207+    for (const pos of coefPos) {208+      p = polyMul(p, polyAdd(new Uint8Array([1]), new Uint8Array([gfPow(2, pos), 0])));209+    }210+    return p;211+  })();212+213+  const syndRev = new Uint8Array(synd.length - 1);214+  for (let i = 0; i < syndRev.length; i++) syndRev[i] = synd[synd.length - 1 - i];215+  let omega: Uint8Array = polyMul(syndRev, errLocPoly);216+  omega = omega.slice(omega.length - errLocPoly.length + 1);217+218+  // Formal derivative of errLocPoly (skip even-degree terms)219+  const errLocDeriv = new Uint8Array(errLocPoly.length - 1);220+  for (let i = 0; i < errLocDeriv.length; i++) {221+    const deg = errLocPoly.length - 1 - i;222+    if (deg % 2 === 1) errLocDeriv[i] = errLocPoly[i];223+  }224+225+  const out = new Uint8Array(received);226+  for (let i = 0; i < coefPos.length; i++) {227+    const pos = coefPos[i];228+    const xi = gfPow(2, pos);229+    const xiInv = gfInv(xi);230+    const num = polyEval(omega, xiInv);231+    const den = polyEval(errLocDeriv, xiInv);232+    if (den === 0) throw new Error("rs: denominator zero in Forney");233+    const magnitude = gfMul(xi, gfDiv(num, den));234+    out[received.length - 1 - pos] ^= magnitude;235+  }236+  return out;237+}238+239+function decodeBlock(received: Uint8Array, parity: number): Uint8Array | null {240+  const synd = calcSyndromes(received, parity);241+  let hasError = false;242+  for (let i = 1; i < synd.length; i++) {243+    if (synd[i] !== 0) { hasError = true; break; }244+  }245+  if (!hasError) return received.slice(0, received.length - parity);246+247+  const errLoc = findErrorLocator(synd, parity);248+  if (!errLoc) return null;249+250+  const positions = findErrors(errLoc, received.length);251+  if (!positions) return null;252+253+  let corrected: Uint8Array;254+  try {255+    corrected = correctErrors(received, synd, positions);256+  } catch {257+    return null;258+  }259+260+  // Verify by recomputing syndromes261+  const synd2 = calcSyndromes(corrected, parity);262+  for (let i = 1; i < synd2.length; i++) {263+    if (synd2[i] !== 0) return null;264+  }265+  return corrected.slice(0, received.length - parity);266+}267+268+// ── Public block-mode API ───────────────────────────────────────────────────269+270+/**271+ * Encode `data` with Reed-Solomon. Splits into BLOCK_DATA-byte data blocks272+ * (zero-padded), appends BLOCK_PARITY parity bytes per block, returns the273+ * concatenation. `_advisoryParityBytes` is accepted for spec compatibility274+ * and ignored — the actual parity volume is fixed by the block geometry.275+ */276+export function reedSolomonEncode(277+  data: Uint8Array,278+  _advisoryParityBytes?: number,279+): Uint8Array {280+  if (data.length === 0) return new Uint8Array(0);281+  const blocks = Math.ceil(data.length / BLOCK_DATA);282+  const out = new Uint8Array(blocks * BLOCK_TOTAL);283+  for (let b = 0; b < blocks; b++) {284+    const start = b * BLOCK_DATA;285+    const slice = data.slice(start, start + BLOCK_DATA);286+    const padded = slice.length === BLOCK_DATA287+      ? slice288+      : (() => { const p = new Uint8Array(BLOCK_DATA); p.set(slice); return p; })();289+    const enc = encodeBlock(padded, BLOCK_PARITY);290+    out.set(enc, b * BLOCK_TOTAL);291+  }292+  return out;293+}294+295+/**296+ * Decode an RS-encoded stream. Splits into BLOCK_TOTAL-byte received blocks297+ * and decodes each. `dataBytes` is the original payload length; the result298+ * is trimmed accordingly. Throws if any block fails (too many errors).299+ *300+ * `_parityBytes` is accepted for spec compatibility; the per-block parity301+ * count is fixed by the block geometry.302+ */303+export function reedSolomonDecode(304+  received: Uint8Array,305+  dataBytes: number,306+  _parityBytes?: number,307+): Uint8Array {308+  if (dataBytes === 0) return new Uint8Array(0);309+  if (received.length % BLOCK_TOTAL !== 0) {310+    throw new Error(311+      `rs: received length ${received.length} not a multiple of ${BLOCK_TOTAL}`,312+    );313+  }314+  const blocks = received.length / BLOCK_TOTAL;315+  const out = new Uint8Array(blocks * BLOCK_DATA);316+  for (let b = 0; b < blocks; b++) {317+    const slice = received.slice(b * BLOCK_TOTAL, (b + 1) * BLOCK_TOTAL);318+    const decoded = decodeBlock(slice, BLOCK_PARITY);319+    if (!decoded) throw new Error(`rs: block ${b} uncorrectable`);320+    out.set(decoded, b * BLOCK_DATA);321+  }322+  return out.slice(0, dataBytes);323+}324+325+/** Constants the Mayan layer needs for capacity math. */326+export const RS_BLOCK_DATA = BLOCK_DATA;327+export const RS_BLOCK_PARITY = BLOCK_PARITY;328+export const RS_BLOCK_TOTAL = BLOCK_TOTAL;