Project Euler 143 / 三角形のトリチェリ点

問題概略
余弦定理→整数のパラメータ表示
整数の組を数える
- FindCycle で長さ 3 の閉路を探す
- 隣接リストを使う
最終結果

問題概略

ABC をすべての内角が 120 度未満の三角形とする。
X を三角形の内点として $\mathrm{XA}=p$ , $\mathrm{XC}=q$ , $\mathrm{XB}=r$ とおく。
フェルマーは「 $p+q+r$ を最小にする X を探す方法はあるか？」とトリチェリに問題を出した。
正三角形 AOB, BNC, AMC が ABC の各辺に構成できるならば，AOB, BNC, AMC に外接する 3 つの円が三角形の内部の 1 点 T で交わることをトリチェリは示した。
さらに，フェルマー点（トリチェリ点，等角中心）と呼ばれる点 T が $p+q+r$ を最小にすることも示した。
$p+q+r$ が最小のとき $\mathrm{AN}=\mathrm{BM}=\mathrm{CO}=p+q+r$ であり，AN, BM, CO が点 T を通ることも示せる。
$p+q+r$ が最小で $a$ , $b$ , $c$ , $p$ , $q$ , $r$ がすべて自然数のとき，三角形 ABC をトリチェリ三角形と呼ぶ。
たとえば $a=399$ , $b=455$ , $c=511$ は $p+q+r=784$ のトリチェリ三角形である。
トリチェリ三角形について，異なる値をとる $p+q+r\ (\leqq 120{,}000)$ の総和を求めよ。
https://projecteuler.net/problem=143

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

drive.google.com

余弦定理→整数のパラメータ表示

与えられた図において点 T のまわりの角はすべて 120 度です。
$\triangle\mathrm{TAB}$ , $\triangle\mathrm{TBC}$ , $\triangle\mathrm{TCA}$ に余弦定理を使います。

$q^2+qr+r^2=a^2,\, p^2+pq+q^2=b^2,\, r^2+rp+p^2=c^2$

これをみたす自然数はピタゴラス数の一般項 $\left(k(m^2-n^2),\, 2kmn,\, k(m^2+n^2)\right)$ と同じ要領でパラメータ表示することができます。

$q^2+qr+r^2=a^2$ を $a^2$ で割って $q/a=x$ , $r/a=y$ とおきます。 $x$ , $y$ は有理数です。

$x^2+xy+y^2=1$

この楕円上の格子点として $(\pm 1,\, 0)$ と $(0,\, \pm 1)$ がすぐみつかります。

$(-1,\, 0)$ を通る直線として $y=t(x+1)$ を考えて楕円との交点を求めると， $(-1,\, 0)$ でない方の座標は次のようになります。

$(x,\, y)=\left(\frac{1-t^2}{1+t+t^2},\, \frac{2t+t^2}{1+t+t^2}\right)$

$t=n/m$ （ $m$ , $n$ は互いに素な自然数）とおきます。

$(x,\, y)=\left(\frac{q}{a},\, \frac{r}{a}\right) =\left(\frac{m^2-n^2}{m^2+mn+n^2},\, \frac{2mn+n^2}{m^2+mn+n^2}\right)$

$(q,\, r,\, a)=(m^2-n^2,\, 2mn+n^2,\, m^2+mn+n^2)$ がみつかりました。

次はこれらが互いに素な自然数になる条件を求めます。

「 $m$ と $n$ は互いに素」と $q>0$ からわかる「 $m>n$ 」の他に条件がないか mathematica で調べてみました。

  In[]:= ClearAll["Global`*"];
    q[m_, n_] := m^2 - n^2;
    r[m_, n_] := 2 m*n + n^2;
    isBad[{m_, n_}] := 
      m > n && CoprimeQ[m, n] && ! CoprimeQ[q[m, n], r[m, n]];
    Select[Tuples[Range@10, 2], isBad]
  
  Out[]= {{4, 1}, {5, 2}, {7, 1}, {7, 4}, {8, 5}, {10, 1}, {10, 7}}

$m-n$ が 3 の倍数だとまずいみたいですね。これを証明します。

$q$ , $r$ の最大公約数を $\langle q,\, r \rangle$ のようにあらわします。

$\langle q,\, r\rangle=\langle m^2-n^2,\, 2mn+n^2\rangle =\langle m^2+2mn,\, 2mn+n^2\rangle=\langle m(m+2n),\, n(2m+n)\rangle$

$m$ と $n$ は互いに素なので「 $m$ と $n$ , $2m+n$ は互いに素」と「 $n$ と $m$ , $m+2n$ は互いに素」がいえます。
上式の最右辺の値は $\langle m+2n,\, 2m+n \rangle$ と同じです。

$\langle q,\, r\rangle=\langle m+2n,\, 2m+n \rangle$

$M=m+2n$ , $N=2m+n$ とおくと $3m=-M+2N$ , $3n=2M-N$ がいえます。
$M$ と $N$ の最大公約数は $3m$ と $3n$ の公約数です。1 か 3 ですね。

この 3 を否定すればいいので， $q$ と $r$ が互いに素な条件は $M=m+2n$ , $N=2m+n$ の少なくとも一方が 3 で割りきれないことです。
casework により「 $m\not \equiv n\pmod{3}$ 」がわかります。

まとめましょう。

$q^2+qr+r^2$ が平方数になるのは $q=k(m^2-n^2)$ , $r=k(2mn+n^2)$ のとき
$m$ と $n$ は互いに素な自然数で $m>n$ かつ $m\not \equiv n\pmod{3}$ をみたす
$k$ は自然数

整数の組を数える

FindCycle で長さ 3 の閉路を探す

上の条件をみたす 2 つの数を辺で結んでグラフを作ります。

「 $p$ と $q$ 」「 $q$ と $r$ 」「 $r$ と $p$ 」が平方数を与えるとき， $p$ , $q$ , $r$ は長さ 3 の閉路を作ります。

mathematica の組み込み関数 FindCycle でこれを探して， $p+q+r\leqq 12\cdot 10^4$ をみたすものの和をとると解けるはず……と思いましたが実際には時間がかかりすぎてダメですね。
「 $12\cdot 10^3$ 以下」は 15 秒程度で答えが出ますが，「 $12\cdot 10^4$ 以下」だと固まってしまいます。

  In[]:= AbsoluteTiming[
    lim = 12*10^3;
    g = {};
    For[m = 2, m^2 + 2 m < lim, m++,
     For[n = 1, n < m, n++,
      If[! CoprimeQ[m, n] || Divisible[m - n, 3], Continue];
      p = m^2 - n^2;
      q = 2*m*n + n^2;
      Do[AppendTo[g, k*p <-> k*q], {k, 1, Quotient[lim, p + q]}]]];
    ans = Total@Select[DeleteDuplicates[
        Total /@ VertexList /@ FindCycle[g, {3}, All]], # < lim &]]
   
  Out[]= {15.2798, 251752}

遅い原因は FindCycle と AppendTo です。

FindCycle ではデフォルトで閉路を 1 つ探すようになっている。閉路を全部探すのは現実的ではない
リストに新しい要素を追加する AppendTo は遅い

stackoverflow.com

隣接リストを使う

グラフを無向辺の集合として作るのをやめて，隣接リストで管理します。

隣接リストは小さい方から大きい方に辺を張って作ります。
$p$ , $q$ , $r$ $(p< q< r)$ がループを作っていることの確認は「 $p\to q$ , $p\to r$ 」→「 $q\to r$ 」の 2 段階に分けて行います。

$p$ の隣接リストから大きさ 2 の部分集合 $x,\, y\}\ (x< y)$ を作る。（辺 $p\to x$ , $p\to y$ がある）
$x$ の隣接リストに $y$ が含まれていれば OK（辺 $x\to y$ がある）

他にもいくつか工夫しました。

ループさせない

mathematica はループ処理が苦手です。 $(m,\, n)$ のリストを Table で一気に作りました。
要素数は約 15,000 です。この変更で計算時間はほとんど変わりませんでした。

  lim = 12*10^4;
  mnTable = Flatten[Table[{m, n}, {m, 2, Sqrt[1 + lim] - 1},
     {n, 1, Min[m, Quotient[lim - m^2, 2 m]]}], 1];
  mnTable = 
   Select[mnTable, 
    CoprimeQ[#[[1]], #[[2]]] && ! Divisible[#[[1]] - #[[2]], 3] &];

$m$ , $n$ の範囲は $p+q=(m^2-n^2)+(2mn+n^2)=m^2+2mn$ が $\mathrm{lim}=12{,}000$ を下回るように決めました。

$n=1$ のときを考えると $m^2+2m< \mathrm{lim}$ より

$2\leqq m<\sqrt{1+\mathrm{lim}}-1$

$n$ の範囲も $m^2+2mn< \mathrm{lim}$ から決めました。

$1\leqq n<\min\left\{n,\, \frac{\mathrm{lim}-m^2}{2m}\right\}$

隣接リストを FixedArray で作る

隣接リストは mathematica version 12.1 で導入された FixedArray で作りました。

  adjacencyList = CreateDataStructure["FixedArray", {}, lim];

公式いわく「要素は素早く更新できる」らしいです。
たしかに adjacencyList = ConstantArray[ {}, lim] のように普通の多次元配列を使うのと比べるとちょっと速かったですね。

隣接リストへの要素の追加はまとめて行う

隣接リストに頂点を 1 つ 1 つ追加するかわりに $(p,\, q)$ の倍数のリスト $\left\{(p,\,q ),\, (2p,\,2q),\, \cdots\right\}$ をまとめて追加するようにしました。

  makeList[{m_, n_}] := Module[{p, q, lst},
    {p, q} = MinMax[{m^2 - n^2, 2 m*n + n^2}];
    lst = Times[{p, q}, #] & /@ Range@Quotient[lim, p + q];
    (adjacencyList["Part", #[[1]]] 
      = {adjacencyList["Part", #[[1]]], #[[2]]}) & /@ lst];

  makeList /@ mnTable;

最小の $(p,\, q)$ は $3,\, 5)$ で，これから作られる $(kp,\, kq)$ は $120{,}000/(3+5)=15{,}000$ 個あります。
これらの登録が 1 回で終わるようにしたわけです。

要素数 1 の隣接リストは調べない

$p$ の隣接リストの要素数が 1 だったら「 $p\to q$ , $p\to r$ 」をみたす $q$ , $r$ が存在しないので，こういう $p$ は除外して考えました。

ちなみに要素数 1 の隣接リストは約 13,000 個で，要素数 2 個以上のものは約 22,000 個でした。

  vertex = Flatten@Position[adjacencyList, {Repeated[_, {2, Infinity}]}];

最終結果

以上をふまえて計算したところ約 2.5 秒でした。

  In[]:= ClearAll["Global`*"];
  AbsoluteTiming[
   lim = 12*10^4;
   mnTable = Flatten[Table[{m, n}, {m, 2, Sqrt[1 + lim] - 1},
      {n, 1, Min[m, Quotient[Quotient[lim, m] - m, 2]]}], 1];
   mnTable = 
    Select[mnTable, 
     CoprimeQ[#[[1]], #[[2]]] && ! Divisible[#[[1]] - #[[2]], 3] &];
   adjacencyList = CreateDataStructure["FixedArray", {}, lim];
   makeList[{m_, n_}] := Module[{p, q, lst},
     {p, q} = MinMax[{m^2 - n^2, 2 m*n + n^2}];
     lst = Times[{p, q}, #] & /@ Range@Quotient[lim, p + q];
     (adjacencyList[
          "Part", #[[1]]] = {adjacencyList["Part", #[[1]]], #[[2]]}) & /@
       lst];
   makeList /@ mnTable;
   adjacencyList = Union@Flatten[#, Infinity] & /@ Normal@adjacencyList;
   vertex = 
    Flatten@Position[adjacencyList, {Repeated[_, {2, Infinity}]}];
   calc[p_] := Module[{nums},
     nums = Select[Subsets[adjacencyList[[p]], {2}],
       MemberQ[adjacencyList[[#[[1]]]], #[[2]]] &];
     p + Total /@ nums];
   ans = Total[
     Select[DeleteDuplicates@Flatten[calc /@ vertex], # < lim &]];]
  
  Out[]= {2.51049, Null}

砂場で遊ぼう

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