大きな数の組み合わせを1000000007で割った余りを計算する方法-ぐぐりらにっき

こんにちは、ぐぐりら(<a href="https://twitter.com/guglilac">@guglilac</a>)です。

大きな数の組み合わせの$10^9+7$で割ったあまりを求めたくなることってありますよね(唐突)
$M=10^9+7$とします

$x={}_n C _r \ mod\  M$

の、$n$とか$r$が大きい場合。

${}_n C _r=\frac{n!}{r!(n-r)!}$

であり、$n$とか$r$が大きい数の場合にはこれを計算したあとに$M$で割ったあまりを求めようとしてもオーバーフローしてしまうためにそもそも組み合わせ自体が計算できません。

でもどうにかして$x$を求めたい。

そういう場面で使える方法についての話。

## とりあえずやってみる

普通に分母がなくて、なにか大きい数の階乗のmodを計算したい場合は簡単で、

```c
int a=1;
for (size_t i = 1; i <= N; i++) {
        a = (a * i) % M;
    }
```
                         
みたいに一個ずつかけたら$M$で割ることを繰り返せばオーバーフローしないで済みます。
これができるのは、合同式が積については計算可能であることが背景にあります。

つまり
mod $M$に対して
$a \equiv b$、$c \equiv d$であれば$ac \equiv bd$が成り立つことが裏にあります。
これにより、全て掛け合わせたあとにあまりを出さなくても、階乗の途中でmodを計算してもいい、ということになります。

今回難しいのは、組み合わせの式には分母にも数字が来ることです。

分母の数字も同様にmodを計算して、その数でさっき求めた分子の合同式での値を割ればいいのでは、と思われるかもしれません。

ですが、これは誤りです。合同式は割り算に対しては上手く計算できません。

## 割り算を掛け算に変換

合同式は割り算が苦手なので、掛け算に直してしまう方法が解決策となります。

$M=10^9+7$は素数なので、フェルマーの小定理を用いることでこれが実現できます。

>フェルマーの小定理
>
> $p$を素数とすると、以下の合同式がなりたつ。
>
> $a^{p-1}\equiv 1\  mod\  p$

これを変形すると
$a^{p-2}\equiv \frac{1}{a}\  mod\  p$
となり、合同式では$a$で割る操作を$a^{p-2}$をかける操作に変換することができることがわかりました。

(もちろん$p$が素数である条件が必要ですが。)

これを用いると、今回の問題では$r!$や$(n-r)!$を合同式上では掛け算に変換することができます。

一見すると、
$(r!)^{M-2}\  mod\ M$などを求めて、分子の方の合同式とかけたくなりますが、$r!$が大きいのでここもばらしましょう。

$r!$を分解します。
$r!=\prod_{i=1}^{r}i$
この$i$について、$(i)^{M-2}\  mod\ M$
を求めて、後で掛け合せればいいです。割り算を掛け算に直せればあとは余りをその都度計算していけばいいので、分割は細かくやっても大丈夫です。

## $(i)^{M-2}\  mod\ M$を計算する部分

一般に、$a^{b}\ mod\ p$を計算することを考えます。

普通に$a$をかけて$p$で割った余りで更新する、という処理を$b$回行うことでももちろん計算できますが、
これは$O(b)$になります。

$O(\log b)$のアルゴリズムがあります。

$a^b$の$b$が偶数だった場合、$a^{b/2}\ mod\ p$を計算して2乗することで求めることができます。
$a^{b/2}\ mod\ p$の部分が再帰になっています。

$a^b$の$b$が奇数だった場合、$a^{b-1}\ mod\ p$を計算し、それに$a$をかけてmodをとります。
$b$が奇数なので$a^{b-1}\ mod\ p$の呼び出しはbが偶数となって再帰呼び出しになります。

終了条件はb=0のときでreturn 1にすればおっけい。

二回に一回は$b$が偶数になり、$b$が偶数の場合は$b$が半分になるので、計算量は$O(\log b)$です。

再帰関数calcをサンプルとして載せておきます。

```c
typedef long long ll;

ll calc(ll a, ll b, ll p) {
    if (b == 0) {
        return 1;
    } else if (b % 2 == 0) {
        ll d = calc(a, b / 2, p);
        return (d * d) % p;
    } else {
        return (a * calc(a, b - 1, p)) % p;
    }
}

```

## おわりに
競プロで使えそうなネタでした。

強い人は知ってるんだろうな〜と思いながら、初心者の自分は備忘録として残しておきました。

大きな数の組み合わせを1000000007で割った余りを計算する方法

0 件のコメント:

コメントを投稿

PROFILE

このブログを検索

POPULAR POSTS

SPONSORED LINK

RECENT POSTS

TWEETS

ARCHIVES

LABEL

CONTACT