[JavaScript] Programming Polynomial and Galois Field and Implementing Reed-Solomon code

bchcode-example.js

import {PF} from "./pf.js";
import {GF} from "./gf.js";
import {GF2n} from "./gf2n.js";
import {BCHCode} from "./bchcode.js";
import {Polynomial, PolynomialUtils} from "./polynomial.js";

{
  console.log("[BCHCode with GF]");

  const gf = GF(PF(2), 4, [1, 1, 0, 0, 1]);
  const N = 15, t = 2, c = 1;
  const bch = BCHCode(gf, N, t, c);
  const msg = [1,0,1,0,0,0,0];
  console.log("msg", msg);
  let coded = bch.encode(msg);
  console.log("coded", coded);
  coded[14] = +!coded[14];
  coded[12] = +!coded[12];
  console.log("error coded", coded);
  
  const fixed = bch.decode(coded);
  console.log("fixed", fixed);
  console.log();
}

{
  console.log("[BCHCode with GF2n]");
  const gf = GF2n(4, 0b10011);
  const N = 15, t = 2, c = 1;
  const bch = BCHCode(gf, N, t, c);
  console.log("gen", bch.gen.toString(2));
  console.log("K", bch.K);
  
  const msg = 0b0000101;
  console.log("msg", msg.toString(2).padStart(bch.K, "0"));
  let coded = bch.encode(msg);
  console.log("coded", coded.toString(2).padStart(N, "0"));
  coded = coded ^ (1 << 14);
  coded = coded ^ (1 << 12);
  console.log("error coded", coded.toString(2).padStart(N, "0"));
  
  const fixed = bch.decode(coded);
  console.log("fixed", fixed.toString(2).padStart(bch.K, "0"));
  console.log();
}

{
  console.log("[BCHCode for format info on QR code]");
  // https://www.thonky.com/qr-code-tutorial/format-version-information
  const gf = GF2n(4, 0b10011);
  const N = 15, t = 3, c = 1;
  const bch = BCHCode(gf, N, t, c);
  console.log("gen", bch.gen.toString(2)); // 0b10100110111
  console.log("K", bch.K); // 5

  const msg = 0b001100; //erro L & mask 4
  let coded = bch.encode(msg); // then masked XOR with 0b101010000010010 on QR
  console.log("coded", coded.toString(2).padStart(N, "0"));
  coded = coded ^ (1 << 14);
  coded = coded ^ (1 << 12);
  coded = coded ^ (1 << 5);
  console.log("error coded", coded.toString(2).padStart(N, "0"));

  const fixed = bch.decode(coded);
  console.log("fixed", fixed.toString(2).padStart(bch.K, "0"));
  console.log();
}

bchcode.js

import {range, uniq, FieldUtils} from "./field-utils.js";
import {Polynomial, PolynomialUtils} from "./polynomial.js";
import {GFUtils} from "./gf.js";

export const BCHCode = (gf, N, t, c) => {
  const gfutils = GFUtils(gf), gfPoly = Polynomial(gf);
  const futils = FieldUtils(gf);
  const pfPoly = gf.poly;

  const pfPoly2gfPoly = (pfp, n) =>
        gfPoly.fromArray(pfPoly.toArray(pfp, n).map(gf.fromSF));
  const gfPoly2pfPoly = (gfp, n) =>
        pfPoly.fromArray(gfPoly.toArray(gfp, n).map(gf.toSF));
  
  // encode
  const bchGenerator = () => {
    const polys = uniq(range(2 * t).map(
      k => gfutils.findIrreducible(gf.alpha(c + k))), pfPoly.eql);
    return pfPoly.prod(polys);
  };
  const gen = bchGenerator();
  const d = pfPoly.order(gen), K = N - d;
  
  const encode = msg => {
    const m = pfPoly.carry(msg, d);
    const r = pfPoly.mod(m, gen);
    return pfPoly.sub(m, r);
  };

  // decode
  const bchSyndromes = rx =>
        range(2 * t).map(k => gfPoly.apply(rx, gf.alpha(c + k)));
  const hasError = synds => !synds.every(gf.isZero);
  
  const bchLx = synds => {
    const leqs0 = range(t).map(i => range(t + 1).map(j => synds[i + j]));
    const v = futils.rank(leqs0);
    const leqs = v === t ? leqs0 :
          range(v).map(i => range(v + 1).map(j => synds[i + j]));
    const rlx = futils.solve(leqs);
    return [gf.one(), ...rlx.reverse()];
  };
  const bchPositions = lx => range(N).filter(
    k => gf.isZero(gfPoly.apply(lx, gf.alpha(-k))));
  
  const bchOx = (sx, lx) =>
        gfPoly.mod(gfPoly.mul(sx, lx), gfPoly.carry(gfPoly.one(), t));
  const bchDLx = lx => gfPoly.diff(lx);

  const bchErrors = (positions, ox, dlx) => positions.map(k => {
    const akinv = gf.alpha(-k);
    const oAkinv = gfPoly.apply(ox, akinv);
    const dlAkinv = gfPoly.apply(dlx, akinv);
    const ak1b = gf.alpha(k * (1 - c));
    return gf.neg(gf.mul(ak1b, gf.div(oAkinv, dlAkinv)));
  });
  
  const bchErrorsPF2 = positions => gfPoly.sum(positions.map(
    k => gfPoly.carry(gfPoly.one(), k)));
  
  const bchEx = (positions, errors) =>
        gfPoly.sum(positions.map((k, i) => gfPoly.monomial(errors[i], k)));
  
  const decode = code => {
    const rx = pfPoly2gfPoly(code, N);
    const synds = bchSyndromes(rx);
    if (!hasError(synds)) return gfPoly2pfPoly(rx.slice(d), K);
    
    //const lx = bchLx(synds);
    const lx = PolynomialUtils(gfPoly).findLinearRecurrence(synds);
    const positions = bchPositions(lx);
    if (positions.length === 0) {
      throw Error(`Cannot recover: error count > ${t}`);
    }

    if (pfPoly.K.size() === 2) {
      const ex = bchErrorsPF2(positions);
      const cx = gfPoly.sub(rx, ex);
      return gfPoly2pfPoly(cx.slice(d), K);
    } else {
      const sx = synds;
      const ox = bchOx(sx, lx);
      const dlx = bchDLx(lx);
      const errors = bchErrors(positions, ox, dlx);
      
      const ex = bchEx(positions, errors);
      const cx = gfPoly.sub(rx, ex);
      return gfPoly2pfPoly(cx.slice(d), K);
    }
  };
  
  return {K, d, gen, encode, decode};
};

berlekamp-massey-example.js

import {Q} from "./field-utils.js";
import {Polynomial, PolynomialUtils} from "./polynomial.js";

{
  const poly = Polynomial(Q());
  const polyutils = PolynomialUtils(poly);

  // fib sequence
  const sx = [1, 1, 2, 3, 5, 8].map(poly.K.fromNum);
  //console.log(sx);
  console.log(polyutils.findLinearRecurrence(sx));
  // => [1,-1,-1] as s[n] - s[n-1] - s[n-2] = 0 
}

field-utils.js

// Array utils
export const range = n => [...Array(n).keys()];
export const uniq = (a, eql = (e1, e2) => e1 === e2) =>
      a.filter((e1, i) => a.findIndex(e2 => eql(e1, e2)) === i);
export const transpose = a => range(a[0].length).map(k => a.map(e => e[k]));

// Field suppliment: generate isZero, sub, div, sum, prod
export const Field = K => {
  const {eql, zero, one, add, mul, neg, inv} = K;
  const sub = (a, b) => add(a, neg(b));
  const div = (a, b) => mul(a, inv(b));
  const sum = es => es.reduce((r, e) => add(r, e), zero());
  const prod = es => es.reduce((r, e) => mul(r, e), one());
  const isZero = e => eql(e, zero());
  return {sub, div, sum, prod, isZero, ...K};
};

// Example Field Q: Rational number 
export const Q = () => {
  const gcd = (a, b) => a === 0 ? b : gcd(b % a, a);
  const normalize = ([e0, e1]) => {
    const f = gcd(e1, e0);
    return [e0 / f, e1 / f];
  };

  const eql = (a, b) => a[0] * b[1] === a[1] * b[0];
  const zero = () => [0, 1];
  const one = () => [1, 1];
  const neg = e => [-e[0], e[1]];
  const inv = e => normalize([e[1], e[0]]);
  const add = (a, b) => normalize([a[0] * b[1] + b[0] * a[1], a[1] * b[1]]);
  const mul = (a, b) => normalize([a[0] * b[0], a[1] * b[1]]);

  const times = (e, k) => normalize([e[0] * k, e[1]]);
  const pow = (e, k) => normalize([e[0] ** k, e[1] ** k]);
  const toStr = e => e[1] === 1 ? `${e[0]}` : `${e[0]}/${e[1]}`;
  const toNum = e => Math.round(e[0] / e[1]);
  const fromNum = n => [n, 1];

  return Field({
    eql, zero, one, neg, inv, add, mul, times, pow, toStr, toNum, fromNum});
};

// Field utilities
export const FieldUtils = K => {
  const {isZero, zero, one, neg, inv, sub, mul, div} = K;
  // linear equations solver on the field f
  const pivoting = (r, u, i) => {
    for (let v = u + 1; v < r.length; v++) {
      if (!isZero(r[v][i])) {
        [r[u], r[v]] = [r[v], r[u]];
        break;
      }
    }
  };
  
  const upperTriangle = r => {
    const n = r[0].length - 1, m = r.length;
    for (let i = 0, u = 0; i < n && u < m; i++) {
      if (isZero(r[u][i])) pivoting(r, u, i);
      if (isZero(r[u][i])) continue;
      for (let v = u + 1; v < m; v++) {
        const c = div(r[v][i], r[u][i]);
        for (let j = i; j < n + 1; j++) {
          r[v][j] = sub(r[v][j], mul(c, r[u][j]));
        }
      }
      u++;
    }
  };

  const rank = lineqs => {
    const r = lineqs.map(lineq => lineq.slice());
    upperTriangle(r);
    return r.filter(eq => !eq.slice(0, -1).every(isZero)).length;    
  };
  
  const solve = lineqs => {
    const r = lineqs.map(lineq => lineq.slice());
    upperTriangle(r);
    const n = r[0].length - 1, m = r.length;
    for (let i = n - 1, u = m - 1; i >= 0 && u >= 0; i--) {
      while (u > 0 && isZero(r[u][i])) u--;
      if (isZero(r[u][i])) continue;
      r[u][n] = div(r[u][n], r[u][i]);
      r[u][i] = one();
      for (let v = u - 1; v >= 0; v--) {
        r[v][n] = sub(r[v][n], mul(r[v][i], r[u][n]));
        r[v][i] = zero();
      }
    }
    return r.slice(0, n).map(l => neg(l[n]));
  };
  return {rank, solve};
};

find-polynomial-example.js

import {PF} from "./pf.js";
import {GF, GFUtils} from "./gf.js";

const outputGrouping = (p, n, f) => {
  const gf = GF(PF(p), n, f), gfUtils = GFUtils(gf);
  console.log(`[GF(${p}^${n}) with ${gf.poly.toStr(f)}=0]`);
  const exposList = gfUtils.exponentGroups();
  const elemsList = exposList.map(expos => expos.map(expo => gf.alpha(expo)));
  for (const elems of elemsList) {
    const poly = gfUtils.findIrreducible(elems[0]);
    console.assert(elems.every(e => gf.eql(poly, gfUtils.findIrreducible(e))),
                   "outputGrouping");
    console.log(`Satisfied ${gf.poly.toStr(poly, "x")}=0`);
    console.log(elems.map(e => `- ${gf.toStr(e)}`).join("\n"));
    console.log();
  }
  console.log();
};

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

gf-calc-table.js

import {range} from "./field-utils.js";
import {PolynomialUtils} from "./polynomial.js";
import {PF} from "./pf.js";
import {GF, GFUtils} from "./gf.js";

// Calculation table of GF add/mul
const calcTable = (elems, op) => elems.map(e1 => elems.map(e2 => op(e1, e2)));
const toNumTable = f =>table => table.map(l => l.map(e => f.toNum(e)));

const permutateNumTable = (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 = f => (a, b) => a.every((al, i) => al.every(
  (av, j) => f.eql(av, b[i][j])));

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 = (f, table, mark) => {
  const elems = f.elems().map(e => f.toStr(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 => `$${f.toStr(e).padEnd(len)}$`));
  const heads = [mark, ...strList];
  const elements = strList.map((e, i) => [`**${e}**`, ...strTable[i]]);
  return mdTable(heads, elements);
};

export const outputGF = (p, n, f, GFctor = GF) => {
  const gf = GFctor(PF(p), n, f);
  const gfUtils = GFUtils(gf);
  const polyUtils = PolynomialUtils(gf.poly);
  
  console.assert(gf.isPrimitive(), "f is not a primitive polynomial");
  
  // generate GF(p^n) system
  const fs = polyUtils.irreducibles(n);
  const es = gf.elems();
  const addTable = calcTable(es, gf.add);
  const mulTable = calcTable(es, gf.mul);
  const pows = gf.pows();
  
  const auto = gfUtils.generatorAutomorphism();
  const autoNum = auto.map(e => gf.toNum(e));

  // check automorphism
  const addNum = toNumTable(gf)(addTable);
  const mulNum = toNumTable(gf)(mulTable);
  console.assert(eqlTable(mulNum, permutateNumTable(mulNum, autoNum)));
  console.assert(eqlTable(addNum, permutateNumTable(addNum, autoNum)));

 // output
  console.log();
  console.log(
    `### GF(${gf.size()} = ${p}^${n}) with $${gf.toStr(f, "a")}=0$`);
  console.log();
  console.log(`Elements of GF($${p}^${n}$)`, "\n");
  console.log("-", es.map(e => `$${gf.toStr(e)}$`).join(", "));
  console.log();
  console.log(`Power of $a$ mod $${gf.toStr(f)}$`);
  console.log();
  console.log(mdTable(["order", "number"], pows.map(
    (e, i) => [`$a^{${i}}$`, `$${gf.toStr(e)}$`])));
  console.log();
  console.log(`Calculation Tables of GF($${p}^${n}$)`);
  console.log();
  console.log(mdCalcTable(gf, addTable, "+"));
  console.log();
  console.log(mdCalcTable(gf, mulTable, "*"));
  console.log();

  console.log(`Generator of automorphism group of GF($${p}^${n}$)`);
  console.log();
  console.log(mdTable(["from", "to"], auto.map(
    (m, s) => [`$${gf.toStr(es[s])}$`, `$${gf.toStr(m)}$`])));
  console.log();

  console.log(`${n}-order irreducible polynomials`);
  console.log();
  for (const f of fs) {
    const gf = GFctor(PF(p), n, f);
    const pe = gf.isPrimitive() ? "(primitive element)" : "";
    console.log("-", `$${gf.poly.toStr(f)} = 0$`, pe);
  }
  console.log();
};

gf-table-example.js

import {outputGF} from "./gf-calc-table.js";

//import {PF} from "./pf.js";
//import {Polynomial, PolynomialUtils} from "./polynomial.js";
//import {GF} from "./gf.js";

// 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(PolynomialUtils(Polynomial(PF(p))).irreducibles(n).filter(f => GF(p, n, f).isPrimitive()));
  //outputGF(p, n, f);
}
{
  const p = 2, n = 6, f = [1, 1, 0, 0, 0, 0, 1];
  //console.log(PolynomialUtils(Polynomial(PF(p))).irreducibles(n).filter(f => GF(p, n, f).isPrimitive()));
  //outputGF(p, n, f);
}

gf.js

import {range, uniq, transpose, Field, FieldUtils} from "./field-utils.js";
import {PF} from "./pf.js";
import {Polynomial} from "./polynomial.js";

// GF(p^n)
export const GF = (K, n, f, name="a") => {
  const pn = K.size() ** n, pn1 = pn - 1;
  const poly = Polynomial(K);
  const p2gf = pe => poly.toArray(pe, n);
  const modpn1 = k => (pn1 + k % pn1) % pn1;

  const {eql, add, times, neg, toNum} = poly;
  const zero = () => p2gf(poly.zero());
  const one = () => p2gf(poly.one());
  const a = p2gf(poly.carry(poly.one(), 1));
  const isZero = e => eql(zero(), e);

  const mul0 = (a, b) => poly.mod(poly.mul(a, b), f);
  const pow0 = (e, k) => k === 0 ? one() : mul0(e, pow0(e, k - 1));
  
  const powList = Object.freeze(
    range(pn1).map(k => Object.freeze(pow0(a, k))));
  const exponent = e => powList.findIndex(pe => eql(e, pe));
  const mul = (a, b) => isZero(a) || isZero(b) ? zero() :
        powList[modpn1(exponent(a) + exponent(b))];
  const pow = (e, k) => isZero(e) ? zero() : k == 1 ? e :
        powList[modpn1(exponent(e) * k)];
  const alpha = (k = 1) => pow(a, k);
  const inv = e => powList[modpn1(-exponent(e))];
  const pows = () => powList;

  const size = () => pn;
  const toStr = e => poly.toStr(e, name);
  const fromNum = num => p2gf(poly.fromNum(num));
  const elems = () => poly.elems(n - 1);
  const isPrimitive = () => uniq(powList, eql).length === powList.length;
  const fromSF = e => p2gf(poly.monomial(e, 0));
  const toSF = e => e[0];
  
  return Field({
    K, poly, n, f, eql, zero, one, alpha, isZero,
    add, times, neg, exponent, mul, pow, inv,
    pows, toStr, size, toNum, fromNum, elems, isPrimitive, fromSF, toSF,
  });
};

// GF Utilities: judge primitive polynomial, irreducible formula of a GF elem 
export const GFUtils = gf => {
  const {poly, n, eql, pow, exponent, elems} = gf;
  const futils = FieldUtils(poly.K);
  const p = gf.K.size(), pn1 = gf.size() - 1;

  const cyclicOrder = e => {
    const k = exponent(e);
    for (let i = 1; i <= n; i++) {
      if (k * (p ** i) % pn1 === k) return i;
    }
    throw Error("never reached");
  };
  
  const findIrreducible = e => {
    const d = cyclicOrder(e);
    // gfeq: [1, e, e^2, ..., e^d] for c0 + c1e +...+ c(d-1)e^(d-1) + e^d = 0
    const gfeq = range(d + 1).map(k => poly.toArray(pow(e, k), n));
    // coefEqs: split gfeq with each coeddicients in e^0,e^1,...,e^d
    const coefEqs = transpose(gfeq);
    // modDim: [c0, c1, ...c(d-1)] satisfied above equation
    console.assert(futils.rank(coefEqs) === d, "findIrreducible");
    const cs = futils.solve(coefEqs);
    const modDim = poly.fromArray(cs);
    return poly.add(modDim, poly.carry(poly.one(), d));
  };

  const exponentGroups = () => {
    const used = new Set();
    const gs = [];
    for (const expo of range(pn1)) {
      if (used.has(expo)) continue;
      const g = uniq(range(n).map(k => (expo * (p ** k)) % pn1));
      g.forEach(expo => used.add(expo));
      gs.push(g);
    }
    return gs;
  };

  const generatorAutomorphism = () => elems().map(e => pow(e, p));
  
  return {findIrreducible, exponentGroups, generatorAutomorphism};
};

gf2n-table-example.js

import {outputGF} from "./gf-calc-table.js";
import {GF2n, PF2Polynomial} from "./gf2n.js";

//import {GFUtils} from "./gf.js";
//import {PolynomialUtils} from "./polynomial.js";

const outputGF2n = (n, f) => {
  outputGF(2, n, f, (p, n, f) => GF2n(n, f));
};

// examples
{
  const n = 2, f = 0b111;
  outputGF2n(n, f);
}
{
  const p = 2, n = 3, f = 0b1011;
  outputGF2n(n, f);
}
{
  const n = 4, f = 0b10011;
  outputGF2n(n, f);
}
{
  const n = 6, f = 0b1000011;
  //console.log(GF2n(n, f).isPrimitive());
  //console.log(PolynomialUtils(PF2Polynomial()).irreducibles(n).filter(f => GF2n(n, f)).map(e => e.toString(2)));
  //outputGF2n(n, f);
}

gf2n.js

import {Field, range, uniq} from "./field-utils.js";
import {PF} from "./pf.js";

const msb32 = n => 31 - Math.clz32(n);
const popcnt32 = n => {
  n = (n & 0x55555555) + ((n >>> 1) & 0x55555555);
  n = (n & 0x33333333) + ((n >>> 2) & 0x33333333);
  n = (n & 0x0f0f0f0f) + ((n >>> 4) & 0x0f0f0f0f);
  n = (n & 0x00ff00ff) + ((n >>> 8) & 0x00ff00ff);
  return (n & 0x0000ffff) + ((n >>> 16) & 0x0000ffff);
};

//GF(2) coefficient Polynomial as bits
export const PF2Polynomial = () => {
  const K = PF(2);
  const eql = (a, b) => a === b;
  const zero = () => 0;
  const one = () => 1;
  const add = (a, b) => a ^ b;
  const neg = e => e;
  const sub = (a, b) => add(a, b);
  const sum = es => es.reduce((r, e) => add(r, e), 0);
  const scale = (e, c) => (c & 1) ? e : 0;
  const times = (e, k) => (k & 1) ? e : 0;
  const carry = (e, k) => e << k;
  const mul = (a, b) => {
    let r = 0;
    for (; b > 0; b >>>= 1, a <<= 1) {
      if (b & 1) r ^= a;
    }
    return r;
  };
  const order = e => Math.max(msb32(e), 0);
  const mod = (a, b) => {
    const mb = msb32(b);
    for (let i = msb32(a); i >= mb; i--) {
      if (a & (1 << i)) a ^= b << (i - mb);
    }
    return a;
  };
  const prod = es => es.reduce((r, e) => mul(r, e), one());
  const diff = e => (e & 0xaaaaaaaa) >>> 1;
  const coef = (e, k) => e & (1 << k);
  const monomial = (c, k) => carry(scale(one(), c), k);
  const apply = (e, v) => (v & 1) ? popcnt32(e) & 1 : 0;
  const toStr = (e, name = "x") => {
    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 `${name}${pow}`;
    }).filter(e => e);
    const s2 = s1.length > 1 && s1[0] === "0" ? s1.slice(1) : s1;
    return s2.reverse().join("+");
  };
  const toArray = (e, len) => {
    const r = Array(len).fill(0);
    for (let i = 0; i < len; i++) r[i] = (e >>> i) & 1;
    return r;
  };
  const fromArray = bs => {
    let e = 0;
    for (let i = 0; i < bs.length; i++) e |= (bs[i] & 1) << i;
    return e;
  };
  const toNum = e => e;
  const fromNum = e => e;
  const elems = k => range(2 ** (k + 1));
  return {
    K, eql, zero, one, add, scale, neg, sub, sum, times, carry, mul, mod, prod,
    order, diff, coef, monomial, apply,
    toStr, toArray, fromArray, toNum, fromNum, elems,
  };
};

// GF(2^n) as bits
export const GF2n = (n, f, name="a") => {
  const poly = PF2Polynomial();
  const {K, eql, zero, one, add, times, neg, toNum, fromNum} = poly;
  const pn = 2 ** n, pn1 = pn - 1;
  const modpn1 = n => (pn1 + n % pn1) % pn1;
  
  const isZero = e => e === 0;
  const mul0 = (a, b) => poly.mod(poly.mul(a, b), f);
  const pow0 = (e, k) => k === 0 ? one() : mul0(e, pow0(e, k - 1));
  
  const powList = Object.freeze(range(pn1).map(k => pow0(2, k)));
  const expoList = range(pn).map(k => k === 0 ? NaN : powList.indexOf(k));
  const exponent = e => expoList[e];
  const mul = (a, b) => isZero(a) || isZero(b) ? zero() :
        powList[modpn1(exponent(a) + exponent(b))];
  const pow = (e, k) => isZero(e) ? zero() : k === 1 ? e :
        powList[modpn1(exponent(e) * k)];
  const alpha = (k = 1) => pow(2, k);
  const inv = e => powList[modpn1(-exponent(e))];
  const pows = () => powList;
  
  const size = () => pn;
  const elems = () => poly.elems(n - 1);
  const isPrimitive = () => uniq(powList, eql).length === powList.length;
  const toStr = e => poly.toStr(e, name);
  const fromSF = e => e % 2;
  const toSF = e => e % 2;
  
  return Field({
    K, poly, n, f, eql, zero, one, alpha, add, times, neg,
    exponent, mul, pow, inv,
    pows, toStr, size, toNum, fromNum, elems, isPrimitive, fromSF, toSF,
  });
};

gf4-2-example.js

import {PF} from "./pf.js";
import {Polynomial, PolynomialUtils} from "./polynomial.js";
import {GF, GFUtils} from "./gf.js";
import {GF2n} from "./gf2n.js";

{
  // GF(16) as GF(4^2)
  console.log("[GF(4^2) with GF(PF(2), 2)]");
  const sf = GF(PF(2), 2, [1, 1, 1]);
  const gf = GF(sf, 2, [[0, 1], [1, 0], [1, 0]], "b");
  console.log("GF(4^2) elements:", gf.elems().map(gf.toStr), "\n");

  // pow list
  console.log("GF(4^2) power of b");
  for (let i = 0; i < gf.size() - 1; i++) {
    console.log(`- b^{${i}}: ${gf.toStr(gf.alpha(i))}`);
  }
  console.log();

  // mul
  const A = [[0,1],[0,1]], B = [[1,0],[1,0]];
  const C = gf.mul(A, B);
  console.log(`mul example: ${gf.toStr(A)} * ${gf.toStr(B)} = ${gf.toStr(C)}\n`);

  console.log("2-order Irreducible polynomials on GF(4) coef");
  const poly = Polynomial(sf);
  const polyutils = PolynomialUtils(poly);
  const irreducibles = polyutils.irreducibles(2);
  for (const irr of irreducibles) {
    console.log(`- ${poly.toStr(irr)} = 0`);
  }
  console.log();
  
  const gfutils = GFUtils(gf);
  console.log(`${gf.toStr(A)} is a solution of ${
    gf.poly.toStr(gfutils.findIrreducible(A))}`);
  console.log();
}

{
  // GF(16) as GF(4^2) with GF2n
  console.log("[GF(4^2) with GF2n(2)]");
  const sf = GF2n(2, 0b111);
  const gf = GF(sf, 2, [0b10, 0b01, 0b01], "b");
  console.log("GF(4^2) elements:", gf.elems().map(gf.toStr), "\n");

  // pow list
  console.log("GF(4^2) power of b");
  for (let i = 0; i < gf.size() - 1; i++) {
    console.log(`- b^{${i}}: ${gf.toStr(gf.alpha(i))}`);
  }
  console.log();

  // mul
  const A = [0b10, 0b10], B = [0b01, 0b01];
  const C = gf.mul(A, B);
  console.log(`mul example: ${gf.toStr(A)} * ${gf.toStr(B)} = ${gf.toStr(C)}\n`);

  console.log("2-order Irreducible polynomials on GF(4) coef");
  const poly = Polynomial(sf);
  const polyutils = PolynomialUtils(poly);
  const irreducibles = polyutils.irreducibles(2);
  for (const irr of irreducibles) {
    console.log(`- ${poly.toStr(irr)} = 0`);
  }
  console.log();

  const gfutils = GFUtils(gf);
  console.log(`${gf.toStr(A)} is a solution of ${
    gf.poly.toStr(gfutils.findIrreducible(A))}`);
  console.log();
}

package.json

{
  "type": "module", "engines": {"node": ">=13.2"},
  "scripts": {
    "e:gf-table": "node gf-table-example.js",
    "e:gf2n-table": "node gf2n-table-example.js",
    "e:find-polynomial": "node find-polynomial-example.js",
    "e:rscode": "node rscode-example.js",
    "examples": "npm run e:gf-table && npm run e:gf2n-table && npm run e:find-polynomial && npm run e:rscode"
  }
}

pf.js

import {range, uniq, Field} from "./field-utils.js";

// number utils
const egcd = (a, b) => {
  if (a === 0) return [b, 0, 1];    // b = 0*0 + b*1
  const r = b % a, q = (b - r) / a; // r = b - a*q
  const [g, s, t] = egcd(r, a);     // g = r*s + a*t;
  return [g, t - q * s, s];         // g = (b-a*q)*s + a*t = a*(t-q*s) + b*t
};
const primes = max => {
  const ps = [];
  let ns = range(max).slice(2);
  while (ns.length > 0) {
    const p = ns[0];
    ps.push(p);
    ns = ns.filter(n => n % p);
  }
  return ps;
};
const primeFactors = n => {
  const ps = [];
  for (const p of primes(n)) {
    if (n % p === 0) ps.push(p);
  }
  return ps;
};


// GF(p)
export const PF = (p, pe) => {
  const p1 = p - 1;
  const modp = n => (p + n % p) % p;
  const modp1 = n => (p1 + n % p1) % p1;

  const eql = (a, b) => modp(a) === modp(b);
  const zero = () => 0;
  const one = () => 1;
  const add = (a, b) => modp(a + b);
  const mul = (a, b) => modp(a * b);
  const neg = e => modp(-e);
  const inv = e => modp(egcd(e, p)[1]);
  
  const times = (e, k) => modp(e * k);
  const pow = (e, k) => {
    k = modp1(k);
    return k === 0 ? one() : k === 1 ? e : mul(e, pow(e, k - 1));
  };
  const toStr = e => `${e}`;
  const toNum = e => e;
  const fromNum = n => n;

  // find primitive element for GF(p)
  const isPrimitiveElement = (e, pfs = primeFactors(p1)) => {
    for (const pf of pfs) if (pow(e, p1 / pf) === 1) return false;
    return true;
  };
  const primitiveElement = p => {
    const pfs = primeFactors(p1);
    for (let e = 2; e < p1; e++) if (isPrimitiveElement(e, pfs)) return e;
    throw Error("never reached");
  };
  const a = pe ? pe : p === 2 ? 1 : p === 3 ? 2 : primitiveElement(p);
  const alpha = (k = 1) => pow(a, k); 
  
  const size = () => p;
  const isPrimitive = () => isPrimitiveElement(a);
  const elems = () => range(p);
  
  return Field({
    p, eql, zero, alpha, one, add, mul, neg, inv, times, pow,
    toStr, toNum, fromNum, size, elems});
};

polynomial.js

import {range, uniq} from "./field-utils.js";

// Polynomial[F]
export const Polynomial = K => {
  const zeros = n => Array(n).fill(K.zero());

  const eql = (a, b) => {
    if (a.length <  b.length) [a, b] = [b, a];
    return a.every((c, k) => k < b.length ? K.eql(c, b[k]) : K.isZero(c));
  };
  const zero = () => [K.zero()];
  const one = () => [K.one()];
  const add = (a, b) => {
    if (a.length < b.length) [a, b] = [b, a];
    return a.map((c, k) => k < b.length ? K.add(c, b[k]) : c);
  };
  const neg = e => e.map(c => K.neg(c));
  const sub = (a, b) => add(a, neg(b));
  const sum = es => es.reduce((r, e) => add(r, e), zero());
  const scale = (e, c) => e.map(ci => K.mul(ci, c));
  const carry = (e, k) => [...zeros(k), ...e];
  const mul = (a, b) => sum(b.map((bk, k) => carry(scale(a, bk), k)));
  
  const order = e => {
    for (let i = e.length - 1; i > 0; i--) if (!K.isZero(e[i])) return i;
    return 0;
  };
  const mod = (a, b) => {
    const an = order(a), bn = order(b);
    console.assert(bn > 0, "poly.mod");
    if (an < bn) return a.slice(0, bn);
    const f = K.div(a[an], b[bn]);
    return mod(sub(a, carry(scale(b, f), an - bn)).slice(0, an), b);
  };
  const prod = es => es.reduce((r, e) => mul(r, e), one());
  const times = (e, k) => e.map(c => K.times(c, k));
  const diff = e => e.map((c, k) => K.times(c, k)).slice(1);

  const coef = (e, k) => e[k] || K.zero();
  const monomial = (c, k) => carry(scale(one(), c), k);
  const apply = (e, v) => K.sum(e.map((c, k) => K.mul(c, K.pow(v, k))));
  
  const toStr = (e, name = "x") => {
    const s1 = e.map((c, k) => {
      if (k === 0) {
        const coef = K.toStr(c);
        return coef.includes("+") ? `(${coef})` : coef;
      }
      if (K.isZero(c)) return "";
      const coef = K.eql(c, K.one()) ? "" : K.toStr(c);
      const cstr = coef.includes("+") ? `(${coef})` : coef;
      const pow = k === 1 ? "" : `^{${k}}`;
      return `${cstr}${name}${pow}`;
    }).filter(e => e);
    const s2 = s1.length > 1 && s1[0] === K.toStr(K.zero()) ? s1.slice(1) : s1;
    return s2.reverse().join("+");
  };
  const toArray = (e, len) =>
        e.length < len ? [...e, ...zeros(len - e.length)] : e.slice(0, len);
  const fromArray = cs => cs;
  const R = {
    K, eql, zero, one, add, scale, carry, neg, sub, sum, mul, mod, prod, times,
    order, diff, coef, monomial, apply, toStr, toArray, fromArray,
  };
  
  if (K.size && K.elems) {
    const toNum = e => e.reduceRight((r, c) => r * K.size() + K.toNum(c), 0);
    const fromNum = n => {
      const e = [];
      while (n > 0) {
        const r = n % K.size();
        e.push(K.fromNum(r));
        n = (n - r) / K.size();
      }
      return e;
    };
    const elems = k => {
      if (k < 0) return [[]];
      const base = elems(k - 1);
      return K.elems().flatMap(c => base.map(b => [...b, c]));
    };
    return {...R, toNum, fromNum, elems};
  } else return R;
};

// Polynomial utilities
export const PolynomialUtils = poly => {
  const {K, eql, one, add, sub, scale, carry, mul, coef, elems} = poly;
  // listing irreducible polynomials
  const kPolynoms = k => elems(k - 1).map(e => add(e, carry(one(), k)));
  const reducibles = n => {
    const rs = [];
    for (let i = 1, nh = n >> 1; i <= nh; i++) {
      const l = kPolynoms(i), r = kPolynoms(n - i);
      for (const lp of l) for (const rp of r) rs.push(mul(lp, rp));
    }
    return uniq(rs, eql);
  };
  const irreducibles = n => {
    const reds = reducibles(n);
    return kPolynoms(n).filter(f => !reds.find(r => eql(f, r)));
  };
  
  const findLinearRecurrence = s => {
    // Berlekamp-Massey algorithm
    let cx = one(), cl = 1, bx = one(), bl = 1, b = K.one(), m = 0;
    for (let i = 0; i < s.length; i++) {
      const d = K.sum(range(cl).map(k => K.mul(coef(cx, k), s[i - k])));
      m++;
      if (K.isZero(d)) continue;
      const tx = cx, tl = cl;
      cx = sub(cx, scale(carry(bx, m), K.div(d, b)));
      cl = Math.max(cl, bl + m);
      if (cl > tl) [bx, bl, b, m] = [tx, tl, d, 0];
    }
    return cx;
  };
  
  return {irreducibles, findLinearRecurrence};
};

rscode-example.js

import {PF} from "./pf.js";
import {GF} from "./gf.js";
import {GF2n} from "./gf2n.js";
import {RSCode} from "./rscode.js";

{
  console.log("[Reed-Solomon error detective code]: N=26,K=7,b=0 on GF(2^8)");
  // 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(GF2n(n, f), N, K, b);
  const msg = [
    0b10000000, 0b01000100, 0b10000101, 0b10100111, 0b01001001, 0b10100111,
    0b10001011, 0b01101100, 0b00000000, 0b11101100, 0b00010001, 0b11101100,
    0b00010001, 0b11101100, 0b00010001, 0b11101100, 0b00010001, 0b11101100,
    0b00010001,
  ];
  console.log("msg", msg);
  
  const coded = encode(msg);
  console.log("coded", coded);
  
  // error injection
  coded[24] ^= 0xff;
  //coded[16] ^= 0x02;
  coded[11] ^= 0xac;
  console.log("error coded", coded);

  const fixed = decode(coded);
  console.log("fixed", fixed);
  
  console.log("Same as original", fixed.every((c, i) => c === msg[i]));
}

{
  console.log("[Reed-Solomon error detective code]: N=7,K=k,b=1 on GF(929)");
  // Reed Solomon code on Prime field 929
  const p = 929;
  const N = 7, K = 3, b = 1, debug = false;
  const {encode, decode} = RSCode(PF(p), N, K, b, debug);

  const msg = [3, 2, 1];
  console.log("msg", msg);
  
  const coded = encode(msg);
  console.log("codec", coded);

  coded[0] = (coded[0] + 10) % p;
  coded[5] = (coded[5] + 20) % p;
  console.log("error coded", coded);
  
  const fixed = decode(coded);
  console.log("fixed", fixed);
  
  console.log("Same as original", fixed.every((c, i) => c === msg[i]));  
}

{
  console.log("[Reed-Solomon error detective code]: N=26,K=7,b=0 on GF(2^8)");
  // Slow GF version
  // GF(2^8) System for QR Code
  const p = 2, n = 8, f = [1, 0, 1, 1, 1, 0, 0, 0, 1]; //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(GF(PF(p), n, f), N, K, b);
  const msg = [
    0b10000000, 0b01000100, 0b10000101, 0b10100111, 0b01001001, 0b10100111,
    0b10001011, 0b01101100, 0b00000000, 0b11101100, 0b00010001, 0b11101100,
    0b00010001, 0b11101100, 0b00010001, 0b11101100, 0b00010001, 0b11101100,
    0b00010001,
  ];
  console.log("msg", msg);
  
  const coded = encode(msg);
  console.log("coded", coded);
  
  // error injection
  coded[24] ^= 0xff;
  //coded[16] ^= 0x02;
  coded[11] ^= 0xac;
  console.log("error coded", coded);

  const fixed = decode(coded);
  console.log("fixed", fixed);
  
  console.log("Same as original", fixed.every((c, i) => c === msg[i]));
}

rscode.js

// Reed-Solomon code
import {range, FieldUtils} from "./field-utils.js";
import {Polynomial, PolynomialUtils} from "./polynomial.js";

export const RSCode = (FF, N, K, b, debug = false) => {
  const d = N - K, t = d >>> 1;
  const poly = Polynomial(FF), futils = FieldUtils(FF);
  // explict conversion between FF-element array and FF-polynomial
  const a2p = poly.fromArray, rev2p = a => a2p([...a].reverse());
  const p2a = poly.toArray, p2rev = (e, len) => [...p2a(e, len)].reverse();

  // gen = (x-a^b)*(x-a^(b+1)*...*(x-a^(b+d-1)) as extracted poly
  const rsGenerator = () => poly.prod(range(d).map(k => poly.add(
    poly.monomial(FF.one(), 1),
    poly.monomial(FF.neg(FF.alpha(b + k)), 0))));
  const gen = rsGenerator();
  if (debug) console.log("gen", gen);
  
  const encode = msg => {
    const base = poly.carry(rev2p(msg.map(FF.fromNum)), d);
    const coded = poly.sub(base, poly.mod(base, gen));
    if (debug) console.assert(poly.mod(coded, gen).every(FF.isZero));
    return p2rev(coded, N).map(FF.toNum);
  };

  const rsSyndromes = cx => range(d).map(
    k => poly.apply(cx, FF.alpha(b + k)));
  const hasError = synds => !synds.every(si => FF.isZero(si));
  
  const rsLx = (synds) => {
    const leqs0 = range(t).map(i => range(t + 1).map(j => synds[i + j]));
    const v = futils.rank(leqs0);
    const leqs = v === t ? leqs0 :
          range(v).map(i => range(v + 1).map(j => synds[i + j]));
    const rlx = futils.solve(leqs);
    return a2p([FF.one(), ...rlx.reverse()]);
  };
  const rsPositions = lx => range(N).filter(
    k => FF.isZero(poly.apply(lx, FF.alpha(-k))));

  const rsOx = (sx, lx) =>
        poly.mod(poly.mul(sx, lx), poly.carry(poly.one(), d));
  const rsDLx = lx => poly.diff(lx);
  const rsErrors = (positions, ox, dlx) => positions.map(k => {
    const akinv = FF.alpha(-k);
    const oAkinv = poly.apply(ox, akinv);
    const dlAkinv = poly.apply(dlx, akinv);
    const ak1b = FF.alpha(k * (1 - b));
    return FF.neg(FF.mul(ak1b, FF.div(oAkinv, dlAkinv)));
  });

  const rsEx = (positions, errors) =>
        poly.sum(positions.map((k, i) => poly.monomial(errors[i], k)));
  
  const decode = code => {
    const rx = rev2p(code.map(FF.fromNum));
    const synds = rsSyndromes(rx);
    if (debug) console.log("sx", synds);
    if (!hasError(synds)) return code.slice(0, K);
    
    const lx = rsLx(synds);
    //const lx = PolynomialUtils(poly).findLinearRecurrence(synds);
    if (debug) console.log("lx", lx);
    const positions = rsPositions(lx);
    if (debug) console.log("positions", positions);
    if (positions.length === 0) {
      throw Error(`Cannot recover: error count > ${t}`);
    }
    
    const sx = a2p(synds);
    const ox = rsOx(sx, lx);
    const dlx = rsDLx(lx);
    const errors = rsErrors(positions, ox, dlx);
    if (debug) console.log("errors", errors);
    const ex = rsEx(positions, errors);
    const cx = poly.sub(rx, ex);
    return p2rev(cx, N).slice(0, K).map(FF.toNum);
  };
  return {encode, decode};
};