PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 421 / n^15+1 の素因数

n^15+1 の形の数は n > 1 のすべての整数 n において合成数である。正の整数 n と m に対し,m を超えない n^15+1 の異なる素因数の和を s(n, m) とする。

  • 2^15+1 = 3×3×11×331
  • s(2, 10) = 3
  • s(2,1000) = 3+11+331 = 345

同様に

  • 10^15+1 = 7×11×13×211×241×2161×9091
  • s(10, 100) = 31
  • s(10,1000) = 483

1 ≤ n ≤ 10^11 における Σs(n, 10^8) を求めよ。

Problem 421 - Project Euler

そのまま書く

定義通り n^15+1 を素因数分解して,素数の和をとる関数を作りました。
s(10^10, 10^6) を求めるだけで7秒近くかかります。Σをとるまでもなく,これはダメですね。


素数を数える

s(n, m) の中身は素数の和なので,各素数が何回登場するか数えます。
たとえば 2 は n^15+1 が偶数のときに登場します。

7の場合を考えてみましょう。 n^{15} \equiv -1\pmod{7} の解は

PowerModList[-1, 1/15, 7]

で求められて  n\equiv 3, 5, 6\pmod{7} です。

 n\equiv 3\pmod{7} のとき n=7k+3 とおけて

 0\leqq k\leqq \dfrac{10^{11}-3}{7}

この形の数は  \dfrac{10^{11}-3}{7}+1 回登場します。以上を一般化して計算しました。計算時間は22分。


PowerModList の計算回数を減らす

10^11回も PowerModList を計算すると時間がかかるのは当たり前なので,この回数を減らすことを考えます。まず規則性を調べました。

f:id:variee:20170708193546p:plain

とる値は 1, 3, 5, 15 の4つだけです。これは乗法群  \mathbb{Z}/p\mathbb{Z} の位数が p-1 だからです。PowerModList の長さは

 \mathrm{GCD}(15,\, p-1)=1,\, 3,\, 5,\, 15

となります。「どんな場合にこうなるか?」「PowerModList の中身は?」を考えると長くなってしまいますが,1 の場合は単純で,PowerModList は自明な解 p-1 しかもちません。

GCD(15, p-1)=1 のときは PowerModList を作らないコードを書いたら,計算時間は15分になりました。


原始根を使う

 n^{15} \equiv -1\pmod{p} は原始根を使うと解けるそうです。

math.stackexchange.com

この方法で計算すると200秒くらいになります。1分を切るにはさらに3倍以上速くしないといけないとは……