任意次元空間における三角形の外心の簡単な求め方

目的

AtCoderに外心を求める問題があって解けなかったので，怒りに任せて任意次元空間における三角形の外心の簡単な求め方と実装を書きます．

方針

三角形の3頂点を$A, B, C$とします．求める点をPとすると，とある実数$x, y$を用いて

\vec{AP} = x \vec{AB} + y \vec{AC}

と書けます．

三角形の外心は、3つの辺の垂直二等分線が交わる点である

ので(例えばWikipedia参照)，次の2つの方程式が成り立ちます．

\left(\frac{\vec{AB}}{2} - \vec{AP}\right) \cdot \vec{AB} = 0 \\
\left(\frac{\vec{AC}}{2} - \vec{AP}\right) \cdot \vec{AC} = 0

先の式を代入して解くと

x \| \vec{AB} \|^2 + y \left( \vec{AB} \cdot \vec{AC} \right) = \frac{1}{2} \| \vec{AB} \|^2 \\
x \left( \vec{AB} \cdot \vec{AC} \right) + y \| \vec{AC} \|^2 = \frac{1}{2} \| \vec{AC} \|^2 \\

行列表示すると

\left(x, y\right) \begin{pmatrix}
\| \vec{AB} \|^2 & \vec{AB} \cdot \vec{AC} \\
\vec{AB} \cdot \vec{AC} & \| \vec{AC} \|^2
\end{pmatrix}
= \frac{1}{2} \left(\| \vec{AB} \|^2, \| \vec{AC} \|^2\right)

これを解くと（逆行列かけると楽です）

\begin{pmatrix}
x \\
y
\end{pmatrix}
=
\frac{1}{2}
\frac{1}{\|\vec{AB}\|^2 \|\vec{AC}\|^2 - \left( \vec{AB} \cdot \vec{AC} \right)^2}
\begin{pmatrix}
\| \vec{AB} \|^2 \| \vec{AC} \|^2 - \| \vec{AC} \|^2 \left( \vec{AB} \cdot \vec{AC} \right) \\
\| \vec{AB} \|^2 \| \vec{AC} \|^2 - \| \vec{AB} \|^2 \left( \vec{AB} \cdot \vec{AC} \right) 
\end{pmatrix}

と，内積計算だけで解けます．

実装

vector<double> circumcenter(vector<double> P1, vector<double> P2, vector<double> P3){
      double AB2 = 0.0;
      double AC2 = 0.0;
      double ABAC = 0.0;
      double EPS = 1.0E-12;
      vector<double> ans(D, 0.0);
      for(ll i=0; i<D; ++i){
        double AB = P2[i]-P1[i];
        double AC = P3[i]-P1[i];
        AB2 += AB*AB;
        AC2 += AC*AC;
        ABAC += AB*AC;
      }

      double denominator = (AB2*AC2-ABAC*ABAC);

      if(denominator < EPS){
        fprintf(stderr, "Circumcenter Error\n");
        return ans;
      }

      double x = 0.5*(AB2*AC2-AC2*ABAC)/denominator;
      double y = 0.5*(AB2*AC2-AB2*ABAC)/denominator;

      for(ll i=0; i<D; ++i){
        double AB = P2[i]-P1[i];
        double AC = P3[i]-P1[i];
        ans[i] = x*AB+y*AC+P1[i];
      }
      return ans;
}

分母が0の時（か0に近い時）は3点$A, B, C$が一直線上にある時です．例外処理しましょう．