Project Euler 102 / 三角形の包含

問題概略
ベクトルの外積
面積
組み込み関数RegionMember

問題概略

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 通りの方法で解きました。

ベクトルの外積
面積
組み込み関数 RegionMember

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

drive.google.com

ベクトルの外積

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

外積その１

点 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 でも可。計算時間はほぼ同じ。