問題概略
原点を中心とする半径 の円の内部にある格子点 の集合を であらわす。
\begin{align*}
I_r=\{(x,\, y)\ |\ x^2+y^2< r^2,\, \text{$x$と$y$は整数}\}
\end{align*}は , , , (複号任意)の 9 点からなる。
の点を頂点とする三角形のうち原点を内部に含むものは 8 個ある。
そのうち 2 つを下図に示す。残りは回転で得られる。の点を頂点とする三角形のうち原点を内部に含むものは 360 個あり, では 10600 個ある。
の点を頂点とする三角形のうち原点を内部に含むものの個数を求めよ。
解説の pdf も作りました。きれいなレイアウトで読みたい方はこちらをどうぞ。
グループ分けして数える
の円内部の格子点は約 35000個。
ここから 3 個選んで作った三角形が原点を内部に含むかどうか調べるのは現実的ではありません。
「2 頂点を固定したとき,もう 1 つの頂点が何個あるか」を考えます。
頂点 A, B を固定したとき,もう 1 つの頂点 C は図の網掛け部(境界は含まない)にあります。
直線 OA, OB 上にある他の格子点 A', B' に対しても は原点を内部に含みます。
図のように線分上と領域内の格子点数を ~ とおきます。
直線を固定したとき,条件をみたす三角形の個数は次のようになります。
条件をみたす三角形の総数は次のようになります。 倍しているのはダブルカウント対策です。
\begin{align*}
\mathrm{ans}=\sum (abc\cdot 2+abd\cdot 2)\cdot \frac{1}{3}=\frac{1}{3}\sum ab(2c+2d)
\end{align*}
格子点の総数を とおきます。
より です。
\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*}
Σの簡約
, は適当なインデックスを使って , のように表されます。
この記号を使って の式を整理します。
は簡単です。
の表において表全体の和から対角成分の和を引いたものの半分です。
\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*}
これは をイメージするとわかりやすいでしょう。
も同様です。
\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*}
を利用すると次のように変形できます。
\begin{align*}
\sum ab(a+b)=\sum f_i \sum {f_i}^2-\sum {f_i}^3
\end{align*}
, , とおきます。
です。
\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 で解きました。
まず格子点を求めます。
対称性を考えて , のものだけ求めました。
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] &];
のリストができたら などを求めて最終計算です。
コード全体は次のようになりました。計算時間は 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}