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 * 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 */2324const FIELD_SIZE = 256;25const PRIMITIVE_POLY = 0x11d;2627const BLOCK_DATA = 178;28const BLOCK_PARITY = 77;29const BLOCK_TOTAL = BLOCK_DATA + BLOCK_PARITY; // 2553031// ── GF(2^8) tables ──────────────────────────────────────────────────────────3233const EXP_TABLE = new Uint8Array(FIELD_SIZE * 2);34const LOG_TABLE = new Uint8Array(FIELD_SIZE);3536(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})();4748function 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}5253function 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}5859function gfPow(x: number, power: number): number {60  if (x === 0) return 0;61  return EXP_TABLE[(LOG_TABLE[x] * power) % 255];62}6364function 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}6869// ── Polynomial helpers (coefficients high→low; out[0] = x^(n-1)) ────────────7071function 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}7677function 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}8485function 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}9697function 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}102103function 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}110111// Cache generators per parity length — only one is used in practice.112const generatorCache = new Map<number, Uint8Array>();113function 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}121122// ── Block encode/decode ─────────────────────────────────────────────────────123124function 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}141142function 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}150151function findErrorLocator(synd: Uint8Array, parity: number): Uint8Array | null {152  let errLoc: Uint8Array = new Uint8Array([1]);153  let oldLoc: Uint8Array = new Uint8Array([1]);154155  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  }170171  // Strip leading zeros172  let lead = 0;173  while (lead < errLoc.length && errLoc[lead] === 0) lead++;174  errLoc = errLoc.slice(lead);175176  const errs = errLoc.length - 1;177  if (errs * 2 > parity) return null; // too many errors178  return errLoc;179}180181function 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}197198function 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  })();212213  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);217218  // 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  }224225  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}238239function 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);246247  const errLoc = findErrorLocator(synd, parity);248  if (!errLoc) return null;249250  const positions = findErrors(errLoc, received.length);251  if (!positions) return null;252253  let corrected: Uint8Array;254  try {255    corrected = correctErrors(received, synd, positions);256  } catch {257    return null;258  }259260  // 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}267268// ── Public block-mode API ───────────────────────────────────────────────────269270/**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 */276export 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}294295/**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 */303export 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}324325/** Constants the Mayan layer needs for capacity math. */326export const RS_BLOCK_DATA = BLOCK_DATA;327export const RS_BLOCK_PARITY = BLOCK_PARITY;328export const RS_BLOCK_TOTAL = BLOCK_TOTAL;