Number Theory

// GCD: Euclidean algorithm
// gcd(a, b) = gcd(b, a % b) — reduces faster than subtraction
const gcd = (a, b) => b === 0 ? a : gcd(b, a % b);

// LCM: use gcd to avoid overflow
const lcm = (a, b) => (a / gcd(a, b)) * b;

// Extended Euclidean: finds x, y such that ax + by = gcd(a, b)
// Used for modular inverse when gcd = 1
function extGcd(a, b) {
    if (b === 0) return [a, 1, 0];
    const [g, x1, y1] = extGcd(b, a % b);
    return [g, y1, x1 - Math.floor(a/b) * y1];
}
// modular inverse of a mod m (when gcd(a,m) = 1):
// extGcd(a, m)[1] mod m

// GCD of array: reduce pairwise
const gcdArray = (arr) => arr.reduce((a, b) => gcd(a, b));

// Key properties:
// gcd(a, b) = gcd(b, a % b)
// gcd(a, 0) = a
// lcm(a, b) = a * b / gcd(a, b)
// gcd divides both a and b, and is the LARGEST such number

// Sieve: O(n log log n)
function sieve(n) {
    const isPrime = new Uint8Array(n+1).fill(1);
    isPrime[0] = isPrime[1] = 0;
    for (let p = 2; p*p <= n; p++)
        if (isPrime[p])
            for (let m = p*p; m <= n; m += p) isPrime[m] = 0;
    return isPrime;
}

// Smallest Prime Factor sieve — enables O(log n) factorization
function spfSieve(n) {
    const spf = Array.from({length:n+1}, (_, i) => i);
    for (let p = 2; p*p <= n; p++)
        if (spf[p] === p) // p is prime
            for (let m = p*p; m <= n; m += p)
                if (spf[m] === m) spf[m] = p;
    return spf;
}

// Factorize using SPF: O(log n) per number
function factorize(n, spf) {
    const factors = [];
    while (n > 1) { factors.push(spf[n]); n = Math.floor(n / spf[n]); }
    return factors;
}

// Factorize single number: O(sqrt n)
function primeFactors(n) {
    const factors = [];
    for (let p = 2; p*p <= n; p++)
        while (n % p === 0) { factors.push(p); n = Math.floor(n/p); }
    if (n > 1) factors.push(n);
    return factors;
}

const MOD = 1_000_000_007n;

// Fast power mod: a^b mod m in O(log b)
function powMod(base, exp, mod) {
    base = BigInt(base) % mod;
    let result = 1n;
    exp = BigInt(exp);
    while (exp > 0n) {
        if (exp & 1n) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1n;
    }
    return result;
}

// Modular inverse using Fermat's little theorem (mod must be prime):
// inv(a) = a^(mod-2) mod mod
const modInverse = (a) => powMod(a, MOD - 2n, MOD);

// Precompute factorials + inverse factorials for nCr mod p
const MAX = 200001;
const fact = new Array(MAX);
const inv_fact = new Array(MAX);
fact[0] = 1n;
for (let i=1;i<MAX;i++) fact[i] = fact[i-1] * BigInt(i) % MOD;
inv_fact[MAX-1] = powMod(fact[MAX-1], MOD-2n, MOD);
for (let i=MAX-2;i>=0;i--) inv_fact[i] = inv_fact[i+1] * BigInt(i+1) % MOD;

const nCr = (n, r) => {
    if (r < 0 || r > n) return 0n;
    return fact[n] * inv_fact[r] % MOD * inv_fact[n-r] % MOD;
};

// Totient of a single number: O(sqrt n)
function totient(n) {
    let result = n;
    for (let p = 2; p*p <= n; p++) {
        if (n % p === 0) {
            while (n % p === 0) n = Math.floor(n/p);
            result -= Math.floor(result / p); // result *= (1 - 1/p)
        }
    }
    if (n > 1) result -= Math.floor(result / n); // n is a prime factor
    return result;
}

// Totient sieve: compute φ(i) for all i in [1, n] — O(n log log n)
function totientSieve(n) {
    const phi = Array.from({length:n+1}, (_, i) => i);
    for (let p = 2; p <= n; p++) {
        if (phi[p] === p) { // p is prime
            for (let m = p; m <= n; m += p)
                phi[m] -= Math.floor(phi[m] / p);
        }
    }
    return phi;
}