二項係数の剰余(modcomb)
目次
1. nCk n~107程度
目的
二項係数nCkを素数pで割った際の剰余nCk mod pを求める。
制約
p:素数
1≦k≦n≦107
n<p
オーダー
前処理:O(N) 計算:O(1)
コード
const int MAX = 510000; const int MOD = 1000000007; long long fac[MAX], finv[MAX], inv[MAX]; // テーブルを作る前処理 void modcombinit() { fac[0] = fac[1] = 1; finv[0] = finv[1] = 1; inv[1] = 1; for (int i = 2; i < MAX; i++){ fac[i] = fac[i - 1] * i % MOD; inv[i] = MOD - inv[MOD%i] * (MOD / i) % MOD; finv[i] = finv[i - 1] * inv[i] % MOD; } } // 二項係数計算 long long modcomb(long long n, long long k){ if (n < k) return 0; if (n < 0 || k < 0) return 0; return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD; }
コードテスト
#include <bits/stdc++.h> using namespace std; const int MAX = 510000; const long long MOD = pow(10, 9) + 7; long long fac[MAX], finv[MAX], inv[MAX]; void modcombinit() { // テーブルを作る前処理 fac[0] = fac[1] = 1; finv[0] = finv[1] = 1; inv[1] = 1; for (int i = 2; i < MAX; i++) { fac[i] = fac[i - 1] * i % MOD; inv[i] = MOD - inv[MOD % i] * (MOD / i) % MOD; finv[i] = finv[i - 1] * inv[i] % MOD; } } long long modcomb(long long n, long long k) { //二項係数計算 if (n < k) return 0; if (n < 0 || k < 0) return 0; return fac[n] * (finv[k] * finv[n - k] % MOD) % MOD; } int main() { ios::sync_with_stdio(false); cin.tie(0); long long A,B; cin >> A>>B; modcombinit(); cout << modcomb(A, B) << endl; }
1000 2 499500
出典
2. n固定,k~105程度
3. n~106程度 ACL使用ver
目的
ACLが出たのでrefine
注意点
MAX=107は手元環境で5秒程度でしたが、AtCoder上のコードテストでは1秒程度でした。 適宜調節してください。
コード
using mint = static_modint<1000000007>; const ll MAX = 1000100; //10^6 mint fac[MAX], finv[MAX], inv[MAX]; void COMinit() { fac[0] = fac[1] = 1; finv[0] = finv[1] = 1; inv[1] = 1; for (int i = 2; i < MAX; i++) { fac[i] = fac[i - 1] * i; inv[i] = ((mint)i).inv(); finv[i] = finv[i - 1] * inv[i] ; } } ll COM(int n, int k) { mint tmp; ll ans; if (n < k) return 0; if (n < 0 || k < 0) return 0; tmp = finv[k] * finv[n - k]; tmp *= fac[n]; ans = (ll)(tmp.val()); return ans; }
コードテスト
#include <atcoder/all> #include <bits/stdc++.h> using namespace std; using namespace atcoder; typedef long long int ll; using mint = static_modint<1000000007>; // using mint=static_modint<998244353>; const ll MAX = 1000100; //10^6 mint fac[MAX], finv[MAX], inv[MAX]; void COMinit() { fac[0] = fac[1] = 1; finv[0] = finv[1] = 1; inv[1] = 1; for (int i = 2; i < MAX; i++) { fac[i] = fac[i - 1] * i; inv[i] = ((mint)i).inv(); finv[i] = finv[i - 1] * inv[i] ; } } ll COM(int n, int k) { mint tmp; ll ans; if (n < k) return 0; if (n < 0 || k < 0) return 0; tmp = finv[k] * finv[n - k]; tmp *= fac[n]; ans = (ll)(tmp.val()); return ans; } int main() { COMinit(); // 計算例 cout << COM(1000, 500) << endl; }
159835829