Figure 1: full ellipsoid outlines
Figure 2: hidden outlines computed using
simplistic version of Ferrari's algorithm
Figure 3: hidden outlines computed using
the combined algorithm described in this article
Don Herbison-Evans
donherbisonevans@yahoo.com
Technical Report TR94-487
Basser Department of Computer Science
University of Sydney
Australia
(updated 22 July 2005)
ABSTRACT
A number of problems in computer graphics reduce to finding approximate real roots of quartic and cubic equations in one unknown. Various solution techniques are discussed briefly. The algorithms for analytic solution are discussed at length.
Methods are presented for controlling round-off error and overflow in the analytic solution of such equations. The solution of a quartic requires the solution of a subsidiary cubic equation. The cubics derived in the algorithms by Descartes, by Neumark, by Ferrari, by Brown, Yacoub & Fraidenraich, and by Christianson are examined.
An operation count of the resulting algorithms is presented, and statistics presented after performing a wide ranging comparison of their rounding-error behaviours.
INTRODUCTION
This article addresses the problems of solving quartic and cubic equations in computer graphics. These equations are interesting as they can be solved by analytic algorithms, and in principle need no iterative techniques.
Quartic equations need to be solved when ray tracing 4th degree surfaces e.g., a torus. Quartics also need to be solved in a number of problems involving quadric surfaces. Quadric surfaces (i.e. ellipsoids, paraboloids, hyperboloids, cones) are useful in computer graphics for generating objects with curved surfaces (Badler, 1979). Fewer primitives are required than with planar surfaces to approximate a curved surface to a given accuracy (Herbison-Evans, 1982).
Bicubic surfaces may also be used for the composition of curved objects. They have the advantage of being able to incorporate recurves: lines of inflection. There is a problem, however, when drawing the outlines of bicubics in the calculation of hidden arcs. The visibility of an outline can change where its projection intersects that of another outline. The intersection can be found as the simultaneous solution of the two projected outlines. For bicubic surfaces, these outlines are cubics, and the simultaneous solution of two of these is a sextic which can only be solved by iterative techniques. For quadric surfaces, the projected outlines are quadratic. The simultaneous solution of two of these leads to a quartic equation.
The need to solve cubic equations in computer graphics arises in the solution of the quartic equations mentioned above. Also, a number of problems which involve the use of cubic splines require the solution of cubic equations.
One simplifying feature of the computer graphics problem is that often, only the real roots (if there are any) are required. The full solution of the quartic in the complex domain (Nonweiler, 1967) is then an unnecessary use of computing resources.
Another simplification in the graphics problem is that displays have only a limited resolution, so that only a limited number of accurate digits in the solution of a cubic or quartic may be required. A resolution of 1 in 1,000,000 should in principle be achievable using single precision floating point (32 bit) arithmetic, which would be more than adequate for most current displays.
ITERATIVE TECHNIQUES
The roots of quartic and cubic equations can be obtained by iterative techniques. These techniques can be useful in animation where scenes change little from one frame to the next. Then the roots for the equations in one frame are good starting points for the solution of the equations in the next frame. There are two problems with this approach.
One is storage. For a scene composed of 'n' quadric surfaces, storage may be required for 4n(n - 1) roots between frames. A compromise is to store pointers to those pairs of quadrics which give no roots. This trivial idea can be used to halve the number of computations within a given frame, for if quadric 'Q1' has no intersection with quadric 'Q2', then 'Q2' will not intersect 'Q1'.
The other problem is more serious: it is the problem of deciding when the number of roots changes. There appears to be no simple way to find the number of roots of a cubic or quartic. The most well-known algorithm for finding the number of real roots, the Sturm sequence, involves approximately as much computation as solving the equations directly by radicals (Ralston, 1965). Without information about the number of roots, iteration where a root has disappeared can waste a lot of computer time, and searching for new roots that may have appeared becomes difficult.
Even when a root has been found, deflation of the polynomial to the next lower degree is prone to severe round-off exaggeration (Conte and de Boor, 1980).
Thus there may be an advantage in examining the techniques available for obtaining the real roots of quartics and cubics analytically.
QUARTIC EQUATIONS
Quartics are the highest degree polynomials which can be solved analytically in general by the method of radicals i.e.: operating on the coefficients with a sequence of operators from the set: sum, difference, product, quotient, and the extraction of an integral order root. An algorithm for doing this was first published by Cardano in the 16th century (Cardano, 1545). A number of other algorithms have subsequently been published. The question arises: which algorithm is for best to use on a computer to finding the real roots in terms of speed and stability for computer graphics.
Very little attention appears to have been given to a comparison of the algorithms. They have differing properties with regard to overflow and the exaggeration of round-off errors. Where a picture results from the computation, any errors may be rather obvious. Figures 1, 2, and 3 show a computer bug composed of ellipsoids with full outlines, incorrect hidden outlines, and correct hidden outlines, respectively. In computer animation, the flashing of incorrectly calculated hidden arcs is most disturbing.
Figure 2: hidden outlines computed using
simplistic version of Ferrari's algorithm
Figure 3: hidden outlines computed using
the combined algorithm described in this article
Many algorithms use the idea of first solving a particular cubic equation, the coefficients of which are derived from those of the quartic. A root of the cubic is then used to factorize the quartic into quadratics, which may then be solved. The algorithms may then be classified according to the way the coefficients of the quartic are combined to form the coefficients of the subsidiary cubic equation. For a general quartic equation of the form:
x4 + ax3 + bx2 + cx + d = 0
the subsidiary cubic can be one of the forms:
Christianson-Brown (Christianson, 1991)
y3 + (4a2b - 4b2 - 4ac + 16d - 3a4/4)y2/(a3 - 4ab + 8c) + (3a2/16 - b/2)y + (ab/16 - c/8 + a3/64) = 0
Descartes-Euler-Cardano (Strong, 1859)
y3 + (2b - 3a2/4)y2 + (3a4/16 - a2b + ac + b2 - 4d)y + (abc - a6/64 + a4b/8 - a3c/4 - a2b2/4 - c2) = 0
Ferrari-Lagrange (Turnbull, 1947)
y3 + by2 + (ac - 4d)y + (a2d + c2 - 4bd) = 0
Neumark (Neumark, 1965)
y3 - 2by2 + (ac + b2 - 4d)y + (a2d - abc + c2) = 0
Yacoub-Fraidenraich-Brown (Yacoub & Fraidenraich, 2003)
(a3 - 4ab + 8c)y3 + (a2b - 4b2 + 2ac + 16d)y2 + (a2c - 4bc + 8ad)y + (a2d - c2) = 0
The casual user of the literature may be confused by variations in the presentation of quartic and cubic equations. Sometimes, the highest degree term has a non-unit coefficient. Sometimes the coefficients are labelled from the lowest degree term to the highest. Sometimes numerical factors of 3, 4 and 6 are included. There are also a number of trivial changes to the cubic caused by the following:
if
y3 + py2 + qy + r = 0
then
z3 - pz2 + qz - r = 0 for z = - y
and
z3 + 2pz2 + 4qz + 8r = 0 for z = 2y
Brown (Brown, 1967) originally published a method of solving palindromic quartics. Christianson (Christianson, 1991) showed how to apply this to a general quartic equation, using the subsidiary cubic equation coefficients:
p = (4a2b - 4b2 - 4ac + 16d - 3a4/4)/(a3 - 4ab + 8c)
q = 3a2/16 - b/2
r = ab/16 - c/8 + a3/64
___________________________________________
| | a+ | a- |
| |_______________|_______________|
| | b+ | b- | b+ | b- |
|_________|_______|_______|_______|_______|
| | d+ | | q | r | q |
| c+ |____|_______|_______|_______|_______|
| | d- | | pq | r | q |
|____|____|_______|_______|_______|_______|
| | d+ | r | q | | q |
| c- |____|_______|_______|_______|_______|
| | d- | r | q | p | pq |
|____|____|_______|_______|_______|_______|
Table 1
The coefficients of the subsidiary cubic for
Christianson and Brown's algorithm that can be stably
computed from the coefficients of the quartic.
Unfortunately, no combinations of signs of the quartic coefficients give a stable computation of the cubic coefficients.
Having solved the cubic, any solution 'y' may be used to calculate 'k' using either
k4 = y4 + y2e2 + ye1 + e0
or
k2 = y2 + e2/2 + e1/4y
where
e0 = d - ac/4 + a2b/16 - 3a4/256
e1 = c - ab/2 + a3/8
e2 = b - 3a2/8
The use of the equation for k4 is only of benefit if either 'y=0' or else 'y' and 'e1' are negative and 'e0' and 'e2' are positive.
The value of 'k' is used to calculate the quantities :
g = 4y
k3
h = (6y2 + e2)k2
These are used to form the quadratic in 'Z' :
Z2 + gZ/k4 + (h/k4 - 2) = 0
The roots of this, 'Z1' and 'Z2', are then use form the quadratics in 'z' :
z2 - Z1*z + 1 = 0
and
z2 - Z2*z + 1 = 0
the roots of which are each used to form the corresponding root of the original quartic :
x = y - kz - a/4
This algorithm uses a subsidiary cubic with coefficients :
p = 2b - 3a2/4
q = 3a4/16 - a2b + ac + b2 - 4d
r = abc - a6/64 + a4b/8 - a3c/4
- a2b2/4 - c2
There is only one combination of quartic coefficients for which the evaluation of the subsidiary cubic coefficients is stable:
___________________________________________ | | a+ | a- | | |_______________|_______________| | | b+ | b- | b+ | b- | |_________|_______|_______|_______|_______| | | d+ | | p r | | p | | c+ |____|_______|_______|_______|_______| | | d- | | p q r | | p | |____|____|_______|_______|_______|_______| | | d+ | | p | | p | | c- |____|_______|_______|_______|_______| | | d- | | p | | p q | |____|____|_______|_______|_______|_______| Table 2 For Descartes' algorithm, the coefficients of the subsidiary cubic that can be stably computed from the coefficients of the quartic.
However, if 'a' and 'c' are negated, the combination (a-, b-, c-, d-) becomes stable at the trivial cost of having to negate the resulting quartic solutions.
In this algorithm, the calculation of the cubic coefficients involves significantly more operations than Ferrari's or Neumark's algorithms. Also, the high power of 'a' in the coefficients makes this algorithm prone to loss of precision and also overflow.
In this algorithm, if the greatest root of the cubic, 'y', is negative, the quartic has no real roots. Otherwise, the solutions 'x' of the quartic are obtained from the quadratics:
x2 + kx + g = 0
and
x2 - kx + h = 0
where
k = y(1/2)
g = (y + e2 + e1/k)
h = (y + e2 - e1/k)
and
e1 = c + a3/8 - ab/2
e2 = b - 3a2d/8
There appears to be no way of making the evaluation of 'e1', 'e2', 'g' and 'h' stable. Some quantities are bound to be subtracted, leading to possible loss of precision.
Ferrari's algorithm has the following coefficients for the subsidiary cubic :
p = b
q = ac - 4d
r = a2d + c2 - 4bd
These have two combinations of signs of 'a', 'b', 'c' and 'd' for which the derivation of all of the coefficients of the cubic, 'p', 'q', and 'r' are stable:
___________________________________________
| | a+ | a- |
| |_______________|_______________|
| | b+ | b- | b+ | b- |
|_________|_______|_______|_______|_______|
| | d+ | p | p r | p q | p q r |
| c+ |____|_______|_______|_______|_______|
| | d- | p q | p q | p | p |
|____|____|_______|_______|_______|_______|
| | d+ | p q | p q r | p | p r |
| c- |____|_______|_______|_______|_______|
| | d- | p | p | p q | p q |
|____|____|_______|_______|_______|_______|
Table 3
For Ferrari's algorithm, the coefficients
of the subsidiary cubic that can be stably
computed from the coefficients of the quartic.
The coefficients of the subsequent quadratics depend on two intermediate quantities, 'e' and 'f', where
e2 = a2/4 - b - y
f2 = y2/4 - d
ef = (ay/4 + c/2)
Using Ferrari's method, the quadratic equations giving the solutions to the quartic are :
x2 + Gx + H = 0
and
x2 + gx + h = 0
where
G = a/2 + e
g = a/2 - e
H = -y/2 + f
h = -y/2 - f
The signs of each of the quartic coefficients 'a', 'b', 'c', 'd' and the root of the cubic 'y', may be positive or negative, giving 32 possible combinations of signs. Of these, only 12 can be clearly solved in a stable fashion for 'e' and 'f' by the choice of 2 out of the 3 equations involving them. In the remaining 20 cases, the most stable choices are unclear. This is shown in tables 2 and 3. Of the 12 stable cases, 2 are from the stable cases for the calculation of 'p', 'q' and 'r'. In the other 10 cases, the value of y may be unreliable.
_________________________________________________________________
| | a+ | a- |
| y+ |__________________________|__________________________|
| | b+ | b- | b+ | b- |
|_________|____________|_____________|____________|_____________|
| | | | | | |
| | d+ | ef | ef | | |
| c+ |____|____________|_____________|____________|_____________|
| | | 2 | 2 | 2 | 2 |
| | d- | ef f | ef f | f | f |
|____|____|____________|_____________|____________|_____________|
| | | | | | |
| | d+ | | | ef | ef |
| c- |____|____________|_____________|____________|_____________|
| | | 2 | 2 | 2 | 2 |
| | d- | f | f | ef f | ef f |
|____|____|____________|_____________|____________|_____________|
Table 4
For Ferrari's algorithm, the intermediate
quantities that can be stably computed from
the coefficients of the quartic and a
positive root of the cubic.
________________________________________________________________
| | a+ | a- |
| y- |__________________________|__________________________|
| | b+ | b- | b+ | b- |
|_________|____________|_____________|____________|_____________|
| | | | 2 | | 2 |
| | d+ | | e | ef | ef e |
| c+ |____|____________|_____________|____________|_____________|
| | | 2 | 2 2 | 2 | 2 2 |
| | d- | f | e f | ef f | ef e f |
|____|____|____________|_____________|____________|_____________|
| | | | 2 | | 2 |
| | d+ | ef | ef e | | e |
| c- |____|____________|_____________|____________|_____________|
| | | 2 | 2 2 | 2 | 2 2 |
| | d- | ef f | ef e f | f | e f |
|____|____|____________|_____________|____________|_____________|
Table 5
For Ferrari's algorithm, the intermediate
quantities that can be stably computed from
the coefficients of the quartic and a
negative root of the cubic.
If 'a' and 'e' are the same sign, and 'b' and 'y' are the same sign, then 'g' may be more accurately computed using:
g = (b + y)/G
If 'a' and 'e' are opposite signs then 'G' can be more accurately computed from 'g' in a similar fashion.
If 'y' and 'f' are the same sign, then 'H' may be more accurately computed using:
H = d/h
If 'y' and 'f' are opposite in sign, then 'h' can be computed similarly from H more accurately.
The solution of the quadratic equations requires the evaluation of the discriminants:
g2 - 4h and G2 - 4H
Unless 'h' and 'H' are negative, one or both of these evaluations will be unstable. Unfortunately, positive 'h' and 'H' values are incompatible with the 2 stable cases for the evaluation of 'p', 'q', 'r', 'e' and 'f', so there is no combination of coefficients for which Ferrari's algorithm can be made entirely stable.
NEUMARK
The algorithm of Neumark unfortunately has no combinations of signs of the quartic coefficients for which there is a stable computation of all of the cubic coefficients :
p = 2b
q = ac + b2 - 4d
r = a2d - abc + c2
___________________________________________ | | a+ | a- | | |_______________|_______________| | | b+ | b- | b+ | b- | |_________|_______|_______|_______|_______| | | d+ | p | p r | p r | p | | c+ |____|_______|_______|_______|_______| | | d- | p q | p q | p | p | |____|____|_______|_______|_______|_______| | | d+ | p r | p | p | p r | | c- |____|_______|_______|_______|_______| | | d- | p | p | p q | p q | |____|____|_______|_______|_______|_______| Table 6 For Neumark's algorithm, the coefficients of the subsidiary cubic that can be stably computed from the coefficients of the quartic.
Again, the solutions 'x' of the quartic are obtained from the quadratics:
x2 + Gx + H = 0
and
x2 + gx + h = 0
where:
G = (a + (a2 - 4y)(1/2))/2
g = (a - (a2 - 4y)(1/2))/2
and
H = (b-y)/2 + (a(b-y) - 2c)/(2(a2 - 4y)(1/2))
h = (b-y)/2 - (a(b-y) - 2c)/(2(a2 - 4y)(1/2))
In Nemark's algorithm, the coefficients of the quadratic equations can be computed in a stable fashion from the solution of the cubic. Any cancellation due to the additions and subtractions can be eliminated by writing :
G = g1 + g2
g = g1 - g2
H = h1 + h2
h = h1 - h2
where
g1 = a/2
g2 = ((a2 - 4y)(1/2))/2
h1= (b - y)/2
h2 = (a(b-y)/2 - c)/(a2 - 4y)(1/2)
and using the identities
g.G = y
h.H = d
Thus if 'g1' and 'g2' are the same sign, 'G' will be accurate but 'g' will lose significant digits by cancellation. Then the value of 'g' can be better obtained using:
g = y/G
If 'g1' and 'g2' are of opposite signs, then 'g' will be accurate, and 'G' better obtained using:
G = y/g
Similarly, 'h' and 'H' can be obtained without cancellation from 'h1', 'h2' and 'd'.
The computation of 'g2' and 'h2' can be made more stable under some circumstances using the alternative formulation:
h2 = (((b - y)2 - 4d)(1/2))/2
Furthermore
g2 = (a - c)/((b - y)2 - 4d)(1/2)
Thus 'g2' and 'h2' can both be computed either using
m = (b - y)2 - 4d
or using
n = a2 - 4y
If 'y' is negative, 'n' should be used. If 'y' is positive and 'b' and 'd' are negative, 'm' should be used. Thus 7 of the 32 sign combinations give stable results with this algorithm. These are shown in tables 7 and 8.
____________________________________________________ | | a+ | a- | | y+ |_____________________|__________________| | | b+ | b- | b+ | b- | |_________|_______|_____________|_______|__________| | | d+ | g1 | g1 h1 | g1 | g1 h1 | | c+ |____|_______|_____________|_______|__________| | | d- | g1 | g1 g2 h1 h2 | g1 | g1 h1 h2 | |____|____|_______|_____________|_______|__________| | | d+ | g1 | g1 h1 | g1 | g1 h1 | | c- |____|_______|_____________|_______|__________| | | d- | g1 | g1 h1 h2 | g1 | g1 h1 h2 | |____|____|_______|_____________|_______|__________| Table 7 For Neumark's algorithm, the intermediate quantities that can be stably computed from the coefficients of the quartic and a positive root of the cubic.
_______________________________________________________ | | a+ | a- | | y- |_____________________|_____________________| | | b+ | b- | b+ | b- | |_________|_____________|_______|_____________|_______| | | d+ | g1 g2 h1 | g1 g2 | g1 g2 h1 h2 | g1 g2 | | c+ |____|_____________|_______|_____________|_______| | | d- | g1 g2 h1 h2 | g1 g2 | g1 g2 h1 h2 | g1 g2 | |____|____|_____________|_______|_____________|_______| | | d+ | g1 g2 h1 h2 | g1 g2 | g1 g2 h1 | g1 g2 | | c- |____|_____________|_______|_____________|_______| | | d- | g1 g2 h1 h2 | g1 g2 | g1 g2 h1 h2 | g1 g2 | |____|____|_____________|_______|_____________|_______| Table 8 For Neumark's algorithm, the intermediate quantities that can be stably computed from the coefficients of the quartic and a negative root of the cubic.
For other cases, a rough guide to which expression to use can be found by assessing the errors of each of these expressions by summing the moduli of the addends:
e(m) = b2 + 2|by| + y2 + 4|d|
e(n) = a2 + 4|y|
Thus, if
|m|.e(n) > |n|.e(m)
then 'm' should be used, otherwise 'n' is more accurate.
YACOUB - FRAIDENRAICH - BROWN
Brown (Brown, 1967) originally showed how to solve palindromic quartics, and Yacoub and Fraidenraich (Yacoub & Fraidenraich, 2003) showed a different way of applying this to a general quartic equation using, except for special cases, the subsidiary cubic equation coefficients:
p = P/U, q = Q/U and r = R/U with
P = a2b - 4b2 + 2ac + 16d
Q = a2c - 4bc + 8ad
R = a2d - c2
U = a3 - 4ab + 8c
___________________________________________
| | a+ | a- |
| |_______________|_______________|
| | b+ | b- | b+ | b- |
|_________|_______|_______|_______|_______|
| | d+ | | QU | | |
| c+ |____|_______|_______|_______|_______|
| | d- | R | RU | R | PQR |
|____|____|_______|_______|_______|_______|
| | d+ | | | U | Q |
| c- |____|_______|_______|_______|_______|
| | d- | R | PQR | RU | R |
|____|____|_______|_______|_______|_______|
Table 9
The coefficients of the subsidiary cubic for
Yacoub, Fraidenraich and Brown's algorithm that can be stably
computed from the coefficients of the quartic.
Unfortunately, no combinations of signs of the quartic coefficients give a stable computation of the cubic coefficients.
Having solved the cubic, any solution 'y' may be used to calculate 'k' using
k = a + 4y
and this is used to calculate the quantities :
e = (a3 - 4c - 2ab + 6a2y - 16by)/k
f2 = (a3 + 8c - 4ab)/k
and then :
g2 = 2(e + f*k)
h2 = 2(e - f*k)
where the positive square root may be taken for each of 'f','g', and 'h'. The quartic roots are calculated from these using :
x1 = (-a - f - g)/4
x2 = (-a - f + g)/4
x3 = (-a + f - h)/4
x4 = (-a + f + h)/4
If there are three real roots of the cubic, then choosing one with the same sign as 'a' will make the calculation of 'k' more accurate. However there is no guarantee that this alone produces the best results, and all three cubic roots may need to be tried to find the most accurate answers.
THE CUBIC
Following Tartaglia, let the cubic equation be
y3 + py2 + qy + r = 0
The solution has been published elsewhere in many texts, but here is expressed (Littlewood, 1950) using:
u = q - p2/3
v = r - pq/3 + 2p3/27
and the discriminant:
j = 4(u/3)3 + v2
If 'j' is positive then there is one root 'y' to the cubic, which may be found using:
y = ((w - v)/2)(1/3) - (u/3)(2/(w - v))(1/3) - p/3
where
w = j(1/2)
As 'w' is chosen here to be positive, this formulation is suitable if 'v' is negative. However, the calculation in this form can lose accuracy if 'v' is positive. This problem can be overcome by the rationalization:
(w - v)/2 = (w2 - v2)/(2(w + v)) = 2(u/3)3/(w + v)
giving the alternative formulation of the root:
y = (u/3)(2/(w + v))(1/3) - ((w + v)/2)(1/3) - p/3
A computational problem with this algorithm is overflow while calculating 'w', for
O(j) = O(p6) + O(q3) + O(r2)
OVERFLOW
If the cubic is the subsidiary cubic of a quartic, the different algorithms have differing overflow behaviours. The most obviously affected quantity is 'j'. In terms of the individual coefficients :
Christianson: O(j) = O(a6) + O(b6) + O(c2) + O(c-6) + O(d6)
Descartes: O(j) = O(a12) + O(b6) + O(c4) + O(d3)
Ferrari: O(j) = O(a3) + O(b6) + O(c4) + O(d3)
Neumark: O(j) = O(a6) + O(b6) + O(c4) + O(d3)
Yacoub & Fraidenraich: O(j) = O(a6) + O(a-6) + O(b6) + O(c4) + O(d6)
Before evaluating the terms of 'w', it is wise to test 'p', 'q' and 'r' against the appropriate root of the maximum number represented on the machine ('M'). The values of 'u' and 'v' should similarly be tested. In the event that some value is too large, various approximations may be employed: e.g.
if |p| > 27M(1/3), then y ~= -p
if |v| > M(1/2), then y ~= v(1/3)
if |u| > 3M(1/3)/4, then y ~= 4(1/3)u/3
If the discriminant 'j' is negative, then there are 3 real roots to the cubic. These real roots of the cubic may then be obtained via parameters 's', 't' and 'k':
s = (-u/3)(1/2)
t = -v/(2s3)
k = arccos(t)/3
giving
y1 = 2s.cos(k) - p/3
y2 = s(-cos(k) + 3(1/2)sin(k)) - p/3
y3 = s(-cos(k) - 3(1/2)sin(k)) - p/3
Note that if the discriminant 'j' is negative, then 'u' must also be negative, guaranteeing a real value for 's'. This value may be taken as positive without loss of generality. Also, 'k' will lie in the range 0 to 60 degrees, so that cos(k) and sin(k) are both positive. Thus we have:
y1 >= y2 >= y3.
If the cubic is a subsidiary of a quartic, then either 'y1' or 'y3' may be the most useful root. For example, in Neumark's algorithm b = 2p , so although 'y1' may be the largest root, it may not be positive. Then if 'b' and 'd' are both negative, it would be advantageous to use the most negative root 'y3'.
The functions sine and cosine of arccos(t)/3 may be tabulated to speed the calculation (Herbison-Evans, 1983). Sufficient accuracy (1 in 107) can be obtained with a table of 200 entries with linear interpolation, requiring 4 multiplications, 8 additions and 2 tests for each function. When t is near its extremes, the asymptotic forms may be useful:
for t -> 0 :
sin(arccos(t)/3) ~= (1 - t/3(1/2))/2 + O(t2)
cos(arccos(t)/3) ~= 3(1/2))/2 + t/6 + O(t2)
for t -> 1 :
sin(arccos(t)/3) ~= (2(1 - t))(1/2)/3 + O((1-t)2)
cos(arccos(t)/3) ~= 1 - (1-t)/9 + O((1-t)2)
If the discriminant 'j' is expanded in terms of the coefficients of the cubic, it has 10 terms. Two pairs of terms cancel and another pair coalesce, leaving 5 independent terms. In principle, any pair of subsets of these may cancel catastrophically, leaving an incorrect value or even an incorrect sign for the discriminant. This problem can be alleviated by calculating the 5 terms separately, and then combining them in increasing order of magnitude (Wilkinson, 1963). When solving quartics, the discriminant can be expanded in terms of the quartic coefficients directly. This gives fifteen terms, which ideally can be sorted by modulus, and combined in increasing order.
NO REAL ROOTS
The absence of real roots to a quartic becomes apparent in some of the algorithms examined here after the subsidiary cubic has been solved, when the arguements of square roots of intermediate quantities are found to be negative. Alternatively, the extrema of the quartic can be found by solving the cubic:
4x3 + 3ax2 + 2bx + c = 0
and substituting each root back into the quartic to find its values at its extrema. If all these values are positive, then the quartic has no real roots.
More economical methods for discovering that there are no real roots may be available.
CONCLUSIONS
There have been many algorithms proposed for solving quartic and cubic equations, but most have been proposed with aims of elegance, generality or simplicity rather than error minimisation or overflow avoidance.
The operation counts of the best combination of stabilized algorithms for non-special cases are summarized in table 10:
| cubic | worst 13 |
worst 15 |
worst 3 |
worst 19 |
| rest of quartic | worst 16 |
worst 34 |
worst 2 |
worst 39 |
| 2 x quadratic | ||||
| totals | worst 35 |
worst 59 |
worst 9 |
worst 64 |
Table 10 Operation counts for a best combination of stabilized algorithms.
A computer program was written to perform a comparison of the stabilities of the five algorithms for the solution of quartic equations. Quartics were examined which had all combinations and permutations of coefficients from the set:
108, 104, 1, 10(-4), 10(-8), -108, -104, -1, -10-4, -10-8
Of the 10,000 equations, all five algorithms agreed on the number of real roots in only 5,947 cases. For these quartics, all five agreed that 726 had no real roots. Of the remaining equations, Neumark's algorithm had the least worst error in 2,080 quartics, Ferrari's in 863, Yacoub & Fraidenraich's 544, Descartes' 42, and Christianson's 13. For the other quartics, there were problems in trying to compare the algorithms.
These occurred when the subsidiary cubic produced 3 roots, and the use of these different roots produced a different number of roots for the quartic. This problem was apparent for all the algorithms in a significant number of equations. It is unclear how to compare their accuracies in these situations. The statistics of these findings are presented in Table 11.
| algorithm | ||
| Christianson | ||
| Descartes | ||
| Ferrari | ||
| Neumark | ||
| Yacoub |
Table 11 Comparison of cases where different cubic roots produced different numbers of quartic roots in a benchmark test of 10,000 quartics.
A check on the accuracy of the roots can be done at the cost of more computation. Each root may be substituted back into the original equation and the residual calculated. This can then be substituted into the derivative to give an estimate of the error of the root or used for a Reguli-Falsi or better still a Newton-Raphson correction.
A further comment may be useful here concerning the language used to implement these algorithm. Compilers for the language C often perform double precision operations on single precision variables ('float'), converting back to single precision for storage. Thus there might be little speed advantage in using 'float' variables compared with using 'double' for these algorithms. Fortran compilers may not do this. For example, using a VAX8600, the time taken to solve the 10,000 different quartics was 6 seconds, for Fortran single precision (using f77), 15 seconds for C single precision (using cc), and 16 seconds for C using double precision.
It may have been observed that the 2 stable cases for the computation of the cubic coefficients for Ferrari's algorithm are different from those of Descartes' algorithm, so that the cubic coefficients may be computed in a stable fashion for 4 of the possible 16 sign combinations of the quartic coefficients. Further work may be able to improve this situation. It would be good to find more algorithms or variations on the five listed here which allowed other combinations of coefficient signs to be handled in a stable fashion. It is far from clear why there are five and only five algorithms so far discovered. Are there more? Are there an infinite number or algorithms?
ACKNOWLEDGEMENTS
Thanks are due to the Basser Department of Computer Science at the University of Sydney (Australia), the Department of Computer Science at the University of Waterloo (Canada), and the Department of Software Engineering at the University of Technology, Sydney (Australia) where much of this work was done. Thanks are also due to Charles Prineas on whose initial investigations this work was based. Thanks are also due to the late Alan Tritter for discussions, Zoe Kaszas for the initial preparation of this paper, Alan Paeth for formatting this paper for publication, and to Professor John Bennett for his continual encouragement.
REFERENCES
Badler, N.I. & Smoliar, S.W., Digital Representations of Human Movement, ACM Computing Surveys, Vol. 11, No. 1, pp. 24-27 (1979)
Brown, K.S., Reducing Quartics to Cubics, http://mathpages.com/home/kmath296.htm (1967)
Cardano, G., Ars Magna, Nurmberg (1545)
Christianson, B., Solving Quartics Using Palindromes, Mathematical Gazette, Vol. 75, pp. 327-328 (1991)
Conte, S.D. & de Boor, C., Elementary Numerical Analysis, McGraw-Hill, p. 117 (1980)
Herbison-Evans, D., Real Time Animation of Human Figure Drawings with Hidden Lines Omitted, I.E.E. Computer Graphics and Applications, Vol. 2, No. 9, pp. 27-33 (1982)
Herbison-Evans, D., Caterpillars and the Inaccurate Solution of Cubic and Quartic Equations, Australian Computer Science Communications, Vol. 5, No. 1, pp. 80-9-(1983)
Littlewood, D.E., A University Algebra, Heineman, London, p. 173 (1950)
Nonweiler, T.R.F., Roots of Low Order Polynomial Equations, Collected Algorithms of the ACM, (C2), Algorithm 326 (1967)
Neumark, S., Solution of Cubic and Quartic Equations, Pergamon Press, Oxford (1965)
Ralston, A., A First Course in Numerical Analysis, McGraw-Hill, p. 351 (1965)
Strong, T., Elementary and Higher Algebra, Pratt and Oakley, p. 469 (1859)
Turnbull, H.W., Theory of Equations, Oliver and Boyd, London, Fourth Edition, p. 130 (1947)
Wilkinson, J.J., Rounding Errors in Algebraic Processes, Prentice-Hall, p. 17 (1963)
Yacoub, M. D. & Fraidenraich G., A New Simple Solution of the General Quartic Equation, submitted for publication (2003)