砂場で遊ぼう

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

トップ > project euler > Project Euler 184 / 原点を含む三角形

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

project euler mathematica

問題概略

原点を中心とする半径 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}