0 はじめに
トーラスの描画の関係で実数係数四次方程式の実数解を求める必要が生じたので、まとめたいと思います。Inigo Quilezさんのウェブサイトに載っていたトーラスを直線との交点で描画する方法:
intersectors
に載っていた変数変換を参考にしつつ、基本的にはWikipediaのサイトから、
四次方程式 三次方程式
を見て自分なりに再構成しました。
1 正規形
方程式は次の形をしているとします。
x^4 + 4k_3x^3 + 4k_2x^2 + 8k_1x + 4k_0 = 0.
これはトーラスと直線の交点を求めようとして4次方程式を作るとこんな感じで整理されるゆえに係数がこんなことになっているのですが、実際この方がうまくいくからというのもあります。
2 最初の変数変換と複二次式
まず$x = y - k_3$と変換すると、次のように$y^3$の係数が消えます。
y^4 + 6c_2y^2 + 4c_1y + 3c_0 = 0.
ただし、
\begin{align}
c_2 &= \frac{1}{3}(2k_2 - 3{k_3}^2), \\[5pt]
c_1 &= 2({k_3}^3 - k_2k_3 + k_1), \\[5pt]
c_0 &= \frac{1}{3}(-3{k_3}^4 + 4{k_3}^2k_2 - 8k_1k_3 + 4k_0).
\end{align}
これは解の平均を取ることでその和が0になるように平行移動することに相当します。ここでもし $c_1=0$ であるなら、これは
y^4 + 6c_2y^2 + 3c_0 = 0
となり、$y^2$ についての二次方程式(複二次式)になるので普通に解けます。すなわち、
d = 9{c_2}^2 - 3c_0
とおくと、
- $d < 0$ のとき・・実数解なし。
- $d \geq 0$ のとき・・
y^2 = -3c_2 \pm \sqrt{d}
となるので、4つの解:
\sqrt{-3c_2+\sqrt{d}}, ~~-\sqrt{-3c_2+\sqrt{d}}, ~~\sqrt{-3c_2-\sqrt{d}}, ~~-\sqrt{-3c_2-\sqrt{d}}
のうち根号内が0以上のものが解となります(根号内が0になる場合は適宜退化する)。こうして0個~4個の実数解が得られます。
以降の議論では $c_1 \ne 0$ であるとします。
3 三次分解方程式
次に、フェラリの方法を使います。これはWikiなどでも紹介されているものです。次の方程式:
4u^3 + 12c_2u^2 + (9{c_2}^2 - 3c_0)u - {c_1}^2 = 0
が正の数の解 $u$ を持つならば、これにより $y$ についての方程式を、
\left( y^2 + 3c_2 + 2u \right)^2 - u\left( 2y - \frac{c_1}{u} \right)^2 = 0
の形にできるのです(単純に展開して確かめられます)。正の数の解 $u$ は常に存在します。なぜなら ${c_1}^2 > 0$ なので、左辺は $u=0$ で負となり、また十分大きい $u$ について正なので、中間値の定理により少なくともひとつの正の解を持ちます。それは最も大きい実数解を取れば確実なので、この後のステップでそれを求めることにします。
4 三次分解方程式の正の解
さて、まず $u = v - c_2$ と変換します。すると $v^2$ の係数は消えて、
4v^3 - 3Qv = R
とできます。この方程式の最大の実数解 $v$ を求めればいいです。ここに、
Q = c_0 + {c_2}^2,~~~~R = {c_1}^2 + {c_2}^3 - 3c_0c_2
です。まず $Q$ や $R$ が0の場合は最大の実数解が楽に得られるので割愛します。以降では $Q, R$ は0でないとします。$R$ の符号(正なら1, 負なら-1)を $sgn(R)$ と書くことにします。
次のようにおきます。
\Delta = R^2 - Q^3, ~~~~h = \frac{R}{|Q|\sqrt{|Q|}}
- $\Delta < 0$ のとき。
この場合、$|h|<1$ となり、また $Q>0$ となります。なので、
v = \sqrt{Q}w, ~~~~4w^3 - 3w = h
と変換すると、3倍角の公式:
4{\cos}^3\theta - 3\cos\theta = \cos{3\theta}
を用いて、 $h=\cos{3\theta}$ を満たす $\theta~(0 < \theta < \pi /3)$ に対して $\cos\theta$ が要求を満たします。実際、この場合解というのはこのような $\theta$ に対して
w=\cos\theta,~~\cos(\theta + \frac{2\pi}{3}),~~\cos(\theta + \frac{4\pi}{3})
の3つですが、これらはそれぞれ $1 >w> 1/2,~-1/2 >w> -1,~1/2 >w> -1/2$ と異なる区間に属し、そのうち最大のものが $\cos\theta$ というわけです。
2. $\Delta \geq 0$ のとき
この場合は $|h| \geq 1$ となります。$Q$ の符号により分けることになります。
2.1 $Q > 0$ のとき
$R$の符号を考慮して $v = sgn(R)\sqrt{Q}w$ と変換します。符号は3乗しても変わらないので、 $sgn(R)R = |R|$ により、
4w^3 - 3w = |h|
となります。そこで双曲線関数:
\cosh x = \frac{e^x + e^{-x}}{2}, ~~~~ \sinh x = \frac{e^x - e^{-x}}{2}
の出番です。これらは、
4{\cosh}^3x - 3\cosh x = \cosh 3x,~~~~4{\sinh}^3 x + 3\sinh x = \sinh 3x
という3倍角の公式を満たします(ただの計算です)。そこで、一つ目を用いて、
\cosh 3x = |h|, ~~~~ x = \frac{1}{3}\log(|h| \pm \sqrt{h^2 - 1})
を得ます。二つの $x$ は正負が逆なので $\cosh$ をかませると同じ値になります。この $\cosh x$が唯一の実数解 $w$ になります。
2.2 $Q < 0$のとき
このときは $v = \sqrt{-Q}w$ と変換することにします。これにより、
4w^3 + 3w = h
となります。そこで、
\sinh 3x = h,~~~~ x = \frac{1}{3}\log( h + \sqrt{h^2 + 1} )
とおけば求める唯一の解 $w$ は $\sinh x$ として得られます。
場合によっては $w$ が負になることもあり得ますが、議論により最終的に得られる $u$ は必ず正の数になります。
5 仕上げ
正の解 $u$ が得られたので、変形:
\left( y^2 + 3c_2 + 2u \right)^2 - u\left( 2y - \frac{c_1}{u} \right)^2 = 0
を用いて実数解を求めることにします。これは二つに分解して、
y^2 - 2\sqrt{u}~y + (3c_2 + 2u + \frac{c_1}{\sqrt{u}}) = 0, \\[15pt]
y^2 + 2\sqrt{u}~y + (3c_2 + 2u - \frac{c_1}{\sqrt{u}}) = 0~
とすればそれぞれ解くことで合計0個~4個の実数解が然るべく得られます。めでたし。
6 さいごに
トーラスの描画はレイマーチングの方が楽だと思います。