三维空间中与三个通用球面相切的8个平面方程的紧凑求解算法咨询

嗨，我完全懂你的感受——这问题乍一看好像挺琐碎，但真要展开硬算的话，结果肯定特别冗长，谁都想找个“巧劲”，不用花好几个小时在没必要的复杂计算上对吧？

我一直在找一种紧凑的算法，用来计算三维空间中与三个通用球面相切的8个平面的方程。

先说说我已经有的基础：我找到了二维版本的解法（也就是求与两个圆相切的4条直线），对应的方程组和推导参数如下：

二维方程组：
$$\begin{cases} (x-x_1)^2+(y-y_1)2=r_1^2\(x-x_2)2+(y-y_2)^2=r_22\end{cases}\qquad{\color{gray}{\left[\begin{matrix}d_x=x_2-x_1\qquad d_y=y_2-y_1\d=\sqrt{d_x^2+d_y2} \end{matrix}\right.}}$$
切线方程形式：$r: y-y_0=m(x-x_0)$

对应的参数表格：
$$\begin{array}{|c||c|c|}\hline&\text{Internal}&\text{External}\\hline(x_0,y_0)&\left(\dfrac{x_2 r_1+x_1 r_2}{r_1+r_2},\dfrac{y_2 r_1+y_1 r_2}{r_1+r_2}\right)&\left(\dfrac{x_2 r_1-x_1 r_2}{r_1-r_2},\dfrac{y_2 r_1-y_1 r_2}{r_1-r_2}\right)\m&\dfrac{d_xd_y\pm(r_1+r_2)\sqrt{d^2-(r_1+r_2)2}}{d_x^2-(r_1+r_2)2}&\dfrac{d_xd_y\pm(r_1-r_2)\sqrt{d^2-(r_1-r_2)2}}{d_x^2-(r_1-r_2)2}\\hline\end{array}$$

现在我想把这个思路推广到三维的情况，对应的三个球面方程组是：
$$\begin{cases} (x-x_1)^2+(y-y_1)2+(z-z_1)^2=r_12\(x-x_2)^2+(y-y_2)2+(z-z_2)^2=r_22\(x-x_3)^2+(y-y_3)2+(z-z_3)^2=r_32\end{cases}$$

这里有个关键观察：切平面一共有8个，这让我想到结果可能和一种“二进制编码”式的符号组合有关，也就是类似 $\pm a_1\pm a_2\pm a_3$ 的形式——结合二维的公式来看，我猜测这里的$a_i$应该就是各个球面的半径。

另外，我还考虑用向量形式来简化这个问题，把球面方程写成：
$$\begin{cases}
\Vert\mathbf{x}-\mathbf{x}_1\Vert=r_1\
\Vert\mathbf{x}-\mathbf{x}_2\Vert=r_2\
\Vert\mathbf{x}-\mathbf{x}_3\Vert=r_3
\end{cases}$$
而切平面的方程则用点积的形式来表达：
$$\langle \mathbf{x}-\mathbf{x}_0,\mathbf{a}\rangle=0$$
其中 $\mathbf{x}:=(x,y,z)$，$\mathbf{x}_0:=(x_0,y_0,z_0)$，$\mathbf{a}:=(a,b,c)$ 分别是空间点向量、切平面上的定点向量和切平面的法向量。

备注：内容来源于stack exchange，提问作者Math Attack