PEをMathematicaで

Project Eulerに挑戦してみよう

算チャレ第1021回を解いてみた

ある4ケタの整数があり,次のような性質を満たしています。

(性質)各ケタの数字を並べかえて作ることのできる整数のうち最大のものから最小のものを引くと,もとの整数になる。

例えば1345という数は「1, 3, 4, 5」の4つの数字で出来ていて,これらを並べ替えてできる最大の整数は5431,最小の整数は1345で,その差は
5431-1345=4086
となり,条件を満たしませんね。

では,上記の性質を満たす4ケタの整数を求めてください。

http://www.sansu.org/used-html/index1021.html


ひさしぶりに算数チャレンジさんを見てみたら,Mathematica 向きの問題があったので解いてみました。
www.sansu.org

IntegerDigits で桁ごとに区切って Sort するだけなので,一瞬で答えが出ます。6174でした。

これで終わってはつまらないので,検索範囲を広げてみました。10^8以下で探すと,条件をみたす数は3桁に1個,4桁に1個,6桁に2個,8桁に2個あることがわかります。

Project Euler 512 / べき乗のトーシェント数の和

オイラーのトーシェント関数を φ(n) とする。

 f(n)=\left(\sum_{i=1}^n \varphi (n^i)\right)\pmod{n+1}

 g(n)=\sum_{i=1}^n f(i)

とすると g(100) = 2007 となる。

g(5*10^8) を求めよ。

Problem 512 - Project Euler


 \varphi(n^i)=n^{i-1} \varphi(n) を使って f(n) の式を整理します。

 \displaystyle\sum_{i=1}^n \varphi (n^i)=\varphi(n)\displaystyle\sum_{i=1}^n n^{i-1}\equiv \varphi(n)\displaystyle\sum_{i=1}^n (-1)^{i-1} \pmod{n+1}

これは n が偶数のとき 0,奇数のとき φ(n) です。奇数について φ(n) の和を求めます。5分かかりました。

Project Euler 451 / モジュラ逆数

数 15 について考えよう。15と互いに素となる15以下の正数は8個ある。1, 2, 4, 7, 8, 11, 13, 14 である。それらの数の15を法とするモジュラ逆数は 1, 8, 4, 13, 2, 11, 7, 14 である。

  • 1*1 mod 15 = 1
  • 2*8 = 16 mod 15 = 1
  • 4*4 = 16 mod 15 = 1
  • 7*13 = 91 mod 15 = 1
  • 11*11 = 121 mod 15 = 1
  • 14*14 = 196 mod 15 = 1

m の法 n に対するモジュラ逆数が m 自身となるような n-1 より小さい最大の正数 m を I(n) としよう.
I(15)=11, I(100)=51, I(7)=1

3 ≤ n ≤ 2*10^7 における ΣI(n) を求めよ。

Problem 451 - Project Euler


「m の法 n に対するモジュラ逆数が m 自身」を合同式で表します。

 m^2\equiv 1\pmod{n}\quad \therefore (m+1)(m-1)\equiv 0\pmod{n}

n を素因数分解してこれを処理すればいいことになります。m+1, m-1 両方の約数になりうる素数は2だけです。

ほとんどのプラミング言語ではこのように解くと思いますが,Mathematica はモジュラ逆数を直接計算できます。100のモジュラ逆数を求めてみましょう。

2番目に大きい 51 が I(100) です。因数分解した合同式から m=n-1 が必ずモジュラ逆数になることがわかるので,一般に I(n) はこの PowerModList の後ろから2番目の数です。その和を求めればできあがりです。

Project Euler 258 / ラグ付フィボナッチ数列

数列を以下のように定義する。

  • 0 ≤ k ≤ 1999 に対して gk = 1
  • k ≥ 2000 に対して gk = g(k-2000) + g(k-1999)

k = 10^18 に対して gk mod 20092010 を求めよ。

Problem 258 - Project Euler


10^18 という大きな添字,20092010という大きな法をどうするかがポイントです。添字は特性方程式で処理し,法は中国式剰余定理で処理しました。

■10^18 の処理

特性方程式 x^(2000)-x-1=0 を使って x^(10^18) を次数下げします。普通に

f[x_]:=PolynomialRemainder[x^(10^18), x^2000 - x - 1];
f[1]

とすると,かなり時間がかかった後で「再帰の回数が限界を超えた」エラーが出てしまいます。ここでしばらくハマりましたが,回避法がわかりました。

  • 係数がかなり大きくなるので,Modulus を指定してあらかじめ割っておく
  • 代入は「/. x->1」を使った ReplaceAll で行う
■20092010 の処理

法の 20092010 を素因数分解します。

 20092010=2\cdot 5\cdot 859\cdot 2339

これらを PolynomialRemainder の Modulus に指定して,最後に ChineseRemainder でまとめます。似たような計算を4回行うことになりますが,意外に時間はかかりませんでした。

PolynomialRemainder は素数しか法に指定できないので,20092010 が各素数を1個ずつしか含まなくて助かりました。

Project Euler 268 / 少なくとも4つの異なる100未満の素数を因数に持つ数

少なくとも 4 つの異なる 100 未満の素数で割り切れる 1000 未満の正の整数は 23 ある。

少なくとも 4 つの異なる 100 未満の素数で割り切れる 10^16 未満の正の整数はいくつあるか答えよ。

Problem 268 - Project Euler


まずは単純な場合を考えました。たとえば「2, 3, 5 のうち少なくとも2つで割り切れるような n 以下の自然数の個数を求めよ」だったら? 素数を1個しか含まない数を数えて全体から引くのが自然ですが,この解法は今考えている問題には向かないので普通に数えます。

f(x)=[n/x]([ ]はガウス記号)とおいて,ベン図を考えると

 f(2\cdot 3)+f(3\cdot 5)+f(5\cdot 2)-2f(2\cdot 3\cdot 5)

これを一般化します。

上式の最後の項の係数が1ではないのが,この問題一番の難所でした。k(≧4)個の数を考えるとき,ここには
C(k,4) - C(k,5) + C(k,6) + ... + (-1)^(k-4)C(k,k)
が来ます。このままでも計算できますが,C(k-1,3) とまとめられることに気づくと速いです。

また,100未満の素数は25個ありますが,これら全部を相手にする必要はありません。素数を14個以上かけると 10^16 を超えてしまうので,最大13個の素数で割り切れる数を考えれば十分です。

あと,細かいことですが「割り算してガウスをとる」よりも「商を求める」の方が速いようです。これだけで時間が約半分になりました。