イマドキの方程式の求め方 三点円弧の中心座標 in 3D （修正）
と言うタイトルでさっき記事を上げたのだが、見事に間違ってました。今回は三平面の交点が三点円弧の中心になることを利用して内積を使って計算。こちらの方が式としては簡単になる。
では改めて。
Q. p1(x1,y1,z1), p2(x2,y2,z2), p3(x3,y3,z3)を通る三次元上の三点円弧の中心 c(Cx,Cy,Cz) を求める。但しp1, p2, p3は同一直線上に無く、かつ異なる点であるとする。
A. p1, p2, p3を含む平面と、線分p1p2および線分p2p3の中点を通り、それぞれの線分を法線とする二つの平面の交点が三点円弧の中心である。
三点から二本の線分を法線ベクトルを得る。
n1 = p1p2 = (x2 - x1, y2 - y1, z2 - z1)
n2 = p1p2 = (x3 - x2, y3 - y2, z3 - z2)
n1を法線ベクトルとし、その中点を通る平面の式
(x2 - x1)*(x - (x1 + x2)/2) + (y2 - y1)*(y - (y1 + y2)/2) + (z2 - z1)*(z - (z1 + z2)/2) = 0 …(1)
n2を法線ベクトルとし、その中点を通る平面の式
(x3 - x2)*(x - (x2 + x3)/2) + (y3 - y2)*(y - (y2 + y3)/2) + (z3 - z2)*(z - (z2 + z3)/2) = 0 …(2)
3点p1,p2,p3を通る平面の式
*1*(x - x1) + *2*(y - y1) + *3*(z - z1) = 0 …(3)
(1)(2)(3)の交点が円の中心である。これをMathematicaのSolve関数に打ち込んで次の解を得る。（ちょっと整形している…ここでも良く間違うんだよね。）
cx=-*4+x1*(-y2^2+y1*(y2-y3)+y2*y3 +*(z1-z2)*(z2-z3)))*(-y2^2*z1+y3^2*z1+x3^2*(z1-z2)+x1^2*z2+y1^2*z2-y3^2*z2+z1^2*z2-z1*z2^2-x1^2*z3-y1^2*z3+y2^2*z3-z1^2*z3+z2^2*z3+z1*z3^2-z2*z3^2+x2^2*(-z1+z3)) +*(y3*(-z1+z2)+y2*(z1-z3)+y1*(-z2+z3))*(-(1/ 2)*(x3*(y1-y2)+x1*(y2-y3)+x2*(-y1+y3))*(x1^2-x2^2+y1^2-y2^2+z1^2-z2^2) -*(z1-z2)*(x3*y2*z1-x2*y3*z1-x3*y1*z2+x1*y3*z2+x2*y1*z3-x1*y2*z3)))/*5**6 +*(z1-z2)*(x3*(z1-z2)+x1*(z2-z3)+x2*(-z1+z3))) +*(y3*(-z1+z2)+y2*(z1-z3)+y1*(-z2+z3))**7 +*(z1-z2)*(y3*(-z1+z2)+y2*(z1-z3)+y1*(-z2+z3)))))
cy=(x3^2*y1^3-x3^2*y1^2*y2-x3^2*y1*y2^2+x3^2*y2^3+x2^3*x3*(y1-y3)-x1^3*(x2-x3)*(y2-y3)+x3^2*y1*z1^2+y2^3*z1^2+x3^2*y3*z1^2-y2^2*y3*z1^2-y2*y3^2*z1^2+y3^3*z1^2-x3^2*y1*z1*z2-x3^2*y2*z1*z2-y1^2*y2*z1*z2-y1*y2^2*z1*z2-2*x3^2*y3*z1*z2+y1^2*y3*z1*z2+y2^2*y3*z1*z2+y1*y3^2*z1*z2+y2*y3^2*z1*z2-2*y3^3*z1*z2-y2*z1^3*z2+y3*z1^3*z2+y1^3*z2^2+x3^2*y2*z2^2+x3^2*y3*z2^2-y1^2*y3*z2^2-y1*y3^2*z2^2+y3^3*z2^2+y1*z1^2*z2^2+y2*z1^2*z2^2-2*y3*z1^2*z2^2-y1*z1*z2^3+y3*z1*z2^3-x3^2*y1*z1*z3+x3^2*y2*z1*z3+y1^2*y2*z1*z3+y1*y2^2*z1*z3-2*y2^3*z1*z3-y1^2*y3*z1*z3+y2^2*y3*z1*z3-y1*y3^2*z1*z3+y2*y3^2*z1*z3+y2*z1^3*z3-y3*z1^3*z3+x3^2*y1*z2*z3-2*y1^3*z2*z3-x3^2*y2*z2*z3+y1^2*y2*z2*z3+y1*y2^2*z2*z3+y1^2*y3*z2*z3-y2^2*y3*z2*z3+y1*y3^2*z2*z3-y2*y3^2*z2*z3-2*y1*z1^2*z2*z3+y2*z1^2*z2*z3+y3*z1^2*z2*z3+y1*z1*z2^2*z3-2*y2*z1*z2^2*z3+y3*z1*z2^2*z3+y1*z2^3*z3-y3*z2^3*z3+y1^3*z3^2-y1^2*y2*z3^2-y1*y2^2*z3^2+y2^3*z3^2+y1*z1^2*z3^2-2*y2*z1^2*z3^2+y3*z1^2*z3^2+y1*z1*z2*z3^2+y2*z1*z2*z3^2-2*y3*z1*z2*z3^2-2*y1*z2^2*z3^2+y2*z2^2*z3^2+y3*z2^2*z3^2-y1*z1*z3^3+y2*z1*z3^3+y1*z2*z3^3-y2*z2*z3^3+x1^2*(y2^3+x2^2*(y1+y2-2*y3)-y2^2*y3-y2*y3^2+y3^3+x3^2*(y1-2*y2+y3)+x2*x3*(-2*y1+y2+y3)-y2*z1*z2+y3*z1*z2+y1*z2^2+y2*z2^2+y2*z1*z3-y3*z1*z3-2*y1*z2*z3-y2*z2*z3-y3*z2*z3+y1*z3^2+y3*z3^2)+x2^2*(y1^3-y1^2*y3+y3^3+x3^2*(-2*y1+y2+y3)+y2*z1^2+y3*z1*z2-y1*(y3^2 -*(z1-z2)*(z1-z3))-2*y2*z1*z3-y3*z1*z3-y3*z2*z3+y2*z3^2+y3*z3^2)+x2*x3*(-2*y1^3+x3^2*(y1-y2)-y2^2*y3-y2*y3^2+y1^2*(y2+y3)-y2*z1^2-y3*z1^2+2*y3*z1*z2-y3*z2^2+2*y2*z1*z3-y2*z3^2+y1*(y2^2+y3^2-2*z1^2+2*z1*z2+z2^2+2*z1*z3-4*z2*z3+z3^2))+x1*(x2^3*(-y1+y3)+x2^2*x3*(y1-2*y2+y3)+x3*(-2*y2^3+x3^2*(-y1+y2)+y1^2*(y2-y3)+y2^2*y3+y2*y3^2+y2*z1^2-y3*z1^2+2*y2*z1*z2+2*y3*z1*z2-2*y2*z2^2-y3*z2^2+y1*(y2^2-y3^2 -*(z2-z3)^2)-4*y2*z1*z3+2*y2*z2*z3+y2*z3^2)+x2*(x3^2*(y1+y2-2*y3)+y2^2*y3+y2*y3^2-2*y3^3+y1^2*(-y2+y3)-y2*z1^2+y3*z1^2-4*y3*z1*z2+y3*z2^2-y1*(y2^2-y3^2 +*(z2-z3)^2)+2*y2*z1*z3+2*y3*z1*z3+2*y3*z2*z3-y2*z3^2-2*y3*z3^2)))/(2*(x1^2*y2^2-2*x1^2*y2*y3+x1^2*y3^2+y2^2*z1^2-2*y2*y3*z1^2+y3^2*z1^2+x3^2*(y1^2-2*y1*y2+y2^2 +*(z1-z2)^2)-2*y1*y2*z1*z2+2*y1*y3*z1*z2+2*y2*y3*z1*z2-2*y3^2*z1*z2+x1^2*z2^2+y1^2*z2^2-2*y1*y3*z2^2+y3^2*z2^2+x2^2*(y1^2-2*y1*y3+y3^2 +*(z1-z3)^2)-2*x1*x3*(y2^2-y2*y3+y1*(-y2+y3) -*(z1-z2)*(z2-z3))+2*y1*y2*z1*z3-2*y2^2*z1*z3-2*y1*y3*z1*z3+2*y2*y3*z1*z3-2*x1^2*z2*z3-2*y1^2*z2*z3+2*y1*y2*z2*z3+2*y1*y3*z2*z3-2*y2*y3*z2*z3+x1^2*z3^2+y1^2*z3^2-2*y1*y2*z3^2+y2^2*z3^2-2*x2*(x3*(y1^2+y2*y3-y1*(y2+y3) +*(z1-z2)*(z1-z3))+x1*(y1*(y2-y3)-y2*y3+y3^2+z1*z2-z1*z3-z2*z3+z3^2))))
cz=(x3^2*y1^2*z1-x3^2*y1*y2*z1+y1^2*y2^2*z1-y1*y2^3*z1-x3^2*y1*y3*z1+x3^2*y2*y3*z1-2*y1^2*y2*y3*z1+y1*y2^2*y3*z1+y2^3*y3*z1+y1^2*y3^2*z1+y1*y2*y3^2*z1-2*y2^2*y3^2*z1-y1*y3^3*z1+y2*y3^3*z1+x3^2*z1^3+y2^2*z1^3-2*y2*y3*z1^3+y3^2*z1^3-x3^2*y1*y2*z2-y1^3*y2*z2+x3^2*y2^2*z2+y1^2*y2^2*z2+x3^2*y1*y3*z2+y1^3*y3*z2-x3^2*y2*y3*z2+y1^2*y2*y3*z2-2*y1*y2^2*y3*z2-2*y1^2*y3^2*z2+y1*y2*y3^2*z2+y2^2*y3^2*z2+y1*y3^3*z2-y2*y3^3*z2-x3^2*z1^2*z2-y1*y2*z1^2*z2+y1*y3*z1^2*z2+y2*y3*z1^2*z2-y3^2*z1^2*z2-x3^2*z1*z2^2-y1*y2*z1*z2^2+y1*y3*z1*z2^2+y2*y3*z1*z2^2-y3^2*z1*z2^2+x3^2*z2^3+y1^2*z2^3-2*y1*y3*z2^3+y3^2*z2^3+x2^3*x3*(z1-z3)-x1^3*(x2-x3)*(z2-z3)+x3^2*y1^2*z3-2*x3^2*y1*y2*z3+y1^3*y2*z3+x3^2*y2^2*z3-2*y1^2*y2^2*z3+y1*y2^3*z3-y1^3*y3*z3+y1^2*y2*y3*z3+y1*y2^2*y3*z3-y2^3*y3*z3+y1^2*y3^2*z3-2*y1*y2*y3^2*z3+y2^2*y3^2*z3+y1*y2*z1^2*z3-y2^2*z1^2*z3-y1*y3*z1^2*z3+y2*y3*z1^2*z3-y1^2*z2^2*z3+y1*y2*z2^2*z3+y1*y3*z2^2*z3-y2*y3*z2^2*z3+y1*y2*z1*z3^2-y2^2*z1*z3^2-y1*y3*z1*z3^2+y2*y3*z1*z3^2-y1^2*z2*z3^2+y1*y2*z2*z3^2+y1*y3*z2*z3^2-y2*y3*z2*z3^2+y1^2*z3^3-2*y1*y2*z3^3+y2^2*z3^3+x1^2*(y2^2*z1-2*y2*y3*z1+y3^2*z1-y1*y2*z2+y2^2*z2+y1*y3*z2-y2*y3*z2+z2^3+x2^2*(z1+z2-2*z3)+y1*y2*z3-y1*y3*z3-y2*y3*z3+y3^2*z3-z2^2*z3-z2*z3^2+z3^3+x3^2*(z1-2*z2+z3)+x2*x3*(-2*z1+z2+z3))+x2*x3*(y2^2*z1-4*y2*y3*z1+y3^2*z1-2*z1^3+x3^2*(z1-z2)-y3^2*z2+z1^2*z2+z1*z2^2-y2^2*z3+z1^2*z3-z2^2*z3+z1*z3^2-z2*z3^2-y1^2*(2*z1+z2+z3)+2*y1*(y3*(z1+z2)+y2*(z1+z3)))+x2^2*(y2*y3*z1+z1^3+y3^2*z2+y1^2*(z1+z2)-y2*y3*z3+y3^2*z3-z1^2*z3-z1*z3^2+z3^3+x3^2*(-2*z1+z2+z3)-y1*(y2*(z1-z3)+y3*(z1+2*z2+z3)))+x1*(x2^3*(-z1+z3)+x2^2*x3*(z1-2*z2+z3)+x2*(-y3^2*z1-y1^2*z2+2*y1*y3*z2-y3^2*z2-z1^2*z2-z1*z2^2+x3^2*(z1+z2-2*z3)+y1^2*z3+2*y1*y3*z3-2*y3^2*z3+z1^2*z3+z2^2*z3+z1*z3^2+z2*z3^2-2*z3^3+y2^2*(-z1+z3)+2*y2*(-2*y1*z3+y3*(z1+z3)))+x3*(-y3^2*z1+y1^2*z2-4*y1*y3*z2+y3^2*z2+z1^2*z2+z1*z2^2-2*z2^3+x3^2*(-z1+z2)-y1^2*z3-z1^2*z3+z2^2*z3-z1*z3^2+z2*z3^2-y2^2*(z1+2*z2+z3)+2*y2*(y3*(z1+z2)+y1*(z2+z3)))))/(2*(x1^2*y2^2-2*x1^2*y2*y3+x1^2*y3^2+y2^2*z1^2-2*y2*y3*z1^2+y3^2*z1^2+x3^2*(y1^2-2*y1*y2+y2^2 +*(z1-z2)^2)-2*y1*y2*z1*z2+2*y1*y3*z1*z2+2*y2*y3*z1*z2-2*y3^2*z1*z2+x1^2*z2^2+y1^2*z2^2-2*y1*y3*z2^2+y3^2*z2^2+x2^2*(y1^2-2*y1*y3+y3^2 +*(z1-z3)^2)-2*x1*x3*(y2^2-y2*y3+y1*(-y2+y3) -*(z1-z2)*(z2-z3))+2*y1*y2*z1*z3-2*y2^2*z1*z3-2*y1*y3*z1*z3+2*y2*y3*z1*z3-2*x1^2*z2*z3-2*y1^2*z2*z3+2*y1*y2*z2*z3+2*y1*y3*z2*z3-2*y2*y3*z2*z3+x1^2*z3^2+y1^2*z3^2-2*y1*y2*z3^2+y2^2*z3^2-2*x2*(x3*(y1^2+y2*y3-y1*(y2+y3) +*(z1-z2)*(z1-z3))+x1*(y1*(y2-y3)-y2*y3+y3^2+z1*z2-z1*z3-z2*z3+z3^2))))
（簡約化の下処理のために後日ちょくちょく変更）
高速化のためにはできるだけ共通項を括り出して簡約化するものだが、今のところ予定なし。誰かC#あたりのプログラム化して教えて下さい。
但し、誰かがこの記事の内容を使って何らかのトラブルに巻き込まれても当方は一切関知しないので念のため。
ささやかに思うことを述べさせていただければ、数式解法を覚えるのは無論OKなんだけど、パンピーはそれをどのように使うのか、どのような時に適用できるのか、について広く応用力を付ける教育に変遷していくべきじゃないかな。ラズパイのMathematica使ったけど、入力一時間（文法忘れてたので）計算数秒だよ。
後日追記：
んーと、割り算があるから、プログラムするならゼロ割り算の可能性を除外しないとこの「公式」は使えない。MathematicaのSolve関数の解は漸化式なので、実際に値をいきなり入れ込むとこんなことが起こる。具体的には直線P1-P2もしくはP2-P3が各座標軸と平行になる＝傾きがゼロになると発生するだろう。Wolframエンジンならちゃんと解は出るが、普通のプログラミング言語では面倒な場合分けが必要になる。
と言うことで長い上記の式を記憶しても（笑）、あまり有用ではない。二次元空間にマッピングして、さらにそこでゼロ除算の可能性を排除して…とやり始めると実際に値を入れて計算するのは楽だなーと感じる（やったけど）。でも世の中、その楽な道の解説ばかりだったので、この記事を書いた。
では改めて。
Q. p1(x1,y1,z1), p2(x2,y2,z2), p3(x3,y3,z3)を通る三次元上の三点円弧の中心 c(Cx,Cy,Cz) を求める。但しp1, p2, p3は同一直線上に無く、かつ異なる点であるとする。
A. p1, p2, p3を含む平面と、線分p1p2および線分p2p3の中点を通り、それぞれの線分を法線とする二つの平面の交点が三点円弧の中心である。
三点から二本の線分を法線ベクトルを得る。
n1 = p1p2 = (x2 - x1, y2 - y1, z2 - z1)
n2 = p1p2 = (x3 - x2, y3 - y2, z3 - z2)
n1を法線ベクトルとし、その中点を通る平面の式
(x2 - x1)*(x - (x1 + x2)/2) + (y2 - y1)*(y - (y1 + y2)/2) + (z2 - z1)*(z - (z1 + z2)/2) = 0 …(1)
n2を法線ベクトルとし、その中点を通る平面の式
(x3 - x2)*(x - (x2 + x3)/2) + (y3 - y2)*(y - (y2 + y3)/2) + (z3 - z2)*(z - (z2 + z3)/2) = 0 …(2)
3点p1,p2,p3を通る平面の式
*1*(x - x1) + *2*(y - y1) + *3*(z - z1) = 0 …(3)
(1)(2)(3)の交点が円の中心である。これをMathematicaのSolve関数に打ち込んで次の解を得る。（ちょっと整形している…ここでも良く間違うんだよね。）
cx=-*4+x1*(-y2^2+y1*(y2-y3)+y2*y3 +*(z1-z2)*(z2-z3)))*(-y2^2*z1+y3^2*z1+x3^2*(z1-z2)+x1^2*z2+y1^2*z2-y3^2*z2+z1^2*z2-z1*z2^2-x1^2*z3-y1^2*z3+y2^2*z3-z1^2*z3+z2^2*z3+z1*z3^2-z2*z3^2+x2^2*(-z1+z3)) +*(y3*(-z1+z2)+y2*(z1-z3)+y1*(-z2+z3))*(-(1/ 2)*(x3*(y1-y2)+x1*(y2-y3)+x2*(-y1+y3))*(x1^2-x2^2+y1^2-y2^2+z1^2-z2^2) -*(z1-z2)*(x3*y2*z1-x2*y3*z1-x3*y1*z2+x1*y3*z2+x2*y1*z3-x1*y2*z3)))/*5**6 +*(z1-z2)*(x3*(z1-z2)+x1*(z2-z3)+x2*(-z1+z3))) +*(y3*(-z1+z2)+y2*(z1-z3)+y1*(-z2+z3))**7 +*(z1-z2)*(y3*(-z1+z2)+y2*(z1-z3)+y1*(-z2+z3)))))
cy=(x3^2*y1^3-x3^2*y1^2*y2-x3^2*y1*y2^2+x3^2*y2^3+x2^3*x3*(y1-y3)-x1^3*(x2-x3)*(y2-y3)+x3^2*y1*z1^2+y2^3*z1^2+x3^2*y3*z1^2-y2^2*y3*z1^2-y2*y3^2*z1^2+y3^3*z1^2-x3^2*y1*z1*z2-x3^2*y2*z1*z2-y1^2*y2*z1*z2-y1*y2^2*z1*z2-2*x3^2*y3*z1*z2+y1^2*y3*z1*z2+y2^2*y3*z1*z2+y1*y3^2*z1*z2+y2*y3^2*z1*z2-2*y3^3*z1*z2-y2*z1^3*z2+y3*z1^3*z2+y1^3*z2^2+x3^2*y2*z2^2+x3^2*y3*z2^2-y1^2*y3*z2^2-y1*y3^2*z2^2+y3^3*z2^2+y1*z1^2*z2^2+y2*z1^2*z2^2-2*y3*z1^2*z2^2-y1*z1*z2^3+y3*z1*z2^3-x3^2*y1*z1*z3+x3^2*y2*z1*z3+y1^2*y2*z1*z3+y1*y2^2*z1*z3-2*y2^3*z1*z3-y1^2*y3*z1*z3+y2^2*y3*z1*z3-y1*y3^2*z1*z3+y2*y3^2*z1*z3+y2*z1^3*z3-y3*z1^3*z3+x3^2*y1*z2*z3-2*y1^3*z2*z3-x3^2*y2*z2*z3+y1^2*y2*z2*z3+y1*y2^2*z2*z3+y1^2*y3*z2*z3-y2^2*y3*z2*z3+y1*y3^2*z2*z3-y2*y3^2*z2*z3-2*y1*z1^2*z2*z3+y2*z1^2*z2*z3+y3*z1^2*z2*z3+y1*z1*z2^2*z3-2*y2*z1*z2^2*z3+y3*z1*z2^2*z3+y1*z2^3*z3-y3*z2^3*z3+y1^3*z3^2-y1^2*y2*z3^2-y1*y2^2*z3^2+y2^3*z3^2+y1*z1^2*z3^2-2*y2*z1^2*z3^2+y3*z1^2*z3^2+y1*z1*z2*z3^2+y2*z1*z2*z3^2-2*y3*z1*z2*z3^2-2*y1*z2^2*z3^2+y2*z2^2*z3^2+y3*z2^2*z3^2-y1*z1*z3^3+y2*z1*z3^3+y1*z2*z3^3-y2*z2*z3^3+x1^2*(y2^3+x2^2*(y1+y2-2*y3)-y2^2*y3-y2*y3^2+y3^3+x3^2*(y1-2*y2+y3)+x2*x3*(-2*y1+y2+y3)-y2*z1*z2+y3*z1*z2+y1*z2^2+y2*z2^2+y2*z1*z3-y3*z1*z3-2*y1*z2*z3-y2*z2*z3-y3*z2*z3+y1*z3^2+y3*z3^2)+x2^2*(y1^3-y1^2*y3+y3^3+x3^2*(-2*y1+y2+y3)+y2*z1^2+y3*z1*z2-y1*(y3^2 -*(z1-z2)*(z1-z3))-2*y2*z1*z3-y3*z1*z3-y3*z2*z3+y2*z3^2+y3*z3^2)+x2*x3*(-2*y1^3+x3^2*(y1-y2)-y2^2*y3-y2*y3^2+y1^2*(y2+y3)-y2*z1^2-y3*z1^2+2*y3*z1*z2-y3*z2^2+2*y2*z1*z3-y2*z3^2+y1*(y2^2+y3^2-2*z1^2+2*z1*z2+z2^2+2*z1*z3-4*z2*z3+z3^2))+x1*(x2^3*(-y1+y3)+x2^2*x3*(y1-2*y2+y3)+x3*(-2*y2^3+x3^2*(-y1+y2)+y1^2*(y2-y3)+y2^2*y3+y2*y3^2+y2*z1^2-y3*z1^2+2*y2*z1*z2+2*y3*z1*z2-2*y2*z2^2-y3*z2^2+y1*(y2^2-y3^2 -*(z2-z3)^2)-4*y2*z1*z3+2*y2*z2*z3+y2*z3^2)+x2*(x3^2*(y1+y2-2*y3)+y2^2*y3+y2*y3^2-2*y3^3+y1^2*(-y2+y3)-y2*z1^2+y3*z1^2-4*y3*z1*z2+y3*z2^2-y1*(y2^2-y3^2 +*(z2-z3)^2)+2*y2*z1*z3+2*y3*z1*z3+2*y3*z2*z3-y2*z3^2-2*y3*z3^2)))/(2*(x1^2*y2^2-2*x1^2*y2*y3+x1^2*y3^2+y2^2*z1^2-2*y2*y3*z1^2+y3^2*z1^2+x3^2*(y1^2-2*y1*y2+y2^2 +*(z1-z2)^2)-2*y1*y2*z1*z2+2*y1*y3*z1*z2+2*y2*y3*z1*z2-2*y3^2*z1*z2+x1^2*z2^2+y1^2*z2^2-2*y1*y3*z2^2+y3^2*z2^2+x2^2*(y1^2-2*y1*y3+y3^2 +*(z1-z3)^2)-2*x1*x3*(y2^2-y2*y3+y1*(-y2+y3) -*(z1-z2)*(z2-z3))+2*y1*y2*z1*z3-2*y2^2*z1*z3-2*y1*y3*z1*z3+2*y2*y3*z1*z3-2*x1^2*z2*z3-2*y1^2*z2*z3+2*y1*y2*z2*z3+2*y1*y3*z2*z3-2*y2*y3*z2*z3+x1^2*z3^2+y1^2*z3^2-2*y1*y2*z3^2+y2^2*z3^2-2*x2*(x3*(y1^2+y2*y3-y1*(y2+y3) +*(z1-z2)*(z1-z3))+x1*(y1*(y2-y3)-y2*y3+y3^2+z1*z2-z1*z3-z2*z3+z3^2))))
cz=(x3^2*y1^2*z1-x3^2*y1*y2*z1+y1^2*y2^2*z1-y1*y2^3*z1-x3^2*y1*y3*z1+x3^2*y2*y3*z1-2*y1^2*y2*y3*z1+y1*y2^2*y3*z1+y2^3*y3*z1+y1^2*y3^2*z1+y1*y2*y3^2*z1-2*y2^2*y3^2*z1-y1*y3^3*z1+y2*y3^3*z1+x3^2*z1^3+y2^2*z1^3-2*y2*y3*z1^3+y3^2*z1^3-x3^2*y1*y2*z2-y1^3*y2*z2+x3^2*y2^2*z2+y1^2*y2^2*z2+x3^2*y1*y3*z2+y1^3*y3*z2-x3^2*y2*y3*z2+y1^2*y2*y3*z2-2*y1*y2^2*y3*z2-2*y1^2*y3^2*z2+y1*y2*y3^2*z2+y2^2*y3^2*z2+y1*y3^3*z2-y2*y3^3*z2-x3^2*z1^2*z2-y1*y2*z1^2*z2+y1*y3*z1^2*z2+y2*y3*z1^2*z2-y3^2*z1^2*z2-x3^2*z1*z2^2-y1*y2*z1*z2^2+y1*y3*z1*z2^2+y2*y3*z1*z2^2-y3^2*z1*z2^2+x3^2*z2^3+y1^2*z2^3-2*y1*y3*z2^3+y3^2*z2^3+x2^3*x3*(z1-z3)-x1^3*(x2-x3)*(z2-z3)+x3^2*y1^2*z3-2*x3^2*y1*y2*z3+y1^3*y2*z3+x3^2*y2^2*z3-2*y1^2*y2^2*z3+y1*y2^3*z3-y1^3*y3*z3+y1^2*y2*y3*z3+y1*y2^2*y3*z3-y2^3*y3*z3+y1^2*y3^2*z3-2*y1*y2*y3^2*z3+y2^2*y3^2*z3+y1*y2*z1^2*z3-y2^2*z1^2*z3-y1*y3*z1^2*z3+y2*y3*z1^2*z3-y1^2*z2^2*z3+y1*y2*z2^2*z3+y1*y3*z2^2*z3-y2*y3*z2^2*z3+y1*y2*z1*z3^2-y2^2*z1*z3^2-y1*y3*z1*z3^2+y2*y3*z1*z3^2-y1^2*z2*z3^2+y1*y2*z2*z3^2+y1*y3*z2*z3^2-y2*y3*z2*z3^2+y1^2*z3^3-2*y1*y2*z3^3+y2^2*z3^3+x1^2*(y2^2*z1-2*y2*y3*z1+y3^2*z1-y1*y2*z2+y2^2*z2+y1*y3*z2-y2*y3*z2+z2^3+x2^2*(z1+z2-2*z3)+y1*y2*z3-y1*y3*z3-y2*y3*z3+y3^2*z3-z2^2*z3-z2*z3^2+z3^3+x3^2*(z1-2*z2+z3)+x2*x3*(-2*z1+z2+z3))+x2*x3*(y2^2*z1-4*y2*y3*z1+y3^2*z1-2*z1^3+x3^2*(z1-z2)-y3^2*z2+z1^2*z2+z1*z2^2-y2^2*z3+z1^2*z3-z2^2*z3+z1*z3^2-z2*z3^2-y1^2*(2*z1+z2+z3)+2*y1*(y3*(z1+z2)+y2*(z1+z3)))+x2^2*(y2*y3*z1+z1^3+y3^2*z2+y1^2*(z1+z2)-y2*y3*z3+y3^2*z3-z1^2*z3-z1*z3^2+z3^3+x3^2*(-2*z1+z2+z3)-y1*(y2*(z1-z3)+y3*(z1+2*z2+z3)))+x1*(x2^3*(-z1+z3)+x2^2*x3*(z1-2*z2+z3)+x2*(-y3^2*z1-y1^2*z2+2*y1*y3*z2-y3^2*z2-z1^2*z2-z1*z2^2+x3^2*(z1+z2-2*z3)+y1^2*z3+2*y1*y3*z3-2*y3^2*z3+z1^2*z3+z2^2*z3+z1*z3^2+z2*z3^2-2*z3^3+y2^2*(-z1+z3)+2*y2*(-2*y1*z3+y3*(z1+z3)))+x3*(-y3^2*z1+y1^2*z2-4*y1*y3*z2+y3^2*z2+z1^2*z2+z1*z2^2-2*z2^3+x3^2*(-z1+z2)-y1^2*z3-z1^2*z3+z2^2*z3-z1*z3^2+z2*z3^2-y2^2*(z1+2*z2+z3)+2*y2*(y3*(z1+z2)+y1*(z2+z3)))))/(2*(x1^2*y2^2-2*x1^2*y2*y3+x1^2*y3^2+y2^2*z1^2-2*y2*y3*z1^2+y3^2*z1^2+x3^2*(y1^2-2*y1*y2+y2^2 +*(z1-z2)^2)-2*y1*y2*z1*z2+2*y1*y3*z1*z2+2*y2*y3*z1*z2-2*y3^2*z1*z2+x1^2*z2^2+y1^2*z2^2-2*y1*y3*z2^2+y3^2*z2^2+x2^2*(y1^2-2*y1*y3+y3^2 +*(z1-z3)^2)-2*x1*x3*(y2^2-y2*y3+y1*(-y2+y3) -*(z1-z2)*(z2-z3))+2*y1*y2*z1*z3-2*y2^2*z1*z3-2*y1*y3*z1*z3+2*y2*y3*z1*z3-2*x1^2*z2*z3-2*y1^2*z2*z3+2*y1*y2*z2*z3+2*y1*y3*z2*z3-2*y2*y3*z2*z3+x1^2*z3^2+y1^2*z3^2-2*y1*y2*z3^2+y2^2*z3^2-2*x2*(x3*(y1^2+y2*y3-y1*(y2+y3) +*(z1-z2)*(z1-z3))+x1*(y1*(y2-y3)-y2*y3+y3^2+z1*z2-z1*z3-z2*z3+z3^2))))
（簡約化の下処理のために後日ちょくちょく変更）
高速化のためにはできるだけ共通項を括り出して簡約化するものだが、今のところ予定なし。誰かC#あたりのプログラム化して教えて下さい。
但し、誰かがこの記事の内容を使って何らかのトラブルに巻き込まれても当方は一切関知しないので念のため。
ささやかに思うことを述べさせていただければ、数式解法を覚えるのは無論OKなんだけど、パンピーはそれをどのように使うのか、どのような時に適用できるのか、について広く応用力を付ける教育に変遷していくべきじゃないかな。ラズパイのMathematica使ったけど、入力一時間（文法忘れてたので）計算数秒だよ。
後日追記：
んーと、割り算があるから、プログラムするならゼロ割り算の可能性を除外しないとこの「公式」は使えない。MathematicaのSolve関数の解は漸化式なので、実際に値をいきなり入れ込むとこんなことが起こる。具体的には直線P1-P2もしくはP2-P3が各座標軸と平行になる＝傾きがゼロになると発生するだろう。Wolframエンジンならちゃんと解は出るが、普通のプログラミング言語では面倒な場合分けが必要になる。
と言うことで長い上記の式を記憶しても（笑）、あまり有用ではない。二次元空間にマッピングして、さらにそこでゼロ除算の可能性を排除して…とやり始めると実際に値を入れて計算するのは楽だなーと感じる（やったけど）。でも世の中、その楽な道の解説ばかりだったので、この記事を書いた。
*1:y2 - y1)*(z3 - z1) - (y3 - y1)*(z2 - z1
*2:z2 - z1)*(x3 - x1) - (z3 - z1)*(x2 - x1
*3:x2 - x1)*(y3 - y1) - (x3 - x1)*(y2 - y1
*4:-(1/2)*(x3*(y1^2-2*y1*y2+y2^2 +*(z1-z2)^2)-x2*(y1^2+y2*y3-y1*(y2+y3) +*(z1-z2)*(z1-z3
*5:x3*(z1-z2)+x1*(z2-z3)+x2*(-z1+z3
*6:y1-y2)*(x3*(y1-y2)+x1*(y2-y3)+x2*(-y1+y3
*7:x1-x2)*(x3*(y1-y2)+x1*(y2-y3)+x2*(-y1+y3
