mcalc/modules/math.js

/*
Copyright (C) 2025 filifa

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with this program.  If not, see <https://www.gnu.org/licenses/>.
*/
function xgcd(a, b) {
	let [old_r, r] = [a, b];
	let [old_s, s] = [1n, 0n];
	let [old_t, t] = [0n, 1n];

	while (r !== 0n) {
		const quotient = old_r / r;
		[old_r, r] = [r, old_r - quotient * r];
		[old_s, s] = [s, old_s - quotient * s];
		[old_t, t] = [t, old_t - quotient * t];
	}

	return [old_r, old_s, old_t];
}

function modinv(x, modulus) {
	let [r, s, t] = xgcd(x, modulus);
	if (r !== 1n) {
		throw new Error(`no inverse exists - ${x} and ${modulus} are not coprime`);
	}

	if (s < 0n) {
		s += modulus;
	}

	return s;
}

function modpow(base, exponent, modulus) {
	if (exponent < 0n) {
		const p = modpow(base, -exponent, modulus);
		return modinv(p, modulus);
	}

	if (modulus === 1n) {
		return 0n;
	}

	let result = 1n;
	base %= modulus;

	while (exponent > 0n) {
		if (exponent % 2n === 1n) {
			result *= base;
			result %= modulus;
		}

		exponent >>= 1n;
		base *= base;
		base %= modulus;
	}

	return result;
}

function factorTwos(n) {
	let s = 0;
	while (n % 2n === 0n) {
		n /= 2n;
		s++;
	}

	return [n, s];
}

function witness(a, n) {
	const [d, s] = factorTwos(n - 1n);

	let x = modpow(a, d, n);
	let y = null;
	for (let i = 0; i < s; i++) {
		y = modpow(x, 2n, n);
		if (y === 1n && x !== 1n && x !== n - 1n) {
			return true;
		}
		x = y
	}

	return y !== 1n
}

function randbigint() {
	return BigInt(Math.floor(Math.random() * Number.MAX_SAFE_INTEGER))
}

function isprime(n) {
	if (n === 2n) {
		return true;
	} else if (n % 2n === 0n) {
		return false;
	}

	const trials = 10;
	for (let i = 0; i < trials; i++) {
		const a = randbigint() % (n - 1n) + 1n;
		if (witness(a, n)) {
			return false;
		}
	}

	return true;
}

function legendreSymbol(a, p) {
	let r = modpow(a, (p - 1n)/2n, p);
	if (r === p - 1n) {
		r = -1n;
	}

	return r;
}

function quadraticNonResidue(p) {
	// TODO: consider randomizing this
	for (let a = 2n; a < p; a++) {
		if (legendreSymbol(a, p) === -1n) {
			return a;
		}
	}
}

function tonelliShanks(n, p) {
	const [q, s] = factorTwos(p - 1n);
	const z = quadraticNonResidue(p);

	let m = s;
	let c = modpow(z, q, p);
	let t = modpow(n, q, p);
	let r = modpow(n, (q+1n)/2n, p);

	while (true) {
		if (t === 0n) {
			return 0n;
		} else if (t === 1n) {
			return r;
		}

		let k = t;
		let i = null;
		for (i = 1; i < m; i++) {
			k = modpow(k, 2n, p);
			if (k === 1n) {
				break;
			}
		}

		if (i === m) {
			throw new Error("radicand is not a quadratic residue of the modulus");
		}

		const e = BigInt(Math.pow(2, m - i - 1));
		const b = modpow(c, e, p);
		m = i;
		c = modpow(b, 2n, p);
		t = (t * c) % p;
		r = (r * b) % p;
	}
}

function modsqrt(n, modulus) {
	// TODO: add support for prime power modulus (Hensel's lemma)
	if (!isprime(modulus)) {
		throw new Error("modulus must be prime to compute square root");
	}

	n %= modulus;
	if (n < 0n) {
		n += modulus;
	}

	if (n % modulus === 0n) {
		return 0n;
	} else if (modulus === 2n) {
		return n % 2n;
	} else if (legendreSymbol(n, modulus) !== 1n) {
		throw new Error("radicand is not a quadratic residue of the modulus");
	} else if (modulus % 4n === 3n) {
		return modpow(n, (modulus+1n)/4n, modulus);
	}

	return tonelliShanks(n, modulus);
}

function ord(n, modulus) {
	n %= modulus;
	if (n < 0n) {
		n += modulus;
	}

	const [g, s, t] = xgcd(n, modulus);
	if (g !== 1n) {
		throw new Error(`can't compute multiplicative order - ${n} and ${modulus} are not coprime`);
	}

	// NOTE: this is a hard problem, but there are more efficient approaches

	let k = 1n;
	let a = n;
	while (a !== 1n) {
		a *= n;
		a %= modulus;
		k += 1n;
	}

	return k;
}