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

PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 142 / 完全平方数コレクション

x+y, x-y, x+z, x-z, y+z, y-z がすべて平方数となる整数の組 x > y > z > 0で, 最小の x+y+z を求めよ。

Problem 142 - Project Euler


与えられた数を順に a^2, b^2, ..., f^2 とおいて x, y, z について解きます。

  • x=(a^2+b^2)/2=(c^2+d^2)/2
  • y=(a^2-b^2)/2=(e^2+f^2)/2
  • z=(c^2-d^2)/2=(e^2-f^2)/2

整理すると

  • b^2=c^2-e^2
  • d^2=a^2-e^2
  • f^2=a^2-c^2
  • aとb, cとd, eとfは偶奇が等しい
  • a^2>b^2, c^2>d^2, e^2>f^2, a^2>c^2>e^2
  • x+y+z=(a^2+c^2+e^2)/2

はじめの3式を使うと,偶奇と大小の条件が整理できます。

 a^2\equiv c^2+e^2\pmod{2},\, c^2+e^2>a^2

この方法で計算してみたら7分半かかりました。


a, c, e のかわりに a, b, c を使うと,「aとbは偶奇が等しい」が言えて速くなりますが,それでも5分以上かかります。

  • d^2=a^2+b^2-c^2
  • e^2=c^2-b^2
  • f^2=a^2-c^2
  • aとbは偶奇が等しい
  • a^2>c^2>b^2, 2c^2>a^2+b^2
  • x+y+z=c^2+(a^2-b^2)/2

また,この方法では c^2+(a^2-b^2)/2 を求めるのですが,「a は大きいが,a≒b のため全体としては小さな値をとる」ケースが考えられ,最小値を求めているとはいえないと思います。

ピタゴラス数を使うともっと速くなりますが,私の腕では3分くらいにしかなりませんでした。1分の壁は厚いです。

Project Euler 133 / レピュニットの非因数

1のみからなる数をレピュニット数という。長さ k のレピュニット数を R(k) であらわす。
R(6) = 111111

R(10^n) というレピュニット数について考える。

R(10), R(100), R(1000) は 17 では割り切れないが,R(10000) は 17 で割り切れる。さらに,R(10^n) が 19 で割り切られるような n は存在しない。驚くべきことに, R(10^n) の因数となりうる100未満の素数は 11, 17, 41, 73 の4個のみである。

R(10^n) の因数となりえない100000未満の素数の和を求めよ。

Problem 133 - Project Euler


レピュニット数の問題の中では問題129と似ています。

variee.hatenadiary.com

R(10^n) が素数 p で割り切れる条件は

 p\, \left|\, \dfrac{10^{10^n}-1}{9} \right.\ \Leftrightarrow\ 9p\, |\, (10^{10^n}-1)

 \therefore 10^{10^n} \equiv 1\pmod{9p}

p≠3 はわかっているので,法から9を落とせます。

 10^{10^n} \equiv 1\pmod{p}

これをみたす n が存在すれば,その n に対して R(10^n) は p で割り切れます。問題129同様,乗法的位数を求める関数
MultiplicativeOrder を使って言い換えます。

R(10^n) の因数になりえない
= MultiplicativeOrder[10, p] は 10^n の約数ではない
= MultiplicativeOrder[10, p] は 2, 5 以外の素因数をもつ


最初に書いたコードがこちらです。

2, 5 以外の素因数をもつかどうか判定する方法がわからなかったので,2, 5 以外の素数を全部かけたものを用意して,それと互いに素かどうか調べるという力技をやっています。意外に時間はかからなくて0.18秒でした。

また,2と5は MultiplicativeOrder[10, p] が定義できないので,2, 3, 5は別扱いとしています。

素因数の判定部分を MatchQ で書き直したのが次のコードです。だいぶスッキリしましたが,時間はほとんど変わりませんでした。

Project Euler 105 / 特殊和集合:検査

大きさ n の集合 A の要素の和を S(A) で表す。空でなく共通要素を持たないいかなる 2 つの部分集合 B と C に対しても以下の性質が真であれば,A を特殊和集合と呼ぼう。

  1. S(B)≠S(C); つまり,部分集合の和はすべて異なる
  2. B が C より多くの要素を含んでいれば S(B) > S(C) となる

{81, 88, 75, 42, 87, 84, 86, 65} は
65 + 87 + 88 = 75 + 81 + 84
であるため特殊和集合ではないが,
{157, 150, 164, 119, 79, 159, 161, 139, 158} はすべての可能な部分集合の対の組み合わせについて両方のルールをみたしており,また S(A)=1286 である。

7 から 12 の要素を含む 100 個の集合が記された 4K のテキストファイル sets.txt を用いて特殊和集合 A1, A2, ... , Ak をすべて特定し,
S(A1) + S(A2) + ... + S(Ak) を求めよ。


Problem 105 - Project Euler


素朴に考えて解きました。下手にサブリストを作ると Map や Apply の作用の仕方がわからなくなってしまいそうなので,sets.txt の要素を1つずつ取り出してチェックしています。

まず,部分集合を要素数 i でグループ分けして,和の最大値 max(i) と最小値 min(i) を求めます。

条件1のチェックは簡単で, DuplicateFreeQ するだけ。

条件2は max(i) < min(i+1) で調べました。

特に工夫するまでもなく1.4秒で終わって,ホッとしました。

Project Euler 132 / 巨大なレピュニットの素因数

1のみからなる数をレピュニット数という。長さ k のレピュニット数を R(k) であらわす。

R(10) = 1111111111 = 11×41×271×9091
の素因数の和は 9414 である。

R(10^9) の最初の40個の素因数の和を求めよ。

Problem 132 - Project Euler


非常に大きな数が相手ですが,Mathematica に組み込まれている
PowerMod[a, b, m](m を法として a^b を与える関数)はものともしない感じでした。

 p\, \left|\, \dfrac{10^k-1}{9} \right.\ \Leftrightarrow\ 9p\, |\, (10^k-1)

 \therefore 10^k \equiv 1\pmod{9p}

これをみたす素数を探します。最初は適当に範囲を広く取って探しました。0.033秒でした。

次に普通にループさせて40個数えました。p=2, 3, 5 は明らかに不適なので除きましたが,ほとんど影響はないはずです。ほんのすこし速くなって0.024秒でした。


variee.hatenadiary.com

Project Euler 130 / 素数桁レピュニットと同じ性質を持つ合成数

1のみからなる数をレピュニット数という。長さ k のレピュニット数を R(k) であらわす。
R(6) = 111111

GCD(n, 10)=1 となる正整数 n について, 必ず正整数 k が存在し n が R(k) を割り切ることが証明できる。A(n) でそのような最小の k を表す。
A(7) = 6, A(41) = 5

5より大きい素数 p について, A(p) が p-1 を割り切ることが知られている。p = 41 のときには A(41)=5 であり,40は5で割り切れる。

非常に少ないのだが,合成数においても上が成立する場合がある。最初の5つの例は 91, 259, 451, 481, 703 である。

GCD(n, 10)=1 かつ A(n) が n-1 を割り切るような最初の25個の合成数 n の総和を求めよ。

Problem 130 - Project Euler

この前後のレピュニットの問題を解いたあとでは易しい問題です。他の問題用のコードをコピペしてちょっといじるだけで解けます。

条件をみたす25個の数はこのようになっています。