PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 126 / 直方体層

3x2x1 の直方体の表面すべてを覆いつくすのに必要最小限の立方体の数は 22 個である。

f:id:variee:20180608145751g:plain

さらにこの立体に表面を覆いつくすように2層目を追加すると,46 個の立方体が必要である。3層目は 78 個,4層目は 118 個の立方体が必要である。

ところで 5x1x1 の直方体への1層目も 22 個の立方体が必要である。同様に 5x3x1, 7x2x1, 11x1x1 の直方体への1層目もすべて 46 個の立方体である。

何層目かが n 個の立方体からなる直方体の数を C(n) と定義する。
C(22) = 2, C(46) = 4, C(78) = 5, C(118) = 8

154 は C(n) = 10 をみたす最小の n である。C(n)=1000 をみたす最小の n を求めよ。

本家:Problem 126 - Project Euler
日本語訳:Problem 126 - PukiWiki

C(n)を求める

受験算数でよくある表面積の問題のように真正面,真上,真横から見た図を描いて数えます。

 C(1)=2(ab+bc+ca)
 C(n+1)=C(n)+4(ab+bc+ca)+8(n-1)

これを使って一般項を求めました。
 C(n)=4n^2+4(a+b+c-3)n
 \qquad\qquad +8-4(a+b+c)+2(ab+bc+ca)

C(n)=1000 をみたす n を探す

a, b, c と層の数 k についての4重ループを使って,C(n)=1000 をみたす最小の n を探します。

アタリをつけるため,はじめは「a, b, c は100以下で,k は10以下」のような条件で探したのですが,この方法では C(n) の正しい値を求められませんし,時間も非常にかかります。

C(n) の上限を20000として検索するよう直したところ,Mathematica を使う意味を感じない,普通の手続き型っぽいソースとなりました。

計算時間は約60秒。C++でやってみたら1秒でした。通常,C++の100倍以上の時間がかかるので,Mathematicaとしては頑張ったほうかなと思います。

Project Euler 106 / 特殊和集合:メタ検査

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

  1. S(B) ≠ S(C):部分集合の要素の総和は異なる
  2. B が C より多くの要素を含んでいれば S(B) > S(C)

本問題に対しては,与えられた集合は n 個の単調増加する要素を含み,かつ第二のルールをすでにみたしているものとする。

n = 4 の集合から得ることができる 25 個の可能な部分集合の対のうち,1 個の対のみが 同一性(第一のルール)をテストされる必要がある。
同様に n = 7 のときは 966 個の部分集合の対のうち 70 個のみがテストされる必要がある。

n = 12 に対して得られる 261625 個の部分集合の対のうち,同一性をテストされる必要があるものは何個か。

注意:この問題は Problem 103 と 105 に関連している*1

本家:Problem 106 - Project Euler
日本語訳:Problem 106 - PukiWiki

n=4 のとき

「n = 4 (中略) 25 個の可能な部分集合の対のうち~」を理解するのにかなり時間がかかりました。

A の要素を小さい方から順にa1, a2,..., an とします。この問題で考えているのは部分集合そのものではなく,部分集合の対なので
「25個ある」=「添字の組の選び方が25通りある」
です。

また,この25という数字は問題文中の「第二のルールをすでにみたしているものとする」を無視した場合の全組み合わせ数のことでした。実際に調べなければならない組はもっと少ないです。

  • ルール2 により,要素の個数が異なる集合を比較する必要はない
  • 要素数1の集合も比較する必要がない

結局,n=4 のとき考えなければならないのは要素数2の場合だけです。
二項係数を C(n, k) のように書くと

 C(4,2)\times C(2,2)/2=\text{3 通り}

あります。

  • {1, 2}, {3, 4}
  • {1, 3}, {2, 4}
  • {1, 4}, {2, 3}

a1 < a3, a2 < a4 より1つ目のケースは調べる必要がありません。
2つ目も同様。
3つ目は a1 < a2, a4 > a3 より調べる必要があります。問題文中の「1 個の対のみ」はこれのことです。

一般化

要素数を k として,選んだ添字を

  • {x1, x2, ..., xk} (x1< x2< ... < xk)
  • {y1, y2, ..., yk} (y1< y2< ... < yk)

とします(x1< y1)。

x2< y2, x3 < y3, ..., xk < yk がすべて成立すれば調べる必要がなく,少なくとも1つが不成立なら調べる必要があります。

調べなくてよいケースの方を数えましょう。これは1から 2k までの整数を2段に並べるときに「右の数は左より大きい」「下の数は上より大きい」をみたす並べ方と同じだけあり,カタラン数 Ck で与えられます。

 C_k=\displaystyle \frac{{}_{2k} C_{k}}{k+1}

x1 ~ yk に使う数字の選び方が C(n, 2k) 通りあるので,要素数 k で調べなくてよいケースの個数は次のとおりです。

 \displaystyle \frac{{}_{2k} C_{k}}{k+1}\cdot {}_{n} C_{2k}

これを組み合わせの総数から引きます。要素数が k で調べなければならないものの個数は次のようになります。

 f(n,k)= \displaystyle  {}_{n} C_{k} \cdot {}_{n-k} C_{k}\cdot \frac{1}{2} - \frac{{}_{2k} C_{k}}{k+1}\cdot {}_{n} C_{2k}
 = \displaystyle  \frac{(k-1)\, n!}{2\cdot (k+1)!\, k!\, (n-2k)!}

これの和をとりました。





*1:実際にはほぼ関係ない。

Project Euler 624 / Two heads are better than one

コイン投げで表が連続して2回出るまでに M 回かかるとする。M が n で割り切れる確率を P(n) であらわす。
P(2)=3/5, P(3)=9/31 である。

有理数 a/b と素数p に対して a≡bq (mod p) をみたす最小の自然数 q を Q(a/b, p) であらわす。

たとえば Q(P(2), 109)=Q(3/5, 109)=66 である。
5x66 =330 ≡ 3 (mod 109)
同様に Q(P(3), 109)=46 である。

Q(P(10^18), 10^9+9) を求めよ。

本家:Problem 624 - Project Euler

下準備

M-1 回目までは「裏は連続していいが,表は連続しない」というルールで表と裏が並びます。これは入試数学でもよくある設定で,フィボナッチ数 f[M-1] 通りあります(証明略)。

「M が n で割り切れる」=「M は n の倍数」から M=kn とおけて,P(n) は次のように表せます。

 P(n)=\displaystyle\sum_{k=1}^{\infty} \frac{f_{kn-1}}{2^{kn}}

フィボナッチ数列の一般項を使ってこの和を求め,P(2)~P(10) を求めました。

この形の P(n) で計算できる Q はせいぜい n=100 程度までなので,もっとうまい方法で P(n) を求めないといけません。

行列を使う

上で求めた P(n) は分母,分子とも指数関数です。底は

  • 分子は -1, 1-√5, 1+√5
  • 分母は -1, 4, 1-√5, 1+√5

これらを解にもつ特性方程式を作って,対応する漸化式を求めました。
P(n)=a(n)/b(n) とおきます。

 a_1=1,\, a_2=3,\, a_3=9
 a_{n+3}=a_{n+2}+6 a_{n+1}+4 a_{n}

 b_1=1,\, b_2=5,\, b_3=31,\, b_4=145
 b_{n+4}=5 b_{n+3}+2 b_{n+2}-20 b_{n+1}-16 b_{n}

これを行列に直して P(n) の分母,分子を求めます。

a≡bq (mod p) は

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

でOK。あっさり終わります*1

おまけ

MatrixPowerMod は昔は普通に使えたような気がしますが,現行バージョンのマニュアルには載っていません。stackexchange に助けられました。

mathematica.stackexchange.com

*1:最近の問題なので答えは載せないことにします

Project Euler 128 / 六角形タイルの差分

六角形のタイルを図のように反時計回りに置いていく。

f:id:variee:20180602222929g:plain

タイル n とそれに隣接するタイルのうち,書いてある数字と n の差が素数となるタイルの枚数を PD(n) と定義する。

たとえばタイル8の場合,まわりのタイルの数字と8の差は 12, 29, 11, 6, 1, 13 となるので PD(8) = 3 であり,同様に PD(17) = 2 である。

PD(n) の最大値は3であることが示せる。
PD(n) = 3 をみたす n を昇順に並べた数列の10番目の項は271である。

この数列について2000番目のタイルを求めよ。

本家:Problem 128 - Project Euler
日本語訳:Problem 128 - PukiWiki

実験

与えられた図をもとに調べると最初の6項は 1, 2, 8, 19, 20, 37 で,まわりの数や差は次のようになっています。カッコで囲んだ数は素数です。

n まわりの数 差(ソート後)
1 2,3,4,5,6,7 1,(2),(3),4,(5),6 1,(2),(3),4,(5),6
2 1,3,7,8,9,19 1,1,(5),6,(7),(17) 1,1,(5),6,(7),(17)
8 2,9,19,20,21,37 6,1,(11),12,(13),(29) 1,6,(11),12,(13),(29)
19 2,7,8,18,36,37 (17),12,(11),1,(17),18 1,(11),12,(17),(17),18
20 8,21,37,38,39,61 12,1,(17),18,(19),(41) 1,12,(17),18,(19),(41)
37 8,19,20,36,60,61 (29),18,(19),1,(23),24 1,18,(19),(23),24,(29)

規則性があるのかないのかはっきりしない感じですが,中央のマスの数字をxとして,まわりに x+1 や x-1 があったり,連続する2数(または3数)があったりすることがわかります。

以下,8~19や20~37のような並びを「周期」と呼ぶことにし,n 個目の周期の先端を H(n),末端を T(n) であらわします。8~19について考えます。

まわりのタイルのパターンは4つ

中央のマスの数字をxとして,そのまわりの6枚のタイルに書かれた数字のパターンを分類してみました。

辺の途中

9, 11, 13, 15, 17, 19は
{x-1, x+1, y, y+1, z, z+1} 型です。

  • x-(x-1)=1, (x+1)-x=1 は素数ではない
  • x-y, x-(y+1) のうち素数はたかだか1個
  • z-x, (z+1)-x のうち素数はたかだか1個

素数はたかだか1+1=2個しかありません。

先端,末端以外の頂点

10, 12, 14, 16, 18は
{x-1, x+1, y, z-1, z, z+1} 型です。
x-1, x+1が素数を与えないことは上で見たので,他の数について考えます。

たとえばx=10のときy=3, z=23です。H(1)=2, H(2)=8, H(3)=20から順に1,2,3進んだところにy=3, x=10, z=23があり,一般に
{z-H(n+1)} - {x-H(n)} と {x-H(n)}-{y-H(n-1)} は一致します。H(n)はつねに6の倍数なので,次のことが言えます。

  • y, zの偶奇は等しい
  • zとz+1, z-1の偶奇は異なる

要するに偶奇に注目するとxのまわりの残り4数は {y, z}, {z+1, z-1} の2グループに分類できるので,素数はたかだか2個しかありません。

先端

H(n) =3n(n-1)+2 の真上から反時計回りに
H(n+1), H(n+1)+1, H(n)+1, H(n-1), T(n), T(n+1)
が並びます。H(n), T(n) を6で割った余りが0, 5であることを考えると,素数を与える可能性があるのは次の3箇所です。
H(n+1)+1, T(n), T(n+1)

末端

T(n) の真上から反時計回りに
T(n+1), H(n), H(n-1), T(n-1), T(n)-1, T(n+1)-1
が並びます。先端の場合と同様,素数を与える可能性があるのは3箇所です。
H(n), H(n-1), T(n+1)-1

解法

先端と末端について調べればOK。0.5秒でした。

tail[n] - head[n]などを展開した式(すべて1次式)を使うと0.3秒に縮まりますが,この差は誤差の範囲内ですね。

Project Euler 111 / 重複桁を持つ素数

重複した桁を含む 4 桁の素数を考える。
各桁に1は最大3個含まれ,このような素数は9つある。

1117, 1151, 1171, 1181, 1511, 1811, 2111, 4111, 8111


n 桁の素数を考える。各桁に d が最大 M(n, d) 個含まれるとする。このような素数が N(n, d) 個あり,それらの和が S(n, d) であるとする。


Digit, d M(4, d) N(4, d) S(4, d)
0 2 13 67061
1 3 9 22275
2 3 1 2221
3 3 12 46214
4 3 2 8888
5 3 1 5557
6 3 1 6661
7 3 9 57863
8 3 1 8887
9 3 7 48073

S(4, d) (0 <= d <= 9) の総和は 273700 である。
S(10, d) の総和を求めよ。

本家:Problem 111 - Project Euler
日本語訳:Problem 111 - PukiWiki

n=4の場合

4桁の素数のうち,1を3個含むものは9個あります。これは

 \{1,1,1,x\}\quad (x=0,\, 1,\, \cdots,\, 9)

の permutation が素数になるかどうか調べればわかります。

0を2個含むものが13個あることは {0, 0, x, y} の permutation について調べればわかります。

n=10の場合

以上を一般化すればn=10の場合も解けます。ただ,桁数が増える分,調べる回数も増えます。本当は

「ある数を9個含むものがあるかどうか調べる」
→「なければ8個含むものを調べる」
→「なければ7個含むものを調べる」
→……

のようにして解くべきなのですが,面倒だったので8個までで十分なことを事前に調べた上で答えを出しました。


内訳はこうなっていました。

d M(10,d) N(10,d) S(10,d)
0 8 8 38000000042
1 9 11 12882626601
2 8 39 97447914665
3 9 7 23234122821
4 9 1 4444444447
5 9 1 5555555557
6 9 1 6666666661
7 9 9 59950904793
8 8 32 285769942206
9 9 8 78455389922

n=4 の場合もそうですが,M(n, d)=n-1, n-2 まで考えれば十分なのが不思議です。