前書き
AGC025のB問題にて,膨大なCombination計算を強いられたため,フェルマーの小定理を用いることにした.
フェルマーの小定理
$M$が素数,$k$が任意の自然数で,$M$と$k$が互いに素のとき,
k^{M-1} \equiv 1 \pmod{M}
となる.[証明はこちら]
これを変形して,
k^{M-1} = k \cdot k^{M-2} \equiv 1 \pmod{M}
ゆえに
k^{-1} \equiv k^{M-2} \pmod{M}
が導かれる.
Combination計算への応用
{}_n C _k = \cfrac{n!}{k!\cdot(n-k)!}=n!\cdot (k!)^{-1}\cdot ((n-k)!)^{-1}
と書ける.
上式($k^{-1} \equiv k^{M-2} \pmod{M}$)と組み合わせると,
{}_n C_k\equiv n!\cdot (k!)^{M-2} \cdot ((n-k)!)^{M-2} \pmod{M}
となる.以下ではこれを実装する.
実装
実際にやってみた.以下に${}_n C_k$の計算を行うコードを示す.ただし,$0< n \leqq 300000$であるとする.
#include <iostream>
#include <vector>
using namespace std;
typedef long long ll;
const ll M = 998244353;
vector<ll> fac(300001); //n!(mod M)
vector<ll> ifac(300001); //k!^{M-2} (mod M)
//a,bの範囲的にこれだけ配列を用意していけば十分
ll mpow(ll x, ll n){ //x^n(mod M) ←普通にpow(x,n)では溢れてしまうため,随時mod計算
ll ans = 1;
while(n != 0){
if(n&1) ans = ans*x % M;
x = x*x % M;
n = n >> 1;
}
return ans;
}
ll comb(ll a, ll b){ //aCbをmod計算
if(a == 0 && b == 0)return 1;
if(a < b || a < 0)return 0;
ll tmp = ifac[a-b]* ifac[b] % M;
return tmp * fac[a] % M;
}
int main()
{
ll n,k;
cin >> n >> k;
//大した量ではないので,先にfax[i]とifax[i]を全て計算しておく
fac[0] = 1;
ifac[0] = 1;
for(ll i = 0; i<300000; i++){
fac[i+1] = fac[i]*(i+1) % M; // n!(mod M)
ifac[i+1] = ifac[i]*mpow(i+1, M-2) % M; // k!^{M-2} (mod M) ←累乗にmpowを採用
}
ll ans=0;
ans = comb(n, k)%M;
cout << ans << endl;
return 0;
}
応用
実際の問題に応用してみた.問題はyukicoder No.117を用いる.当然フェルマーの小定理は順列や重複組み合わせにも応用可能.
順列計算への応用
{}_n P _k = \cfrac{n!}{(n-k)!}=n!\cdot ((n-k)!)^{-1}
と書ける.
上式($k^{-1} \equiv k^{M-2} \pmod{M}$)と組み合わせると,
{}_n P_k\equiv n!\cdot ((n-k)!)^{M-2} \pmod{M}
となる.
重複組み合わせ計算への応用
{}_n H _k = \cfrac{(n+k-1)!}{k!\cdot(n-1)!}=(n+k-1)!\cdot (k!)^{-1}\cdot ((n-1)!)^{-1}
と書ける.
上式($k^{-1} \equiv k^{M-2} \pmod{M}$)と組み合わせると,
{}_n H_k\equiv (n+k-1)!\cdot (k!)^{M-2} \cdot ((n-1)!)^{M-2} \pmod{M}
となる.
#include <iostream>
#include <vector>
#include <cmath>
#define FOR(i,a,b) for(int i=(a);i<(b);++i)
using namespace std;
typedef long long ll;
const ll M = pow(10,9)+7;
vector<ll> fac(1000001); //n!(mod M)
vector<ll> ifac(1000001); //k!^{M-2} (mod M)
ll mpow(ll x, ll n){ //x^n(mod M)
ll ans = 1;
while(n != 0){
if(n&1) ans = ans*x % M;
x = x*x % M;
n = n >> 1;
}
return ans;
}
ll comb(ll a, ll b){ //aCb(mod M)
if(a == 0 && b == 0)return 1;
if(a < b || a < 0)return 0;
ll tmp = ifac[a-b]* ifac[b] % M;
return tmp * fac[a] % M;
}
ll perm(ll a, ll b){
if(a < b) return 0;
ll tmp = ifac[a-b] % M;
return tmp * fac[a] % M;
}
ll dupc(ll a,ll b){
if(a == 0 && b == 0)return 1;
if(a < 0 || b == 0)return 0;
ll tmp = ifac[a-1]* ifac[b] % M;
return tmp * fac[a+b-1] % M;
}
int main()
{
ll t;
cin >> t;
vector<char> c(t);
vector<int> a(t),b(t);
FOR(i, 0, t){
scanf(" %c(%d,%d)",&c[i] ,&a[i] ,&b[i]);
}
fac[0] = 1;
ifac[0] = 1;
for(ll i = 0; i<1000000; i++){
fac[i+1] = fac[i]*(i+1) % M; // n!(mod M)
ifac[i+1] = ifac[i]*mpow(i+1, M-2) % M; // k!^{M-2} (mod M)
}
FOR(i, 0, t){
ll ans=0;
if(c[i] == 'P')ans = perm(a[i],b[i]);
if(c[i] == 'C')ans = comb(a[i],b[i]);
if(c[i] == 'H')ans = dupc(a[i],b[i]);
cout << ans << endl;
}
return 0;
}
あとがき
とりあえず,膨大なCombination(組み合わせ)計算を強いられたときは,素数$M$についてmodをとるフェルマーの小定理を用いて計算すべしということが分かった.
参考文献
マツシタのお勉強:http://keita-matsushita.hatenablog.com/entry/2016/12/05/184011
高校数学の美しい物語:https://mathtrain.jp/fermat_petit