[JavaScript] Programming for Galois Field (and applied to encode/decode Reed-Solomon code)

find-polynomial.js

const egcd = (a, b) => {
  if (a === 0) return [b, 0, 1];    // b = 0*0 + b*1
  const r = b % a, d = (b - r) / a; // r = b - a*d
  const [g, s, t] = egcd(r, a);     // g = r*s + a*t
  return [g, t - d * s, s];         // g = (b-a*d)*s + a*t = a*(t-d*s) + b*s
};
//console.log(egcd(13, 5));
//console.log(egcd(5, 13));

// positive mod
const modp = p => n => (p + (n % p)) % p;
// 1/n mod p: n * 1/n mod p = 1
const invp = (n, p) => {
  const [one, m, q] = egcd(n, p);
  console.assert(one === 1 && n * m + p * q === one, "invp");
  return modp(p)(m);
};
//console.log(invp(5, 13));


// solve equations on mod p
const pivoting = (r, i) => {
  for (let j = i + 1; j < r.length; j++) {
    if (r[j][i] !== 0) {
      [r[i], r[j]] = [r[j], r[i]];
      break;
    }
  }
};
// Gaussian linear equations solver on mod p
// - lineqs: list of c1*s1+c2*s1+...+cn*sn++c=0 as [c1, c2, ..., cn, c]
// - returns [s1, s2, ..., sn]
const solvep = p => lineqs => {
  console.assert(lineqs.every(f => f.length - 1 <= lineqs.length), "solvep");
  const mod = modp(p);
  const r = lineqs.map(f => f.slice());
  const n = r[0].length - 1;
  for (let i = 0; i < n; i++) for (let j = i + 1; j < r.length; j++) {
    if (r[i][i] === 0) pivoting(r, i);
    const f = mod(r[j][i] * invp(r[i][i], p));
    for (let k = i; k < n + 1; k++) r[j][k] = mod(r[j][k] - f * r[i][k]);
  }
  for (let i = n - 1; i >= 0; i--) {
    r[i][n] = mod(r[i][n] * invp(r[i][i], p));
    r[i][i] = 1;
    for (let j = i - 1; j >= 0; j--) {
      r[j][n] = mod(r[j][n] - r[j][i] * r[i][n]);
      r[j][i] = 0;
    }
  }
  return r.slice(0, n).map(l => mod(-l[n]));
};

//BEGIN [gfPow from gfpn.js]
const range = n => [...Array(n).keys()];
const zeros = n => Array(n).fill(0);
const onehot = (n, k) => [...zeros(k), 1, ...zeros(n - k - 1)];
const eql = (a, b) => a.every((ai, i) => ai === b[i]);
const uniq = (a, eql = (e1, e2) => e1 === e2) =>
      a.filter((e1, i) => a.findIndex(e2 => eql(e1, e2)) === i);

const polyAdd = p => (a, b) => {
  if (a.length < b.length) [a, b] = [b, a];
  return a.map((c, i) => i < b.length ? modp(p)(c + (b[i])) : c);
};
const polyTimes = p => (a, m) => a.map((c, i) => modp(p)(c * m));
const polyNeg = p => a => polyTimes(p)(a, -1);
const polySub = p => (a, b) => polyAdd(p)(a, polyNeg(p)(b));
const polySum = p => pls => pls.reduce((r, pl) => polyAdd(p)(r, pl), []);
const polyCarrys = (a, k) => [...zeros(k), ...a];
const polyMul = p => (a, b) => {
  const carrys = polyCarrys, times = polyTimes(p), sum = polySum(p);
  const bn = b.length - 1, n = a.length + bn;
  return sum(b.map((c, k) => carrys(times(a, c), k)));
};
const polyMod = p => (a, b) => {
  console.assert(b[b.length - 1] !== 0);
  const sub = polySub(p), carrys = polyCarrys, times = polyTimes(p);
  const ld = a.length - b.length;
  let r = a;
  for (let i = 0; i <= ld; i++) {
    const c = r[r.length - 1 - i];
    r = sub(r, carrys(times(b, c), ld - i));
  }
  return r.slice(0, b.length - 1);
};
const polyToString = (e, val="a") => {
  const s1 = e.map((v, i) => {
    if (i === 0) return `${v}`;
    if (v === 0) return "";
    const vs = v === 1 ? "" : `${v}`;
    const pow = i === 1 ? "" : `^{${i}}`;
    return `${vs}${val}${pow}`;
  }).filter(e => e);
  const s2 = s1.length > 1 && s1[0] === "0" ? s1.slice(1) : s1;
  return s2.reverse().join("+");
};

const gfOne = (p, n) => onehot(n, 0);
const gfA = (p, n) => onehot(n, 1);
const gfAdd = (p, n) => polyAdd(p);
const gfTimes = (p, n) => polyTimes(p);
const gfMul = (p, n) => (f) => (a,b) => polyMod(p)(polyMul(p)(a,b), f);
const gfPow = (p, n) => f => {
  const mul = gfMul(p, n)(f);
  const one = gfOne(p, n);
  const pow = (e, k) => k === 0 ? one : mul(e, pow(e, k - 1));
  return pow;
};
//END [gfPow from gfpn.js]

// transpose([[1, 2], [3, 4], [5, 6]]) => [[1,3,5],[2,4,6]]
const transpose = a => range(a[0].length).map(k => a.map(e => e[k]));

// build coefficients linear equations for element e on GF(p^n) with f(a)=0
// e.g. a=[0, 1] of GF(3^2) with a^2+2a+2=0 (f=[2, 2, 1])
// c0*1 + c1*e + e^2 = 0 =>
//   c0*[1, 0]+c1*[0, 1]+[1, 1]=[0, 0] as [[1, 0, 1], [0, 1, 1]]
const gfCoeffEquations = (p, n) => f => (e, dim = n) =>
      transpose(range(dim + 1).map(k => gfPow(p, n)(f)(e, k)));
//console.log(gfCoeffEquations(2, 4)([1,1,0,0,1])([0,1,1,0], 2));

const gfFindFormula = (p, n) => f => (e, dim = n) =>
      [...solvep(p)(gfCoeffEquations(p, n)(f)(e, dim)), 1, ...zeros(n - dim)];

const gfExpo = (p, n) => f => expo => gfPow(p, n)(f)(gfA(p, n), expo);
const gfExpoGrouping = (p, n) => {
  const m = p ** n - 1;
  const used = new Set();
  const gs = [];
  for (const expo of range(m)) {
    if (used.has(expo)) continue;
    const g = [expo];
    for (let e = expo * p % m; e !== 0 && e !== expo; e = e * p % m) {
      g.push(e);
    }
    g.forEach(e => used.add(e));
    gs.push(g);
  }
  return gs;
};
//console.log(gfExpoGrouping(3, 2));
//console.log(gfExpoGrouping(5, 2));
//console.log(gfExpoGrouping(2, 4));
//console.log(gfExpoGrouping(2, 6));

const outputGrouping = (p, n, f) => {
  console.log(`[GF(${p}^${n}) with ${polyToString(f)}=0]`);
  const es = gfExpoGrouping(p, n);
  const gs = es.map(expos => expos.map(expo => gfExpo(p, n)(f)(expo)));
  for (const g of gs) {
    const ff = gfFindFormula(p, n)(f)(g[0], g.length);
    console.assert(g.every(e => eql(ff, gfFindFormula(p, n)(f)(e, g.length))),
                   "outputGrouping");
    console.log(`Satisfied ${polyToString(ff, "x")}=0`);
    console.log(g.map(e => `- ${polyToString(e)}`).join("\n"), "\n");
  }
};

// example
{
  const p = 2, n = 6, f = [1,1,0,0,0,0,1];
  //const p = 3, n = 2, f = [2,2,1];
  outputGrouping(p, n, f);
}

gf2n.js

// Specialized codes GF(2^n) from gfpn.js
const range = n => [...Array(n).keys()];
const uniq = a => a.filter((e1, i) => a.findIndex(e2 => e1 === e2) === i);
const msb32 = n => 31 - Math.clz32(n);

// polynomial 2^k as bit list
const polynom2s = k => range(2 ** (k + 1));
const poly2Add = (a, b) => a ^ b;
const poly2Times = (a, m) => m % 2 === 0 ? 0 : a;
const poly2Carrys = (a, k) => a << k;
// mul: case of enough bit (numbit >= 2 * n)
const poly2Mul = (a, b) => {
  let r = 0;
  for (let i = 0; b > 0; i++, b >>>= 1) {
    if (b & 1) r ^= a << i;
  }
  return r;
};
const poly2Mod = n => f => e => {
  const mf = msb32(f);
  for (let i = msb32(e); i >= n; i--) {
    if (e & (1 << i)) e ^= f << (i - mf);
  }
  return e;
};
const poly2ToString = (e, val="a") => {
  const s1 = [...e.toString(2).split("")].reverse().map((v, i) => {
    if (i === 0) return v;
    if (v === "0") return "";
    const pow = i === 1 ? "" : `^{${i}}`;
    return `${val}${pow}`;
  }).filter(e => e);
  const s2 = s1.length > 1 && s1[0] === "0" ? s1.slice(1) : s1;
  return s2.reverse().join("+");
};

// find irreducible polynomals
const kPolynom2s = k => polynom2s(k - 1).map(v => (1 << k) | v);
const reducibles = n => {
  const rs = [];
  for (let i = 1, nh = n >>> 1; i <= nh; i++) {
    const l = kPolynom2s(i), r = kPolynom2s(n - i);
    for (const lp of l) for (const rp of r) rs.push(poly2Mul(lp, rp));
  }
  return uniq(rs);
};
//console.log(reducibles(3).map(v => v.toString(2)));
const irreducibles = n => {
  const reds = reducibles(n);
  return kPolynom2s(n).filter(p => !reds.includes(p));
};
//console.log(irreducibles(3).map(v => v.toString(2)));

// GF(2^n)
const gf2 = n => polynom2s(n - 1);
const gf2Add = n => (a, b) => poly2Add(a, b);
const gf2Times = n => (a, m) => poly2Times(a, m);
const gf2Carry = n => f => e => {
  const overflow = 1 << msb32(f);
  e <<= 1;
  return (e & overflow) ? (e ^ f) : e;
};

//console.log(polyMod(3, 0b1011)(0b10000).toString(2));
const gf2Carrys0 = n => f => (e, k) => poly2Mod(n)(f)(e << k);
const gf2Mul0 = n => f => (a, b) => poly2Mod(n)(f)(poly2Mul(a, b));

// mul: case of ferwer bit
const gf2Carrys = n => f => (e, k) => {
  if (msb32(e) + k <= 32) return gf2Carrys0(n)(f)(e, k);
  for (let i = 0; i < k; i++) e = gf2Carry(n)(f)(e);
  return e;
};
const gf2Mul = n => f => (a, b) => {
  if (msb32(a) + msb32(b) <= 32) return gf2Mul0(n)(f)(a, b);
  let r = 0;
  for (let i = 0; b > 0; i++, b >>>= 1) {
    if (b & 1) r ^= gf2Carrys(n)(f)(a, i);
  }
  return r;
};

// pow
const gf2Pow = n => f => (e, k) => {
  const mul = gf2Mul(n)(f);
  let r = 1;
  for (let i = 0; i < k; i++) r = mul(r, e);
  return r;
};

// f(a)=0 is primitive element of GF(2^n) or not
const gf2PowList = n => f => {
  const carry = gf2Carry(n)(f);
  const r = [];
  for (let k = 0, m = 2 ** n - 1, c = 1; k < m; k++, c = carry(c)) r.push(c);
  return r;
};
const gf2IsPrimitive = n => f => {
  const powList = gf2PowList(n)(f);
  return uniq(powList).length === powList.length;
};

// automorphism
const calcTable = (elems, op) => elems.map(e1 => elems.map(e2 => op(e1, e2)));
const permutateTable = (table, moves) => {
  const rev = moves.map((_, i) => moves.indexOf(i));
  return table.map((_, i) => table[rev[i]].map((_, j, l) => moves[l[rev[j]]]));
};
const eqlTable = (a, b) => a.every((al, i) => al.every(
  (av, j) => av === b[i][j]));


const gf2Moves = n => f => {
  const pn = 2 ** n;
  const powList = gf2PowList(n)(f);
  const revList = range(pn).map(i => i === 0 ? null : powList.indexOf(i));
  return range(pn).map(i => i === 0 ? 0 : powList[revList[i] * 2 % (pn - 1)]); 
};

// applications
const mdTable = (heads, elements) => {
  const maxes = heads.map(
    (h, i) => Math.max(h.length, ...elements.map(l => l[i].length)));
  const top = `| ${heads.map((h, i) => h.padEnd(maxes[i])).join(" | ")} |`;
  const guide = `|-${heads.map((_, i) => "-".repeat(maxes[i])).join("-|-")}-|`;
  const lines = elements.map(
    l => `| ${l.map((e, i) => e.padEnd(maxes[i])).join(" | ")} |`);
  return [top, guide, ...lines].join("\n");
};

const mdCalcTable = (es, table, mark) => {
  const elems = es.map(e => poly2ToString(e));
  const len = Math.max(...elems.map(e => e.length));
  const strList = elems.map(e => `$${e.padEnd(len)}$`);
  const strTable = table.map(
    l => l.map(e => `$${poly2ToString(e).padEnd(len)}$`));
  const heads = [mark, ...strList];
  const elements = strList.map((e, i) => [e, ...strTable[i]]);
  return mdTable(heads, elements);
};

const outputGF2 = (n, f) => {
  console.assert(f.length === n + 1 && gf2IsPrimitive(n)(f),
                 "f is not a polynomial for primitive element");

  // generate GF(p^n) system
  const fs = irreducibles(n);
  const es = gf2(n);
  const add = gf2Add(n);
  const mul = gf2Mul(n)(f);
  const addTable = calcTable(es, add);
  const mulTable = calcTable(es, mul);
  const powList = gf2PowList(n)(f);
  const auto = gf2Moves(n)(f);

  // check automorphism
  console.assert(eqlTable(mulTable, permutateTable(mulTable, auto)));
  console.assert(eqlTable(addTable, permutateTable(addTable, auto)));
  // output
  console.log();
  console.log(
    `### GF(${2 ** n} = ${2}^${n}) with $${poly2ToString(f, "a")}=0$`, `\n`);
  console.log(`Elements of GF($${2}^${n}$)`, "\n");
  console.log("-", es.map(e => `$${poly2ToString(e)}$`).join(", "), "\n");

  console.log(`Power of $a$ mod $${poly2ToString(f)}$`, "\n");
  console.log(mdTable(["order", "number"], powList.map(
    (e, i) => [`$a^{${i}}$`, `$${poly2ToString(e)}$`])), "\n");
  
  console.log(`Calculation Tables of GF($${2}^${n}$)`, "\n");
  console.log(mdCalcTable(es, addTable, "+"), "\n");
  console.log(mdCalcTable(es, mulTable, "*"), "\n");

  const n2s = num => poly2ToString(num);
  console.log(`Generator of automorphism group of GF($${2}^${n}$)`, "\n");
  console.log(mdTable(["from", "to"], auto.map(
    (m, s) => [`$${n2s(s)}$`, `$${n2s(m)}$`])), "\n");
    console.log(`${n}-order irreducible polynomials`, "\n");
  for (const f of fs) {
    const pe = gf2IsPrimitive(n)(f) ? "(primitive element)" : "";
    console.log("-", `$${poly2ToString(f, "x")} = 0$`, pe);
  }
  console.log();
};

// examples
{
  const n = 2, f = 0b111;
  outputGF2(n, f);
}
{
  const n = 3, f = 0b1011;
  outputGF2(n, f);
}
{
  const n = 4, f = 0b10011;
  outputGF2(n, f);
}
{
  //const n = 6, f = 0b1000011;
  //console.log(gf2IsPrimitive(n)(f));
  //console.log(irreducibles(6).filter(gf2IsPrimitive(n)).map(e => poly2ToString(e)));
  //outputGF2(n, f);
}

gfpn.js

// Programming for Galois Field

// 0. Array Utilities
const range = n => [...Array(n).keys()];
const zeros = n => Array(n).fill(0);
const onehot = (n, k) => [...zeros(k), 1, ...zeros(n - k - 1)];
const eql = (a, b) => a.every((ai, i) => ai === b[i]);
const uniq = (a, eql = (e1, e2) => e1 === e2) =>
      a.filter((e1, i) => a.findIndex(e2 => eql(e1, e2)) === i);
//console.log(onehot(3, 0));                      // [1, 0, 0]
//console.log(onehot(3, 2));                      // [0, 0, 1]
//console.log(uniq([1, 3, 2, 1, 2, 3, 4, 2, 5])); // [1, 3, 2, 4, 5]

// 1. mod-p polynomial as coefficient array: a+b*x+c*x^2+... as [a, b, c, ...]
const modp = p => a => (p + a % p) % p;
const polynoms = (p, k) => {
  if (k < 0) return [[]];
  const base = polynoms(p, k - 1);
  return range(p).flatMap(k => base.map(e => [...e, k]));
};
const polyAdd = p => (a, b) => {
  if (a.length < b.length) [a, b] = [b, a];
  return a.map((c, i) => i < b.length ? modp(p)(c + (b[i])) : c);
};
const polyTimes = p => (a, m) => a.map((c, i) => modp(p)(c * m));

const polyNeg = p => a => polyTimes(p)(a, -1);
const polySub = p => (a, b) => polyAdd(p)(a, polyNeg(p)(b));
const polySum = p => pls => pls.reduce((r, pl) => polyAdd(p)(r, pl), []);

const polyCarrys = (a, k) => [...zeros(k), ...a];
const polyMul = p => (a, b) => {
  const carrys = polyCarrys, times = polyTimes(p), sum = polySum(p);
  const bn = b.length - 1, n = a.length + bn;
  return sum(b.map((c, k) => carrys(times(a, c), k)));
};
const polyMod = p => (a, b) => {
  console.assert(b[b.length - 1] !== 0);
  const sub = polySub(p), carrys = polyCarrys, times = polyTimes(p);
  const ld = a.length - b.length;
  let r = a;
  for (let i = 0; i <= ld; i++) {
    const c = r[r.length - 1 - i];
    r = sub(r, carrys(times(b, c), ld - i));
  }
  return r.slice(0, b.length - 1);
};
//console.log(polynoms(2, 1));             // [[0, 0], [1, 0], [0, 1], [1, 1]]
//console.log(polyAdd(2)([1, 0], [1, 1])); // [0, 1]
//console.log(polyTimes(3)([2, 0], 2));    // [1, 0]
//console.log(polyMul(2)([1,1], [1,1,1]));// [1, 0, 0, 1]
//console.log(polyDivMod(2)([0,0,0,1], [1,1,1]));//[1,1], [1, 0]

// polynomial notation as MathJAX latex form
const polyToString = (e, val="a") => {
  const s1 = e.map((v, i) => {
    if (i === 0) return `${v}`;
    if (v === 0) return "";
    const vs = v === 1 ? "" : `${v}`;
    const pow = i === 1 ? "" : `^{${i}}`;
    return `${vs}${val}${pow}`;
  }).filter(e => e);
  const s2 = s1.length > 1 && s1[0] === "0" ? s1.slice(1) : s1;
  return s2.reverse().join("+");
};
//console.log(polyToString([0, 0, 0])); //"0"
//console.log(polyToString([0, 1, 0])); //"a"
//console.log(polyToString([0, 1, 2])); //"2a^{2}+a"


// 2. listing irreducible polynomials (equations)
// list all k-order polymolials (as n-order polynomial = (n+1) elements array)
const kPolynoms = (p, k) => polynoms(p, k - 1).map(e => [...e, 1]);
//console.log(kPolynoms(3, 2, 1));        //[[0, 1, 0], [1, 1, 0], [2, 1, 0]]
//console.log(kPolynoms(3, 2, 2).length); //9 

// all reducible n-order polynomials
const reducibles = (p, n) => {
  const mul = polyMul(p);
  const rs = [];
  for (let i = 1, nh = n >>> 1; i <= nh; i++) {
    const l = kPolynoms(p, i), r = kPolynoms(p, n - i);
    for (const lp of l) for (const rp of r) rs.push(mul(lp, rp));
  }
  return uniq(rs, eql);
};
//console.log(reducibles(3, 2).length); //6
//console.log(reducibles(2, 3).length); //6

// all irreducible n-order polynomials
const irreducibles = (p, n) => {
  const reds = reducibles(p, n);
  return kPolynoms(p, n).filter(p => !reds.find(r => eql(p, r)));
};
//console.log(irreducibles(3, 2));        //[[1, 0, 1], [2, 1, 1], [2, 2, 1]]
//console.log(irreducibles(2, 3));        //[[1, 1, 0, 1], [1, 0, 1, 1]]
//console.log(irreducibles(2, 4).length); //3
//console.log(irreducibles(2, 8).length); //30


// 3. Garois Field system
// elements on GF(p^n)
const gf = (p, n) => polynoms(p, n - 1);
//console.log(gf(3, 2)); //[[0,0],[1,0],[2,0],[0,1],[1,1],[2,1],[0,2],[1,2],[2,2]]
//console.log(gf(2, 3)); //[[0,0,0],[1,0,0],[0,1,0],[1,1,0],[0,0,1],[1,0,1],[0,1,1],[1,1,1]]

// basic G(^n) elements: 0, 1, a
const gfZero = (p, n) => zeros(n);
const gfOne = (p, n) => onehot(n, 0);
const gfA = (p, n) => onehot(n, 1);
//console.log(gfA(3, 2)); //[0. 1]
//console.log(gfA(2, 3)); //[0. 1. 0]

// add two GF(p^n) elements: a+b
const gfAdd = (p, n) => polyAdd(p);
// mul GF(p^n) element with integer: a*k
const gfTimes = (p, n) => polyTimes(p);
// mul two GF(p^n) elements: a*b
const gfMul = (p, n) => (f) => (a,b) => polyMod(p)(polyMul(p)(a,b), f);

// 3.1 Judge solutions of irreducible formula is primitive element or not
// power of GF(p^n) element: e^k
const gfPow = (p, n) => f => {
  const mul = gfMul(p, n)(f);
  const one = gfOne(p, n);
  const pow = (e, k) => k === 0 ? one : mul(e, pow(e, k - 1));
  return pow;
};
// list of power of a=[0, 1, 0, 0, ...]: [a^0, a^1, a^2, ..., a^(p^n - 1)]
const gfPowList = (p, n) => f => {
  const pow = gfPow(p, n)(f);
  const a = gfA(p, n);
  return range(p ** n - 1).map(k => pow(a, k));
};
//console.log(gfPowList(2, 3)([1, 1, 0, 1]));
//console.log(gfPowList(3, 2)([2, 2, 1]));

// answer a=[0, 1, 0, ...] of f(a) = 0 is primitive element of GF(p^n) or not
const gfIsPrimitive = (p, n) => f => {
  const powList = gfPowList(p, n)(f);
  return uniq(powList, eql).length === powList.length;
};
//console.log(gfIsPrimitive(3, 2)([2, 2, 1])); // true
//console.log(gfIsPrimitive(3, 2)([1, 0, 1])); // false (irreducible)
//console.log(gfIsPrimitive(3, 2)([1, 1, 1])); // false (reducible)
//console.log(gfIsPrimitive(3, 2)([0, 0, 1])); // false (reducible)
//console.log(gfIsPrimitive(3, 2)([0, 1, 1])); // false (reducible)

// assign GF(p^n) element e to a polynomial poly(x=e)
const gfApply = (p, n) => f => {
  const add = gfAdd(p, n), times = gfTimes(p, n), pow = gfPow(p, n)(f);
  return (poly, e) => poly.reduce(
    (r, c, i) => add(r, times(pow(e, i), c)), zeros(n));
};
//console.log(gfApply(3, 2)([2,2,1])([2, 2, 1], [0, 1])); // [0, 0]
//console.log(gfApply(3, 2)([2,2,1])([1, 0, 1], [2, 2])); // [0, 0]
//console.log(gfApply(3, 2)([2,2,1])([2, 1, 1], [2, 1])); // [0, 0]

// 4. automorphism of GF
// convert polynomial and number
const polyToNumber = (p, n) => e => e.reduce((r, ei, i) => r + ei * p ** i, 0);
const numberToPoly = (p, n) => i => range(n).map(k => {
  const t = i % (p ** (k + 1)), q = p ** k;
  return (t - (t % q)) / q;
});
//console.log(polyToNumber(3, 2)([2, 2]));
//console.log(numberToPoly(2, 3)(6));

// Calculation table of GF add/mul
const calcTable = (elems, op) => elems.map(e1 => elems.map(e2 => op(e1, e2)));
const polyToNumberTable = (p, n) =>
      table => table.map(l => l.map(e => polyToNumber(p, n)(e)));

// permutation for automorphism
const permutateTable = (table, moves) => {
  const rev = moves.map((_, i) => moves.indexOf(i));
  return table.map((_, i) => table[rev[i]].map((_, j, l) => moves[l[rev[j]]]));
};
const eqlTable = (a, b) => a.every((al, i) => al.every(
  (av, j) => av === b[i][j]));

// generator of automophism group of GF(p^n)
const gfMoves = (p, n) => f => {
  const pn = p ** n;
  const powList = gfPowList(p, n)(f).map(polyToNumber(p, n));
  const revList = range(pn).map(i => i === 0 ? null : powList.indexOf(i));
  return range(pn).map(i => i === 0 ? 0 : powList[revList[i] * p % (pn - 1)]);
};
//console.log(gfMoves(3, 2)([2, 2, 1]));
//console.log(gfMoves(2, 3)([1, 1, 0, 1]));

// 5. application
// output as markdown
const mdTable = (heads, elements) => {
  const maxes = heads.map(
    (h, i) => Math.max(h.length, ...elements.map(l => l[i].length)));
  const top = `| ${heads.map((h, i) => h.padEnd(maxes[i])).join(" | ")} |`;
  const guide = `|-${heads.map((_, i) => "-".repeat(maxes[i])).join("-|-")}-|`;
  const lines = elements.map(
    l => `| ${l.map((e, i) => e.padEnd(maxes[i])).join(" | ")} |`);
  return [top, guide, ...lines].join("\n");
};
const mdCalcTable = (es, table, mark) => {
  const elems = es.map(e => polyToString(e));
  const len = Math.max(...elems.map(e => e.length));
  const strList = elems.map(e => `$${e.padEnd(len)}$`);
  const strTable = table.map(
    l => l.map(e => `$${polyToString(e).padEnd(len)}$`));
  const heads = [mark, ...strList];
  const elements = strList.map((e, i) => [e, ...strTable[i]]);
  return mdTable(heads, elements);
};

const outputGF = (p, n, f) => {
  console.assert(f.length === n + 1 && gfIsPrimitive(p, n)(f),
                 "f is not a polynomial for primitive element");

  // generate GF(p^n) system
  const fs = irreducibles(p, n);
  const es = gf(p, n);
  const add = gfAdd(p, n);
  const mul = gfMul(p, n)(f);
  const addTable = calcTable(es, add);
  const mulTable = calcTable(es, mul);
  const powList = gfPowList(p, n)(f);
  const auto = gfMoves(p, n)(f);

  // check automorphism
  const addNum = polyToNumberTable(p, n)(addTable);
  const mulNum = polyToNumberTable(p, n)(mulTable);
  console.assert(eqlTable(mulNum, permutateTable(mulNum, auto)));
  console.assert(eqlTable(addNum, permutateTable(addNum, auto)));

  // output
  console.log();
  console.log(
    `### GF(${p ** n} = ${p}^${n}) with $${polyToString(f, "a")}=0$`, `\n`);
  console.log(`Elements of GF($${p}^${n}$)`, "\n");
  console.log("-", es.map(e => `$${polyToString(e)}$`).join(", "), "\n");

  console.log(`Power of $a$ mod $${polyToString(f)}$`, "\n");
  console.log(mdTable(["order", "number"], powList.map(
    (e, i) => [`$a^{${i}}$`, `$${polyToString(e)}$`])), "\n");
  
  console.log(`Calculation Tables of GF($${p}^${n}$)`, "\n");
  console.log(mdCalcTable(es, addTable, "+"), "\n");
  console.log(mdCalcTable(es, mulTable, "*"), "\n");

  const n2s = num => polyToString(numberToPoly(p, n)(num));
  console.log(`Generator of automorphism group of GF($${p}^${n}$)`, "\n");
  console.log(mdTable(["from", "to"], auto.map(
    (m, s) => [`$${n2s(s)}$`, `$${n2s(m)}$`])), "\n");

  console.log(`${n}-order irreducible polynomials`, "\n");
  for (const f of fs) {
    const pe = gfIsPrimitive(p, n)(f) ? "(primitive element)" : "";
    console.log("-", `$${polyToString(f, "x")} = 0$`, pe);
  }
  console.log();
};


// examples
{
  const p = 2, n = 2, f = [1, 1, 1];
  outputGF(p, n, f);
}

{
  const p = 3, n = 2, f = [2, 2, 1];
  outputGF(p, n, f);
}

{
  const p = 2, n = 3, f = [1, 1, 0, 1];
  outputGF(p, n, f);
}

{
  //const p = 2, n = 4, f = [1, 1, 0, 0, 1];
  //console.log(irreducibles(p, n).filter(gfIsPrimitive(p, n)));
  //outputGF(p, n, f);
}
{
  //const p = 2, n = 6, f = [1, 1, 0, 0, 0, 0, 1];
  //console.log(irreducibles(p, n).filter(gfIsPrimitive(p, n)));
  //outputGF(p, n, f);
}

reed-solomon.js

// functions from gf2n.js
const range = n => [...Array(n).keys()];
const uniq = a => a.filter((e1, i) => a.findIndex(e2 => e1 === e2) === i);
const msb32 = n => 31 - Math.clz32(n);

const polynom2s = k => range(2 ** (k + 1));
const poly2Add = (a, b) => a ^ b;
const poly2Times = (a, m) => m % 2 === 0 ? 0 : a;

const poly2Mul = (a, b) => {
  let r = 0;
  for (let i = 0; b > 0; i++, b >>>= 1) {
    if (b & 1) r ^= a << i;
  }
  return r;
};

const gf2 = n => polynom2s(n - 1);
const gf2Add = n => (a, b) => poly2Add(a, b);
const gf2Times = n => (a, m) => poly2Times(a, m);
const gf2Carry = n => f => e => {
  const overflow = 1 << msb32(f);
  e <<= 1;
  return (e & overflow) ? (e ^ f) : e;
};

// mul: case of enough bit (numbit >= 2 * n)
const poly2Mod = n => f => e => {
  const mf = msb32(f);
  for (let i = msb32(e); i >= n; i--) {
    if (e & (1 << i)) e ^= f << (i - mf);
  }
  return e;
};
//console.log(polyMod(3, 0b1011)(0b10000).toString(2));
const gf2Carrys0 = n => f => (e, k) => poly2Mod(n)(f)(e << k);
const gf2Mul0 = n => f => (a, b) => poly2Mod(n)(f)(poly2Mul(a, b));

// mul: case of ferwer bit
const gf2Carrys = n => f => (e, k) => {
  if (msb32(e) + k <= 32) return gf2Carrys0(n)(f)(e, k);
  for (let i = 0; i < k; i++) e = gf2Carry(n)(f)(e);
  return e;
};
const gf2Mul = n => f => (a, b) => {
  if (msb32(a) + msb32(b) <= 32) return gf2Mul0(n)(f)(a, b);
  let r = 0;
  for (let i = 0; b > 0; i++, b >>>= 1) {
    if (b & 1) r ^= gf2Carrys(n)(f)(a, i);
  }
  return r;
};

// pow
const gf2Pow = n => f => (e, k) => {
  const mul = gf2Mul(n)(f);
  let r = 1;
  for (let i = 0; i < k; i++) r = mul(r, e);
  return r;
};

// f(a)=0 is primitive element of GF(2^n) or not
const gf2PowList = n => f => {
  const carry = gf2Carry(n)(f);
  const r = [];
  for (let k = 0, m = 2 ** n - 1, c = 1; k < m; k++, c = carry(c)) r.push(c);
  return r;
};
// gf2n.js

// Appendix: Reed-Solomon code
// Coefficient is GF(2^n) => reversed array of GF(2^n) elements
{
  console.log("[Reed-Solomon error detective code]");
  // GF(2^8) System for QR Code
  const n = 8, f = 0b100011101; //a^8+a^4+a^3+a^2+1=0
  const zero = 0, one = 1, a = 2, isZero = e => e === zero;
  const add = gf2Add(n), times = gf2Times(n), mul = gf2Mul(n)(f), pow = gf2Pow(n)(f);

  // example from https://kouyama.sci.u-toyama.ac.jp/main/education/2007/infomath/pdf/text/text09.pdf
  const N = 26, K = 19, t = 3, b = 0; // N-K > 2*t
  const msg = [
    0b10000000, 0b01000100, 0b10000101, 0b10100111, 0b01001001, 0b10100111,
    0b10001011, 0b01101100, 0b00000000, 0b11101100, 0b00010001, 0b11101100,
    0b00010001, 0b11101100, 0b00010001, 0b11101100, 0b00010001, 0b11101100,
    0b00010001,
  ];

  // utilities for RS encode
  const combinations = (a, n) => {
    if (n === 0) return [[]];
    return range(a.length - n + 1).flatMap(
      k => combinations(a.slice(k + 1), n - 1).map(s => [a[k], ...s]));
  };
  const zeros = n => Array(n).fill(0);
  const sum = a => a.reduce((r, e) => r + e, 0); // sum for usual numbers

  const neg = e => times(e, -1);
  const sub = (a, b) => add(a, neg(b));
  const adds = (es) => es.reduce((r, e) => add(r, e), zero);
  const muls = (es) => es.reduce((r, e) => mul(r, e), one);

  // 1. Find generator polynomial G(x)
  // reverse order polynomial coefficients for `mod gen`
  const rsGen = (N, K, b) => {
    const ss = range(N - K).map(v => v + b);
    return range(N - K + 1).map(k => {
      const c = adds(combinations(ss, k).map(g => pow(a, sum(g))));
      return times(c, (-1) ** (N - K - k));
    });
  };
  const gen = rsGen(N, K, b);
  console.log("Generator polynomial", gen);

  // 2. RS Encode
  // NOTE: message byte sequence as reversed order: [x^{N-1}, x^{N-2}..., x, 1]
  const deconv = (a, g) => {
    if (a.length < g.length) return [[], a];
    const k = a.length - g.length + 1;
    const d = [], r = a.slice();
    for (let i = 0; i < k; i++) {
      const c = r[i];
      d.push(c);
      for (let j = 0; j < g.length; j++) {
        r[i + j] = sub(r[i + j], mul(g[j], c));
      }
    }
    return [d, r.slice(k)];
  };
  const rsRx = (N, K, msg, gen) => {
    const padded = [...msg, ...zeros(N - K)];
    const [d, r] = deconv(padded, gen);
    return r;
  };
  const rx = rsRx(N, K, msg, gen);
  console.log("R[x]", rx.map(e => e.toString(2).padStart(8, "0")));
  const code = [...msg, ...rx.map(c => neg(c))];
  // NOTE: code is multiple of G(x): code(x) become 0 when applied x=a^(b+0) or ... or a^(b^+N-K-1)

  // 2.1. Prepare data for decoding
  // For decoding, code array as normal polynomial order: [1,x,x^2,...,x^{N-1}]
  const coder0 = code.slice().reverse(); 

  // Inject example errors (upto t=3 errors)
  const coder = coder0.slice();
  coder[2] = add(coder[2], 0xff);
  //coder[10] = add(coder[10], 0x02);
  coder[15] = add(coder[15], 0xac);
  console.log("No-Error code", coder0);
  console.log("Error injected", coder);
  console.log("Error values", coder.map((c, i) => sub(coder0[i], c)));
  
  // 3. RS Decode
  // 3.1. create Syndrome
  const rsSyndromes = (N, K, b, coder) => range(N - K).map(
    k => adds(coder.map((e, i) => mul(e, pow(a, i * (b + k))))));
  // syndromes are all 0 when no error
  const syndromes = rsSyndromes(N, K, b, coder);
  console.log("S[x]", syndromes);
  //console.log(range(N).map(k => adds(coder0.map((e, i) => mul(e, pow(a, i * (b + k)))))));
  
  // 3.2. Find error positions
  const rsEquations = (t, syndromes) => range(t).map(i => range(t + 1).map(j => syndromes[i + j]));
  const sequations = rsEquations(t, syndromes);
  console.log("Equations from S[x]", sequations);
  
  const pn1 = 2 ** n - 1;
  const modpn1 = v => (pn1 + v % pn1) % pn1;
  const powList = [...gf2PowList(n)(f), zero];
  const revList = powList.map((_, i) => powList.indexOf(i));
  const inv = e => powList[modpn1(-revList[e])];

  // solve equations on GF(2^n)-coeff equations
  const pivoting = (m, i) => {
    for (let j = i + 1; j < m.length; j++) {
      if (!isZero(m[j][i])) {
        [m[i], m[j]] = [m[j], m[i]];
        break;
      }
    }
  };
  const solve = m => {
    const r = m.map(l => l.slice());
    const n = r[0].length - 1;
    for (let i = 0; i < n; i++) for (let j = i + 1; j < r.length; j++) {
      if (isZero(r[i][i])) pivoting(r, i);
      const f = mul(r[j][i], inv(r[i][i]));
      for (let k = i; k < n + 1; k++) r[j][k] = sub(r[j][k], mul(f, r[i][k]));
    }
    for (let i = n - 1; i >= 0; i--) {
      r[i][n] = mul(r[i][n], inv(r[i][i]));
      r[i][i] = one;
      for (let j = i - 1; j >= 0; j--) {
        r[j][n] = sub(r[j][n], mul(r[j][i], r[i][n]));
        r[j][i] = zero;
      }
    }
    return r.slice(0, n).map(l => neg(l[n]));
  };
  const rsLx = sequations => {
    // error count is non-zero element of ls
    const ls = solve(sequations);
    return [one, ...ls.slice().reverse()];
  };
  const lx = rsLx(sequations);
  console.log("L[x]", lx); // [1,0,0,a] if 1-error, [1,0,b,a] if 2-error, [1,c,b,a] if 3-error

  const apply = (f, x) => adds(f.map((c, i) => mul(c, pow(x, i))));
  //console.log(range(N).map(k => apply(lx, inv(pow(a, k)))));
  const rsPositions = (N, lx) => range(N).filter(k => isZero(apply(lx, pow(a, modpn1(-k)))));
  const positions = rsPositions(N, lx);
  console.log("Error positions", positions); // same as injected error positions

  // 3.3. Error values for each error positions
  const conv = (a, b) => {
    const an = a.length, bn = b.length;
    const ps = a.map((ai, i) => [...zeros(i), ...b.map(bj => mul(ai, bj)), ...zeros(an - i)]);
    return ps.reduce((r, p) => r.map((e, i) => add(e, p[i])), zeros(an + bn));
  };
  //const rsOx = (t, syndromes, lx) => conv(syndromes, lx).slice(0, 2 * t);
  //const rsOx = (t, syndromes, lx) => conv(syndromes.slice(0, 2 * t), lx);
  const rsOx = (t, syndromes, lx) => conv(syndromes.slice(0, 2 * t), lx).slice(0, 2 * t);
  const ox = rsOx(t, syndromes, lx);
  console.log("O[x]", ox);

  const rsDLx = lx => lx.map((c, i) => times(c, i)).slice(1);
  const dlx = rsDLx(lx);
  console.log("L'[x]", dlx);

  const rsErrors = (b, positions, ox, dlx) => positions.map(k => {
    const akinv = pow(a, modpn1(-k));
    const oakinv = apply(ox, akinv), dlakinv = apply(dlx, akinv);
    const akb = pow(a, k * modpn1(1 - b)); 
    return times(mul(akb, mul(neg(oakinv), inv(dlakinv))),  -1);
  });

  const errors = rsErrors(b, positions, ox, dlx);
  console.log("Error values", errors);

  // 3.4 Apply error positions and these values into code
  const rsFixed = (coder, positions, errors) => coder.map((c, i) => {
    const idx = positions.indexOf(i);
    return idx < 0 ? c : sub(c, errors[idx]);
  });
  const fixed = rsFixed(coder, positions, errors);
  console.log("Fixed code", fixed);
  console.log("Same as original", fixed.every((c, i) => c === coder0[i]));
}

rs-code-note.txt

# Decoding Reed Solomon code

---

## Example properties of Reed Solomon

- N: size of code array, K: size of message array
- m = N - K = 7 : size of parity data array
- t = 7 // 2 = 3 : max errors in code array
- b = 0 : base order of generator polynomial as Pi[i=0..<m]{x - a^(b+i)}
- GN(2^8): Galois-Field of polynomial coefficients (as GF(2^8) array)

---

## Generator polynomial G(x)

G(x) = (x-1)(x-a)(x-a^2)(x-a^3)...(x-a^6)
     = x^6-(1+a+a^2+...a^6)x^5+...+a^21

- G(1)=G(a)=G(a^2)=...G(a^6)=0
- G(a^7)!=0, ..., G(a^254)!=0

---

## Summary of RS Encode

- Message data M[] = [M[0], M[1], ..., M[K-1]] : M{i] as polynomial emenet of GN(2^8) 
- Message polynomial M(x) = M[0]x^(k-1) + M[1]x^(k-2) + ... + M[K-1]
- Parity polynomial P(x) = M(x)*x^m mod G(x) as P[0]x^(m-1) + P[1]x^(m-2) + ... + P[m-1]
- Coded polynomial C(x) = M[0]x^(N-1) + ... + N[K-1]x^m - P[0]x^(m-1) - ... - P[m-1] 

Coded array C[] = [M[0], ..., N[K-1], -P[0], ..., -P[m-1]]

----

## Error injected message

- Valid message V[] = [C{N-1], C[N-2], ..., C[0]] 
- 3 Errors E[] = [0, ..., 0, e0, 0, ..., 0, e1, 0, ..., 0, e2, 0, ...0]
- Recieved data R[] = V[] + E[]
    - R[i0] = V[i0]+e0, R[i1] = V[i1] + e1, R[i2] = V[i2] + e2
    - other R[i] = V[i]

---

## Compute Syndromes

- s0 = R[0]+R[1]+...+R[N-1]
- s1 = R[0]a+R[1]a^2+...+R[N-1]a^(N-1)
- s2 = R[0]a^2+R[1]a^4+...+R[N-1]a^2(N-1)
- ...
- s6 = R[0]a^6+R[1]a^12+...+R[N-1]a^6(N-1)

if no error, s0 = s1 = ... = s6 = 0

- sk = V[0]a^k+...V[N-1]a^k(N-1) = 0

if error, s0 != 0, s1 != 0, ..., s6 != 0

- R[i] = V[i]+E[i]

sk = (V[0]+E[0])a^k + ... + (V[N-1]+E[N-1])a^k(N-1)
   = V[0]a^k+...V[N-1]a^k(N-1) + E[0]a^k+...+E[N-1]a^k(N-1)
   = E[0]a^k+...+E[N-1]a^k(N-1)
   = e0*a^(k*i0) + e1*a^(k*i1) + e2*a^(k*i2)

- s0 = e0+e1+e2
- s1 = e0a^i0 + e1a^i1 + e2a^i2
- s2 = e0a^2i0 + e1a^2i1 + e2a^2i2
- s3 = e0a^3i0 + e1a^3i1 + e2a^3i2
- s4 = ...
- s5 = ...
- s6 = ...

---

## Error position polynomial

l(x) = (1-a^i0*x)(1-a^i1*x)(1-a^i2*x)

- l(a^-i0) = l(a^-i1) = l(a^-i2) = 0


l(x) = 1 - (a^i0+a^i1+a^i2)x + (a^(i0+i1)+a^(i0+i2)+a^(i1+i2))x^2+a^(i0+i1+i2)x^3
     = l0 + l1*x + l2*x^2 + l3*x^3 

- l0 = 1
- l1 = -(a^i0+a^i1+a^i2)
- l2 = a^(i0+i1)+a^(i0+i2)+a^(i1+i2)
- l3 = -a^(i0+i1+i2)

----

## Equations of l1,l2,l3 and syndromes

- s0*l3 + s1*l2 + s2*l1 + s3 = 0
- s1*l3 + s2*l2 + s3*l1 + s4 = 0
- s2*l3 + s3*l2 + s4*l1 + s5 = 0

Find l1, l2, l3 from solving above equations

On each i=0..<N, search i=i0,i1,i2 with l(a^-i) = 0

----

## Forcus on e0 terms of the equations

s0*l3 = -e0*a^(i0+i1+i2) + ...

s1*l2 = e0*{a^i0(a^(i0+i1) + a^(i0+i2) + a^(i1+i2))} + ...
      = e0*{a^(2i0+i1) + a^(2i0+i2) + (i0+i1+i2)} + ...

s2*l1 = e0*a^2i0{a^i0 + a^i1 + a^i2} + ...
      = -e0*{a^3i0 + a^(2i0+i1) + a^(2i0+i2)} + ...

s3 = e0*a^3i0 + ...

- As a result:  s0*l3 + s1*l2 + s2*l3 + s3 = 0

----

## Error value polynomials

s(x) = s0 + s1*x + s2*x^2 + s3*x^3 + s4*x^4 + s5*x^5

o(x) = s(x)*l(x) mod x^6

dl(x) = d/dx{l(x)} = d/dx{(1-a^i0*x)(1-a^i1*x)(1-a^i2*x)}
      = -a^i0(1-a^i1*x)(1-a^i2*x) - a^i1(1-a^i0*x)(1-a^i2*x) - a^i2(1-a^i0*x)(1-a^i1*x)


---

s(x) = e0(1+a^i0*x+(a^i0*x)^2+...+(a^i0*x)^5) + e1(...) + e2(...)

- 1+b+b^2+..+b^5 = (b^6-1)/(b-1)

s(x) = e0((a^i0*x)^6-1)/(a^i0*x-1) + e1((a^i1*x)^6-1)/(a^i1*x-1) + e2((a^i2*x)^6-1)/(a^i2*x-1)

s(x)l(x) = -e0((a^i0*x)^6-1)(1-a^i1*x)(1-a^i2*x) - e1((a^i1*x)^6-1)(1-a^i0*x)(1-a^i2*x) - e2((a^i2*x)^6-1)(1-a^i0*x)(1-a^i1*x)

- ((a^i0*x)^6-1) mod x^6 = -1

o(x) = s(x)l(x) mod x^6 = e0(1-a^i1*x)(1-a^i2*x) + e1(1-a^i0*x)(1-a^i2*x) + e2(1-a^i0*x)(1-a^i1*x)

---

[apply a^-i0]

o(a^-i0) = e0(1-a^(i1-i0))(1-a^(i2-i0)
dl(a^-i0) = -a^i0(1-a^(i1-i0))(1-a^(i2-i0)

o(a^-i0)/dl(a^-i0) = -e0a^-i0

e0 = -a^i0*o(a^-i0)/dl(a^-i0)

=> ek = -a^ik*o(a^-ik)/dl(a^-ik)


---

## Recover valid data

- V[] = R[] - E[]

Valid message as M[] = [V[N-1], ..., V[m-1]]

rs-code.js

// functions from gf2n.js
const range = n => [...Array(n).keys()];
const uniq = a => a.filter((e1, i) => a.findIndex(e2 => e1 === e2) === i);
const msb32 = n => 31 - Math.clz32(n);

const polynom2s = k => range(2 ** (k + 1));
const poly2Add = (a, b) => a ^ b;
const poly2Times = (a, m) => m % 2 === 0 ? 0 : a;

const poly2Mul = (a, b) => {
  let r = 0;
  for (let i = 0; b > 0; i++, b >>>= 1) {
    if (b & 1) r ^= a << i;
  }
  return r;
};

const gf2 = n => polynom2s(n - 1);
const gf2Add = n => (a, b) => poly2Add(a, b);
const gf2Times = n => (a, m) => poly2Times(a, m);
const gf2Carry = n => f => e => {
  const overflow = 1 << msb32(f);
  e <<= 1;
  return (e & overflow) ? (e ^ f) : e;
};

// mul: case of enough bit (numbit >= 2 * n)
const poly2Mod = n => f => e => {
  const mf = msb32(f);
  for (let i = msb32(e); i >= n; i--) {
    if (e & (1 << i)) e ^= f << (i - mf);
  }
  return e;
};
//console.log(polyMod(3, 0b1011)(0b10000).toString(2));
const gf2Carrys0 = n => f => (e, k) => poly2Mod(n)(f)(e << k);
const gf2Mul0 = n => f => (a, b) => poly2Mod(n)(f)(poly2Mul(a, b));

// mul: case of ferwer bit
const gf2Carrys = n => f => (e, k) => {
  if (msb32(e) + k <= 32) return gf2Carrys0(n)(f)(e, k);
  for (let i = 0; i < k; i++) e = gf2Carry(n)(f)(e);
  return e;
};
const gf2Mul = n => f => (a, b) => {
  if (msb32(a) + msb32(b) <= 32) return gf2Mul0(n)(f)(a, b);
  let r = 0;
  for (let i = 0; b > 0; i++, b >>>= 1) {
    if (b & 1) r ^= gf2Carrys(n)(f)(a, i);
  }
  return r;
};

// pow
const gf2Pow = n => f => (e, k) => {
  const mul = gf2Mul(n)(f);
  let r = 1;
  for (let i = 0; i < k; i++) r = mul(r, e);
  return r;
};

// f(a)=0 is primitive element of GF(2^n) or not
const gf2PowList = n => f => {
  const carry = gf2Carry(n)(f);
  const r = [];
  for (let k = 0, m = 2 ** n - 1, c = 1; k < m; k++, c = carry(c)) r.push(c);
  return r;
};
// end from gf2n.js

// Utilities for RS code
const combinations = (a, n) => {
  if (n === 0) return [[]];
  return range(a.length - n + 1).flatMap(
    k => combinations(a.slice(k + 1), n - 1).map(s => [a[k], ...s]));
};
const zeros = n => Array(n).fill(0);
const sum = a => a.reduce((r, e) => r + e, 0); // sum for usual numbers

// Reed-Solomon code 
const rscode = ({N, K, b = 0, n = 8, f = 0b100011101} = {}) => {
  console.assert(K > 0 && N > K);
  const t = (N - K) >>> 1;
  const zero = 0, one = 1, a = 2, isZero = e => e === zero;
  const add = gf2Add(n), times = gf2Times(n), mul = gf2Mul(n)(f), pow = gf2Pow(n)(f);

  const neg = e => times(e, -1);
  const sub = (a, b) => add(a, neg(b));
  const adds = (es) => es.reduce((r, e) => add(r, e), zero);
  const muls = (es) => es.reduce((r, e) => mul(r, e), one);

  // generator
  const rsGen = (N, K, b) => {
    const ss = range(N - K).map(v => v + b);
    return range(N - K + 1).map(k => {
      const c = adds(combinations(ss, k).map(g => pow(a, sum(g))));
      return times(c, (-1) ** (N - K - k));
    });
  };
  const gen = rsGen(N, K, b);

  // encode
  // deconvolution a(x),g(x) to d(x),r(x) as a(x) = d(x)*g(x)+r(x)
  const deconv = (a, g) => {
    if (a.length < g.length) return [[], a];
    const k = a.length - g.length + 1;
    const d = [], r = a.slice();
    for (let i = 0; i < k; i++) {
      const c = r[i];
      d.push(c);
      for (let j = 0; j < g.length; j++) {
        r[i + j] = sub(r[i + j], mul(g[j], c));
      }
    }
    return [d, r.slice(k)];
  };
  // compute parity part
  const rsRx = (N, K, msg, gen) => {
    const padded = [...msg, ...zeros(N - K)];
    const [d, r] = deconv(padded, gen);
    return r;
  };
  const rsEncode = msg => {
    console.assert(msg.length === K, `msg.length should be ${K}`);
    const rx = rsRx(N, K, msg, gen);
    return [...msg, ...rx.map(c => neg(c))];
  };
  
  // decode
  // syndromes
  const rsSyndromes = (N, K, b, coder) => range(N - K).map(
    k => adds(coder.map((e, i) => mul(e, pow(a, i * (b + k))))));

  // error positions
  // gaussion elimination solver on GF(2^n) coefficient equations
  const pn1 = 2 ** n - 1;
  const modpn1 = v => (pn1 + v % pn1) % pn1;
  const powList = [...gf2PowList(n)(f), zero];
  const revList = powList.map((_, i) => powList.indexOf(i));
  const inv = e => powList[modpn1(-revList[e])];
  const pivoting = (m, i) => {
    for (let j = i + 1; j < m.length; j++) {
      if (!isZero(m[j][i])) {
        [m[i], m[j]] = [m[j], m[i]];
        break;
      }
    }
  };
  const solve = m => {
    const r = m.map(l => l.slice());
    const n = r[0].length - 1;
    for (let i = 0; i < n; i++) for (let j = i + 1; j < r.length; j++) {
      if (isZero(r[i][i])) pivoting(r, i);
      const f = mul(r[j][i], inv(r[i][i]));
      for (let k = i; k < n + 1; k++) r[j][k] = sub(r[j][k], mul(f, r[i][k]));
    }
    for (let i = n - 1; i >= 0; i--) {
      r[i][n] = mul(r[i][n], inv(r[i][i]));
      r[i][i] = one;
      for (let j = i - 1; j >= 0; j--) {
        r[j][n] = sub(r[j][n], mul(r[j][i], r[i][n]));
        r[j][i] = zero;
      }
    }
    return r.slice(0, n).map(l => neg(l[n]));
  };

  // error polision polynomial l(x)
  const rsEquations = (t, syndromes) => range(t).map(i => range(t + 1).map(j => syndromes[i + j]));
  const rsLx = sequations => {
    // error count is non-zero element of ls
    const ls = solve(sequations);
    return [one, ...ls.slice().reverse()];
  };
  const apply = (f, x) => adds(f.map((c, i) => mul(c, pow(x, i))));
  const rsPositions = (N, lx) => range(N).filter(k => isZero(apply(lx, pow(a, modpn1(-k)))));

  // convolution on GF(2^n) as polynomial multiplication
  const conv = (a, b) => {
    const an = a.length, bn = b.length;
    const ps = a.map((ai, i) => [...zeros(i), ...b.map(bj => mul(ai, bj)), ...zeros(an - i)]);
    return ps.reduce((r, p) => r.map((e, i) => add(e, p[i])), zeros(an + bn));
  };
  // polynomials for error values
  const rsOx = (t, syndromes, lx) => conv(syndromes.slice(0, 2 * t), lx).slice(0, 2 * t);
  const rsDLx = lx => lx.map((c, i) => times(c, i)).slice(1);
  const rsErrors = (b, positions, ox, dlx) => positions.map(k => {
    const akinv = pow(a, modpn1(-k));
    const oakinv = apply(ox, akinv), dlakinv = apply(dlx, akinv);
    const akb = pow(a, k * modpn1(1 - b)); 
    return times(mul(akb, mul(neg(oakinv), inv(dlakinv))), -1);
  });

  // apply errors
  const rsFixed = (coder, positions, errors) => coder.map((c, i) => {
    const idx = positions.indexOf(i);
    return idx < 0 ? c : sub(c, errors[idx]);
  });

  const rsDecode = code => {
    console.assert(code.length === N, `ode.length should be ${N}`);
    const coder = code.slice().reverse();
    const syndromes = rsSyndromes(N, K, b, coder);
    const sequations = rsEquations(t, syndromes);
    const lx = rsLx(sequations);
    const positions = rsPositions(N, lx);
    const ox = rsOx(t, syndromes, lx);
    const dlx = rsDLx(lx);
    const errors = rsErrors(b, positions, ox, dlx);
    const fixed = rsFixed(coder, positions, errors);
    return fixed.reverse().slice(0, K);
  };

  return {encode: rsEncode, decode: rsDecode};
};
 

// Example
{
  console.log("[Reed-Solomon error detective code]");
  // GF(2^8) System for QR Code
  const n = 8, f = 0b100011101; //a^8+a^4+a^3+a^2+1=0
  // example from https://kouyama.sci.u-toyama.ac.jp/main/education/2007/infomath/pdf/text/text09.pdf
  const N = 26, K = 19, b = 0;
  const {encode, decode} = rscode({N, K, b, n, f});
  const msg = [
    0b10000000, 0b01000100, 0b10000101, 0b10100111, 0b01001001, 0b10100111,
    0b10001011, 0b01101100, 0b00000000, 0b11101100, 0b00010001, 0b11101100,
    0b00010001, 0b11101100, 0b00010001, 0b11101100, 0b00010001, 0b11101100,
    0b00010001,
  ];

  // encode
  const code = encode(msg);
  // error injection
  code[24] ^= 0xff;
  //code[16] ^= 0x02;
  code[11] ^= 0xac;

  // decode with error fixed
  const fixed = decode(code);
  //console.log(fixed.map(e => e.toString(2).padStart(8, "0")));
  console.log("Same as original", fixed.every((c, i) => c === msg[i]));
}

refined as https://gist.github.com/bellbind/ba4b307fea8251f7f40aeb7e9349814b