ABC110 D - Factorization メモ - sashimingの競プロメモ

逆元の mod のとり方を知ったので自分用にメモ。

問題はこちらから

大まかな解き方

問題を解く上での考え方は単純。 $M$ の素因数の種類の数を $k$ として、

$M = {p_1}^{b_1} \times {p_2}^{b_2} \times \cdots \times {p_k}^{b_k}$

と素因数分解されるとします。

それぞれの素因数 $p_i$ が $b_i$ 個あるので、これらを $N$ 個のグループに分配します（０個のグループがあってもよい）。この分け方の総数は ${}_{N+{b_i}-1}\mathrm{C}_{b_i}(={}_{N}\mathrm{H}_{b_i})$ です（１列に並んだ $b_i$ 個の ◯ を $N-1$ 本の区切りで分けると考えればよいです、下図参照）。

この値を $i = 1, 2, \cdots , k$ のそれぞれについて全てかけ合わせたものを 1e9+7 (= 1000000007)で割った余りが答えです。
要するに

$\displaystyle \left( \prod_{i=1}^k {}_{N+{b_i}-1}\mathrm{C}_{b_i} \right) \bmod 10^{9}+7$

が答えです。

f:id:sashiming:20190202191504p:plain:w500

問題点

${}_{N+{b_i}-1}\mathrm{C}_{b_i}$ を愚直に求めて mod をとろうとしても、値が long long でも扱えないほどクソでかくなるときが出てきます。 ${}_{10^{5}}\mathrm{C}_{10}$ の時点でムリです。これを解決するにはどうすればいいでしょうか？

解法

${\displaystyle {}_{n}\mathrm{C}_{r} = \frac{n!}{r!(n-r)!} }$ の定義に除算が入っているので、通常では計算途中に mod 計算をすることができません。しかし、次のようにすれば除算でも mod をとることが出来ます。

ここからは ${\frac{b}{a}}$ の mod をとることについて考えましょう。

${\frac{b}{a} = b \times a^{-1}}$ なので、 ${\frac{1}{a}}$ の mod をとることは $a^{-1}$ の mod をとることと同義です。（ちなみに $a^{-1}$ のことを $a$ の逆元と呼びます）

ここで、フェルマーの小定理より

$\displaystyle a^{p-1} \equiv 1 \pmod{p} \quad (a,p \,$ は互いに素 $)$

この両辺に $a^{-1}$ をかけて、

$\displaystyle a^{p-2} \equiv a^{-1} \pmod{p} \quad (a,p \,$ は互いに素 $)$

（p = 1e9+7 であり、1e9+7 は素数なので、a が p の倍数でない限り a, p は互いに素になります。今回の場合、制約上 a が p の倍数になることはないためフェルマーの小定理が適用できます）

以上より、 $b \div a$ の答えを $p$ で割った余りは、 $b \times a^{p-2}$ を $p$ で割った余りに等しいです。

${}_{n}\mathrm{C}_{r}$ の定義は次のように変形されます。

${\displaystyle {}_{n} \mathrm{C} _{r} = \frac{n!}{r!(n-r)!} = n! \times {r!}^{-1} \times {(n-r)!}^{-1} }$

剰余を求めたいので、 $a^{p-2} \equiv a^{-1} \pmod{p}$ を利用して

${\displaystyle {}_{n} \mathrm{C} _{r} \equiv n! \times {r!}^{-1} \times {(n-r)!}^{-1} \equiv n! \times {r!}^{p-2} \times {(n-r)!}^{p-2} \pmod{p} }$

コード解説

ACコードは以下の通りです。（そこそこ長いです。もっと洗練できるかも）

#include <bits/stdc++.h>
using namespace std;
#define int long long // 簡略化のため longlong を int に
const long long MOD = 1e9+7;

int fact[150000], revfact[150000]; // fact[i]はiまでの階乗、revfact[i]はiまでの階乗の逆元
queue<int> prime_count;

int powmod(int a, int p){
    int res = 1, tmp = a;
    while(p != 0){
        if(p % 2) res = (res*tmp) % MOD;
        tmp = (tmp*tmp) % MOD;
        p /= 2;
    }
    return res;
}

void factmod(){
    fact[0] = revfact[0] = 1;
    for(int i = 1; i < 150000; i++){
        fact[i] = (fact[i-1] * i) % MOD;
        revfact[i] = (revfact[i-1] * powmod(i, MOD-2)) % MOD;
    }
}

int nCrmod(int n, int r){ // nCr % MOD
    return (((fact[n] * revfact[r]) % MOD) * revfact[n-r]) % MOD;
}

signed main(){ // ここの signed を int にするとCEします
    int N, M; cin >> N >> M;

    for(int i = 2; i <= sqrt(M); i++){ // 素因数分解して指数を queue に突っ込む
        if(M % i == 0){
            int cnt = 0;
            while(M % i == 0){
                M /= i;
                cnt++;
            }
            prime_count.push(cnt);
        }
        if(M == 1) break;
    }
    if(M != 1) prime_count.push(1); // Mが素数の場合はqueueに1だけが入ることになります

    factmod();

    int ans = 1;
    while(!prime_count.empty()){
        int bi = prime_count.front();
        ans *= nCrmod(N+bi-1, bi); ans %= MOD;
        prime_count.pop();
    }
    cout << ans << endl;

    return 0;
}

一部分ずつ解説していきます。

int powmod(int a, int p){
    int res = 1, tmp = a;
    while(p != 0){
        if(p % 2)
            res = (res*tmp) % MOD;
        tmp = (tmp*tmp) % MOD;
        p /= 2;
    }
    return res;
}

$a^{p}$ (mod 1e9+7) を計算します。
今回の場合、 $p = 10^{9}+5$ なので、通常の累乗の求め方（計算量 $O(p)$ ）ではTLEしてしまいます。

そこで、繰り返し二乗法という手法を使います。

例として $187^{82}$ について考えることにしましょう。 $187^{82}$ は次のように変形されます。

$\begin{align} 187^{82} & = 187^{64} \times 187^{16} \times 187^{2} \\ & = 187^{2^{6}} \times 187^{2^{4}} \times 187^{2^{1}} \end{align}$

次に、指数である $82$ を２進数で表すと $1010010$ です。ここで、ビットが立っているところは $187^{82}$ を分解した時のそれぞれの指数に対応しています。
したがって、 $187^{82}$ は $187^{2^{1}} , 187^{2^{2}} , \cdots , 187^{2^{6}}$ をそれぞれ計算すると求められます。
（僕の説明はわかりにくいのでこちらを参考にすると良いと思います）

繰り返し二乗法を使うと、 $a^{p}$ が $O(\log{p})$ と高速に求められます。

void factmod(){
    fact[0] = revfact[0] = 1;
    for(int i = 1; i < 150000; i++){
        fact[i] = (fact[i-1] * i) % MOD;
        revfact[i] = (revfact[i-1] * powmod(i, MOD-2)) % MOD;
    }
}

fact[i] は i までの階乗 (mod 1e9+7)の配列、revfact[i] は i までの階乗の逆元 (mod 1e9+7) の配列です。制約上、150000 ぐらいまで計算しておけば十分です。十分すぎるくらいです。

そのままの階乗の計算は特に言うことはありませんが、 $a$ の階乗の逆元の計算は

$\begin{align} \displaystyle ({a!})^{-1} & = \frac{1}{a!} \\ &= \frac{1}{(a-1)!} \times \frac{1}{a} \\ &= (a-1)!^{-1} \times a^{-1} \end{align}$

で求められます。

int nCrmod(int n, int r){
    return (((fact[n] * revfact[r]) % MOD) * revfact[n-r]) % MOD;
}

nCr (mod 1e9+7) を求めています。nCr の定義の通りに計算しています。

${\displaystyle {}_{n} \mathrm{C} _{r} = n! \times {r!}^{-1} \times {(n-r)!}^{-1} }$

for(int i = 2; i <= sqrt(M); i++){
    if(M % i == 0){
        int cnt = 0;
        while(M % i == 0){
            M /= i;
            cnt++;
        }
        prime_count.push(cnt);
    }
    if(M == 1) break;
}
if(M != 1) prime_count.push(1); // Mが素数の場合はqueueに1だけが入ることになります

$M$ を素因数分解し、指数のみを queue に突っ込んでいます（指数しか計算に利用しないためです。素因数分解をライブラリ化する場合は map で素因数を key にしておくとよさそう）。

$M$ の素因数分解では、 $2$ から $\sqrt{\mathstrut M}$ までを調べれば十分です（ $M$ が合成数の場合、 $\sqrt{\mathstrut M}$ より大きい素因数は存在しないため）。 $M$ が素数かそうでないかは、for文を抜けた後の $M$ が $1$ であるかどうかで判断できます。