AtCoder Beginner Contest 110: D - Factorizatiom

競技プログラミング AtCoder 整数問題 素数 重複組合せ 二項係数 拡張ユークリッドの互除法 逆元

問題

問題文

https://atcoder.jp/contests/abc110/tasks/abc110_d

問題概要

正整数 N, M が与えられる.

a_1 \times a_2 \times ... \times a_N = M となるような正整数からなる長さ N の数列 a が何通りあるかを 109+7 で割った剰余の形で答えよ.

制約

  • 1 \leqq N \leqq 10^{5}

  • 1 \leqq M \leqq 10^{9}

解答例

指針

  • 重複組合せ

解説

まず M がある素数 p の冪乗の形で表される場合を考える.

a_1 \times a_2 \times ... \times a_N = p^{k}

ここで, a_i = p^{c_i} とおくと, 以下のように書ける.

p^{c_1} \times p^{c_2} \times ... \times p^{c_N} = p^{k}

したがって,

c_1 + c_2 + ... + c_N = k

という方程式を満たす非負整数で構成される数列 c が何通りかを考えればいいことが分かる.

上記の整数方程式の解の個数は重複組合せの問題として解くことができる.

区別しないボールが k 個, 仕切りが N-1 個あると仮定すると非負整数解の個数は, _{N+k-1} C_k となる.

例えば, N = 3, k = 4 のとき, ボールを O しきりを | で表すとすると,

OOOO|| を並び替えることで方程式を表現できる.

ex.)
OOOO|| => 4 + 0 + 0
OO|O|O => 2 + 1 + 1
|O|OOO => 0 + 1 + 3

次に M が素数 p の冪乗の形ではない場合を考える.

a_1 \times a_2 \times ... \times a_N = p_1^{e_1} \times p_2^{e_2} \times ... \times p_k^{e_k}

この場合は,

a_1 \times a_2 \times ... \times a_N = p_1^{e_1}

a_1 \times a_2 \times ... \times a_N = p_2^{e_2}

...

a_1 \times a_2 \times ... \times a_N = p_k^{e_k}

それぞれの場合の掛け算で求めることができる.

実装

素因数分解

整数 N の素因数分解は小さい素数から sqrt(N) まで試し割りしていけばいい.

小さい方から割っていくと合成数で割る心配はない.

例えば 2 で試し割りしたあとは, 4 や 6 では割り切れなくなる.

  • C++ による実装例
vector<pair<ll, ll> > prime_factorize(ll n) {
    vector<pair<ll, ll> > ret;
    // 小さい素数から順に試し割りしていく
    for (ll i = 2; i*i <= n; i++) {
        if (n % i != 0) { continue; }
        int e = 0;
        while (n % i == 0) {
            e++;
            n /= i;
        }
        ret.push_back(make_pair(i, e));
    }
    // 2 ~ sqrt(N) で試し割りして残ったものは必ず素数
    if (n != 1) {
        ret.push_back(make_pair(n, 1));
    }
    return ret;
}
  • Python3 による実装例
import math

def prime_factorize(n):
    ret = []
    # 小さい素数から順に試し割りしていく
    for i in range(2, int(math.sqrt(n))+1):
        if n % i != 0:
            continue
        e = 0
        while n % i == 0:
            e += 1
            n //= i
        ret += [(i, e)]

    if n != 1:
        ret += [(n, 1)]

    return ret

nCr mod p の計算方法

愚直な方針

_{n} C_r = \frac{n!}{(n-r)!r!} であるからそのまま計算したくなるが, その場合算術オーバーフローとなる.

パスカルの三角形

1 \leqq r \leqq n \leqq 2000 程度であれば, パスカルの三角形で二項係数 _{n} C_r を求めることができる. 時間計算量は O(N^{2})

しかし今回はこの方法だと TLE となる.

  • C++ による実装例
const long long MOD = 1000000007;
const int MAX_V = 1000;
long long table[MAX_V][MAX_V];

void cal() {
    table[0][0] = 1;
    for (int i = 1; i < MAX_V; i++) {
        table[i][0] = 1;
        for (int j = 1; j < MAX_V; j++) {
            table[i][j] = (table[i-1][j-1] + table[i-1][j]) % MOD;
        }
    }
}
  • Python3 による実装例
def calc(MAX_V, MOD):
    table = [ [0]*MAX_V for _ in range(MAX_V) ]
    table[0][0] = 1
    for i in range(1, MAX_V):
        table[i][0] = 1
        for j in range(1, MAX_V):
            table[i][j] = (table[i-1][j-1] + table[i-1][j]) % MOD

    return table

逆元を利用した計算

mod の世界では逆元 (モジュラ逆数) というものがあり, 割り算を掛け算として処理することができる. 掛け算として処理することができれば計算途中で剰余をとっても最終的な結果に影響しないので嬉しい. 逆元の計算は RSA暗号の p, q から秘密鍵 d を求めるときに使うことでお馴染みの拡張ユークリッドの互除法を使う.

フェルマーの小定理を利用して逆元を計算することもできるらしい.

この方法は mod が素数のときしか使えない.

時間計算量は O(nlogn)

  • C++ による実装例
typedef long long ll;

ll extgcd(ll a, ll b, ll &x, ll &y) {
    ll d = a;
    if (b != 0) {
        d = extgcd(b, a%b, y, x);
        y -= (a/b) * x;
    } else {
        x = 1; y = 0;
    }
    return d;
}

ll modinv(ll a, ll m) {
    ll x, y;
    extgcd(a, m, x, y);
    return (m + x % m) % m;
}

ll mod_comb(ll n, ll r, ll mod) {
    ll ans_mul = 1, ans_div = 1;

    for (int i = 0; i < r; i++) {
        ans_mul *= (n-i);
        ans_div *= (i+1);
        ans_mul %= mod;
        ans_div %= mod;
    }
    return ans_mul * modinv(ans_div, mod) % mod;
}
  • Python3 による実装例
def extgcd(a, b):
    x,y, u,v = 0,1, 1,0
    while a != 0:
        q, r = b//a, b%a
        m, n = x-u*q, y-v*q
        b,a, x,y, u,v = a,r, u,v, m,n
        g = b
    return x, y, g


def modinv(a, m):
    x, y, g = extgcd(a, m)
    if g != 1:
        print ("[+]Inverse does not exist.")
    else:
        return ((m + x) % m) % m


def calc_comb(n, r, mod):
    ans_mul, ans_div = 1, 1

    for i in range(r):
        ans_mul *= (n-i)
        ans_div *= (i+1)
        ans_mul %= mod
        ans_div %= mod

    return ans_mul * modinv(ans_div, mod) % mod

以上より, この問題を解くことができた.

  • C++ による実装例
// submission: https://atcoder.jp/contests/abc110/submissions/3904238
// Language: C++14 (GCC 5.4.1)
// Time: 1 ms
// Memoly: 256 KB

#include <iostream>
#include <vector>

using namespace std;

typedef long long ll;

vector<pair<ll, ll> > prime_factorize(ll n) {
    vector<pair<ll, ll> > ret;
    // 小さい素数から順に試し割りしていく
    for (ll i = 2; i*i <= n; i++) {
        if (n % i != 0) { continue; }
        int e = 0;
        while (n % i == 0) {
            e++;
            n /= i;
        }
        ret.push_back(make_pair(i, e));
    }
    // 2 ~ sqrt(N) で試し割りして残ったものは必ず素数
    if (n != 1) {
        ret.push_back(make_pair(n, 1));
    }
    return ret;
}

ll extgcd(ll a, ll b, ll &x, ll &y) {
    ll d = a;
    if (b != 0) {
        d = extgcd(b, a%b, y, x);
        y -= (a/b) * x;
    } else {
        x = 1; y = 0;
    }
    return d;
}

ll modinv(ll a, ll m) {
    ll x, y;
    extgcd(a, m, x, y);
    return (m + x % m) % m;
}

ll mod_comb(ll n, ll r, ll mod) {
    ll ans_mul = 1, ans_div = 1;

    for (int i = 0; i < r; i++) {
        ans_mul *= (n-i);
        ans_div *= (i+1);
        ans_mul %= mod;
        ans_div %= mod;
    }
    return ans_mul * modinv(ans_div, mod) % mod;
}

int main() {
    ll N, M;
    const ll MOD = ll(1e9+7);

    cin >> N >> M;

    vector<pair<ll, ll> > v = prime_factorize(M);

    ll ans = 1;

    for (int i = 0; i < v.size(); i++) {
        ans *= mod_comb(N+v[i].second-1, v[i].second, MOD);
        ans %= MOD;
    }
    cout << ans << endl;

    return 0;
}
  • Python3 による実装例
# submission: https://atcoder.jp/contests/abc110/submissions/3904263
# Language: Python3 (3.4.3)
# Time: 22 ms
# Memoly: 3188 KB

import math

def prime_factorize(n):
    ret = []
    # 小さい素数から順に試し割りしていく
    for i in range(2, int(math.sqrt(n))+1):
        if n % i != 0:
            continue
        e = 0
        while n % i == 0:
            e += 1
            n //= i
        ret += [(i, e)]

    if n != 1:
        ret += [(n, 1)]

    return ret


def extgcd(a, b):
    x,y, u,v = 0,1, 1,0
    while a != 0:
        q, r = b//a, b%a
        m, n = x-u*q, y-v*q
        b,a, x,y, u,v = a,r, u,v, m,n
        g = b
    return x, y, g


def modinv(a, m):
    x, y, g = extgcd(a, m)
    if g != 1:
        print ("[+]Inverse does not exist.")
    else:
        return ((m + x) % m) % m


def mod_comb(n, r, mod):
    ans_mul, ans_div = 1, 1

    for i in range(r):
        ans_mul *= (n-i)
        ans_div *= (i+1)
        ans_mul %= mod
        ans_div %= mod

    return ans_mul * modinv(ans_div, mod) % mod


def main():
    N, M = map(int, input().split())
    v = prime_factorize(M)
    MOD = int(10**9+7)
    ans = 1

    for i in v:
        ans *= mod_comb(N+i[1]-1, i[1], MOD)
        ans %= MOD
    print (ans)

main()

参考文献