Solving Cubics and Quartics.pdf

(279 KB) Pobierz
Solving Cubics and Quartics
Solving Cubics and Quartics
http://www-staff.it.uts.edu.au/~don/pubs/solving.html
Solving Quartics and Cubics for Graphics
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
1 z 16
2008-08-18 22:29
 
Solving Cubics and Quartics
http://www-staff.it.uts.edu.au/~don/pubs/solving.html
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 1: full ellipsoid outlines
2 z 16
2008-08-18 22:29
 
18637552.003.png
Solving Cubics and Quartics
http://www-staff.it.uts.edu.au/~don/pubs/solving.html
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:
x 4 + ax 3 + bx 2 + cx + d = 0
the subsidiary cubic can be one of the forms:
Christianson-Brown (Christianson, 1991)
y 3 + (4a 2 b - 4b 2 - 4ac + 16d - 3a 4 /4)y 2 /(a 3 - 4ab + 8c) + (3a 2 /16 - b/2)y + (ab/16 - c/8 + a 3 /64) = 0
Descartes-Euler-Cardano (Strong, 1859)
y 3 + (2b - 3a 2 /4)y 2 + (3a 4 /16 - a 2 b + ac + b 2 - 4d)y + (abc - a 6 /64 + a 4 b/8 - a 3 c/4 - a 2 b 2 /4 - c 2 ) = 0
Ferrari-Lagrange (Turnbull, 1947)
y 3 + by 2 + (ac - 4d)y + (a 2 d + c 2 - 4bd) = 0
Neumark (Neumark, 1965)
y 3 - 2by 2 + (ac + b 2 - 4d)y + (a 2 d - abc + c 2 ) = 0
Yacoub-Fraidenraich-Brown (Yacoub & Fraidenraich, 2003)
(a 3 - 4ab + 8c)y 3 + (a 2 b - 4b 2 + 2ac + 16d)y 2 + (a 2 c - 4bc + 8ad)y + (a 2 d - c 2 ) = 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:
3 z 16
2008-08-18 22:29
 
18637552.001.png 18637552.002.png
Solving Cubics and Quartics
http://www-staff.it.uts.edu.au/~don/pubs/solving.html
if
y 3 + py 2 + qy + r = 0
then
z 3 - pz 2 + qz - r = 0 for z = - y
and
z 3 + 2pz 2 + 4qz + 8r = 0 for z = 2y
CHRISTIANSON - BROWN
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 = (4a 2 b - 4b 2 - 4ac + 16d - 3a 4 /4)/(a 3 - 4ab + 8c)
q = 3a 2 /16 - b/2
r = ab/16 - c/8 + a 3 /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
k 4 = y 4 + y 2 e 2 + ye 1 + e 0
or
k 2 = y 2 + e 2 /2 + e 1 /4y
where
e 0 = d - ac/4 + a 2 b/16 - 3a 4 /256
e 1 = c - ab/2 + a 3 /8
e 2 = b - 3a 2 /8
The use of the equation for k 4 is only of benefit if either 'y=0' or else 'y' and 'e 1 ' are negative and 'e 0 ' and
4 z 16
2008-08-18 22:29
 
Solving Cubics and Quartics
http://www-staff.it.uts.edu.au/~don/pubs/solving.html
'e 2 ' are positive.
The value of 'k' is used to calculate the quantities :
g = 4y k 3
h = (6y 2 + e 2 )k 2
These are used to form the quadratic in 'Z' :
Z 2 + gZ/k 4 + (h/k 4 - 2) = 0
The roots of this, 'Z 1 ' and 'Z 2 ', are then use form the quadratics in 'z' :
z 2 - Z 1 *z + 1 = 0
and
z 2 - Z 2 *z + 1 = 0
the roots of which are each used to form the corresponding root of the original quartic :
x = y - kz - a/4
DESCARTES - EULER - CARDANO
This algorithm uses a subsidiary cubic with coefficients :
p = 2b - 3a 2 /4
q = 3a 4 /16 - a 2 b + ac + b 2 - 4d
r = abc - a 6 /64 + a 4 b/8 - a 3 c/4 - a 2 b 2 /4 - c 2
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.
5 z 16
2008-08-18 22:29
 
Zgłoś jeśli naruszono regulamin