Galois fields of non-prime order in Javascript

galois_arithmetic.mjs

/**
 * A nonnegative integer.
 * @typedef {number} NaturalNumber
 */

/**
 * Datatype representing a polynomial as an array of coefficients.
 * e.g. "5X^3 + X + 2" is represented by [2, 1, , 5].
 * e.g. "0" is represented by [].
 * Must have no trailing zeroes. Array slots containing 0 are equivalent to empty (sparse) array slots.
 * @typedef {Array.<NaturalNumber>} Polynomial
 */


export class GaloisField {
	#ring;
	#irr;

	/**
	 * @param {NaturalNumber} p - The characteristic of the field.
	 * @param {NaturalNumber} [k=1] - The degree of the field.
	 * @param {Polynomial} [irr] - The irreducible polynomial used to construct the field, when k !== 1.
	 * @returns {GaloisField}
	 */
	constructor(p, k=1, irr=null) {
		if (typeof p === 'number' && p > Number.MAX_SAFE_INTEGER)
			throw new RangeError("characteristic too large for Number");
		else if (typeof p === 'bigint')
			throw new TypeError("Bignum support not implemented yet");

		if (k === 1 && irr === null)
			// Shim that will allow the k>1 codepaths to handle the k===1 case
			irr = [1];
		else if (k !== 1 && (irr === null || irr.length !== k+1))
			throw new Error("prime power fields need an irreducible polynomial of degree k");
		else if (k === 1 && irr !== null)
			throw new Error("prime order fields do not use an irreducible polynomial");

		this.#ring = new PolyRing(p);
		this.#irr = irr;
	}

	/**
	 * @param {NaturalNumber} x - An integer representative of a polynomial. Must be a natural number less than the order of the field.
	 * @returns {Polynomial}
	 */
	element(x) {
		if (x < 0 || x >= this.order)
			throw new RangeError();

		return this.#ring.element(x);
	}

	/**
	 * @param {Polynomial} a - An element of the field.
	 * @returns {NaturalNumber}
	 */
	representative(a) {
		return this.#ring.representative(a);
	}

	/**
	 * Test whether a Polynomial is an element of the field.
	 * Returns true if so;
	 * returns false if it is a well-formed Polynomial but not in the field;
	 * throws TypeError if it is not a well-formed polynomial.
	 * @param {Polynomial} a
	 * @returns {boolean}
	 */
	element_ok(a) {
		if (!this.#ring.element_ok(a))
			return false;

		if (a.length > this.degree)
			return false;

		return true;
	}

	/**
	 * Addition in the field.
	 * @param {Polynomial} a - Addend. Must be an element of the field.
	 * @param {Polynomial} b - Addend. Must be an element of the field.
	 * @returns {Polynomial}
	 */
	add(a, b) {
		// No need to reduce because addition cannot increase the degree
		return this.#ring.add(a, b);
	}

	/**
	 * Subtraction in the field.
	 * @param {Polynomial} a - Minuend. Must be an element of the field.
	 * @param {Polynomial} b - Subtrahend. Must be an element of the field.
	 * @returns {Polynomial}
	 */
	sub(a, b) {
		// No need to reduce because subtraction cannot increase the degree
		return this.#ring.sub(a, b);
	}

	/**
	 * Additive reciprocal in the field.
	 * @param {Polynomial} a - The polynomial to find the reciprocal of. Must be an element of the field.
	 * @returns {Polynomial}
	 */
	neg(a) {
		return this.#ring.neg(a);
	}

	/**
	 * Multiplication in the field.
	 * @param {Polynomial} a - Multiplicand. Must be an element of the field.
	 * @param {Polynomial} b - Multiplicand. Must be an element of the field.
	 * @returns {Polynomial}
	 */
	mul(a, b) {
		// Need to reduce because multiplication usually increases the degree
		return this.#reduce(this.#ring.mul(a, b));
	}

	/**
	 * Multiplication in the field.
	 * @param {Polynomial} a - Dividend. Must be an element of the field.
	 * @param {Polynomial} b - Divisor. Must be a nonzero element of the field.
	 * @returns {Polynomial}
	 */
	div(a, b) {
		// TODO: Do we really need to reduce here?
		return this.#reduce(this.#ring.mul(a, this.inv(b)));
	}

	/**
	 * Multiplicative reciprocal in the field.
	 * @param {Polynomial} a - The polynomial to find the reciprocal of. Must be a nonzero element of the field.
	 * @returns {Polynomial}
	 */
	inv(a) {
		if (a.length === 0)
			throw new RangeError("Polynomial division by zero");

		// TODO: Do we really need to reduce here?
		return this.#reduce(this.#ring.egcd(a, this.#irr)[1]);
	}

	/**
	 * @param {Polynomial} a - Polynomial. Must be an element of the underlying integer ring.
	 * @returns {Polynomial} - An element of the field.
	 */
	#reduce(a) {
		// Divide the value by our irreducible polynomial, and return the remainder.
		return this.#ring.divmod(a, this.#irr)[1];
	}

	/**
	 * @type {NaturalNumber}
	 */
	get characteristic() {
		return this.#ring.characteristic;
	}

	/**
	 * @type {NaturalNumber}
	 */
	get degree() {
		return this.#irr.length - 1;
	}

	/**
	 * @type {NaturalNumber}
	 */
	get order() {
		return this.characteristic ** this.degree;
	}
}


class PolyRing {
	#characteristic;

	/**
	 * @param {NaturalNumber} p - The order of the ring.
	 * @returns {PolyRing}
	 */
	constructor(p) {
		if (typeof p === 'bigint')
			throw new TypeError("Bignum support not implemented yet");

		this.#characteristic = p;
	}

	/**
	 * @param {NaturalNumber} x - An integer representative of a polynomial. Must be a natural number.
	 * @returns {Polynomial}
	 */
	element(x) {
		var result = Array.from(x.toString(this.#characteristic), (n) => +n);
		result.reverse();
		trimTrailingZeroes(result);
		return result;
	}

	/**
	 * @param {Polynomial} a - An element of the ring.
	 * @param {NaturalNumber}
	 */
	representative(a) {
		var result = 0,
		    p = this.characteristic;
		for ( let i = a.length - 1 ; i > -1 ; i-- ) {
			result *= p;
			result += (a[i] ?? 0);
		}
		return result;
	}

	/**
	 * Test whether a Polynomial is an element of the ring.
	 * Returns true if so;
	 * returns false if it is a well-formed Polynomial but not in the ring;
	 * throws TypeError if it is not a well-formed polynomial.
	 * @param {Polynomial} a
	 * @returns {boolean}
	 */
	element_ok(a) {
		if (!(
			Object.getPrototypeOf(a) === Array
			&& !a.every((x) =>
				(typeof x === 'number')
				&& (x <= Number.MAX_SAFE_INTEGER && (x%1) === 0)
				&& (x >= 0)
			)
		))
			throw new TypeError("Not an Array of natural numbers");

		if (
			a.length > 0 && (
				!((a.length - 1) in a)
				|| a[a.length - 1] === 0
			)
		)
			throw new TypeError("Has trailing zeroes or trailing empty slots");

		const p = this.characteristic;
		if (!a.every((x) => (x >= 0 && x < p)))
			// Has some coefficient out of range.
			return false;

		// Array of natural numbers
		// with no trailing zeroes
		// and every coefficient in-range for this ring.
		return true;
	}

	/**
	 * Addition in the ring.
	 * @param {Polynomial} a - Addend. Must be an element of the ring.
	 * @param {Polynomial} b - Addend. Must be an element of the ring.
	 * @returns {Polynomial}
	 */
	add(a, b) {
		const m = this.#characteristic;

		var c = a.slice();
		for (let [i, n] of Array_entriesSparse(b)) {
			c[i] = ((c[i] ?? 0) + n) % m;
		}
		trimTrailingZeroes(c);
		return c;
	}


	/**
	 * Subtraction in the ring.
	 * @param {Polynomial} a - Minuend. Must be an element of the ring.
	 * @param {Polynomial} b - Subtrahend. Must be an element of the ring.
	 * @returns {Polynomial}
	 */
	sub(a, b) {
		const m = this.#characteristic;

		var c = a.slice();
		for (let [i, n] of Array_entriesSparse(b)) {
			c[i] = ((c[i] ?? 0) - n + m) % m;
		}
		trimTrailingZeroes(c);
		return c;
	}

	/**
	 * Additive reciprocal in the ring.
	 * @param {Polynomial} a - The polynomial to find the reciprocal of. Must be an element of the ring.
	 * @returns {Polynomial}
	 */
	neg(a) {
		const m = this.#characteristic;

		return a.map((n) => (m - n));
	}

	/**
	 * Multiplication in the ring.
	 * @param {Polynomial} a - Multiplicand. Must be an element of the ring.
	 * @param {Polynomial} b - Multiplicand. Must be an element of the ring.
	 * @returns {Polynomial}
	 */
	mul(a, b) {
		const m = this.#characteristic;

		var c = [];
		for (let [aExp, aCoeff] of Array_entriesSparse(a)) {
			for (let [bExp, bCoeff] of Array_entriesSparse(b)) {
				let i = aExp + bExp,
				    n = (aCoeff * bCoeff) % m;
				c[i] = ((c[i] ?? 0) + n) % m;
			}
		}
		trimTrailingZeroes(c);
		return c;
	}

	/**
	 * Euclidean division in the ring. Returns the quotient and remainder, respectively.
	 * @param {Polynomial} a - Dividend. Must be an element of the ring.
	 * @param {Polynomial} b - Divisor. Must be an element of the ring.
	 * @returns {[Polynomial, Polynomial]}
	 */
	divmod(a, b) {
		const m = this.#characteristic;

		var q = [],
		    r = a.slice(),
		    bLeadExp = b.length - 1,
		    bLeadCoeff = b[bLeadExp];
		for (let [aiExp, aiCoeff] of Array_entriesSparseReversed(r)) {
			if (aiExp < bLeadExp)
				break;
			let qiExp = aiExp - bLeadExp,
			    qiCoeff = (aiCoeff * modInv(bLeadCoeff, m)) % m;

			// 1. Record this term in the quotient
			q[qiExp] = ((q[qiExp] ?? 0) + qiCoeff) % m;

			// 2. Subtract the product of this term and the divisor from the remainder
			for (let [biExp, biCoeff] of Array_entriesSparse(b)) {
				let i = biExp + qiExp,
				    n = (biCoeff * qiCoeff) % m;
				r[i] = ((r[i] ?? 0) - n + m) % m;
			}
		}
		trimTrailingZeroes(q, 1);
		trimTrailingZeroes(r, 1);
		return [q, r];
	}

	/**
	 * Return elements [g, x, y] such that this.representative(this.add(this.mul(a, x), this.mul(b, y))) === this.representative(g)
	 * @param {Polynomial}
	 * @param {Polynomial}
	 * @returns {[Polynomial,Polynomial,Polynomial]}
	 */
	egcd(a, b) {
		const divmod = this.divmod.bind(this),
		      sub = this.sub.bind(this),
		      mul = this.mul.bind(this);

		// h/t https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm?oldid=1294792213#Description
		var r1 = a.slice(),
		    s1 = this.element(1),
		    t1 = this.element(0);
		for (
			let r2 = b.slice(),
			    s2 = this.element(0),
			    t2 = this.element(1),
			    q2, r3, s3, t3;
			r2.length !== 0 ;
			[r1, s1, t1, r2, s2, t2] = [r2, s2, t2, r3, s3, t3]
		) {
			[q2, r3] = divmod(r1, r2);
			s3 = sub(s1, mul(q2, s2));
			t3 = sub(t1, mul(q2, t2));
		}
		return [r1, s1, t1];
	}

	/**
	 * @type {NaturalNumber}
	 */
	get characteristic() {
		return this.#characteristic;
	}
}


function* Array_entriesSparse(a) {
	for ( let end = a.length, i = 0 ; i < end ; i++ )
		if (i in a) // Skip empty slots
			yield [i, a[i]];
}


function* Array_entriesSparseReversed(a) {
	for ( let i = a.length - 1 ; i > -1 ; i-- ) {
		if (i in a) // Skip empty slots
			yield [i, a[i]];
	}
}


function trimTrailingZeroes(a) {
	var newLength = a.length;
	for ( let i = newLength - 1 ; i > -1 ; i-- ) {
		if (i in a && a[i] !== 0)
			break;
		// Empty slot, or slot containing 0; remove it
		newLength = i;
	}
	a.length = newLength;
}


/**
 * Modular inverse.
 * @param {NaturalNumber} a - The number to find the multiplicative reciprocal of.
 * @param {NaturalNumber} m - The modulus.
 * @returns {NaturalNumber}
 */
function modInv(a, m) {
	var a_ = a, b = 1;
	for (
		let m_ = m,
		    c = 0,
		    q, r;
		m_ !== 0;
		[a_, b, c, m_] = [m_, c, b - q * c, r]
	) {
		[q, r] = [Math.floor(a_ / m_), a_ % m_];
	}
	if (a_ !== 1)
		throw new RangeError("Not invertible");
	if (b < 0)
		b += m;
	return b;
}