Fast Exponentiation

Recursive Approach

function fastExponentiation(base, exponent, modulus) {
    if (exponent === 0) return 1;
    
    const halfExp = fastExponentiation(base, Math.floor(exponent / 2), modulus);
    let result = (halfExp * halfExp) % modulus;
    
    if (exponent % 2 === 1) {
        result = (result * base) % modulus;
    }
    
    return result;
}

Iterative Approach

function fastExponentiation(base, exponent, modulus) {
    let result = 1;
    
    while (exponent > 0) {
        if (exponent % 2 === 1) {
            result = (result * base) % modulus;
        }
        
        base = (base * base) % modulus;
        exponent = Math.floor(exponent / 2);
    }
    
    return result;
}

Matrix Fast Exponentiation

Use with the Matrix class implementation provided below.

function matrixFastExponentiation(matrix, power) {
    let result = matrix;
    power--;
    
    while (power > 0) {
        if (power % 2 === 1) {
            result = multiplyMatrices(result, matrix);
        }
        
        matrix = multiplyMatrices(matrix, matrix);
        power = Math.floor(power / 2);
    }
    
    return result;
}

Euclidean Algorithm

Greatest Common Divisor (GCD)

function gcd(a, b) {
    return b === 0 ? a : gcd(b, a % b);
}

Least Common Multiple (LCM)

function lcm(a, b) {
    return (a * b) / gcd(a, b);
}

Extended Euclidean Algroithm

function extendedEuclidean(a, b, coefficients) {
    if (b === 0) {
        coefficients.x = 1;
        coefficients.y = 0;
        return a;
    }
    
    const tempCoefficients = {x: 0, y: 0};
    const result = extendedEuclidean(b, a % b, tempCoefficients);
    
    coefficients.x = tempCoefficients.y;
    coefficients.y = tempCoefficients.x - Math.floor(a / b) * tempCoefficients.y;
    
    return result;
}

Multiplicative Inverse

Using Fermat's Little Theorem

function multiplicativeInverse(num, prime) {
    return fastExponentiation(num, prime - 2, prime);
}

Inverses of All Elements in an Array

Formula: \(a_i^{-1} = \prod_{j=1}^{i-1}a_j \times \left( \prod_{j=1}^{i}a_j \right)^{-1} \)

function computeArrayInverses(array, size, modulus) {
    const prefixProduct = new Array(size + 1).fill(1);
    
    // Compute prefix products
    for (let i = 1; i <= size; i++) {
        prefixProduct[i] = (prefixProduct[i - 1] * array[i]) % modulus;
    }
    
    // Compute inverse of the full product
    const totalInverse = fastExponentiation(prefixProduct[size], modulus - 2, modulus);
    
    const inverses = new Array(size + 1);
    let currentInverse = totalInverse;
    
    // Compute individual inverses
    for (let i = size; i >= 1; i--) {
        inverses[i] = (currentInverse * prefixProduct[i - 1]) % modulus;
        currentInverse = (currentInverse * array[i]) % modulus;
    }
    
    return inverses;
}

Permutations and Combinations

Comibnation Using Multiplicative Inverses

function precomputeFactorials(maxN, modulus) {
    const factorial = new Array(maxN + 1).fill(1);
    const inverseFactorial = new Array(maxN + 1).fill(1);
    
    for (let i = 1; i <= maxN; i++) {
        factorial[i] = (factorial[i - 1] * i) % modulus;
    }
    
    inverseFactorial[maxN] = fastExponentiation(factorial[maxN], modulus - 2, modulus);
    
    for (let i = maxN - 1; i >= 0; i--) {
        inverseFactorial[i] = (inverseFactorial[i + 1] * (i + 1)) % modulus;
    }
    
    return { factorial, inverseFactorial };
}

function combination(n, k, factorial, inverseFactorial, modulus) {
    if (k < 0 || k > n) return 0;
    
    return (factorial[n] * inverseFactorial[n - k] % modulus) * 
           inverseFactorial[k] % modulus;
}

Gaussian Elimination

Gauss-Jordan Elimination

const EPSILON = 1e-6;

function gaussJordanElimination(matrix, n) {
    let currentRow = 0;
    
    for (let col = 0; col < n; col++) {
        // Find pivot row
        let pivotRow = currentRow;
        while (pivotRow < n && Math.abs(matrix[pivotRow][col]) < EPSILON) {
            pivotRow++;
        }
        
        if (pivotRow === n) continue; // No pivot found in this column
        
        // Swap current row with pivot row
        [matrix[currentRow], matrix[pivotRow]] = [matrix[pivotRow], matrix[currentRow]];
        
        // Normalize the pivot row
        const pivotValue = matrix[currentRow][col];
        for (let j = col; j <= n; j++) {
            matrix[currentRow][j] /= pivotValue;
        }
        
        // Eliminate this variable from all other rows
        for (let i = 0; i < n; i++) {
            if (i === currentRow) continue;
            
            const factor = matrix[i][col];
            for (let j = col; j <= n; j++) {
                matrix[i][j] -= factor * matrix[currentRow][j];
            }
        }
        
        currentRow++;
    }
    
    // Check for solutions
    if (currentRow < n) {
        for (let i = currentRow; i < n; i++) {
            if (Math.abs(matrix[i][n]) > EPSILON) {
                return -1; // No solution
            }
        }
        return 0; // Infinite solutions
    }
    
    return 1; // Unique solution
}

Matrix Operations

class Matrix {
    constructor(rows = 0, cols = 0, modulus = 1000000007) {
        this.rows = rows;
        this.cols = cols;
        this.modulus = modulus;
        this.data = Array(rows).fill().map(() => Array(cols).fill(0));
    }
    
    multiply(other) {
        if (this.cols !== other.rows) {
            throw new Error("Matrix dimensions incompatible for multiplication");
        }
        
        const result = new Matrix(this.rows, other.cols, this.modulus);
        
        for (let i = 0; i < this.rows; i++) {
            for (let j = 0; j < other.cols; j++) {
                for (let k = 0; k < this.cols; k++) {
                    result.data[i][j] = (result.data[i][j] + 
                                         this.data[i][k] * other.data[k][j]) % this.modulus;
                }
            }
        }
        
        return result;
    }
    
    static identity(size, modulus) {
        const matrix = new Matrix(size, size, modulus);
        for (let i = 0; i < size; i++) {
            matrix.data[i][i] = 1;
        }
        return matrix;
    }
}

Extendde Chinese Remainder Theorem

function extendedChineseRemainderTheorem(remainders, moduli) {
    let currentRemainder = remainders[0];
    let currentModulus = moduli[0];
    let hasNoSolution = false;
    
    for (let i = 1; i < remainders.length; i++) {
        const a = remainders[i];
        const m = moduli[i];
        
        // Solve: currentRemainder ≡ x (mod currentModulus)
        //         a ≡ x (mod m)
        const coefficients = {x: 0, y: 0};
        const gcdValue = extendedEuclidean(currentModulus, -m, coefficients);
        
        if ((a - currentRemainder) % gcdValue !== 0) {
            hasNoSolution = true;
            break;
        }
        
        const k = (a - currentRemainder) / gcdValue;
        const lcmValue = (currentModulus * m) / gcdValue;
        
        const newRemainder = (currentRemainder + coefficients.y * k * m) % lcmValue;
        currentRemainder = (newRemainder + lcmValue) % lcmValue;
        currentModulus = lcmValue;
    }
    
    if (hasNoSolution) {
        return -1; // No solution exists
    }
    
    return (currentRemainder + currentModulus) % currentModulus;
}