100兆の階乗の25兆桁めの数字は？【python】

概要

「100兆の階乗の25兆桁目の数字を求めよ(表記は10進法）」という問題をpythonで解いてみました。（元ネタは「100兆の階乗の右から数えて25兆番目にある数字は偶数であるか奇数であるか」という有名な問題です)

方針

100兆の階乗を愚直に求めると計算資源が足りませんので、ある程度効率的に計算する工夫をします。

一般的に$n!$の$k$桁めの数字を求める問題を考えます。まず、$n！$の右側に0が何個連続するかを考えます。これは$n!$を$2$または$5$で割り切れる回数をそれぞれ$c2$,$c5$とすると、$\min(c2,c5)$です。今回は明らかに$c2>c5$なので、$0$が$c5$個連続するとわかります。
したがって$n!$の$k$桁目の数字を求めるには$n!/10^{c5}$の$k-c5$桁目の数を求めればよいです。これは$k-c5$がある程度小さければ、効率的に求めることができます。

やり方としてはまず、階乗を構成する数字のうち、$5$の倍数を除く奇数を全てかけ合わせた値を$M=10^{k-c5}$を法として求めます。$M$を法としたときは、$M$ごとに値がループするので、M以下の数の積についてだけ求めて、それを$n//M$倍し、さらに中途半端に余る$n$ mod $M$について計算した積を乗じればよいです($//$は余りを切り捨てたわり算の結果を求める演算とします)。計算量は$O(M)$です。
次に、飛ばしていた偶数部分($2\times 4\times 6\times ...)$の処理について考えますが、これは各要素を2で割ると、$(n//2)!$になるので、上記の$5$の倍数を除く奇数をすべて掛け合わせた積を求める操作を再度行うことができます。この処理を再帰的に繰り返すことで、求めたかった値が得られます。また5の倍数の部分についても似たような処理でMを法とした値を求めることができます。合わせてループの回数は$O((\log n)^2)$です。

上記の計算が終わったら、計算を後回しにしていた2の処理をします。と言っても$2^{c2-c5}$を乗じるだけです。下記コードで$O((\log n)^2 + M)$の計算量で解けます。

答えが7となったのですが、あっているでしょうか？どうやって検算したらよいかわかりません。

コード

# n!のk桁めの数字を出力する
def calc(n,k):
  #n!は5で何回割り切れるか
  c5=0
  n_tmp=n
  while n_tmp>0:
    n_tmp//=5
    c5+=n_tmp
  #もしk<c5なら0を出力して終了
  if k<=c5:
    print(0)
    return
  #n!は2で何回割り切れるか
  c2=0
  n_tmp = n
  while n_tmp>0:
    n_tmp//=2
    c2+=n_tmp
  #n!//(5**c5)//(2**c2)をmod(10**(k-c5)）で求める
  #計算量がO(10**(k-c5))程度になるので、a-c5が大きすぎる（例えば>9)とこの方法は使えない。
  M=10**(k-c5)
  #1!~M!をmod M で求める
  d = [1]*M
  r=1
  for i in range(1,M):
    if i%5==0 or i%2==0:
      d[i]=d[i-1]
      continue
    r*=i
    r%=M
    d[i]=r
  d[0]=d[M-1]
  #nをループのたびに5で割りながら、5の倍数を除く奇数部分の積を求める操作を繰り返す
  r,n_tmp = 1,n
  while n_tmp>0:
    n_tmp2 = n_tmp
    while n_tmp2>0:
      if n_tmp2>=M:
        r*=pow(d[0],n_tmp2//M,M)
      if n_tmp2%M != 0:
        r*=d[n_tmp2%M]
      r%=M
      n_tmp2//=2
    n_tmp//=5
  #2を乗じる
  r*=pow(2,c2-c5,M)
  print(str(r%M+M)[1])

calc(10**14,25*(10**12))