Add BigInt utility functions for number theory operations (#2205)

Co-authored-by: GCHQDeveloper581 <63102987+GCHQDeveloper581@users.noreply.github.com> (initial tests)

src/core/config/Categories.json

@@ -222,6 +222,9 @@
             "Subtract",
             "Multiply",
             "Divide",
+            "Modular Exponentiation",
+            "Modular Inverse",
+            "Extended GCD",
             "Mean",
             "Median",
             "Standard Deviation",

src/core/lib/BigIntUtils.mjs

@@ -0,0 +1,73 @@
+/**
+ * @author p-leriche [philip.leriche@cantab.net]
+ * @copyright Crown Copyright 2025
+ * @license Apache-2.0
+ */
+
+import OperationError from "../errors/OperationError.mjs";
+
+/**
+ * Number theory utilities used by cryptographic operations.
+ *
+ * Currently provides:
+ *   - parseBigInt
+ *   - Extended Euclidean Algorithm
+ *   - Modular Exponentiation
+ *
+ * Additional algorithms may be added as required.
+ */
+
+/**
+ * parseBigInt helper operation
+ */
+export function parseBigInt(value, param) {
+    const v = (value ?? "").trim();
+    if (/^0x[0-9a-f]+$/i.test(v)) return BigInt(v);
+    if (/^[+-]?[0-9]+$/.test(v)) return BigInt(v);
+    throw new OperationError(param + " must be decimal or hex (0x...)");
+}
+
+/**
+ * Extended Euclidean Algorithm
+ *
+ * Returns [g, x, y] such that:
+ *   a*x + b*y = g = gcd(a, b)
+ *
+ *   (Uses an iterative algorithm to avoid possible stack overflow)
+ */
+export function egcd(a, b) {
+    let oldR = a, r = b;
+    let oldS = 1n, s = 0n;
+    let oldT = 0n, t = 1n;
+
+    while (r !== 0n) {
+        const quotient = oldR / r;
+
+        [oldR, r] = [r, oldR - quotient * r];
+        [oldS, s] = [s, oldS - quotient * s];
+        [oldT, t] = [t, oldT - quotient * t];
+    }
+
+    // oldR is the gcd
+    // oldS and oldT are the Bézout coefficients
+    return [oldR, oldS, oldT];
+}
+
+/**
+ * Modular exponentiation
+ */
+export function modPow(base, exponent, modulus) {
+    let result = 1n;
+    base %= modulus;
+
+    while (exponent > 0n) {
+        if (exponent & 1n) {
+            result = (result * base) % modulus;
+        }
+        base = (base * base) % modulus;
+        exponent >>= 1n;
+    }
+
+    return result;
+}
+