PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 183 / 分割した積の最大値

正整数 N に対して P=(N/k)^k (k >=2) の最大値を M(N) であらわす。

 M(11)=\left(\dfrac{11}{4}\right)^4=\dfrac{14641}{256}=57.19140625

M(N)が無限小数のとき D(N)=N とし,M(N) が有限小数のとき D(N)=-N とする。
5 ≦ N ≦ 100の とき,ΣD(N) = 2438

5 ≦ N ≦ 10000 のとき,ΣD(N) を求めよ。

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


受験数学っぽい問題でした。f(x)=(N/x)^x は x=N/e で最大になります。
k1=floor(N/e), k2=k1+1 として

 M(N)=\max\{f(k_1),\, f(k_2)\}

対応する k を kk とおきます。N/kk の分母が2と5以外の素因数をもつかどうか調べれば D(N) を決定できます。

Project Euler 159 / 因数分解の数字根和

合成数は多くの異なった方法で因数分解することができる。たとえば 24 は 7 通りに因数分解される。

  • 24 = 2x2x2x3
  • 24 = 2x3x4
  • 24 = 2x2x6
  • 24 = 4x6
  • 24 = 3x8
  • 24 = 2x12
  • 24 = 24

ある数の各桁の数字を足し合わせることを和が 10 未満になるまで繰り返したときに得られる数を数字根 (digital root, DR) と呼ぶ。467 の数字根は 8 である。
4+6+7=17 → 1+7=8
dr(467)=8

それぞれの因数の数字根の和を数字根和 (Digital Root Sum , DRS) と呼ぶことにする。以下に 24 の DRS を示す。

因数分解Digital Root Sum
2x2x2x39
2x3x49
2x2x610
4x610
3x811
2x125
246

24 の数字根和の最大値は 11 である。mdrs(n) を n の数字根和の最大値と定義する。
mdrs(24) = 11

1 < n < 1,000,000 について ∑mdrs(n) を求めよ。

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


数字根は

  • n が9で割り切れるなら9
  • 割り切れないなら 9 で割った余り

となります。

mdrs(n)は再帰で求めます。mdrs(n)=f(n)として

 f(24)=\max\{f(2)+f(12),\, f(3)+f(8),\, f(4)+f(6),\, dr(24)\}
 f(12)=\max\{f(2)+f(6),\, f(3)+f(4),\, dr(12)\}

などが言えます。下から順に求めていけばOK。

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:最近の問題なので答えは載せないことにします