PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 59 / XOR暗号解読

各文字はそれぞれ一意のコードに割り当てられている。よく使われる標準としてASCII (American Standard Code for Information Interchange) がある。ASCIIでは大文字 A = 65,アスタリスク (*) = 42,小文字 k = 107 のようにコードが割り当てられている。

モダンな暗号化の方法としてテキストファイルの各バイトをASCIIに変換し,秘密鍵から計算された値とXORを取る手法がある。XOR関数の良い点は暗号化に用いたのと同じ暗号化鍵でXORを取ると平文を復号できる点である。

  • 65 XOR 42 = 107
  • 107 XOR 42 = 65

破られない暗号化のためには鍵は平文と同じ長さのランダムなバイト列でなければならない。ユーザーは暗号文と暗号化鍵を別々の場所に保存する必要がある。もし一方が失われると暗号文を復号することは不可能になる。

悲しいかな,この手法はほとんどのユーザーにとって非現実的である。そこで,鍵の変わりにパスワードを用いる手法が用いられる。パスワードが平文より短ければ(よくあることだが)パスワードは鍵として繰り返し用いられる。この手法では安全性を保つために十分長いパスワードを用いる必要があるが,記憶するためにはある程度短くないといけない。

この問題での課題は簡単になっている。暗号化鍵は3文字の小文字である。cipher1.txt は暗号化されたASCIIのコードを含んでいる。また,平文はよく用いられる英単語を含んでいる。この暗号文を復号し,平文のASCIIでの値の和を求めよ。

Problem 59 - Project Euler

下調べ

手始めに a=97 から z=122 までのアルファベット一文字で復号した場合の文字コードを調べました。以下のコードで myList0 は問題文中の cipher1.txt です。

文字コード32は「スペース」,127は「DEL」です。意外なことに文字ごとの違いはほとんどありませんでした。

次に「aaa」で復号してみました。これは明らかに駄目です。スペースやピリオドによる区切りがないのと記号が多すぎるのとで,文字化けにしか見えません。

ASCII文字コード - IT用語辞典

本番

小文字とスペースは10点,大文字と数字は1点として最大値を与える組み合わせを探しました。

はじめは@などの記号があったら大幅減点するコードを書いたのですが,絶対に出てこないとは言い切れないのでやめました。たとえば本文中にメールアドレスが書かれていれば@が出てくるのは当たり前ですよね。

単純なルールの割にうまくいって,パスワードは「god」,その文字コードの和は314*1でした。

後始末

cipher1.txt を「god」で復号した結果は聖書の一部です。

各文字の出現頻度はこのようになります。最も多いのはスペースで,e の約2倍ありました。

文字頻度表 によると,CNN, ABC, VOA の各ニュースの各分野からランダムに4万1600文字を抜き出したときの出現頻度順は次のようになっていたそうです。大体の傾向は似ていますが,ところどころ違っています。

e, a, t, i, o, s, n, r, h, l, d, c, u, m, p, f, g, y, w, b, v, k, j, x, q, z

*1:正解者掲示板に「円周率」との指摘あり。

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倍以上速くしないといけないとは……

Project Euler 160 / 階乗の下位桁

任意の N に対して f(N) を N! の末尾の 0 の前にある最後の5桁とする。

  • 9! = 362880, f(9)=36288
  • 10! = 3628800, f(10)=36288
  • 20! = 2432902008176640000, f(20)=17664

f(1,000,000,000,000) を求めよ。

Problem 160 - Project Euler

下準備

2と5でどんどん割っていって,mod 10^5 で下5桁を取り出すのが基本的な流れです。まず,実験用の基本的な関数を用意します。

n が p で何回割り切れるかは IntegerExponent[n,p] で求められます。

n! が p で何回割り切れるかは次の g[n,p] で求められ,f[n] も定義できます。

f[n] の計算に中国式剰余定理も使えます。

  1. mod 2^5 と mod 5^5 を求める
  2. 中国式剰余定理でまとめ,10^5で割った余りを求める

n! を10でどんどん割っていくと2は普通5個以上残るので,mod 2^5 は基本的に0です。以下,mod 5^5 のみ考えます。


周期性を調べる

f[10^5] から f[25*10^5] まで一万きざみで求めてみると

 f(5n)=f(n)

が予想できます。

f:id:variee:20170702013941p:plain

もう少しくわしく調べると

 n=2500k\, \Rightarrow\, f(5n)=f(n)

が成り立っています。証明は後回しにして,これで答えを出してしまいましょう。10^12 を5でどんどん割っていって 2,560,000 に落としてから f を計算しました。


周期性の証明

g[n, 5]=g[n] と略記します。

 g(n)=\left[\dfrac{n}{5}\right]+\left[\dfrac{n}{5^2}\right]+\left[\dfrac{n}{5^3}\right]+\cdots

から  g(5n)=n+g(n) が言えます。

 f(5n)=mod\left[\dfrac{(5n)!}{10^{g(5n)}},\, 5^5\right]=mod\left[\dfrac{(5n)!}{n!\cdot 10^n}\cdot \dfrac{n!}{10^{g(n)}},\, 5^5\right]

=mod\left[\dfrac{(5n)!}{n!\cdot 10^n},\, 5^5\right]\cdot mod\left[\dfrac{n!}{10^{g(n)}},\, 5^5\right]

=mod\left[\dfrac{(5n)!}{n!\cdot 10^n},\, 5^5\right]\cdot f(n)

 \varphi(5^5)=4\cdot 5^4=2500 から  2^{2500}\equiv 1\pmod{5^5} が言えるので,n が 2500 の倍数のとき

 2^n \equiv 1\pmod{5^5}

が成り立ちます。 f(5n)=f(n) を示すのに

 \dfrac{(5n)!}{n!\cdot 5^n}\equiv 1\pmod{5^5}

を言えば十分です。

左辺は 5n 以下で 5^5 と互いに素な自然数の積です。wilson の定理の拡張版を使います。

 \displaystyle\prod_{k=1,\, gcd(k,m)=1}^m k\equiv -1 \pmod m

Wilson's theorem - Wikipedia

この定理から次が言えます。

 \displaystyle\prod_{k=1,\, gcd(k,5^5)=1}^{5^5} k\equiv -1 \pmod{5^5}

n が 2500=4*5^4 の倍数なので 5n は 5^5 の倍数であること,余りの周期性を使います。

 \dfrac{(5n)!}{n!\cdot 5^n}\equiv\displaystyle \left(\prod_{k=1,\, gcd(k,5^5)=1}^{5^5} k\right)^{\frac{5n}{5^5}}\equiv (-1)^{\frac{5n}{5^5}}=1 \pmod{5^5}

最後のところで  \frac{5n}{5^5} が偶数であることを使いました。

これで証明完了です。計算時間は1秒未満でしたが,証明には数日かかりました。

Project Euler 272 / モジュラー立方数2

正の整数 n に対し,C(n) を 1< x< n で x3≡1 mod n をみたす整数 x の数と定義する。

n=91 のとき,x は 8 つの値を取りうる。9, 16, 22, 29, 53, 74, 79, 81 であり,C(91)=8 である。

C(n)=242 をみたす正整数 n≤10^11 の合計を求めよ。

Problem 272 - Project Euler

検算用の式を作る1

PowerModList を使うと C(n) が作れます。まずは検算用に絶対に間違えようのないコードを書きました。


C(91)=8 の吟味

中国式剰余定理を使うと C(91)=8 になる理由がわかります。

 91=7\times 13,\, C(7)=3,\, C(13)=3

C(x) は C(pq)=C(p)C(q) (p, q は異なる素数)をみたします。3*3=9 個の整数から自明な1を除いて

C(91)=3\times 3-1=8

となるわけです。もう少し詳しく調べると次のことがわかります。

  • p が3で割って余り1の素数のとき, C(p^n)=3\quad (n\geqq 1)
  • p が3で割って余り2の素数のとき, C(p^n)=1\quad (n\geqq 1)
  •  C(3)=1, C(3^n)=3\quad (n\geqq 2)

3で割って余り1の素数を「よい素数」,余り2の素数を「悪い素数」と呼ぶことにしましょう。
242=3^5-1 なので,条件をみたす数は次のいずれかの形をしています。

  • (5種類のよい素数の積)x(悪い素数の積)
  • (5種類のよい素数の積)x(悪い素数の積)x3
  • (4種類のよい素数の積)x(悪い素数の積)x(2個以上の3の積)

条件をみたす数の最小値は

3^2\times 7\times 13\times 19\times 31=482391

また,よい素数は9268個ありました。

検算用の式を作る2

よい素数と9に注目して検算用の式その2を作りました。できればこれで brute force したかったのですが,10^10以上は「数が大きすぎて扱えない」エラーが出てしまいます。


本番その1

次の手順で解いたところ,約55分かかりました。

  1. 9とよい素数の和集合(myList)を作る
  2. myList の元を5個かけあわせた数 p[a,b,c,d,e] を作る
  3. p[a,b,c,d,e] の倍数で C(n)=242 となるものの和を求める


本番その2

p[a,b,c,d,e] の倍数について調べるのに「かけてから素因数分解」は効率が悪いと思ったので,「かける相手をあらかじめ素因数分解しておいて使いまわす」に変えてみましたが,時間は変わりませんでした。

包除原理を使うとよさそうなのですが,Mathematica っぽい書き方をするとメモリ不足でエラーになってしまうのが難点です。

Project Euler 224 / だいたい直角三角形2

三辺(a≤b≤c)がすべて整数の三角形で,三辺が次の式をみたしているものを"わずかに鈍角(barely obtuse)"と呼ぶ。

 a^2+b^2=c^2-1

周辺の長さが 75,000,000 以下のわずかに鈍角な三角形はいくつあるか。

Problem 224 - Project Euler

第223問と同じ方法で解きました。

variee.hatenadiary.com


「方程式の整数解を求める」で 282分。

「素因数分解する」で 18 分。

第223問と比べると検索範囲の広さは3倍ですが,a^2+1 は a^2-1 にくらべて約数が少ないのか計算時間は 23分→18分と短くなりました。

正解者掲示板を見たところ,速い人は大体1つの解から次の解を作る漸化式を使っているようでした。ただ,「本当にすべてを尽くしているのか?」「その証明は?」と考えると,自分的には素因数分解による解法が現実的に思えます。