大きな整数の階乗の剰余を求める

math.factorialを使う

pythonでは整数Nの階乗の剰余は以下のようにmath.factorialを使えば簡単に求まりますが、Nが大きくなると非常に大きな数を一旦生成するのであまり効率が良くありません。(計算量は$O(\frac{1}{2}n^2)$)

from math import factorial
M = 10**9+7
N = 10**6
print(N, factorial(N)%M)
# 1000000 641102369

Nの階乗の素因数分解を求める

そこでまずNの階乗の素因数分解を行います。考え方は、例えば10の階乗の場合、2の倍数が5個、$2^2=4$の倍数が2個、$2^3=8$の倍数が1個、なので素因数2は計8個あることが分かります。これをN以下のすべての素数について行うとNの階乗の素因数分解が求まります。

from sympy import primerange
def fact2fct(n):   # nの階乗の素因数分解
  ret = dict()
  for p in primerange(2,n+1):
    ret[p], pp = 0, p
    while pp <= n:
      ret[p] += n // pp
      pp *= p
  return ret

N = 10
print(N, fact2fct(N))
#10 {2: 8, 3: 4, 5: 2, 7: 1}

各素因数に対してpow関数で剰余計算を行う

素因数分解の形になればpow関数が使えます。(計算量はほぼ$O(n)$)

def factmod(fct,M):
  ret = 1
  for p in fct.keys():
    ret *= pow(p,fct[p], M)
    ret %= M
  return ret

M = 10**9+7
n = 10**6
stime = time.time()
print(n, factmod(fact2fct(n), M))
#1000000 641102369

(応用) １からNまでの数の最小公倍数の素因数分解

この考え方を使うと１からNまでの数の最小公倍数の素因数分解も同様に実装できます。

from sympy import primerange
def lcmall(n):   # 1からnのすべての数の最小公倍数の素因数分解
  ret = dict()
  for p in primerange(2,n+1):
    pp, i = p, 0
    while pp <= n:
      pp *= p
      i += 1
    ret[p] = i
N = 10
print(N, lcmall(N))
# 10 {2: 3, 3: 2, 5: 1, 7: 1}

（開発環境：Google Colab）

この考え方は以下のProject Eulerの問題を解くのに役に立ちます。

• Problem 429: Sum of squares of unitary divisors.
• Problem 772: Balanceable k-bounded partitions.