The Bisectrix of the Tetrahedron: An Application of Symbolic Computation

Mathematica's algebraic manipulative functions are used to find an 8-term polynomial defining a surface.

by Garry Helzer

The bisectrix of a solid is the envelope of the family of planes that bisect the solid in the sense that they divide it into two equal volumes. The bisectrix of the tetrahedron has been studied by Gunther and Wilker [Gunther and Wilker 1992], who coined the name, and Beyer and Swartz [Beyer and Swartz 1993]. The only published figures of the tetrahedral bisectrix appear in [Beyer and Swartz 1993]. In this article we succeed in obtaining an 88 term polynomial of degree 12 defining a portion of the tetrahedral bisectrix and use this polynomial to correct some inaccuracies in the figures given in [Beyer and Swartz 1993]. 
In detail, the bisectrix of a tetrahedron may be assembled from seven separate pieces. Four of the pieces are triangular shapes cut from a cubic surface which, up to an affine coordinate change, has an equation of the form x  y  z = "const." The four cubic pieces are affinely equivalent to the piece shown in Figure 1. (Two figures are affinely equivalent if one can be mapped to another by an affine transformation—a linear transformation followed by translation. Congruence is more restrictive than similarity which in turn is more restrictive than affine equivalence and the latter is more restrictive than projective equivalence. All triangles are affinely equivalent as are all tetrahedra. Parabolas are not affinely equivalent to hyperbolas or ellipses whereas all conics are projectively equivalent.) 


[FIGURE 1] 
Figure 1 belongs to a particular equilateral tetrahedron used to draw all the figures in this article (to experiment with other tetrahedra, see the electronic supplement). 
The other three pieces of the bisectrix are affinely equivalent to the shape of Figure 2. The pieces are saddle shaped and the piece of Figure 2 lies quite close to the ordinary quadratic saddle z = x  y. 
A previously unsuspected aspect of these saddles is the existence of narrow strips at the edges of the shape in Figure 2 (colored red) that are folded back along a sharp crease. The strips are widest midway between the corners and decrease to points at the corners. Of the four edges of the saddle shape, two opposite edges are folded up and two are folded down. Only the nearest folded up edge is visible in Figure 2. The others are obscured by the rest of the surface. Unlike the quadratic saddle, this saddle is not a ruled surface, but the line segments joining the corners do lie entirely within the surface. 
To assemble the bisectrix the three saddle-shaped pieces are joined along the straight lines connecting their corners. The resulting figure has eight cells, four of which have the folded edges protruding inward (Figure 3). The boundaries of each triplet of folded edges exactly meet the boundaries of one of the four cubic pieces. Figure 4 shows the frontmost cubic piece in place. The other three cubic pieces are obscured by the rest of the figure. The pieces join together smoothly at the edges of the triangular pieces. The surfaces have order of contact two there. 

[FIGURE 2] 

[FIGURE 3] 
New facts about the bisectrix in this article are the existence of the folds, the smoothness of the contact where the blue and red pieces meet in Figure 4 and the fact that saddle pieces are cut from a surface f(x, y, z) = 0 of degree 12. 
Although the stated results are geometric, the purpose of the article is to demonstrate the use of the developing theory of symbolic computation in solving a non-trivial problem: the determination of a polynomial equation satisfied by a saddle piece. The tools used are resultants, factoring algorithms and, most importantly, Gröbner bases. 

[FIGURE 4]

Definitions and Background

In the terminology of [Gunther and Wilker 1992], the four triangular pieces are called (1,3) pieces because they are the envelopes of the set of planes that separate a fixed vertex from the other three vertices. The three saddle pieces are called (2,2) pieces as they are the envelope of the planes separating a fixed pair of vertices from the other pair. 
Beyer and Swartz, using Macsyma on a Vax computer, attempted to find an algebraic equation defining a surface containing a piece congruent to a (2,2) patch, but the computing resources were insufficient for the task. They estimated that the (2,2) pieces satisfied a polynomial equation of total degree lying between 5 and 150. We show how such a polynomial may be computed using Mathematica and a desktop computer (All computations but one were carried out on a Power Macintosh 7100/66 equipped with 24 Megabytes of RAM). The (2,2) pieces turn out to satisfy a polynomial equation of degree 12. Using [Cox, Little, and O'Shea 1992] as a reference, the computations involved are accessible at an advanced undergraduate level.

Triangles and (1,3) pieces

The case of a triangle is shown in Figure 5. For a triangle the bisectrix consists of three segments of hyperbolas. The segments are tangent where they meet. Each segment is the envelope of those bisecting lines which separate one vertex from the other two. The isolated vertex is the center of the hyperbola and the adjacent sides lie on the asymptotes. 

 
Bisectors of a Triangle and the Bisectrix. 
These statements are easy to check if we exploit the relationship between areas and determinants. If the columns of a 2×2 determinant are considered to be vectors in the plane (tails at the origin) then the determinant is the signed area of the parallelogram that has the vectors as adjacent sides. The area is positive if the angle through the parallelogram from the first column to the second is counterclockwise, otherwise it is negative. The area of the triangle with adjacent sides consisting of these vectors is half the determinant, with the same sign convention. 

 

A bisecting line. 
Now consider Figure 6. Let A, B, C be vectors, implemented in Mathematica as two-element lists, specifying the vertices of a triangle. Define parameters p and q by

\(TraditionalForm
    \`P = A + p \((B - A)\) \(, \ \ \ \ \) Q = A + q \((C - A)\) . \)

Then, as p and q vary from zero to one, the points P and Q move from A to B and C respectively. If the line P  Q is to bisect the the area of △  A  B  C then, letting det (X, Y) be the determinant with column vectors X and Y,

({{2 "Area" (A, P, Q),  = det (P - A, Q - A)}, {,  = det (p (B - A), q (C - A))}, {,  = p q det (B - A, C - A)}, {,  = 2 p q "Area" (A, B, C)}})

So the the bisection condition is p  q = 1/2 or q = 1/2  p, where  .5 ≤ p ≤ 1 since 0 ≤ q ≤ 1. 
Now we use a result due to Gunther and Wilker[Gunther and Wilker 1992]: The hyperplanes bisecting the volume of an n-dimensional figure touch its bisectrix at the centroid of their section with the figure. The centroid of a segment is its midpoint so the point

\(TraditionalForm
    \`\(1\/2\) \((P + Q)\) = A + \(p\/2\) \((B - A)\) + \(1\/\(4  p\)\) 
      \((C - A)\)\)

lies on the bisectrix. Choosing any values for the two element lists A, B, C the above formula may be used—exactly as written—as an argument to ParametricPlot to draw a bisectrix segment (See the electronic supplement.) The formula represents an affine transformation which carries the hyperbola 2  p  q = 1 to another hyperbola. The rest of the statements about the bisectrix of a triangle readily follow. 
We have done the triangle in some detail because it is typical and relatively short. For a tetrahedron, add one more vertex, D, and one more parameter r. The (signed) volume of a tetrahedron is one-sixth the corresponding determinant. The condition that the plane bisect the tetrahedron is r = 1/ (2  p  q) and the centroid of a triangle with vertices P, Q, R is (P + Q + R) /3. Choosing \(TraditionalForm
      \`A = \((\(-3\), 3, \(-3\))\), B = \((3, 3, 3)\), C = 
        \((\(-3\), \(-3\), 3)\), D = \((3, \(-3\), \(-3\))\)\) the Mathematica expressions are 
Figure 1 shows this (1,3) piece. 
Although parametric equations for a (1,3) piece are quite simple, plotting is complicated by the fact that the parameter domain is not a rectangle. The domain consists of the portion of the p  q-plane within the rectangle \(TraditionalForm
      \`0 ≤ p ≤ 1, 0 ≤ q ≤ 1\) and above the hyperbola 2  p  q = 1. The pieces may be plotted either by triangulating the domain (see [Wickham-Jones 1994] for details) or by using the package [Helzer].

Resultants and Gröbner bases

The analysis of the (2,2) pieces is more involved. Both Gunther and Wilker [Gunther and Wilker 1992] and Beyer and Swartz [Beyer and Swartz 1993] give parameterizations of the (2,2) patches. Gunther and Wilker give a parameterization by two independent parameters using Hypergeometric functions. Beyer and Swartz provide an algebraic parameterization using three algebraically dependent parameters. 
Beyer and Swartz describe two attempts to find a polynomial equation for a (2,2) patch. One approach uses resultants. 
Given parametric formulas x = (Subscript[f, 1]) (t), y = (Subscript[f, 2]) (t) for a curve in the plane it is sometimes possible to find an implicit formula f (x, y) = 0 for the curve by eliminating the parameter t. If the functions involved are polynomials the parameter may be eliminated using a function, actually a large determinant, called the resultant. If \(TraditionalForm
      \`p \((x, y, t)\) = x - \(f\_1\) \((t)\), q \((x, y, t)\) = y - 
        \(f\_2\) \((t)\)\) then Resultant[p,q,t], if it is not zero, is a polynomial, f (x, y) involving only x and y and the original curve satisfies f (x, y) = 0. The polynomial f (x, y) often has high degree and the set of solutions to f (x, y) = 0 may include points not on the original curve. A discussion of resultants may be found in [Cox, Little, and O'Shea 1992]. 
From their parameterization, Beyer and Swartz find four polynomial equations of the form

\(TraditionalForm\`\(f\_1\)(x, p, q, r) = 0, \ \ \ 
    \(f\_2\)(y, p, q, r) = 0, 
\(f\_3\)(z, p, q, r) = 0, \ \ \ \ 
    \(f\_4\)(p, q, r) = 0\)

where p, q, and r are the parameters. The first three equations are of degree six and the last is of degree four. Using resultants to eliminate r from pairs of these equations gives three equations involving the parameters p and q only. Eliminating q from appropriate pairs of the new equations gives two equations involving the parameter p only. The resultant of these two equations with respect to p is an equation involving only x, y, and z and satisfied by a (2,2) patch. However the last resultant is computationally intractable and the computation was not carried to completion. Beyer and Swartz estimate the degree of their last resultant to be 150. The computational difficulty can be appreciated by noting that the general polynomial of degree n in m variables has Binomial[m+n,m] terms. For n = 150 and m = 3 this is gives potentially 585,276 terms for the last resultant, and the coefficients need not be small. 
In a second approach, C. de Boor used interpolation techniques to find a lower bound on the degree. The conclusion was that the degree is at least five. Thus Beyer and Swartz conclude that the degree of the (2,2) surface patches is between five and 150. We show below that the degree is twelve. 
Another approach to the elimination of parameters in polynomial equations is via Gröbner bases. In fact Gröbner
bases are at the heart of the algorithms used by Mathematica's Solve and Eliminate functions. 
Gröbner bases were introduced in 1965 by Bruno Buchberger [Buchberger 1965] and since then computational commutative algebra has developed into an active area of research. Introductory developments of Gröbner basis theory may be found in [Adams and Loustaunau 1994], [Becker and Weispfenning 1993], and [Cox,David, Little, John, and O'Shea 1992]. The article [Helzer 1995] is a Mathematica oriented overview. 
Given a set of polynomial equations, (Subscript[p, 1]) = 0, (Subscript[p, 2]) = 0, …, (Subscript[p, r]) = 0 in the variables (Subscript[x, 1]), (Subscript[x, 2]), …, (Subscript[x, n]), Buchberger's algorithm produces an equivalent set of equations, called a Gröbner basis, with special properties. The definition of equivalence is in terms of ideals. The ideal associated to a set of equations is the (infinite) set of all polynomials that can be expressed as sums of multiples (by polynomials) of the left hand sides (Subscript[p, 1]), (Subscript[p, 2]), …, (Subscript[p, n]). The set of polynomials is said to be a basis for the ideal or to generate the ideal. Two sets of polynomial equations are equivalent if the polynomials generate the same ideal. This means that each set is obtained from the other by polynomial manipulations. 
The Gröbner bases produced by the built in functions provided in Mathematica version 2 all have the following property. (The version 3.0 functions are more general.) An ordering is chosen for the variables, say \(TraditionalForm
      \`x\_1 ≺ \(x\_2\) 
        \( ≺ ⋯ ≺ \) x\_n\) and a basis G = {(Subscript[g, 1]), (Subscript[g, 2]), …, (Subscript[g, t])} is computed. If the ideal contains any polynomials that depend only upon the first k variables, then (at least) one of the (Subscript[g, i]) will depend only upon the first k variables. 
This applies to parameter elimination as follows. Starting with Beyer and Swartz's four polynomials F = {(Subscript[f, 1]) (x, p, q, r), (Subscript[f, 2]) (y, p, q, r), (Subscript[f, 3]) (z, p, q, r), (Subscript[f, 4]) (p, q, r)} men- tioned above, compute a Gröbner basis with respect to any ordering that puts x, y, and z ahead of p, q, and r, say \(TraditionalForm
      \`x ≺ y ≺ z ≺ p ≺ q ≺ 
        r\). The Gröbner basis computed from F with respect to this ordering will contain at least one polynomial involving only x, y, and z. The parameters p, q and r have been eliminated. The Mathematica expression is 
Unfortunately this calculation is computationally intractable, at least on small machines, and we must use other methods. 
Gröbner bases have another property that we will use below. If the set of polynomials G is a Gröbner basis and f is any polynomial, then, by a generalization of the long division algorithm, one may jointly divide f by the polynomials of G producing a unique remainder called the normal form and written (Subscript[N, G]) (f).Two important properties of the normal form are

◼ The polynomial f is in the ideal generated by G if and only if (Subscript[N, G]) (f) = 0.

◼ If a, b are constants and f, g are polynomials, then

\(TraditionalForm
    \`\(N\_G\) \((a  f + b  g)\) = \(\((a  N)\)\_G\) \((f)\) + 
      \(\((b  N)\)\_G\) \((g)\) . \)

Thus the polynomial

\(TraditionalForm
    \`f \((x\_1, x\_2, …, x\_n)\) = 
      ∑\(a\_\(α\_1, α\_2, …, α\_n\)\) 
      \(x\_1\%α\_1\) \(x\_2\%α\_2\) 
      ⋯  x\_n\%α\_n\)

is in the ideal generated by G if and only if

\(TraditionalForm
    \`0 = \(\(N\_G\)(f(x\_1, x\_2, …, x\_n))\) 
        \( = \ \ \ \ 
          ∑\(a\_\(α\_1, α\_2, …, α\_n\)\) 
              \(\(N\_G\)(
                  \(x\_1\%α\_1\) \(x\_2\%α\_2\) 
                    ⋯  x\_n\%α\_n) . \)\)\)

That is, finding polynomials in the ideal is equivalent to finding linear relations among normal forms of power products. In Mathematica normal forms are computed using the AlgebraicRules function. (In version 3.0 Al- gebraicRules has been superceded by Polynomi- alReduce which computes (Subscript[N, G]) (f) directly. See [Lichtblau 1995].) The AlgebraicRules function computes a Gröbner basis but returns it as a special data structure called an AlgebraicRules object. If arules is an Al- gebraicRules object and f is a polynomial then (Subscript[N, G]) (f) is f /. arules. The arguments of AlgebraicRules are the same as those of GroebnerBasis except that, instead of a list of polynomials in the first argument, AlgebraicRules requires the corresponding list of equations. If F is a valid first argument for GroebnerBa- sis then #==0&/@F is a valid first argument for AlgebraicRules.

Rules for a determinant

In the above computation for a triangle the areas were expressed in terms of determinants but no determinant was actually computed. Using rules for determinants, one determinant was expressed as a multiple of another. We wish to automate this process for 3×3 determinants. To do so we use the rule set 
Given any three expressions X, Y, Z that are expressed as simple linear combinations of w[1], w[2], w[3] the above rules will reduce d[X, Y, Z] to a multiple of d[w[1], w[2], w[3]].

A parameterization for a (2,2) piece

Figure 7 is a schematic of a tetrahedron, A  B  C  D, and a plane intersecting edges of the tetrahedron at P, Q, R, S. 

 

A plane separating two edges. 
As before, introduce a parameters 0 ≤ p, q, r, s ≤ 1 specifying the positions of P, Q, R, S respectively. Let p and r be zero at A and q and s be zero at D. Thus 
where A, B, C, D are unassigned symbols. Define vectors (Subscript[W, 1]) = B - A, (Subscript[W, 2]) = C - A, (Subscript[W, 3]) = D - A. Then the volume of the tetrahedron is (1/6) det((Subscript[W, 1]), (Subscript[W, 2]), (Subscript[W, 3])). To convert expressions involving A, B, C, D to expressions involving (Subscript[W, 1]), (Subscript[W, 2]), (Subscript[W, 3]) and compute the resulting multiples of det((Subscript[W, 1]), (Subscript[W, 2]), (Subscript[W, 3])) use the definitions: 
The function WForm rewrites its argument as a simple linear combination of w[1], w[2], w[3]. 
Now we have tools to derive an algebraic parametric expression for a (2,2) piece of the bisectrix. The first constraint on the four parameters p, q, r, s is that the four points they define lie in a plane. Four points define a tetrahedron. The four points lie in a plane if the tetrahedron is degenerate—has a volume of zero. Thus we require that the polynomial 
be zero. The second constraint is that the plane bisect the tetrahedron. The portion of the tetrahedron in Figure 7 above the plane may be decomposed into three tetrahedra: C  R  S  Q, C  R  S  B, and B  S  P  R. The sum of these three volumes, taken with the appropriate signs, should be half the volume of the tetrahedron. Thus the polynomial 
should be zero. These two constraints insure that P, Q, R, S lie in a bisecting plane. Next we need the centroid of the quadrangle P  Q  R  S. 
Unlike segments and triangles, the centroid of a quadrangle is not the vector average of its vertices (see [Coxeter 1989] for a discussion of this point). We can get a system of linear equations defining the centroid as follows. The diagonal P  Q divides the quadrangle into two triangles. The centroid of the quadrangle (which is the same as the center of gravity of thin sheet with the same shape) will lie on the segment joining the centroids of the two triangles. It will also lie on the segment joining the centroids of the two triangles obtained by using the diagonal S  R to divide the quadrangle. So introduce two new parameters, t and u, by defining 
The components of the vector equation Segment1-Segment2==0 are three scalar equations for the two unknown parameters t, u. The right hand sides of these equations are extracted by the expressions 
We now have five equations in the six variables p, q, r, s, t, u. To obtain parametric expressions for a particular (2,2) piece, choose particular values for the vertices A, B, C, D and use Segment1 or Segment2 to define x, y, z. A choice that works well is 
We now have 11 polynomials in the 9 variables p, q, r, s, t, u, x, y, z which, when equated to zero, define x, y, z as functions of six algebraically dependent parameters. The Gröbner basis of the ideal generated by these polynomials with respect to the ordering \(TraditionalForm
      \`z ≺ y ≺ x ≺ u ≺ t ≺ 
        s ≺ r ≺ q ≺ p\) contains a polynomial of smallest degree involving only x, y, z. This polynomial defines a surface containing the (2,2) piece. However, the computation is intractable, at least on small machines. But if the order is reversed, the corresponding Gröbner basis is relatively easy to compute. The expression that computes it is 
This basis consists of 14 polynomials forming a `triangular' system that allows us to compute the variables successively, starting with a given p, q, both between zero and one. 
Switching to equations, the first basis element is of the form

\(TraditionalForm
    \`a \((p, q)\) r\^2 + b \((p, q)\) r + c \((p, q)\) = 0. \)

This allows r to be computed given p, q. Analysis shows that, for most p, q between zero and one, exactly one root of the quadratic lies between zero and one. (For some values of p and q the equation degenerates to 0 = 0.) 
The next three equations are of the form

a (p, q, r) s + b (p, q, r) = 0.

If an r has been determined from the first equation, at least one of these equations will give a unique value for s. If the first equation is degenerate, these equations give values for both r and s. There is one exception to the previous statement. If p = q = (1/2), then the first four equations reduce to r + s = 1—for any choice of r, 0 ≤ r ≤ 1 one may find a corresponding s, 0 ≤ s ≤ 1. This gives the surprising fact (also noticed in [Gunther and Wilker 1992]): any plane cutting opposite edges of a tetrahedron at their midpoints divides the tetrahedron into two equal volumes. 
The next six equations are all of the form

a (p, q, r, s) t + b (p, q, r, s) = 0

and, in aggregate, determine a unique t from any values of p, q, r, s computed using the previous equations. 
The eleventh equation may be written u = a (p, q, r, s, t). It is the only equation in which u appears and it may be discarded, along with the variable u which is no longer needed. 
The last three equations may be written in the form

\(TraditionalForm
    \`x = a \((p, q, t)\) \(, \ \ \) y = b \((r, s, t)\) \(, \ \ \) z = c 
      \((p, q, r, s, t)\)\)

and uniquely determine x, y, z when p, q are chosen in appropriate ranges and r, s, t are computed from the preceding equations. 
TheBasis may be used to plot a (2,2) piece parametrically, but the result is not very satisfactory because of a discontinuity at p = q = (1/2) and a singularity of the piece itself near the boundary. We do not pursue the subject but proceed directly to the computation of the degree of the surface.

The degree of the Surface

We can search for polynomials of the form f(x, y, z) in the ideal generated by our Gröbner basis by computing the normal forms of the power products of the form (x^a) (y^b) (z^c) and looking for linear relations among them. Each such linear relation gives a polynomial in the ideal. An approach would be to compute some normal forms and check for relations. If no relations are found add some more normal forms and try again. Continue until a relation is found (see [Faugère, Gianni, Lazard, and Mora 1993]). The linear systems to be solved grow in size and resources are exhausted before a relation is found. But a variant process succeeds. First use linear equations to find low degree polynomials with some of the parameters eliminated, then eliminate the rest with resultants. 
The Gröbner basis computed above contains the polynomial

y - 1 + 2  r + 2  t - 2  r  t - 2  s  t .

Thus the ideal contains quadratic polynomials free of p and q. Using the technique described below a second quadratic free if p and q may be found. The resultant of these two quadratics with respect to the variable r is a cubic. So the ideal contains a cubic free of p, q and r. There are no more cubics (except for constant multiples) but the computation below produces a large number of quartics free of p, q and r. 
First to compute normal forms we need an Alge- braicRules object. 
We also need a utility to generate power products. Given a list of variables and an integer n the function PPSet will return the list of power products, in the given variables, of total degree less than or equal to n. 
The procedure is

◼ Generate a set of power products (Subscript[p, 1]), (Subscript[p, 2]), …, (Subscript[p, s]). 
◼ Compute the normal forms \(TraditionalForm
      \`\(N\_G\) \((p\_1)\), \(N\_G\) \((p\_2)\), …, \(N\_G\) 
        \((p\_s)\)\). 
◼ Generate symbolic coefficients (Subscript[a, 1]), (Subscript[a, 2]), …, (Subscript[a, s]), form the big polynomial

\(TraditionalForm
    \`\(a\_1\) \(N\_G\) \((p\_1)\) + \(a\_2\) \(N\_G\) \((p\_2)\) + 
      ⋯ + \(a\_s\) \(N\_G\) \((p\_s)\)\)

and collect like power products. 
◼ Equate the coefficients of the power products to zero obtaining a homogeneous system of linear equations in the unknown (Subscript[a, i]). 
◼ Solve the system. For each solution (Subscript[a, 1]), (Subscript[a, 2]), …, (Subscript[a, s]), the polynomial \(TraditionalForm
      \`\(a\_1\) p\_1 + \(a\_2\) p\_2 + ⋯ + \(a\_s\) p\_s\) is in the ideal.

The linear system is homogeneous and must be solved in exact arithmetic. It is a sparse system and Solve is a better choice of solver than LinearSolve. If there is a single solution it is the trivial solution (all (Subscript[a, i]) = 0), otherwise there are an infinity of solutions with some of the (Subscript[a, i]) arbitrary and the rest written in terms of the arbitrary ones. A particular solution may be obtained by setting one of the arbitrary coefficients to one and the rest to zero. A linearly independent set of solutions may be obtained by doing this for each arbitrary coefficient in turn. This is the purpose of the replacement rule involving the identity matrix below. 
There is a technical problem with step (3) above. When collecting like power products we want the coefficient of p in (Subscript[a, 1]) p + ((Subscript[a, 3]) + (Subscript[a, 5])) p  q to be (Subscript[a, 1]). But Mathematica's Coefficient function considers the coefficient to be (Subscript[a, 1]) + ((Subscript[a, 3]) + (Subscript[a, 5])) q. To get around this problem we replace each power product (p^a) (q^b) (r^c) (s^d) (t^e) by the single variable (Subscript[x, a  b  c  d  e]), or X[a,b,c,d,e] in the expressions below. The vector (a, b, c, d, e) is called the multidegree of the power product and is computed by the function Mdeg below. The replacement rule subrule changes a single power product and the function subscript applies subrule to all the power products in an Expanded polynomial. 
The following sequence of Mathematica expressions computes a linearly independent set of quartics in the ideal and free of the parameters p, q and r. 
The list lowdegpolys contains 21 polynomials, one cubic and 20 quartics. The linear system generated consists of 128 equations in 126 unknowns. The polynomial BIGpoly has 3,804 terms. The computation takes about 8 minutes wall clock time on a Power Macintosh 7100/66 (24 Megabytes RAM). 
Computing resultants for two pairs of these quartics gives two polynomials in the ideal free of p, q, r and s and computing the resultant of these two gives a polynomial in the ideal free of all parameters. 
The polynomial res5 has 2,959 terms but decomposes as a product of 10 factors. Only the last of these factors lies in the ideal. It may be computed as 
The polynomial satisfied by the (2,2) piece of the bisectrix has degree 12 and 88 terms. It is shown in Figure 8.

\(TraditionalForm
    \`l27\ \(x\ \) y - \(108\ \) x\^3\ y + \(162\ \) x\^5\ y - 
      \(108\ \) x\^7\ y + \(27\ \) x\^9\ y - \(108\ \) x\ y\^3 - 
      \(2357\ \) x\^3\ y\^3 - \(1617\ \) x\^5\ y\^3 - \(15\ \) x\^7\ y\^3 + 
      x\^9\ y\^3 + \(162\ \) x\ y\^5 - \(1617\ \) x\^3\ y\^5 - 
      \(96\ \) x\^5\ y\^5 - x\^7\ y\^5 - \(108\ \) x\ y\^7 - 
      \(15\ \) x\^3\ y\^7 - x\^5\ y\^7 + \(27\ \) x\ y\^9 + x\^3\ y\^9 - 
      \(27\ \) z + \(108\ \) x\^2\ z - \(162\ \) x\^4\ z + 
      \(108\ \) x\^6\ z - \(27\ \) x\^8\ z + \(108\ \) y\^2\ z + 
      \(5133\ \) x\^2\ y\^2\ z + \(6297\ \) x\^4\ y\^2\ z + 
      \(759\ \) x\^6\ y\^2\ z - \(9\ \) x\^8\ y\^2\ z - \(162\ \) y\^4\ z + 
      \(6297\ \) x\^2\ y\^4\ z + \(1720\ \) x\^4\ y\^4\ z + 
      \(33\ \) x\^6\ y\^4\ z + \(108\ \) y\^6\ z + \(759\ \) x\^2\ y\^6\ z + 
      \(33\ \) x\^4\ y\^6\ z - \(27\ \) y\^8\ z - \(9\ \) x\^2\ y\^8\ z - 
      \(3255\ \) x\ \(y\ \) z\^2 - \(7395\ \) x\^3\ \(y\ \) z\^2 - 
      \(1557\ \) x\^5\ \(y\ \) z\^2 - \(81\ \) x\^7\ \(y\ \) z\^2 - 
      \(7395\ \) x\ y\^3\ z\^2 - \(2517\ \) x\^3\ y\^3\ z\^2 - 
      \(357\ \) x\^5\ y\^3\ z\^2 - \(3\ \) x\^7\ y\^3\ z\^2 - 
      \(1557\ \) x\ y\^5\ z\^2 - \(357\ \) x\^3\ y\^5\ z\^2 - 
      \(2\ \) x\^5\ y\^5\ z\^2 - \(81\ \) x\ y\^7\ z\^2 - 
      \(3\ \) x\^3\ y\^7\ z\^2 + \(487\ \) z\^3 + \(2691\ \) x\^2\ z\^3 + 
      \(837\ \) x\^4\ z\^3 + \(81\ \) x\^6\ z\^3 + \(2691\ \) y\^2\ z\^3 - 
      \(2787\ \) x\^2\ y\^2\ z\^3 - \(411\ \) x\^4\ y\^2\ z\^3 + 
      \(27\ \) x\^6\ y\^2\ z\^3 + \(837\ \) y\^4\ z\^3 - 
      \(411\ \) x\^2\ y\^4\ z\^3 - \(6\ \) x\^4\ y\^4\ z\^3 + 
      \(81\ \) y\^6\ z\^3 + \(27\ \) x\^2\ y\^6\ z\^3 + 
      \(5829\ \) x\ \(y\ \) z\^4 + \(1818\ \) x\^3\ \(y\ \) z\^4 + 
      \(81\ \) x\^5\ \(y\ \) z\^4 + \(1818\ \) x\ y\^3\ z\^4 + 
      \(399\ \) x\^3\ y\^3\ z\^4 + \(3\ \) x\^5\ y\^3\ z\^4 + 
      \(81\ \) x\ y\^5\ z\^4 + \(3\ \) x\^3\ y\^5\ z\^4 - \(2205\ \) z\^5 - 
      \(1026\ \) x\^2\ z\^5 - \(81\ \) x\^4\ z\^5 - \(1026\ \) y\^2\ z\^5 - 
      \(375\ \) x\^2\ y\^2\ z\^5 - \(27\ \) x\^4\ y\^2\ z\^5 - 
      \(81\ \) y\^4\ z\^5 - \(27\ \) x\^2\ y\^4\ z\^5 - 
      \(153\ \) x\ \(y\ \) z\^6 - \(27\ \) x\^3\ \(y\ \) z\^6 - 
      \(27\ \) x\ y\^3\ z\^6 - x\^3\ y\^3\ z\^6 + \(81\ \) z\^7 + 
      \(27\ \) x\^2\ z\^7 + \(27\ \) y\^2\ z\^7 + 
      \(9\ \) x\^2\ y\^2\ z\^7\)

\(TraditionalForm\`\(
FIGURE\ \(8 : \)\)\)

The polynomial bisectrix. 
The polynomial for any other (2,2) patch of any tetrahedron may be found by an affine change of coordinates. It should be checked that the polynomial bisectrix does indeed lie in the ideal. This is the case if and only if the normal form bisectrix/.NFrules is zero. This is indeed the case but it is the most difficult computation described here, consuming about one and a half hours cpu time on a Sun 360 workstation with 36 Megabytes of RAM.

Drawing a (2,2) Patch and Checking the Contact

The difficulty in accurately drawing a (2,2) patch is that the arcs forming the outer edge of the figure consist of sharp cusp-like singular points where the surface is folded back (Figure 2). Mathematica displays a surface as a collection of polygons. Unless vertices of the appropriate polygons fall directly on the arcs of singular points, the fold will not be apparent. Some sort of smooth bend will probably appear instead. Given parametric equations for the
surface it will not usually be apparent which parameter values correspond to the singular points. However, if we have an implicit equation for the surface in addition to parametric equations we can write down equations that give the set of points in the parameter domain that map to the singular points on the surface. Suppose that \(TraditionalForm
      \`x = x \((a, b)\), y = y \((a, b)\), z = z \((a, b)\)\) are parametric equations for the surface and that f (x, y, z) = 0 is an equation of the surface. The singular points are the points at which the gradient vector ((Subscript[f, x]), (Subscript[f, y]), (Subscript[f, z])) is the zero vector. Thus the solutions of

({{(Subscript[f, x])(x(a, b), y(a, b), z(a, b)),  = 0}, {(Subscript[f, y])(x(a, b), y(a, b), z(a, b)),  = 0}, {(Subscript[f, z])(x(a, b), y(a, b), z(a, b)),  = 0}})

are the points in the parameter domain that map to the singular points of the surface. For Figure 2, f (x, y, z) was taken to be the polynomial of Figure 8. Since the (2,2) patch used is roughly horizontal only the third equation above was needed. For parametric equations the hypergeometric parameterization from [Gunther and Wilker 1992] was used. The domain of the hypergeometicparameterization is the region of the first quadrant bounded by the axes and the interior of a certain ellipse together with symmetrically placed regions in the other quadrants. The Mathematica implementation of the parameterization is somewhat long since the formulas are indeterminate forms at the unit points of the axes. The parameterization and domain boundary are given in the appendix. The package SurfaceClip.m (see [Helzer]) was used to separately plot the regions inside and outside the arcs of singular points and the results merged with Show. The other (2,2) patches visible in Figure 3 were produced by applying homogeneous transformations from [Wickham-Jones 1994]. 
To check the contact between a (1,3) piece and a (2,2) piece we can use the parameterization param13 above and the polynomial bisectrix of Figure 8. If one surface is given parametrically by "\!\(\*StyleBox[\"\\\"x\\\"\", \"TB\"]\)"(p, q) and a second surface is given implicitly by f("\!\(\*StyleBox[\"\\\"x\\\"\", \"TB\"]\)")=0 (here "\!\(\*StyleBox[\"\\\"x\\\"\", \"TB\"]\)"=(x, y, z)), then the surfaces have contact of order 2 at a common point "\!\(\*StyleBox[\"\\\"x\\\"\", \"TB\"]\)"((Subscript[p, 0]), (Subscript[q, 0])) if all partial derivatives of f("\!\(\*StyleBox[\"\\\"x\\\"\", \"TB\"]\)"(p, q)) of order less than three are zero at ((Subscript[p, 0]), (Subscript[q, 0])) but some third order derivative is not zero. The particular (1,3) piece and (2,2) piece at hand meet along an arc which is the image of the curve segment \(TraditionalForm
      \`2 \( p  q\) = 1,  .5 ≤ p, q ≤ 1\) in the parameter domain. For these pieces the composition f("\!\(\*StyleBox[\"\\\"x\\\"\", \"TB\"]\)"(p, q)) is 
and the first partial derivative with respect to p along the arc is 
which is zero. The vanishing of the other partial derivatives can be checked the same way.

Appendix: The Gunther-Wilker parameterization

The hypergeometric parameterization of the `horizontal' (2,2) piece for the tetrahedron with vertices (-3, 3, -3), (3, 3, 3), (-3, -3, 3), (3, -3, -3) is given by {x[a,b], y[a,b], z[a,b]} defined below. 
The domain of definition of the parameterization is that portion of the a, b plane where the function ClipEdge[a,b] is nonpositive. 
The zero set of ClipEdge consists of four elliptical arcs joining the unit points on the axes. The parameters corresponding to the singular points on the (2,2) piece lie quite close to \(TraditionalForm
      \`| a | 
          \(+\( | \(b | \)\)\) = 1\).


References 
Adams, W. and Loustaunau, P., An Introduction to Gröbner Bases, American Mathematical Society, Providence, 1994. 
Becker, T. and Weispfenning, V., Gröbner Bases, A Computational Approach to Commutative Algebra, Springer-Verlag, New York, Berlin, London, Paris, Tokyo, Hong Kong, 1993. 
Beyer, W. A., and Swartz, B., Bisectors of Triangles and Tetrahedra, The American Mathematical Monthly, 100, 626–640, (1993). 
Buchberger, B., Ein Algorithmus zum Affinden der Basiselementen des Restklassenringes nach einem null dimensionalen Polynomial ideal, Ph.D. Thesis, Inst. University of Innsbruck, Innsbruck, Austria, 1965. 
Buchberger, B., Gröbner Bases: An Algorithmic Method in Polynomial Ideal Theory, Recent Trends in Multidimensional Systems Theory, 184–232. D. Reidel Dordrecht Ed. N. K. Bose 
Buchberger, B., Collins, G. E. and Kutzler, B., Annual Review of Computer Science, Algebraic Methods in Geometric Reasoning, 1988, 3, 86–119. 
Cox, D. Little, J. and O'Shea, D., Ideals, Varieties, and Algorithms, Springer-Verlag, New York, Berlin, London, Paris, Tokyo, Hong Kong, 1992. 
Coxeter, H. S. M., Introduction to Geometry, John Wiley and Sons, inc, New York, 1989. 
Faugère, J. C., Gianni, P., Lazard, D., and Mora, T., Efficient Computation of Zero-dimensional Gröbner Bases by Change of Ordering J. Symbolic Computation, 1993, 16, 329–344. 
Gunther, G. and Wilker, J. B. The Bisectrix of a Tetrahedron Mathematika, 1992, 93–103, 39. 
Helzer, G., Gröbner Bases, The Mathematica Journal, 1995, 67–73, 5, No. 1. 
Helzer, G., SurfaceClip.m, A Mathematica package. 
Lichtblau, D., Gröbner Bases in Mathematica 3.0, The Mathematica Journal, 1996, 81–88, 6 No. 4. 
Wickham-Jones, T.,  Mathematica Graphics, Techniques and Applications (and accompanying software), Springer-Verlag/TELOS, New York, 1994. 


About the Author 
Garry Helzer
Department of Mathematics
University of Maryland
College Park, MD 20742
g_helzer@math.umd.edu