大きな数の組み合わせを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)乗の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; } } ``` ### おわりに 競プロで使えそうなネタでした。 強い人は知ってるんだろうな〜と思いながら、初心者の自分は備忘録として残しておきました。

この記事をシェアする

• 
  Twitter
• 
  Facebook
• 
  Google+
• B!
  はてブ
• 
  Pocket
• 
  Feedly