二項係数 (nCr) の計算方法
AtCoder Grand Contest 025 にて を使う問題が出題されました。いい機会だと思ったので の高速な実装をしてみました。
使用言語: C++
1. なぜ逆元 が必要なのか
なので、 が分かれば で が求まります。よって事前に階乗を計算しておけばよいことが分かります。
しかし残念ながら
とは限りません。
例: のとき
です。一方、
なので は計算できません。
よって実際には も事前に計算しておく必要があります。
2. 階乗 と 逆元 の計算
であることを使えば は全体で で求まります。
一方、
であるので が分かれば は全体で で求まります。しかし を求めるのは少々やっかいです。まず思いつくのは逆元の定義
を利用して愚直に全探索する方法ですが、計算量は になってしまいます。実は次に紹介するフェルマーの小定理と繰り返し二乗法を利用すると を で求めることができます。
3. フェルマーの小定理
フェルマーの小定理は を素数としたとき の倍数でない任意の整数 に対して
が成立するというものです。簡単に証明を書いておきます。
(証明)
と は互いに素なので
となる は存在しない。よって はすべて異なり のいずれかと同じになる。それぞれの要素の積をとると
となるので定理が従う。
(証明終)
この定理を利用すると
であることが分かります。
4. 繰り返し二乗法
繰り返し二乗法はその名の通り二乗を繰り返して を計算するアルゴリズムです。 の二進法表現を
とすると、
と表せます。
を使って をそれぞれ求めることができるので、計算量は となります。
例: のとき
繰り返し二乗法を利用すると は で求まります。
5. 実装
以上より が が で前計算できることが分かりました。よって は で求まります。
#include <bits/stdc++.h> using namespace std; typedef long long ll; #define REP(i, n) for(int i=0; i<(n); ++i) #define FOR(i, a, b) for(int i=(a); i<(b); ++i) #define FORR(i, a, b) for(int i=(b)-1; i>=(a); --i) constexpr ll MOD=1000000007ll; // constexpr ll MOD=998244353ll; constexpr int MAX=510000; ll fact[MAX], fact_inv[MAX]; // 繰り返し二乗法 ll power(ll a, ll b){ ll res=1; while(b>0){ if(b&1) res=res*a%MOD; a=a*a%MOD; b>>=1; } return res; } ll comb(ll n, ll r){ return (fact[n]*fact_inv[r])%MOD*fact_inv[n-r]%MOD; } int main(){ ios::sync_with_stdio(false); cin.tie(0); int n=100000; // cin>>n; fact[0]=1; // 階乗の計算 REP(i, n) fact[i+1]=fact[i]*(i+1)%MOD; fact_inv[n]=power(fact[n], MOD-2); // 逆元の計算 FORR(i, 0, n) fact_inv[i]=fact_inv[i+1]*(i+1)%MOD; // nCrの計算 cout<<"6C3: "<<comb(6, 3)<<'\n'; cout<<"666C390: "<<comb(666, 390)<<'\n'; return 0; }