砂場で遊ぼう

c++/python/mathematica などの練習帳

Project Euler 184 / 原点を含む三角形

問題概略

原点を中心とする半径  r の円の内部にある格子点  (x,\, y) の集合を  I_r であらわす。
\begin{align*}
I_r=\{(x,\, y)\ |\ x^2+y^2< r^2,\, \text{$x$と$y$は整数}\}
\end{align*}

 I_2 (0,\, 0),  (\pm 1,\, 0),  (\pm 1,\, \pm 1),  (0,\, \pm 1)(複号任意)の 9 点からなる。
 I_2 の点を頂点とする三角形のうち原点を内部に含むものは 8 個ある。
そのうち 2 つを下図に示す。残りは回転で得られる。

f:id:variee:20220117000005g:plain

 I_3 の点を頂点とする三角形のうち原点を内部に含むものは 360 個あり, I_5 では 10600 個ある。

 I_{105} の点を頂点とする三角形のうち原点を内部に含むものの個数を求めよ。

https://projecteuler.net/problem=184

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

グループ分けして数える

 r=105 の円内部の格子点は約 35000個。
ここから 3 個選んで作った三角形が原点を内部に含むかどうか調べるのは現実的ではありません。
「2 頂点を固定したとき,もう 1 つの頂点が何個あるか」を考えます。

頂点 A, B を固定したとき,もう 1 つの頂点 C は図の網掛け部(境界は含まない)にあります。
直線 OA, OB 上にある他の格子点 A', B' に対しても  \triangle\mathrm{A'B'C} は原点を内部に含みます。
図のように線分上と領域内の格子点数を  a d とおきます。

f:id:variee:20220117000248p:plain

直線を固定したとき,条件をみたす三角形の個数は次のようになります。

f:id:variee:20220117000300p:plain

条件をみたす三角形の総数は次のようになります。 1/3 倍しているのはダブルカウント対策です。
\begin{align*}
\mathrm{ans}=\sum (abc\cdot 2+abd\cdot 2)\cdot \frac{1}{3}=\frac{1}{3}\sum ab(2c+2d)
\end{align*}

格子点の総数を  S とおきます。
 S=2(a+b+c+d) より  2(c+d)=S-2a-2b です。
\begin{align*}
\mathrm{ans}=\frac{1}{3}\sum ab(S-2a-2b)=\frac{1}{3}\left\{S\sum ab-2\sum ab(a+b)\right\}
\end{align*}

Σの簡約

 a,  b は適当なインデックスを使って  a=f_i,  b=f_j のように表されます。
この記号を使って  \sum の式を整理します。

 \sum ab は簡単です。
 f_i\cdot f_j の表において表全体の和から対角成分の和を引いたものの半分です。
\begin{align*}
\sum ab=\sum_{i < j} f_i f_j=\left\{\left(\sum f_i\right)^2-\sum {f_i}^2\right\}\cdot\frac{1}{2}
\end{align*}

これは  (x+y+\cdots)^2=x^2+y^2+\cdots+2(xy+\cdots をイメージするとわかりやすいでしょう。

 \sum ab(a+b) も同様です。
\begin{align*}
\sum ab(a+b) =\sum_{i< j} f_if_j(f_i+f_j)=\sum_{i< j} ({f_i}^2f_j+f_i{f_j}^2)
\end{align*}

 (x+y+\cdots)(x^2+y^2+\cdots)=x^3+y^3+\cdots+x^2y+xy^2+\cdots を利用すると次のように変形できます。
\begin{align*}
\sum ab(a+b)=\sum f_i \sum {f_i}^2-\sum {f_i}^3
\end{align*}

 \sum f_i=S_1,  \sum {f_i}^2=S_2,  \sum {f_i}^3=S_3 とおきます。
 S=2\sum f_i=2S_1 です。
\begin{align*}
\mathrm{ans} &=\frac{1}{3}\left\{2S_1({S_1}^2-S_2)\cdot\frac{1}{2}-2(S_1S_2-S_3)\right\}\\
&=\frac{1}{3}({S_1}^3-3S_1S_2+2S_3)
\end{align*}

mathematicaによる実装

mathematica で解きました。

まず格子点を求めます。
対称性を考えて  1\leqq x\leqq r-1,  0\leqq y\leqq \sqrt{r^2-1-x^2} のものだけ求めました。

  r = 105;
  latticePoint = Flatten[Table[{x, y}, {x, 1, r - 1},
     {y, 0, Floor[Sqrt[r^2 - 1 - x^2]]}], 1];

次に  f_i/2 のリストを作ります。
格子点を偏角でグループ分けして,各グループの要素数を求めました。

偏角の計算法はいろいろありますが, (x,\, y)\to x+iy のように一旦複素数に直すのが速いです。

  lst = Length /@ GatherBy[latticePoint, N@Arg[#[[1]] + #[[2]]*I] &];

 f_i/2 のリストができたら  S_1 などを求めて最終計算です。

コード全体は次のようになりました。計算時間は 0.008 秒です。

  In[]:= ClearAll["Global`*"];
  AbsoluteTiming[
   r = 105;
   latticePoint = Flatten[Table[{x, y}, {x, 1, r - 1},
      {y, 0, Floor[Sqrt[r^2 - 1 - x^2]]}], 1];
   lst = Length /@ GatherBy[latticePoint, N@Arg[#[[1]] + #[[2]]*I] &];
   S1 = Total[lst]*2;
   S2 = Total[lst^2]*2;
   S3 = Total[lst^3]*2;
   ans = (S1^3 - 3*S1*S2 + 2*S3)/3;]
  
  Out[]= {0.008169, Null}

Project Euler 227 / The Chase

問題概略

"The Chase" とは 2 つのサイコロと偶数人のプレイヤーによって行われるゲームである。

各プレイヤーはテーブルのまわりに座っている。
対面にいるある 2 人のプレイヤーがそれぞれ 1 つのサイコロをもっているところからゲームを開始する。
各ターンにサイコロをもっている 2 人のプレイヤーがそれぞれサイコロを振る。

  • 1 を出した場合,プレイヤーは自分のサイコロを左隣の人に渡す
  • 6 を出した場合は右隣の人に渡す
  • 1, 6 以外の目を出した場合はサイコロは移動しない

サイコロの移動が終わったときに 1 人のプレイヤーが 2 つのサイコロをもっていた場合,そのプレイヤーは負けとなりゲームは終了する。

100 人のプレイヤーでゲームを行ったとき,ゲームが終了するまでのターン数の期待値はいくつか? 四捨五入して小数点以下 10 桁まで答えよ。

https://projecteuler.net/problem=227

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

漸化式の立式

プレイヤーの座席に時計まわりに 1, 2, …, 100 と番号をふって,これを座標として扱います。

2 つのサイコロの距離(座標の差)ははじめ 50 で,これが 0 になるまでのターン数の期待値が答えです。

サイコロの距離が  n の状態からゲームが終了するまでのターン数の期待値を  e(n) とおいて,サイコロの座標の差が  k 変化する確率を  \mathrm{prob}(k) とおきます。

 \mathrm{prob}(k) の値は簡単に計算できます。

\begin{align*}
\mathrm{prob}(-2)=\frac{1}{36},\,
\mathrm{prob}(-1)=\frac{2}{9},\,
\mathrm{prob}(0)=\frac{1}{2},\,
\mathrm{prob}(1)=\frac{2}{9},\,
\mathrm{prob}(2)=\frac{1}{36}
\end{align*}

 e(n) については大まかには  e(n+k)\mathrm{prob}(k) の和です。

\begin{align*}
e(0)=0,\, e(n)=1+\sum_{k=-2}^2 e(n+k)\mathrm{prob}(k)\quad (n\geqq 1)
\end{align*}

実際にはプレイヤーが円状に並んでいるので少し修正が必要です。

たとえば  n=50 のときに  k=1 の目が出たら  n=50+1=51 になるのではなく, n=100-(50+1)=49 になります。

 n の状態で  k の目が出たときのサイコロの距離を  \mathrm{pos}(n+k) であらわします。

\begin{align*}
&\mathrm{pos}(n+k)=\min\{|n+k|,\, 100-(n+k)\}\\
&e(0)=0,\, e(n)=1+\sum_{k=-2}^2 e\left(\mathrm{pos}(n+k)\right)\mathrm{prob}(k)\quad (n\geqq 1)
\end{align*}

6人のとき

100 人の場合を考える前に例として 6 人のときを考えます。 n 6/2=3 以下です。

\begin{align*}
e(1) &=1+\sum_{k=-2}^2 e\left(\mathrm{pos}(1+k)\right)\mathrm{prob}(k)\\
&=1+e(\mathrm{pos}(-1))\mathrm{prob}(-2)+e(\mathrm{pos}(0))\mathrm{prob}(-1)
+e(\mathrm{pos}(1))\mathrm{prob}(0)\\
&\qquad +e(\mathrm{pos}(2))\mathrm{prob}(1)+e(\mathrm{pos}(3))\mathrm{prob}(2)\\
&=1+e(1)\cdot\frac{1}{36}+e(0)\cdot\frac{2}{9}+e(1)\cdot\frac{1}{2}
+e(2)\cdot\frac{2}{9}+e(3)\cdot\frac{1}{36}\\
&=1+\frac{19}{36}e(1)+\frac{2}{9}e(2)+\frac{1}{36}e(3)
\end{align*}

\begin{align*}
\therefore \frac{17}{36}e(1)-\frac{2}{9}e(2)-\frac{1}{36}e(3)=1
\end{align*}

 e(2),  e(3) についても同様です。

\begin{align*}
-\frac{2}{9}e(1)+\frac{17}{36}e(2)-\frac{2}{9}e(3)=1,\, -\frac{1}{18}e(1)-\frac{4}{9}e(1)e(2)+\frac{1}{1}e(3)=1
\end{align*}

これらを行列に直して解きます。

\begin{align*}
&\begin{pmatrix}
\frac{17}{36} & -\frac{2}{9} & -\frac{1}{36} \\[3pt]
-\frac{2}{9} & \frac{17}{36} & -\frac{2}{9} \\[3pt]
-\frac{1}{18} & -\frac{4}{9} & \frac{1}{2}
\end{pmatrix}
\begin{pmatrix}
e(1)\\
e(2)\\
e(3)
\end{pmatrix}
=\begin{pmatrix}
1\\
1\\
1
\end{pmatrix}\\
&\therefore \begin{pmatrix}
e(1)\\
e(2)\\
e(3)
\end{pmatrix}
=\begin{pmatrix}
\frac{17}{36} & -\frac{2}{9} & -\frac{1}{36} \\[3pt]
-\frac{2}{9} & \frac{17}{36} & -\frac{2}{9} \\[3pt]
-\frac{1}{18} & -\frac{4}{9} & \frac{1}{2}
\end{pmatrix}^{-1}
\begin{pmatrix}
1\\
1\\
1
\end{pmatrix}=
\begin{pmatrix}
\frac{419}{44}\\[3pt]
\frac{152}{1}\\[3pt]
\frac{675}{44}
\end{pmatrix}
\end{align*}

 e(3) を小数に直したものが答えになります。

mathematicaによる解答

mathematica で解きました。計算時間は約 0.02 秒です。

    In[]:= ClearAll["Global`*"];
    AbsoluteTiming[
     nmax = 100;
     pos[n_] := Min[Abs@n, nmax - Abs@n];
     prob[-2] = 1/36;
     prob[-1] = 2/9;
     prob[0] = 1/2;
     prob[1] = 2/9;
     prob[2] = 1/36;
     mat = IdentityMatrix[nmax/2];
     For[i = 1, i <= Length@mat, i++, 
      For[j = -2, j <= 2, j++, 
       If[pos[i + j] != 0, mat[[i]][[pos[i + j]]] -= prob@j]]];
     a = Transpose@ConstantArray[1, nmax/2];
     ans = N[Last[(Inverse@mat) . a], 10];]
    
    Out[]= {0.0232107, Null}

Project Euler 102 / 三角形の包含

問題概略

3 つの異なる点が  -1000 \leqq x \leqq 1000,  1000 \leqq y \leqq 1000 かつ三角形となるように平面上にランダムに与えられる。
次の 2 つの三角形を考える。
\begin{align*}
&\rm A(-340,\, 495),\, B(-153,\, -910),\, C(835,\, -947)\\
&\rm X(-175,\, 41),\, Y(-421,\, -714),\, Z(574,\, -645)
\end{align*}

三角形 ABC は原点を内部に含み,三角形 XYZ は含まない。

テキストファイル triangles.txt にランダムな 1000 個の三角形が適当なフォーマットで含まれている。内部に原点を含む三角形の個数を答えよ。
https://projecteuler.net/problem=102

3 通りの方法で解きました。

  1. ベクトルの外積
  2. 面積
  3. 組み込み関数 RegionMember

解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。

drive.google.com

ベクトルの外積

ベクトルの外積の符号を調べます。

外積その1

点 O が三角形 ABC の内部に含まれる条件は  \vec{a} \times \vec{b},\, \vec{b} \times \vec{c},\, \vec{c} \times \vec{a} z 成分の符号が等しいことです。
A, B, C が反時計回りならすべて正,時計回りならすべて負になります。

3 頂点の座標が  \{x_1,\, y_1,\, x_2,\, y_2,\, x_3,\, y_3\} の形で与えられるので

 x_1y_2-y_1x_2,\, x_2y_3-y_2x_3,\, x_3y_1-y_3x_1

の符号(Sign)が等しいものを数えました。

  In[]:= dat = Import["p102_triangles.txt", "CSV"];
  AbsoluteTiming[
   isOK[{a_, b_, c_, d_, e_, f_}] :=
    Sign[a*d - b*c] == Sign[c*f - d*e] == Sign[e*b - f*a];
   ans = Length@Select[dat, isOK];]
  
  Out[]= {0.0043344, Null}

外積その2

普通のプログラミング言語だと「符号が等しい」を「すべて正またはすべて負」で処理すると思います。
mathematica でもやってみました。計算時間はちょっとのびます。

  In[]:= dat = Import["p102_triangles.txt", "CSV"];
  AbsoluteTiming[
   isOK[{a_, b_, c_, d_, e_, f_}] := Module[{x, y, z},
     x = a*d - b*c; y = c*f - d*e; z = e*b - f*a;
     (x > 0 && y > 0 && z > 0) || (x < 0 && y < 0 && z < 0)];
   ans = Length@Select[dat, isOK];]
  
  Out[]= {0.0079219, Null}

外積その3

 \{x_1,\, y_1,\, x_2,\, y_2,\, x_3,\, y_3\} の偶数番目だけ,奇数番目だけを取り出したベクトルを考えると外積の計算が 1 回ですみます。

 (x_1,\, x_2,\, x_3)\times (y_1,\, y_2,\, y_3)=(x_2y_3-y_2x_3,\, x_3y_1-y_3x_1,\, x_1y_2-y_1x_2)

ついでに Cross で外積を計算し,Apply と Equal で値が等しいかどうか調べるように書き直したところ,かえって時間がかかる結果に……

  In[]:= dat = Import["p102_triangles.txt", "CSV"];
  AbsoluteTiming[
   isOK[{a_, b_, c_, d_, e_, f_}] :=
    Apply[Equal, Sign /@ Cross[{a, c, e}, {b, d, f}]];
   ans = Length@Select[dat, isOK];]
  
  Out[]= {0.0253116, Null}

面積

点 O が三角形 ABC の内部に含まれる条件は  \triangle \rm ABC=\triangle OAB+\triangle OBC+\triangle OCA です。点 O が外部にあるときは  \triangle \rm ABC=\triangle OAB+\triangle OBC-\triangle OCA のように右辺にマイナスがあらわれます。

Partition で位置ベクトルを作る → Simplex にわたして三角形(単体)を指定*1→ Area で面積計算,とやりました。

  In[]:= dat = Import["p102_triangles.txt", "CSV"];
  AbsoluteTiming[
   z = {0, 0};
   isOK[lst_] := Module[{a, b, c, s0, s1, s2, s3},
     {a, b, c} = Partition[lst, 2];
     s0 = Area[Simplex[{a, b, c}]]; s1 = Area[Simplex[{z, a, b}]];
     s2 = Area[Simplex[{z, b, c}]]; s3 = Area[Simplex[{z, c, a}]];
     s0 == s1 + s2 + s3];
   ans = Length@Select[dat, isOK];]
  
  Out[]= {0.0366012, Null}

組み込み関数RegionMember

この問題にうってつけの関数がありました。RegionMember[reg, {x, y}] で点  (x,\, y) が領域 reg に含まれるかどうか調べられます。
面積を使う解法と同様に単体を作って,原点を含むかどうか調べました。

  In[]:= dat = Import["p102_triangles.txt", "CSV"];
  AbsoluteTiming[
   isOK[lst_] := RegionMember[Simplex[Partition[lst, 2]], {0, 0}];
   ans = Length@Select[dat, isOK];]
  
  Out[]= {0.0305366, Null}

*1:領域指定は Polygon でも Triangle でも可。計算時間はほぼ同じ。