PEをMathematicaで

Project Eulerに挑戦してみよう

Project Euler 90 / 2つの立方体の数字

立方体の6面それぞれに異なる数字(0から9)が書かれている。2番目の立方体も同様になっている。異なる位置に2つの立方体を隣りあわせることで様々な2桁の数を作ることができる。たとえば平方数である64も作ることができる。

f:id:variee:20170812010252g:plain

両方の立方体の数字をうまく選ぶと100以下のすべての平方数(01, 04, 09, 16, 25, 36, 49, 64, 81)を作ることが可能である。

たとえば {0, 5, 6, 7, 8, 9} を一方の立方体に,{1, 2, 3, 4, 8, 9} をもう一方の立方体に配置すればよい。

しかし,6と9を逆さまに回転することを許すと {0, 5, 6, 7, 8, 9} と {1, 2, 3, 4, 6, 7} のような配列で9つすべての平方数をあらわすことができる。

順番ではなくそれぞれの立方体の数字に着目して配列を区別する。

  • {1, 2, 3, 4, 5, 6} は {3, 6, 4, 1, 2, 5} と同じものとし,
  • {1, 2, 3, 4, 5, 6} は {1, 2, 3, 4, 5, 9} と異なるものとする。

6と9を逆さにすることを許すため,最後の例で区別された両方の配列のかわりに {1, 2, 3, 4, 5, 6, 9} という要素数が7つに拡張された配列を使用して2桁の数を作ることにする。

すべての平方数を表示しうる2つの立方体の異なる配列の組はいくつあるか。

Problem 90 - Project Euler

チェック用の関数を作る

まずは問題文中で与えられた2つのサイコロ a, b が条件をみたすことを確かめる関数を作ります。たとえば01を作れる条件は「aが0を含み,bが1を含む」または「aが1を含み,bが0を含む」です。MemberQ を使います。

無事,チェック成功しました。


6と9の扱い

「6しか含まなかったら9を加え,9しか含まなかったら6を加える」と変換することにしました。


本番

Subsets[Range[0, 9], {6}] でサイコロのもとを作ってチェックしました。無事,1.1秒で終了。

Project Euler 417 / 逆数の循環節その2

単位分数とは分子が1の分数である。

分母が 2 および 5 の素因数だけからなるとき,その単位分数は循環節を持たないと考えられる。そのような単位分数の循環節の長さは 0 であるとしよう。

1/n の循環節の長さを L(n) で表す。3 ≤ n ≤ 1,000,000 における ΣL(n) は 55535191115 である。

3 ≤ n ≤ 100,000,000 における ΣL(n) を求めよ。

Problem 417 - Project Euler


「逆数の循環節その1」はなんと第26問。随分と離れた類題です。
variee.hatenadiary.com

「循環節の長さをどう求めるか」「2と5をどう処理するか」がポイントかと思います。

循環節の長さをどう求めるか

問題26では素直に RealDigits[1/n] で循環小数を求めてその長さを調べましたが,さすがにこの方法では 10^8 は扱えそうにありません。
 10^x \equiv 1\pmod{n} をみたす x の最小値として求めました。

n=10^3 ではほぼ同じスピードですが,n=10^5 では約5倍違います。

2と5をどう処理するか

上の f417[n] は n=2^a 5^b のとき1を返します。「このような n の個数を求めて最後に引く」と「If 文を使って 0 を返すようにする」を試してみましたが,大差ありませんでした。


本番

以上をもとに最終的なコードを書きました。いろいろ工夫したものの1分の壁は厚く,5分かかりました。


英語版 wikipedia を参考に

正解者掲示板で英語版 wikipedia がくわしいと知りました。

Repeating decimal - Wikipedia

その情報をもとに書いてみたのがこちらのコードです。残念ながら,私の腕ではかえって遅くなってしまいました。

Project Euler 389 / 正多面体のサイコロ

1つの偏りのない四面体のサイコロを振り,出た目 T を記録する。
T 個の偏りのない六面体のサイコロを振り, 出た目の合計 C を記録する。
C 個の偏りのない八面体のサイコロを振り, 出た目の合計 O を記録する。
O 個の偏りのない十二面体のサイコロを振り, 出た目の合計 D を記録する。
D 個の偏りのない二十面体のサイコロを振り, 出た目の合計 I を記録する。
I の分散を四捨五入して小数点以下第4位まで求めよ。

Problem 389 - Project Euler

まずCの分散を求める

C の分散をどう求めるかを考えて,それを一般化しました。

C の確率分布

C の確率分布は一応求められます。

 P(C=c)=\displaystyle\sum_t P(T=t)P(C=c | T=t)

この式の P(C=c | T=t) は条件つき確率です。T=t の条件下で C=c となる確率を表します。P(後|前) のように書きます。

同じ方法で I の確率分布も表せますが,こういうダサい方法では解きたくないので Law of Total Variance(Eve's law)を使います。

Law of Total Variance(Eve's law)

 Var(C)=E(Var(C|T))+Var(E(C|T))

この公式は wikipedia にも載っています*1が,こちらのページがわかりやすいです。

www.probabilityformula.org

Var(C) の各項の計算に入りましょう。E(C|T) は T がわかっているときの C の期待値です。六面体のサイコロを T 個振ったときの目の和が C なので,1個だけ振ったときの目の T 倍になります。六面体のサイコロを1回振ったときの目を c とすると

 E(C|T)=E(c)T

E(c) は定数なので V(aX)=a^2 V(X) が使えて,Var(C) の右辺第2項は簡単な形になります。

 Var(E(C|T)) =Var(T) \left\{E(c)\right\}^2

次に Var(C) の右辺第1項を求めます。Var(C|T) はTがわかっているときの C の分散です。Cは六面体のサイコロを T 個振ったときの目の和で,各サイコロの出目は独立なので V(X1+X2)=V(X1)+V(X2) が使えます。

 Var(C|T)=Var(c)T
 \therefore E(Var(C|T))= Var(c) E(T)

以上を Var(C) の式に代入します。

 Var(C)=Var(c) E(T)+Var(T) \left\{E(c)\right\}^2

n 面体のサイコロでは E=(n+1)/2, V=(n^2-1)/12 となることを使うと,Var(C) が計算できます。

 Var(C)=\dfrac{6^2-1}{12}\cdot \dfrac{5}{2}+\dfrac{4^2-1}{12}\cdot \left(\dfrac{7}{2}\right)^2=\dfrac{1085}{48}

期待値も必要

上で導いた式を繰り返し使うには E(T), E(C), ……を求めておく必要があります。Adam's law というらしいのですが,この問題に使うのはちょっと大げさです。

 E(Y)=E(E(Y|X))

(サイコロを1個振ったときの目の平均値)*(サイコロを振る回数の期待値)で求めれば十分でしょう。

E(T)=(1+4)/2=5/2
E(C)=E(c)E(T)=(1+6)/2*5/2=35/4
E(O)=E(o)E(C)=315/8
E(D)=E(d)E(O)=4095/16
E(I)=E(i)E(D)=85995/32

解答

係数を全部決めてしまってハードコーディングすることも考えましたが,一応,漸化式にまとめました。
答えの厳密な値は 2464129395/1024 です。


統計学入門 (基礎統計学?)

統計学入門 (基礎統計学?)

Project Euler 147 / 斜交平行格子内の長方形

3x2 の斜め線が引かれた格子には下図のように全部で37個の異なった長方形が存在する。

f:id:variee:20170728013418g:plain

3x2 より横にも縦にも小さい5個の格子(1x1,2x1,3x1,1x2,2x2)を考えると, それらに存在する長方形の数は以下のようになる。

  • 1x1: 1
  • 2x1: 4
  • 3x1: 8
  • 1x2: 4
  • 2x2: 18

これらに 3x2 の格子の37を加えると全部で72個の長方形が存在する。47x43 以下の格子について異なる長方形の数を求めよ。

Problem 147 - Project Euler


「解いた」のではなく,調べました。code golf にそっくりの問題があります。betseg 氏のコードと Neil,Sherlock9 両氏の解説を読めばわかるでしょう。

codegolf.stackexchange.com

m*n (m≧n) の格子を考えます。
傾いていない長方形は C[m+1, 2]*C[n+1, 2] 個。
傾いた長方形は n*n の部分と残りの部分にわけて数えます。
(n-1)n(4n^2+4n+3)/6 + (m-n)n(4n^2-1)/3 個。

これらの和をとります。

Project Euler 297 / ゼッケンドルフ表記

すべての正整数はフィボナッチ数列の連続しない項の和で一意に表せる。たとえば 100=3+8+89 である。 このような表し方はゼッケンドルフ表記と呼ばれる。

整数 n>0 に対して z(n) を n のゼッケンドルフ表記中の項数とする。

z(5)=1, z(14)=2, z(100)=3

0< n< 10^6 に対して Σz(n)=7894453 である。

0< n< 10^17 に対して Σz(n) を求めよ。

Problem 297 - Project Euler

実験

「Zeckendorf mathematica」でググると,mathematica による実装が何件もヒットします。stackexchange のコードを利用しました。

mathematica.stackexchange.com

Σz(n) も求められますが,10^6までで40秒以上かかります。

z(n) のスピードアップ

n のゼッケンドルフ表記は必ず n 以下の最大のフィボナッチ数(fib(i)とする)を含みます。すべての自然数はゼッケンドルフ表記できるので,z(n) を求めるのに n 以下の最大のフィボナッチ数をどんどん引いていっても大丈夫。証明は省略しますが,「連続しない項の和」の条件も問題ありません。

 z(n)=z(n-fib(i))+1

この方法で 10^6 までの和を求めると約4秒。10倍のスピードになりましたが,10^17を相手にするのはまだ無理です。

Σz(n) のスピードアップ

n=a から n=b までの z(n) の和を S(a, b)で表します。

フィボナッチ数までの和

 S(1,\, fib(i))
 =S(1,\, fib(i-1))+S(fib(i-1)+1,\, fib(i)-1)+z\left(fib(i)\right)

第2項の添字をずらします。
fib(i-1)+1 ≦ n ≦ fib(i)-1=fib(i-1)+fib(i-2)-1 のとき

 z(n)=z(n-fib(i-1))+1

このような n の個数は fib(i-2)-1 個。添字を fib(i-1) ずらすと n の範囲は 1 ~ fib(i-2)-1 に変わります。

 S(fib(i-1)+1,\, fib(i)-1)
 =S(1,\, fib(i-2)-1)+fib(i-2)-1
 =S(1,\, fib(i-2))+fib(i-2)-2

はじめの式にもどして  z\left(fib(i)\right)=1 を使うと

 S(1,\, fib(i))
 =S(1,\, fib(i-1))+S(1,\, fib(i-2))+fib(i-2)-1

非フィボナッチ数までの和

非フィボナッチ数 N 未満の最大のフィボナッチ数を fib(i) とします。

 S(1,\, N)=S(1,\, fib(i))+S\left(fib(i)+1,\, N\right)

第2項の添字をずらします。

 S(1,\, N)=S(1,\, fib(i))+S\left(1,\, N-fib(i)\right)+N-fib(i)

本番

上で導いた2式を両方使うのがスジですが,実は片方だけですみます。
たとえば 10^17-1 のゼッケンドルフ表記は次のようになっていて,N-fib(i) がフィボナッチ数になるのは最後の2のときだけです。これは他の数でも同様です。

「非フィボナッチ数までの和」の右辺第1項の添字を1つずらすと,添字がずっと非フィボナッチ数のまま計算できます。

 S(1,\, fib(i))=S(1,\, fib(i)-1)+z\left(fib(i)\right)
 =S(1,\, fib(i)-1)+1

 \therefore S(1,\, N)
 =S(1,\, fib(i)-1)+S\left(1,\, N-fib(i)\right)+N-fib(i)+1

この式を使うと 0.003 秒で答えが出ます。