pythonではべき乗の剰余演算はpow関数で高速に行えますが。二項係数を求めるsympyのbinomialにはmoduloの指定は出来ないようです。
今回はそのbinomialの剰余演算を簡単に高速化する方法を見つけたので報告します。
Mod(binomial(N, R, evaluate=False), p)
調べてみるとSumpy 1.14.0 documentationのbinomialの項目に以下のような説明がありました。
In many cases, we can also compute binomial coefficients modulo a prime p quickly using Lucas’ Theorem [R217], though we need to include evaluate=False to postpone evaluation:
from sympy import Mod
Mod(binomial(156675, 4433, evaluate=False), 10**5 + 3)
# 28625
Lucasの定理を使っているということなので、以前に二項係数を素数で割った余りを求める（Lucasの定理、中国剰余定理で作った関数lucasと計算時間を比較してみました。
条件は、
- Mod(binomial(N, R, evaluate=False), p) で、
- $N, R = 10^{18}, 10^5 $
- $p$は$10^x \ (x=1,2,...,16)$以上の答えが０にならない素数
比較の結果は以下のグラフのように驚くほど一致しました。
結果的には$p \lt R$であればとても効率的だと思います。
Mod(factorial(N, evaluate=False), p)
同様に階乗を求めるfactorialでも計ってみました。比較対象は大きな整数の階乗の剰余を求めるで作った素因数分解を使ったfactmod(fact2fct(N),p) です
今度は単純に$N$を変化させました。
- $N = 10^x (x =1,2,...9)$
- $p = 10^{11}+3$
factorialでも結果は非常に近いものになりましたので、似たようなアルゴリズムを使っているのかなと思います。
今回はbinomial と factorial を評価しました。今後他の関数でも効果があれば追加したいと思います。
（開発環境：Google Colab）