Help us understand the problem. What is going on with this article?

フェルマーの小定理を用いたCombination(組み合わせ)計算

More than 1 year has passed since last update.

前書き

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$であるとする.

test.cpp
#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} 

となる.

main.cpp
#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

ofutonfuton
初心者の備忘録的Qiita
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
ユーザーは見つかりませんでした