Two point raytracing for reflection off a 3D plane

Lanczos solutions for cubics and quartics

From Lanczos (1956), pages 6-8, 19-22:

3. Cubic equations. Equations of third and fourth order are still solvable by algebraic formulas. However, the numerical computations required by the formulas are usually so involved and time-absorbing that we prefer less cumbersome methods which give the roots in approximation only but still close enough for later refinement.

The solution of a cubic equation (with real coefficients) is particularly convenient since one of the roots must be real. After finding this root, the other two roots follow immediately by solving a quadratic equation.

A general cubic equation can be written in the form

$$f(\xi) = \xi^3 + a \xi^2 + b \xi - c = 0$$

The factor of $\xi^3$ can always be normalized to 1 since we can divide through by the highest coefficient. Moreover, the absolute term can always be made negative because, if it is originally positive, we put $\xi_1 = - \xi$ and operate with $\xi_1$ .

Now it is convenient to introduce a new scale factor which will normalize the absolute term to $-1$ . We put

$$x = \alpha \xi , a_1 = \alpha a , b_1 = \alpha^2 b , c_1 = \alpha^3 c$$

and write the new equation

$$f(x) = x^3 + a_1 x^2 + b_1 x - c_1 = 0$$

If we choose

$$\alpha = 1 / {\root 3 \of c}$$

we obtain

$$c_1 = 1 .$$

Now, since $f(0)$ is negative and $f(\infty )$ is positive, we know that there must be at least one root between $x=0$ and $x=\infty$ . We put $x=1$ and evaluate $f(1)$ . If $f(1)$ is positive, the root must be between 0 and $1$ ; if $f(1)$ is negative, the root must be between $1$ and $\infty$ . Moreover, since

$$x_1 \cdot x_2 \cdot x_3 = 1$$

we know in advance that we cannot have three roots between 0 and $1$ , or $1$ and $\infty$ . Hence if $f(1)>0$ , we know that there must be one and only one real root in the interval $[ 0 , 1 ]$ , while if $f(1) < 0$ , we know that there must be one and only one real root in the interval $[1 , \infty]$ . The latter interval can be changed to the interval $[1,0]$ by the transformation

$$\bar x = \frac{1}{x}$$

which simply means that the coefficients of the equation change their sequence:

$$- c_1 {\bar x}^3 + b_1 {\bar x}^2 + a_1 {\bar x} + 1 = 0$$

Hence we have reduced our problem to the new problem: find the real root of a cubic equation in the range $[ 0 , 1 ]$ . We solve this problem in good approximation by taking advantage of the remarkable properties of the Chebyshev polynomials (cf. VII, 9) which enable us to reduce a higher power to lower powers with a small error. In particular, the third Chebyshev polynomial

$$T_3^* (x) = 32x^3 - 48x^2 +18 x -1$$

normalized to the range $[ 0 , 1 ]$ gives

$$x^3 = \frac{48 x^2 - 18x +1}{32} = 1.5x^2 - 0.5625 x +0.03125$$

with a maximum error of $\pm \frac{1}{32}$ . The original cubic is thus reducible to a quadratic with an error not exceeding $3 \%$ .

We now solve this quadratic, retaining only the root between 0 and $1$ .

$$\vdots \makebox[0.9\textwidth][1pt]{\ }$$

11. Equations of fourth order. Algebraic equations of fourth order with generally complex roots occur frequently in the stability analysis of airplanes and in problems involving servomechanisms. The historical method of solving algebraic equations of fourth order (also called biquadratic or quartic equations) involves the following steps. By a transformation of the form $x + \alpha$ the coefficient of the cubic term is annihilated. Then an auxiliary cubic equation is solved. The roots of the original equation are constructed with the help of the three roots of the auxiliary cubic. Numerically this method is lengthy and cumbersome. The following modification of the traditional procedure yields the four roots of an arbitrary quartic equation with real coefficients on the basis of a quick and numerically convenient scheme.

Every equation of the form

$$x^4 + c_1 x^3 + c_2 x^2 + c_3 x + c_4 = 0$$

can be rewritten as follows:

$$( x^2 + \alpha x + \beta )^2 = (a x + b) ^2$$

If the original $c_i$ are real, the new coefficients are also real. Hence the original equation becomes solvable in the form of the quadratic equation

$$x^2 + \alpha x + \beta \pm ( a x + b ) = 0$$

which has four (generally complex) roots, obtainable by the standard formula. The new coefficients can be determined as follows. We evaluate in succession the following numerical constants:

$$\alpha = {\frac{ c_1 }{2} }, A = c_2 - \alpha^2 , B = c_3 - \alpha A$$

and form the cubic equation

$$\xi^3 + ( 2 A - \alpha^2 ) \xi^2 + (A^2 + 2 B \alpha - 4 c_4 ) \xi - B^2 = 0$$

Since the left side is negative at $\xi = 0$ , a positive real root must exist. We determine this root according to the method of § 3. In order to avoid later corrections, it is advisable to add at this point Newton's correction (cf. § 5), obtaining $\xi$ with great accuracy. The coefficients of the reduced equation are then determined as follows:

$$\begin{array}{ll} {\displaystyle \alpha =} {\textstyle \frac{... ...ft( \alpha - {\frac{ B}{\xi}} \right) } \end{array}\mbox{\ \ .}$$

Two point raytracing for reflection off a 3D plane