読者です 読者をやめる 読者になる 読者になる

PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 407 / 冪等元

0 ≤ a ≤ 5 のときの a に対し a^2 mod 6 を計算すると 0,1,4,3,4,1 となる。

a2 ≡ a (mod 6) をみたす最大の a は 4 となる。
a2 ≡ a (mod n) をみたす a < n の最大値を M(n) としよう。M(6) = 4 である。

1 ≤ n ≤ 10^7 のときの ΣM(n) を求めよ.

Problem 407 - Project Euler


a^2-a = a(a-1) において a と a-1 が互いに素なことがポイントです。

第1案:中国式剰余定理

■n = p(素数)のとき

p | a(a-1) より
p | a または p | a-1

a< n=p より a=1 です。同様に n=p^k のときも a=1 が言えます。

■n = pq (p, q は異なる素数)のとき

pq | a(a-1) において a< n=pq から pq | a も pq | a-1 も基本ありません(例外は a=1 のとき)。

「p | a かつ q | a-1」または「p | a-1 かつ q | a」

これを合同式で表します。

「a ≡ 0 (mod p) かつ a ≡ 1 (mod q)」または
「a ≡ 1 (mod p) かつ a ≡ 0 (mod q)」

これらの解のうち大きい方が M(n) です。これは一般の場合にも言えます。
ChineseRemainder[{0, 1}, {p, q}] のようにして求めました。

結果は約35分。まだまだ。

第2案:逆モジュロ

「a ≡ 0 (mod p) かつ a ≡ 1 (mod q)」の場合を考えます。
第1式から a=px とおいて第2式に代入すると

 px \equiv 1 \pmod{q}

x は q を法とする p の逆モジュロです。a=px に代入して

 a=p\cdot (p^{-1} \mod{q})

もう1つの場合も同様で,M(n) はこれら2数のうち大きい方です。

以上を PowerMod を使って書いたところ,約10分になりました。


第3案:互いに素な約数だけ考える

a と a-1 が互いに素であることから,n の約数 d として意味があるのは d と n/d が互いに素なものだけです。これまでは f[n] の CoprimeQ でチェックしていましたが,f[n] に代入する前に絞り込むことにします。

下のコード中では Times やら Subsets やらネストさせていますが,たとえば
360=2^3 * 3^2 * 5
を次のように変換しています。

  1. FactorInteger で素因数分解
  2. Cases で {a, b} を a^b に直す
  3. Subsets で部分集合にわける
  4. Times を Apply して積を求める

{{2, 3}, {3, 2}, {5, 1}}
→ {8, 9, 5}
→ {{}, {8}, {9}, {5}, {8, 9}, {8, 5}, {9, 5}, {8, 9, 5}}
→ {1, 8, 9, 5, 72, 40, 45, 360}

また,「√n 以下」で抽出しておきながら f[n] で2通りの値を計算するのは無駄なので,この処理はなくしました。結果は93秒。残念ながら1分は切れませんでした。