Code Templates

This template collection provides clean implementations of fundamental mathematical algorithms in multiple programming languages commonly used in coding interviews. For a theoretical overview, check out the Math Algorithms Guide.

GCD (Greatest Common Divisor) - Euclidean Algorithm

JavaScript:

function gcd(a, b) {
    while (b !== 0) {
        const temp = b;
        b = a % b;
        a = temp;
    }
    return a;
}

// Extended GCD for modular inverse
function extendedGcd(a, b) {
    if (b === 0) {
        return { gcd: a, x: 1, y: 0 };
    }
    const { gcd, x: x1, y: y1 } = extendedGcd(b, a % b);
    const x = y1;
    const y = x1 - Math.floor(a / b) * y1;
    return { gcd, x, y };
}

function modInverse(a, m) {
    const { gcd, x } = extendedGcd(a, m);
    if (gcd !== 1) {
        throw new Error('Modular inverse does not exist');
    }
    return (x % m + m) % m;
}

Python:

from math import gcd
from typing import Tuple

def extended_gcd(a: int, b: int) -> Tuple[int, int, int]:
    if b == 0:
        return a, 1, 0
    gcd, x1, y1 = extended_gcd(b, a % b)
    x = y1
    y = x1 - (a // b) * y1
    return gcd, x, y

def mod_inverse(a: int, m: int) -> int:
    gcd, x, _ = extended_gcd(a, m)
    if gcd != 1:
        raise ValueError("Modular inverse does not exist")
    return (x % m + m) % m

def gcd_multiple(*args: int) -> int:
    if len(args) == 0:
        raise ValueError("At least one number is required")
    result = args[0]
    for num in args[1:]:
        result = gcd(result, num)
    return result

Java:

public class MathUtils {
    public static long gcd(long a, long b) {
        while (b != 0) {
            long temp = b;
            b = a % b;
            a = temp;
        }
        return a;
    }

    public static class ExtendedGcd {
        long gcd, x, y;
        ExtendedGcd(long gcd, long x, long y) {
            this.gcd = gcd; this.x = x; this.y = y;
        }
    }

    public static ExtendedGcd extendedGcd(long a, long b) {
        if (b == 0) {
            return new ExtendedGcd(a, 1, 0);
        }
        ExtendedGcd result = extendedGcd(b, a % b);
        long x = result.y;
        long y = result.x - (a / b) * result.y;
        return new ExtendedGcd(result.gcd, x, y);
    }

    public static long modInverse(long a, long m) {
        ExtendedGcd result = extendedGcd(a, m);
        if (result.gcd != 1) {
            throw new ArithmeticException("Modular inverse does not exist");
        }
        return (result.x % m + m) % m;
    }
}

C++:

#include <tuple>
#include <stdexcept>

using std::tuple;

long long gcd(long long a, long long b) {
    while (b != 0) {
        long long temp = b;
        b = a % b;
        a = temp;
    }
    return a;
}

tuple<long long, long long, long long> extended_gcd(long long a, long long b) {
    if (b == 0) {
        return {a, 1LL, 0LL};
    }
    auto [gcd_val, x1, y1] = extended_gcd(b, a % b);
    long long x = y1;
    long long y = x1 - (a / b) * y1;
    return {gcd_val, x, y};
}

long long mod_inverse(long long a, long long m) {
    auto [gcd_val, x, y] = extended_gcd(a, m);
    if (gcd_val != 1) {
        throw std::runtime_error("Modular inverse does not exist");
    }
    return (x % m + m) % m;
}

Go:

func gcd(a, b int64) int64 {
    for b != 0 {
        a, b = b, a%b
    }
    return a
}

type ExtGcdResult struct {
    Gcd, X, Y int64
}

func extendedGcd(a, b int64) ExtGcdResult {
    if b == 0 {
        return ExtGcdResult{Gcd: a, X: 1, Y: 0}
    }
    result := extendedGcd(b, a%b)
    x := result.Y
    y := result.X - (a/b)*result.Y
    return ExtGcdResult{Gcd: result.Gcd, X: x, Y: y}
}

func modInverse(a, m int64) (int64, error) {
    result := extendedGcd(a, m)
    if result.Gcd != 1 {
        return 0, fmt.Errorf("modular inverse does not exist")
    }
    inverse := (result.X%m + m) % m
    return inverse, nil
}

Ruby:

def gcd(a, b)
  while b != 0
    a, b = b, a % b
  end
  a
end

# Extended GCD for modular inverse
def extended_gcd(a, b)
  if b == 0
    return { gcd: a, x: 1, y: 0 }
  end
  result = extended_gcd(b, a % b)
  {
    gcd: result[:gcd],
    x: result[:y],
    y: result[:x] - (a / b) * result[:y]
  }
end

def mod_inverse(a, m)
  result = extended_gcd(a, m)
  if result[:gcd] != 1
    raise 'Modular inverse does not exist'
  end
  (result[:x] % m + m) % m
end

Code Breakdown

The Euclidean GCD algorithm uses a simple but powerful approach:

Key Variables:
- a, b: The two numbers to find GCD for
- temp: Temporary storage during swap
- result: Stores final GCD value

Algorithm Flow:

  graph TD
    Start[Start with a, b]
    Check{b == 0?}
    Yes -->|Yes| Return[Return a as GCD]
    No -->|No| Swap[temp = b]
    No --> Mod[b = a mod b]
    No --> Swap2[a = temp]
    No --> Check
    style Return fill:#6f6,stroke:#333

Critical Sections:
- Initialization: Start with two input numbers
- Loop: While b is not zero, keep taking remainder
- Swap: Exchange values to continue with smaller numbers
- Termination: When remainder is zero, return the divisor

The magic is in how quickly numbers reduce - each iteration significantly shrinks the problem size.

Sieve of Eratosthenes

JavaScript:

function sieve(maxN) {
    const isPrime = new Array(maxN + 1).fill(true);
    isPrime[0] = isPrime[1] = false;

    const primes = [];

    for (let i = 2; i <= maxN; i++) {
        if (isPrime[i]) {
            primes.push(i);
            for (let j = i * i; j <= maxN; j += i) {
                isPrime[j] = false;
            }
        }
    }

    return primes;
}

function smallestPrimeFactor(n) {
    const spf = new Array(n + 1);
    for (let i = 0; i <= n; i++) {
        spf[i] = i;
    }

    for (let i = 2; i * i <= n; i++) {
        if (spf[i] === i) { // i is prime
            for (let j = i * i; j <= n; j += i) {
                if (spf[j] === j) {
                    spf[j] = i;
                }
            }
        }
    }

    return spf;
}

Python:

from typing import List

def sieve(max_n: int) -> List[int]:
    if max_n < 2:
        return []

    is_prime = [True] * (max_n + 1)
    is_prime[0] = is_prime[1] = False

    primes = []
    for i in range(2, max_n + 1):
        if is_prime[i]:
            primes.append(i)
            for j in range(i * i, max_n + 1, i):
                is_prime[j] = False

    return primes

def smallest_prime_factor(n: int) -> List[int]:
    spf = list(range(n + 1))
    for i in range(2, int(n**0.5) + 1):
        if spf[i] == i:  # i is prime
            for j in range(i * i, n + 1, i):
                if spf[j] == j:
                    spf[j] = i
    return spf

def prime_factors(x: int, spf: List[int]) -> List[int]:
    factors = []
    while x > 1:
        factors.append(spf[x])
        x //= spf[x]
    return list(set(factors))  # unique factors only

Java:

import java.util.*;

public class SieveUtils {
    public static List<Integer> sieve(int maxN) {
        if (maxN < 2) return new ArrayList<>();

        boolean[] isPrime = new boolean[maxN + 1];
        Arrays.fill(isPrime, true);
        isPrime[0] = isPrime[1] = false;

        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= maxN; i++) {
            if (isPrime[i]) {
                primes.add(i);
                for (long j = (long)i * i; j <= maxN; j += i) {
                    isPrime[(int)j] = false;
                }
            }
        }
        return primes;
    }

    public static int[] smallestPrimeFactor(int n) {
        int[] spf = new int[n + 1];
        for (int i = 0; i <= n; i++) {
            spf[i] = i;
        }

        for (int i = 2; i * i <= n; i++) {
            if (spf[i] == i) { // i is prime
                for (long j = (long)i * i; j <= n; j += i) {
                    if (spf[(int)j] == j) {
                        spf[(int)j] = i;
                    }
                }
            }
        }

        return spf;
    }
}

C++:

#include <vector>
#include <cmath>

using std::vector;

vector<int> sieve(int max_n) {
    if (max_n < 2) return {};

    vector<bool> is_prime(max_n + 1, true);
    is_prime[0] = is_prime[1] = false;

    vector<int> primes;
    for (long long i = 2; i <= max_n; ++i) {
        if (is_prime[i]) {
            primes.push_back(i);
            for (long long j = i * i; j <= max_n; j += i) {
                is_prime[j] = false;
            }
        }
    }

    return primes;
}

vector<int> smallest_prime_factor(int n) {
    vector<int> spf(n + 1);
    for (int i = 0; i <= n; ++i) {
        spf[i] = i;
    }

    for (long long i = 2; i * i <= n; ++i) {
        if (spf[i] == i) { // i is prime
            for (long long j = i * i; j <= n; j += i) {
                if (spf[j] == j) {
                    spf[j] = i;
                }
            }
        }
    }

    return spf;
}

vector<int> prime_factors(int x, const vector<int>& spf) {
    vector<int> factors;
    while (x > 1) {
        factors.push_back(spf[x]);
        x /= spf[x];
    }
    // Remove duplicates
    sort(factors.begin(), factors.end());
    factors.erase(unique(factors.begin(), factors.end()), factors.end());
    return factors;
}

Go:

func sieve(maxN int) []int {
    if maxN < 2 {
        return []int{}
    }

    isPrime := make([]bool, maxN+1)
    for i := range isPrime {
        isPrime[i] = true
    }
    isPrime[0], isPrime[1] = false, false

    primes := []int{}
    for i := 2; i <= maxN; i++ {
        if isPrime[i] {
            primes = append(primes, i)
            for j := i * i; j <= maxN && j > 0; j += i {
                isPrime[j] = false
            }
        }
    }

    return primes
}

func smallestPrimeFactor(n int) []int {
    spf := make([]int, n+1)
    for i := 0; i <= n; i++ {
        spf[i] = i
    }

    for i := 2; i*i <= n; i++ {
        if spf[i] == i { // i is prime
            for j := i * i; j <= n && j > 0; j += i {
                if spf[j] == j {
                    spf[j] = i
                }
            }
        }
    }

    return spf
}

func primeFactors(x int, spf []int) []int {
    factors := []int{}
    for x > 1 {
        factors = append(factors, spf[x])
        x /= spf[x]
    }

    // Remove duplicates
    seen := make(map[int]bool)
    result := []int{}
    for _, factor := range factors {
        if !seen[factor] {
            seen[factor] = true
            result = append(result, factor)
        }
    }

    return result
}

Ruby:

def sieve(max_n)
  return [] if max_n < 2

  is_prime = Array.new(max_n + 1, true)
  is_prime[0] = is_prime[1] = false

  primes = []
  (2..max_n).each do |i|
    if is_prime[i]
      primes << i
      j = i * i
      while j <= max_n
        is_prime[j] = false
        j += i
      end
    end
  end

  primes
end

def smallest_prime_factor(n)
  spf = (0..n).to_a
  (2..Integer.sqrt(n)).each do |i|
    next unless spf[i] == i # i is prime
    j = i * i
    while j <= n
      spf[j] = i if spf[j] == j
      j += i
    end
  end
  spf
end

def prime_factors(x, spf)
  factors = []
  while x > 1
    factors << spf[x]
    x /= spf[x]
  end
  factors.uniq
end

Code Breakdown

The Sieve algorithm efficiently finds all primes up to a given limit:

Key Variables:
- isPrime: Boolean array tracking if each number is prime
- primes: List collecting discovered prime numbers
- spf: Array storing smallest prime factor for each number

Algorithm Flow:

  graph TD
    Init[Initialize isPrime array]
    Start[Start from 2]
    IsPrime{isPrime i?}
    Yes --> Add[Add to primes list]
    Yes --> Mark[Mark multiples as composite]
    No --> Next[Continue to next]
    Next --> IsPrime
    style Add fill:#6f6,stroke:#333

Critical Sections:
- Initialization: Create array with all values initially prime
- Outer Loop: Iterate from 2 to max_n
- Prime Check: If current number is prime, add to result
- Inner Loop: Mark all multiples of prime as composite (start from i*i for efficiency)
- Smallest Prime Factor: Extended version tracks the smallest prime dividing each number for factorization

The key optimization is starting marking from i*i instead of 2*i - smaller multiples are already marked by previous primes.

Fast Exponentiation

JavaScript:

function modPow(base, exponent, modulus) {
    let result = 1;
    base = base % modulus;

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

    return result;
}

// For very large exponents, use array representation
function superPow(a, b) { // b is array of digits
    const MOD = 1337;
    let result = 1;
    a = a % MOD;

    for (let i = b.length - 1; i >= 0; i--) {
        let digit = b[i];
        let current = 1;

        // Calculate a^digit % MOD
        for (let d = 0; d < digit; d++) {
            current = (current * a) % MOD;
        }

        result = (result * current) % MOD;
        a = modPow(a, 10, MOD);
    }

    return result;
}

Python:

def mod_pow(base: int, exponent: int, modulus: int) -> int:
    if modulus == 1:
        return 0

    result = 1
    base = base % modulus
    while exponent > 0:
        if exponent % 2 == 1:
            result = (result * base) % modulus
        base = (base * base) % modulus
        exponent //= 2
    return result

def super_pow(a: int, b: List[int]) -> int:
    # b is array of digits for exponent
    MOD = 1337
    result = 1
    a = a % MOD

    for digit in reversed(b):
        # Calculate pow(a, digit, MOD)
        pow_digit = 1
        for _ in range(digit):
            pow_digit = (pow_digit * a) % MOD

        result = (result * pow_digit) % MOD
        a = mod_pow(a, 10, MOD)

    return result

def matrix_pow(matrix: List[List[int]], exponent: int, mod: int) -> List[List[int]]:
    n = len(matrix)
    result = [[1 if i == j else 0 for j in range(n)] for i in range(n)]  # identity

    while exponent > 0:
        if exponent % 2 == 1:
            result = matrix_multiply(result, matrix, mod)
        matrix = matrix_multiply(matrix, matrix, mod)
        exponent //= 2

    return result

Java:

public class FastExponentiation {
    public static long modPow(long base, long exponent, long modulus) {
        if (modulus == 1) return 0;

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

    public static int superPow(int a, int[] b) {
        final int MOD = 1337;
        int result = 1;
        a = a % MOD;

        for (int digit : b) {
            int powDigit = 1;
            for (int d = 0; d < digit; d++) {
                powDigit = (powDigit * a) % MOD;
            }
            result = (result * powDigit) % MOD;
            a = (int)modPow(a, 10, MOD);
        }

        return result;
    }
}

C++:

#include <vector>

using std::vector;

long long mod_pow(long long base, long long exponent, long long modulus) {
    if (modulus == 1) return 0;

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

int super_pow(int a, vector<int> b) {
    const int MOD = 1337;
    int result = 1;
    a = a % MOD;

    for (int digit : b) {
        int pow_digit = 1;
        for (int d = 0; d < digit; d++) {
            pow_digit = (pow_digit * a) % MOD;
        }
        result = (result * pow_digit) % MOD;
        a = (int)mod_pow(a, 10, MOD);
    }

    return result;
}

using Matrix = vector<vector<long long>>;
Matrix matrix_multiply(const Matrix& a, const Matrix& b, long long mod) {
    int n = a.size();
    int m = b[0].size();
    int p = a[0].size();
    Matrix result(n, vector<long long>(m, 0));

    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            for (int k = 0; k < p; k++) {
                result[i][j] = (result[i][j] + a[i][k] * b[k][j]) % mod;
            }
        }
    }

    return result;
}

Go:

func modPow(base, exponent, modulus int64) int64 {
    if modulus == 1 {
        return 0
    }

    result := int64(1)
    base = base % modulus
    for exponent > 0 {
        if exponent%2 == 1 {
            result = (result * base) % modulus
        }
        base = (base * base) % modulus
        exponent /= 2
    }

    return result
}

func superPow(a int, b []int) int {
    const MOD = 1337
    result := 1
    a = a % MOD

    for _, digit := range b {
        powDigit := 1
        for d := 0; d < digit; d++ {
            powDigit = (powDigit * a) % MOD
        }
        result = (result * powDigit) % MOD
        a = int(modPow(int64(a), 10, MOD))
    }

    return result
}

type Matrix [][]int64
func matrixMultiply(a, b Matrix, mod int64) Matrix {
    n := len(a)
    m := len(b[0])
    p := len(a[0])
    result := make(Matrix, n)
    for i := range result {
        result[i] = make([]int64, m)
    }

    for i := 0; i < n; i++ {
        for j := 0; j < m; j++ {
            for k := 0; k < p; k++ {
                result[i][j] = (result[i][j] + a[i][k]*b[k][j]) % mod
            }
        }
    }

    return result
}

Ruby:

def mod_pow(base, exponent, modulus)
  return 0 if modulus == 1

  result = 1
  base = base % modulus
  while exponent > 0
    result = (result * base) % modulus if exponent.odd?
    base = (base * base) % modulus
    exponent /= 2
  end
  result
end

def super_pow(a, b)
  # b is array of digits for exponent
  MOD = 1337
  result = 1
  a = a % MOD

  b.reverse.each do |digit|
    pow_digit = 1
    digit.times do
      pow_digit = (pow_digit * a) % MOD
    end
    result = (result * pow_digit) % MOD
    a = mod_pow(a, 10, MOD)
  end

  result
end

def matrix_pow(matrix, exponent, mod)
  n = matrix.size
  result = Array.new(n) { |i| Array.new(n) { |j| i == j ? 1 : 0 } } # identity

  while exponent > 0
    result = matrix_multiply(result, matrix, mod) if exponent.odd?
    matrix = matrix_multiply(matrix, matrix, mod)
    exponent /= 2
  end

  result
end

def matrix_multiply(a, b, mod)
  n = a.size
  m = b[0].size
  p = a[0].size
  result = Array.new(n) { Array.new(m, 0) }

  (0...n).each do |i|
    (0...m).each do |j|
      (0...p).each do |k|
        result[i][j] = (result[i][j] + a[i][k] * b[k][j]) % mod
      end
    end
  end

  result
end

Code Breakdown

Fast exponentiation computes base^exponent mod modulus efficiently using binary representation:

Key Variables:
- base: Current power of base being considered
- exponent: Remaining bits to process
- result: Accumulated result
- modulus: Modulo value for keeping numbers manageable

Algorithm Flow:

  graph TD
    Start[Start with base, exponent]
    Loop{exponent > 0?}
    Yes --> BitCheck{exponent odd?}
    BitCheck -->|Yes| Multiply[result = result * base]
    BitCheck -->|No| Square[base = base * base]
    Multiply --> Square
    Square --> Halve[exponent = exponent / 2]
    Halve --> Loop
    No --> Done[Return result]
    style Done fill:#6f6,stroke:#333

Critical Sections:
- Initialization: Set result to 1, apply modulus to base
- Bit Processing: If current bit of exponent is 1, multiply result by current base
- Squaring: Double the power of base in each step
- Halving: Move to next bit of exponent by integer division
- Modulo: Apply modulus at each multiplication to prevent overflow

The power is in the exponent’s binary representation: we multiply in only when bits are set, achieving O(log n) complexity instead of O(n).

Combinatorics

JavaScript:

function binomialCoefficient(n, k, mod = null) {
    if (k > n) return 0;
    if (k === 0 || k === n) return 1;

    k = Math.min(k, n - k);
    let res = 1;

    for (let i = 1; i <= k; i++) {
        res *= (n - i + 1);
        if (mod) res %= mod;
        res /= i;
    }

    return mod ? Math.round(res) % mod : res;
}

// Precompute factorials for efficiency
const MAX_N = 100000;
const MOD = 1000000007;

const factorial = new Array(MAX_N + 1).fill(0);
factorial[0] = 1;
for (let i = 1; i <= MAX_N; i++) {
    factorial[i] = (factorial[i - 1] * i) % MOD;
}

const invFactorial = new Array(MAX_N + 1).fill(0);
invFactorial[MAX_N] = modPow(factorial[MAX_N], MOD - 2, MOD);
for (let i = MAX_N - 1; i >= 0; i--) {
    invFactorial[i] = (invFactorial[i + 1] * (i + 1)) % MOD;
}

function binomialMod(n, k) {
    if (k > n || k < 0) return 0;
    return factorial[n] * invFactorial[k] % MOD * invFactorial[n - k] % MOD;
}

Python:

from functools import lru_cache
import sys

def binomial_coefficient(n: int, k: int, mod: int = None) -> int:
    if k > n:
        return 0
    if k == 0 or k == n:
        return 1

    k = min(k, n - k)
    res = 1
    for i in range(1, k + 1):
        res *= (n - i + 1)
        if mod:
            res %= mod
        res //= i

    return res % mod if mod else res

# Precomputed factorials
MAX_N = 10**5
MOD = 10**9 + 7

factorial = [0] * (MAX_N + 1)
factorial[0] = 1
for i in range(1, MAX_N + 1):
    factorial[i] = (factorial[i - 1] * i) % MOD

def mod_inverse(x: int, mod: int) -> int:
    return pow(x, mod - 2, mod)

inv_factorial = [0] * (MAX_N + 1)
inv_factorial[MAX_N] = mod_inverse(factorial[MAX_N], MOD)
for i in range(MAX_N - 1, -1, -1):
    inv_factorial[i] = (inv_factorial[i + 1] * (i + 1)) % MOD

def binomial_mod(n: int, k: int) -> int:
    if k > n or k < 0:
        return 0
    return factorial[n] * inv_factorial[k] % MOD * inv_factorial[n - k] % MOD

@lru_cache(maxsize=50000)
def stirling2(n: int, k: int) -> int:
    """Stirling numbers of the second kind"""
    if k > n or k < 0:
        return 0
    if k == 0 or k == n:
        return 1
    return (stirling2(n - 1, k - 1) + k * stirling2(n - 1, k)) % MOD

Java:

public class Combinatorics {
    private static final int MAX_N = 100000;
    private static final long MOD = 1000000007L;
    private static long[] factorial = new long[MAX_N + 1];
    private static long[] invFactorial = new long[MAX_N + 1];

    static {
        factorial[0] = 1;
        for (int i = 1; i <= MAX_N; i++) {
            factorial[i] = factorial[i - 1] * i % MOD;
        }
        invFactorial[MAX_N] = modInverse(factorial[MAX_N], MOD);
        for (int i = MAX_N - 1; i >= 0; i--) {
            invFactorial[i] = invFactorial[i + 1] * (i + 1) % MOD;
        }
    }

    public static long binomialMod(int n, int k) {
        if (k > n || k < 0) return 0;
        return factorial[n] * invFactorial[k] % MOD * invFactorial[n - k] % MOD;
    }

    public static long binomialCoefficient(int n, int k) {
        if (k > n) return 0;
        if (k == 0 || k == n) return 1;

        k = Math.min(k, n - k);
        long res = 1;
        for (int i = 1; i <= k; i++) {
            res = res * (n - i + 1) / i;
        }
        return res;
    }

    private static long modInverse(long a, long m) {
        return modPow(a, m - 2, m);
    }

    public static long modPow(long base, long exp, long mod) {
        long result = 1;
        base %= mod;
        while (exp > 0) {
            if (exp % 2 == 1) {
                result = result * base % mod;
            }
            base = base * base % mod;
            exp /= 2;
        }
        return result;
    }
}

C++:

#include <vector>
#include <cstdint>

using namespace std;

const int64_t MAX_N = 1e5;
const int64_t MOD = 1e9 + 7;

vector<int64_t> factorial;
vector<int64_t> inv_factorial;

void precompute_factorials() {
    factorial.resize(MAX_N + 1);
    factorial[0] = 1;
    for (int64_t i = 1; i <= MAX_N; ++i) {
        factorial[i] = factorial[i - 1] * i % MOD;
    }

    inv_factorial.resize(MAX_N + 1);
    inv_factorial[MAX_N] = mod_pow(factorial[MAX_N], MOD - 2, MOD);
    for (int64_t i = MAX_N - 1; i >= 0; --i) {
        inv_factorial[i] = inv_factorial[i + 1] * (i + 1) % MOD;
    }
}

int64_t binomial_mod(int64_t n, int64_t k) {
    if (k > n || k < 0) return 0;
    return factorial[n] * inv_factorial[k] % MOD * inv_factorial[n - k] % MOD;
}

int64_t binomial_coefficient(int64_t n, int64_t k) {
    if (k > n) return 0;
    if (k == 0 || k == n) return 1;

    k = min(k, n - k);
    long double res = 1.0;
    for (int64_t i = 1; i <= k; ++i) {
        res = res * (n - i + 1) / i;
    }
    return (int64_t)round(res);
}

int64_t catalan_number(int64_t n) {
    return binomial_mod(2 * n, n) * mod_pow(n + 1, MOD - 2, MOD) % MOD;
}

Go:

const maxN = 100000
const mod = 1000000007

var factorial []int64
var invFactorial []int64

func init() {
    factorial = make([]int64, maxN+1)
    factorial[0] = 1
    for i := 1; i <= maxN; i++ {
        factorial[i] = factorial[i-1] * int64(i) % mod
    }

    invFactorial = make([]int64, maxN+1)
    invFactorial[MAX_N] = modPow(factorial[MAX_N], mod-2, mod)
    for i := maxN - 1; i >= 0; i-- {
        invFactorial[i] = invFactorial[i+1] * int64(i+1) % mod
    }
}

func binomialMod(n, k int64) int64 {
    if k > n || k < 0 {
        return 0
    }
    return factorial[n] * invFactorial[k] % mod * invFactorial[n-k] % mod
}

func binomialCoefficient(n, k int64) int64 {
    if k > n {
        return 0
    }
    if k == 0 || k == n {
        return 1
    }

    k = min(k, n-k)
    res := int64(1)
    for i := int64(1); i <= k; i++ {
        res = res * (n - i + 1) / i
    }
    return res
}

func catalanNumber(n int64) int64 {
    return binomialMod(2*n, n) * modPow(n+1, mod-2, mod) % mod
}

func modPow(base, exp, mod int64) int64 {
    result := int64(1)
    base = base % mod
    for exp > 0 {
        if exp%2 == 1 {
            result = result * base % mod
        }
        base = base * base % mod
        exp /= 2
    }
    return result
}

func min(a, b int64) int64 {
    if a < b {
        return a
    }
    return b
}

Ruby:

MAX_N = 100_000
MOD = 1_000_000_007

$factorial = [0] * (MAX_N + 1)
$factorial[0] = 1
(1..MAX_N).each do |i|
  $factorial[i] = ($factorial[i - 1] * i) % MOD
end

def mod_inverse(x, mod)
  pow(x, mod - 2, mod)
end

$inv_factorial = [0] * (MAX_N + 1)
$inv_factorial[MAX_N] = mod_inverse($factorial[MAX_N], MOD)
(MAX_N - 1).downto(0).each do |i|
  $inv_factorial[i] = ($inv_factorial[i + 1] * (i + 1)) % MOD
end

def binomial_mod(n, k)
  return 0 if k > n || k < 0
  $factorial[n] * $inv_factorial[k] % MOD * $inv_factorial[n - k] % MOD
end

def binomial_coefficient(n, k)
  return 0 if k > n
  return 1 if k == 0 || k == n

  k = [k, n - k].min
  res = 1
  (1..k).each do |i|
    res = res * (n - i + 1) / i
  end
  res
end

def catalan_number(n)
  binomial_mod(2 * n, n) * mod_pow(n + 1, MOD - 2, MOD) % MOD
end

def mod_pow(base, exp, mod)
  result = 1
  base = base % mod
  while exp > 0
    result = (result * base) % mod if exp.odd?
    base = (base * base) % mod
    exp /= 2
  end
  result
end

Code Breakdown

Combinatorics functions compute permutations and combinations efficiently:

Key Variables:
- factorial: Precomputed array of factorials
- invFactorial: Precomputed array of modular inverses of factorials
- n, k: Parameters for binomial coefficient C(n,k)
- MOD: Prime modulus for operations (typically 10^9 + 7)

Algorithm Flow:

  graph TD
    Precompute["Precompute factorials"] --> Inverse["Compute inverse factorials"]
    Inverse --> Compute["Calculate C(n,k)"]
    Compute --> Formula["n! * inv(k!) * inv((n-k)!) mod MOD"]
    Formula --> Result["Return result"]
    style Result fill:#6f6,stroke:#333

Critical Sections:
- Precomputation: Calculate factorials once for all values up to MAX_N
- Modular Inverse: Compute inverse of MAX_N! factorial using Fermat’s little theorem (a^(p-2) mod p)
- Downward Pass: Fill inverse factorials from MAX_N down to 0 using relation
- Combination Formula: C(n,k) = factorial[n] * invFactorial[k] * invFactorial[n-k] mod MOD
- Edge Cases: Handle k > n, k = 0, k = n

The key insight is using precomputed modular inverses to compute C(n,k) in O(1) after O(N) preprocessing.

Performance Notes: Always use modular inverses for large combinatorics in competitive programming. Precompute factorials when C(n,k) will be called multiple times.