Combinatorics

C(n,k) = n! / (k! * (n-k)!)

template<typename T>
T C(T n, T k) { // O(k)
    assert(0 <= n && 0 <= k);
    if(n < k) return static_cast<T>(0);
    T ans = static_cast<T>(1);
    for(T i = static_cast<T>(1); i <= k; i++) {
        ans *= (n - k + i);
        ans /= i;
    }
    return ans;
}

int ways = C(5, 2); // Returns 10
// Choosing 2 items from 5: (5*4)/(2*1) = 10

// Requires: numeric_mint_compress or numeric_mint_full
// Requires: math_factorial

template<typename T>
struct Combinatorics {
    int sz;
    Factorial<T> fact;
    
    Combinatorics(int n) : sz(n), fact(n) {}
    
    T C(int n, int k) { // O(1) with preprocessing O(n)
        if(n < k) return static_cast<T>(0);
        return fact[n] / (fact[n-k]*fact[k]);
    }

    T H(int n, int k) {
        // Combinations with repetition
        return C(n+k-1, k);
    }
    
    T P(int n, int k) {
        // Permutations without repetition
        if(n < k) return static_cast<T>(0);
        return fact[n] / fact[n-k];
    }
    
    T P(int n, const vector<int> &m) {
        // Permutations with repetition
        T product = static_cast<T>(1);
        for(int i = 0; i < (int) m.size(); ++i) {
            product *= fact[m[i]];
        }
        return fact[n] / product;
    }

    T catalan(int n) {
        // Catalan numbers
        return C(2*n, n) / static_cast<T>(n+1);
    }
};

template<typename T>
using Comb = Combinatorics<T>;

Comb<Mint> comb(100000);
Mint ways = comb.C(100, 50); // Binomial coefficient
Mint perm = comb.P(10, 3);   // Permutations

const int64_t md = 1e9 + 7;

int64_t fastpow(int64_t a, int64_t b, const int64_t &mod) {
    if (b == 0) return 1;
    int64_t half = fastpow(a, b / 2, mod);
    int64_t answer = half * half % mod;
    if (b & 1) answer = answer * a % mod;
    return answer;
}

int64_t inv(int64_t a, const int64_t &mod) { 
    return fastpow(a, mod - 2, mod); 
}

const int MXF = 1e6;
int64_t fact[MXF];

void build_factorial() {
    fact[0] = 1;
    for (int i = 1; i < MXF; ++i) {
        fact[i] = fact[i - 1] * i % md;
    }
}

int64_t nCr(int n, int r) {
    if (n < r) return 0;
    return fact[n] * inv(fact[n - r] * fact[r] % md, md) % md;
}

// At main():
// build_factorial();

build_factorial();
int64_t result = nCr(100, 50) % md;

P(n,k) = n! / (n-k)!
     = n × (n-1) × (n-2) × ... × (n-k+1)

C(n, k) = C(n, n-k)         // Symmetry
C(n, 0) = C(n, n) = 1        // Edge cases
C(n, k) = C(n-1, k-1) + C(n-1, k)  // Pascal's identity
Σ C(n, k) = 2^n              // Sum over all k from 0 to n

T catalan(int n) {
    return C(2*n, n) / static_cast<T>(n+1);
}

vector<int> kth_permutation(vector<int> perm, int k) {
    int64_t factorial = 1LL;
    int n = (int) perm.size();
    for(int64_t num = 2; num < n; ++num)
        factorial *= num; // (n-1)!
    k--; // k-th to 0-indexed
    vector<int> answer; 
    answer.reserve(n);
    while(true) {
        answer.push_back(perm[k / factorial]);
        perm.erase(perm.begin()+(k/factorial));
        if((int) perm.size() == 0) break;
        k %= factorial;
        factorial /= (int) perm.size();
    }
    return answer;
}

vector<int> kth_permutation(int n, int k, int start=0) {
    vector<int> perm(n);
    iota(perm.begin(), perm.end(), start);
    return kth_permutation(perm, k);
}

vector<int> perm = kth_permutation(4, 10);
// 4! = 24 total permutations
// Returns the 10th permutation of {0,1,2,3}

// Works for 1 <= k <= n <= 20 approximately
// Complexity: worst case O(2^n) approximately
bool next_combination(vector<int> &comb, int n) {
    int k = (int) comb.size();
    for (int i = k - 1; i >= 0; i--) {
        if (comb[i] <= n - k + i) {
            ++comb[i];
            while (++i < k) {
                comb[i] = comb[i - 1] + 1;
            }
            return true;
        }
    }
    return false;
}

void all_combinations(int n, int k) {
    vector<int> comb(k);
    iota(comb.begin(), comb.end(), 1);
    do {
        for (const int &v : comb) {
            cout << v << " ";
        }
        cout << endl;
    } while (next_combination(comb, n));
}

// Generate all C(5,3) = 10 combinations
all_combinations(5, 3);
// Output:
// 1 2 3
// 1 2 4
// 1 2 5
// 1 3 4
// ...

vector<int> v = {1, 2, 3};
do {
    // Process permutation
    for (int x : v) cout << x << " ";
    cout << endl;
} while (next_permutation(v.begin(), v.end()));

// Generates all 3! = 6 permutations

// Select k people from n candidates
int count_teams(int n, int k) {
    Comb<Mint> comb(n);
    return comb.C(n, k)();
}

// With constraints: at least 1 from group A and 1 from group B
int count_mixed_teams(int groupA, int groupB, int k) {
    Comb<Mint> comb(groupA + groupB);
    Mint total = comb.C(groupA + groupB, k);
    Mint onlyA = comb.C(groupA, k);
    Mint onlyB = comb.C(groupB, k);
    return (total - onlyA - onlyB)();
}

// Count paths from (0,0) to (n,m) moving only right and down
int count_paths(int n, int m) {
    // Need to make n+m moves: n down, m right
    return C(n + m, n);
}

// With obstacle at (x,y)
int count_paths_obstacle(int n, int m, int x, int y) {
    int total = C(n+m, n);
    int through_obstacle = C(x+y, x) * C((n-x)+(m-y), n-x);
    return total - through_obstacle;
}

// Distribute k identical items into n bins (stars and bars)
int distribute_items(int n, int k) {
    Comb<Mint> comb(n + k);
    // Using H(n,k) = C(n+k-1, k)
    return comb.H(n, k)();
}

// Count permutations where no element is in its original position
Mint derangements(int n) {
    vector<Mint> D(n+1);
    D[0] = 1;
    D[1] = 0;
    for (int i = 2; i <= n; i++) {
        D[i] = (i-1) * (D[i-1] + D[i-2]);
    }
    return D[n];
}
// D(n) ≈ n!/e for large n

Total = All - Bad1 - Bad2 + Bad1∩Bad2