Implicitization via the Gröbner Walk

Daniel Lichtblau
Wolfram Research, Inc.
100 Trade Center Dr.
Champaign IL USA 61820
danl@wolfram.com

ACA 2007
Approximate Algebraic Computation Session
Oakland University
Rochester MI
July 19-22, 2007

Abstract: The Gröbner walk is a useful method for conversion from a "simple" Gröbner basis to a different one in a desired term order.Various issues along the way include coefficient swell (similar to that seen in the classical Buchberger algorithm), polynomials with many initial elements in the "end-game" phase, and the like.We take as benchmark the implicitization of the 32 bicubic parametric patches in the Newell (Utah) teapot.We will see how the Gröbner walk can be used to implicitize all of them, using in some cases approximate arithmetic and an early abort strategy, with a result that can be certified a posteriori.To the best of my knowledge the four spout patches have never before been amenable to a Gröbner basis approach. We show some other nontrivial examples from the implicitization literature.

Introduction

What is the Gröbner walk?

It is a method to convert from an "easy" Gröbner basis to a desired one that might be more difficult to compute directly.

What is it used for?

Many things. Among them, implicitization of parametrized curves and surfaces. This is the application we consider in this talk.

How does it work?

In brief, one moves between cones in the Gröbner fan, beginning with the initial order and moving along a "line" in weight vector space toward the terminal (desired) term order. Good references are (Collart et al 1997) and (Amrhein et al 1997).

The term orders are regarded as weight matrices with each vector comprised of integers. General theory then says they must be square, have full rank, and have the first nonzero in each column be positive.

One could opt for more generality. We do not.

One starts by computing a Gröbner basis with respect to the start order.

Polynomial "initials" are the sum of those monomials in a given polynomial that are equal with respect to the first weight vector.

Each move is relatively simple because it involves only ideals of initials of polynomials in the current basis, and these are mostly monomials.

Convert the initials ideal to get into the neighboring cone.

Use this to carry along the full polynomial set ("lifting" step).

(Optionally) interreduce.

If not yet in terminal cone, restart by looking for the next neighboring cone.

Complexity issues

There are many ways for this process to bog down. The point of this talk is to discuss some, and ways to make them less intolerable. Here are a few items.

The initial basis may be slow to compute.

One typically uses a degree-reverse-lexicographic term order but of course any valid term order might be used.

Coefficient swell as we go through cones can be a problem.

We might encounter cone traversals where the initials ideal contains polynomials with many terms. This more or less defeats the purpose of the Gröbner walk. It is especially common once we are on the boundary of the final cone.

For purposes of variable elimination (needed e.g. for implicitization) it might be unnecessary and inefficient to compute the full basis with respect to a given elimination ordering.

Interreduction of intermediate bases is not actually necessary.

If not done, then one encounters "false cones" problem wherein ties in new initials do not actually get us always to a new cone.

In practice the false cones problem is MUCH worse for efficiency than the slowness of interreduction. We always interreduce.

Tactics we use or otherwise discuss

Perturbation of initial weight vector

This is important because if not done we are more likely to bog down in intermediate Gröbner basis computations involving initials. Reason: We may make a cone crossing that is not through the interior of a proper face, and in this circumstance the initials may have many more monomials than would otherwise be the case. The point of the  Gröbner walk is to work with ideals generated by small initials.

Handling of final weight vector

This is important in those (very common) cases where the first weight vector for the final cone puts us at a border of many cones, by virtue of being not in the interior of that cone. There are many ways discussed in the literature for dealing with this.

Do nothing and bear the consequences of a slow final conversion.

(Amrhein et al 1997) Heuristic perturbation. One adds decreasing multiples of later weight vectors to the first one. With probability 1 this puts us in the interior of SOME cone. If it turns out not to be the correct cone we make the multiples smaller and repeat.

(Amrhein et al 1997) Other methods such as fractal walk.

(Tran 2000) Deterministic perturbation. One computes a term order that is a priori in the interior of the final cone.

(Fukuda et al 2007) Use a symbolic perturbation.

My synopsis

Heuristic or simulated perturbation are the methods of choice.

Deterministic perturbation can generate weight vectors with HUGE integers, and will hang in exponent arithmetic, if there are more than 5 or so variables present.

Doing nothing will hang for same reason as Buchberger algorithm.

Other methods are more trouble for less gain than heuristic and simulated perturbation.

Elimination order

As we do elimination in the examples herein, we use a destination term order well suited for this.

Weight monomials with elimination variables higher than monomials without them.

Break ties based on degree reverse lexicographic ordering of the entire varaible set.

This term order is discussed in (Cox et al 2007).

Initial order

It is helpful to use a "sensible" initial ordering.

Could use degree reverse lexicographic.

We do this, and pay a cost in the need to compute an initial basis with respect to this term order.

If we begin with a rational parametrization, could use an initial order such that input, with denominators cleared, is already a Gröbner basis. This would have a weight matrix with initial vector (0,...,0,1,1,...,1).

In practice such orders can sometimes cause huge inefficiency in the walk process. We will not use them.

Early abort for elimination orders

In implicitization and other elimination tasks we can check at intermediate stages of the walk to see if a basis happens to be an elimination basis for the polynomials.

This is defined as a basis G for the ideal I such that, if the noneliminated variable set is Y and the field is K, then  𝕂[Y]⋂G generates the elimination ideal K[Y]⋂I .

If so, we may often use it just as we would the basis we were trying for. We simply take those polynomials in the basis that are free of the variables we were eliminating (see next subsection for why this works). For many purposes, e.g. surface implicitization, this is a single polynomial.

This applies even though the intermediate Gröbner basis is NOT a Gröbner basis with respect to the provided termination ordering. Hence it is a bona fide "short cut".

Early version of this called "sudden death" (Amrhein et al 1997).

Recent version called "ideal-specific elimination order" (Tran 2004).

Me, I've always viewed it as "opportunistic intermediate processing".

Algorithmic underpinnings of this early abort strategy

The condition is simple to check algorithmically due to a basic result.

Proposition (Tran 2004): A Gröbner basis G for an ideal I  is an elimination basis for a set of variables X if, for each polynomial g∈G, either the head term of g contains variables in X, or else g is entirely free of variables in X.

Proof: As above, denote the complementary variable set as Y and the field as 𝕂. We refer to polynomials in 𝕂[Y]⋂G as (Subscript[G, Y]). Take f∈𝕂[Y]⋂I. Without loss of generality assume f is reduced with respect to (Subscript[G, Y]) and the given term order. If f≠0 then, since G is a Gröbner basis with respect to that term order, f can be further reduced by some g∈G-(Subscript[G, Y]). But the head term condition shows no such g can reduce f, so f=0. So f∈(Subscript[G, Y]). Hence (Subscript[G, Y]) generates 𝕂[Y]⋂I.□

Approximate arithmetic

In some cases we will use finite precision to "guess" a result. This allows us to avoid intermediate coefficient swell. The result can be verified after the fact. We show this in some examples below.

Examples

I will take as benchmarks several examples from the literature of implicitization of bicubic polynomial patch surface parametrizations (Tran 2007). This by no means captures the full range of examples of interest to Gröbner walk enthusiasts. But it does give a means to gauge developments and strategies for some of the issues raised above. These will have Cartesian coordinates {x,y,z} and surface parameters {s,t}. We will need to eliminate the latter. For further examples I will use an algebraic but nonrational parametrization from previous work of my own (Lichtblau 2006), and a rational nonpolynomial parametrization from (Tran 2004) .

Example 1 (Bicubic parametrization attributed to C. Hoffmann)

elims={s,t};
vars={x,y,z};

HoffmannPolys={3 t (((t-1))^2)+(((s-1))^3)+3 s-x,3 s (((s-1))^2)+(t^3)+3 t-y,-3 s ((s^2)-5 s+5) (t^3)-3 ((s^3)+6 (s^2)-9 s+1) (t^2)+t (6 (s^3)+9 (s^2)-18 s+3)-3 s (s-1)-z};

Example 2 (Tran 2007)

TranPolys={-17-27 s+36 (s^2)-19 (s^3)-45 s (t^2)+81 (s^2) (t^2)-54 (s^3) (t^2)+30 s (t^3)-54 (s^2) (t^3)+36 (s^3) (t^3)+10 x,198 s-369 (s^2) t+246 (s^3) t-198 (t^2)+369 (s^2) (t^2)-246 (s^3) (t^2)+100 y,-57-81 (s^2)+42 (s^3)+99 (t^2)-81 s (t^2)-108 (s^2) (t^2)+90 (s^3) (t^2)-66 (t^3)+54 s (t^3)+72 (s^2) (t^3)-60 (s^3) (t^3)+40 z};

Example 3 (Newell)

This is the 32 patches of Newell's notorious "Utah teapot". I grabbed the data points from a MathSource item, rationalized the coordinates, obtaining parametric polynomial systems. As they are long I only show two.

In past some of these have proven difficult for Gröbner bases and Dixon resultants. Robert Lewis' program Fermat now can handle them using the latter approach, with times ranging from split second to approximately a minute (Lewis 2007).

patches3D={{-7/5+(231*s^2)/125-(56*s^3)/125+(3*t)/16-(99*s^2*t)/400+(3*s^3*t)/50-(39*t^2)/80+(1287*s^2*t^2)/2000-(39*s^3*t^2)/250+t^3/5-(33*s^2*t^3)/125+(8*s^3*t^3)/125+x,(294*s)/125-(63*s^2)/125-(56*s^3)/125-(63*s*t)/200+(27*s^2*t)/400+(3*s^3*t)/50+(819*s*t^2)/1000-(351*s^2*t^2)/2000-(39*s^3*t^2)/250-(42*s*t^3)/125+(9*s^2*t^3)/125+(8*s^3*t^3)/125+y,-12/5-(63*t)/160+(63*t^2)/160+z},{(294*s)/125-(63*s^2)/125-(56*s^3)/125-(63*s*t)/200+(27*s^2*t)/400+(3*s^3*t)/50+(819*s*t^2)/1000-(351*s^2*t^2)/2000-(39*s^3*t^2)/250-(42*s*t^3)/125+(9*s^2*t^3)/125+(8*s^3*t^3)/125+x,7/5-(231*s^2)/125+(56*s^3)/125-(3*t)/16+(99*s^2*t)/400-(3*s^3*t)/50+(39*t^2)/80-(1287*s^2*t^2)/2000+(39*s^3*t^2)/250-t^3/5+(33*s^2*t^3)/125-(8*s^3*t^3)/125+y,-12/5-(63*t)/160+(63*t^2)/160+z},{7/5-(231*s^2)/125+(56*s^3)/125-(3*t)/16+(99*s^2*t)/400-(3*s^3*t)/50+(39*t^2)/80-(1287*s^2*t^2)/2000+(39*s^3*t^2)/250-t^3/5+(33*s^2*t^3)/125-(8*s^3*t^3)/125+x,(-294*s)/125+(63*s^2)/125+(56*s^3)/125+(63*s*t)/200-(27*s^2*t)/400-(3*s^3*t)/50-(819*s*t^2)/1000+(351*s^2*t^2)/2000+(39*s^3*t^2)/250+(42*s*t^3)/125-(9*s^2*t^3)/125-(8*s^3*t^3)/125+y,-12/5-(63*t)/160+(63*t^2)/160+z},{(-294*s)/125+(63*s^2)/125+(56*s^3)/125+(63*s*t)/200-(27*s^2*t)/400-(3*s^3*t)/50-(819*s*t^2)/1000+(351*s^2*t^2)/2000+(39*s^3*t^2)/250+(42*s*t^3)/125-(9*s^2*t^3)/125-(8*s^3*t^3)/125+x,-7/5+(231*s^2)/125-(56*s^3)/125+(3*t)/16-(99*s^2*t)/400+(3*s^3*t)/50-(39*t^2)/80+(1287*s^2*t^2)/2000-(39*s^3*t^2)/250+t^3/5-(33*s^2*t^3)/125+(8*s^3*t^3)/125+y,-12/5-(63*t)/160+(63*t^2)/160+z},{-3/2+(99*s^2)/50-(12*s^3)/25-(3*t)/4+(99*s^2*t)/100-(6*s^3*t)/25+t^3/4-(33*s^2*t^3)/100+(2*s^3*t^3)/25+x,(63*s)/25-(27*s^2)/50-(12*s^3)/25+(63*s*t)/50-(27*s^2*t)/100-(6*s^3*t)/25-(21*s*t^3)/50+(9*s^2*t^3)/100+(2*s^3*t^3)/25+y,-12/5+(63*t)/40-(3*t^3)/40+z},{(63*s)/25-(27*s^2)/50-(12*s^3)/25+(63*s*t)/50-(27*s^2*t)/100-(6*s^3*t)/25-(21*s*t^3)/50+(9*s^2*t^3)/100+(2*s^3*t^3)/25+x,3/2-(99*s^2)/50+(12*s^3)/25+(3*t)/4-(99*s^2*t)/100+(6*s^3*t)/25-t^3/4+(33*s^2*t^3)/100-(2*s^3*t^3)/25+y,-12/5+(63*t)/40-(3*t^3)/40+z},{3/2-(99*s^2)/50+(12*s^3)/25+(3*t)/4-(99*s^2*t)/100+(6*s^3*t)/25-t^3/4+(33*s^2*t^3)/100-(2*s^3*t^3)/25+x,(-63*s)/25+(27*s^2)/50+(12*s^3)/25-(63*s*t)/50+(27*s^2*t)/100+(6*s^3*t)/25+(21*s*t^3)/50-(9*s^2*t^3)/100-(2*s^3*t^3)/25+y,-12/5+(63*t)/40-(3*t^3)/40+z},{(-63*s)/25+(27*s^2)/50+(12*s^3)/25-(63*s*t)/50+(27*s^2*t)/100+(6*s^3*t)/25+(21*s*t^3)/50-(9*s^2*t^3)/100-(2*s^3*t^3)/25+x,-3/2+(99*s^2)/50-(12*s^3)/25-(3*t)/4+(99*s^2*t)/100-(6*s^3*t)/25+t^3/4-(33*s^2*t^3)/100+(2*s^3*t^3)/25+y,-12/5+(63*t)/40-(3*t^3)/40+z},{-2+(66*s^2)/25-(16*s^3)/25+(3*t^2)/2-(99*s^2*t^2)/50+(12*s^3*t^2)/25-t^3+(33*s^2*t^3)/25-(8*s^3*t^3)/25+x,(84*s)/25-(18*s^2)/25-(16*s^3)/25-(63*s*t^2)/25+(27*s^2*t^2)/50+(12*s^3*t^2)/25+(42*s*t^3)/25-(9*s^2*t^3)/25-(8*s^3*t^3)/25+y,-9/10+(27*t)/20-(27*t^2)/40+(3*t^3)/40+z},{(84*s)/25-(18*s^2)/25-(16*s^3)/25-(63*s*t^2)/25+(27*s^2*t^2)/50+(12*s^3*t^2)/25+(42*s*t^3)/25-(9*s^2*t^3)/25-(8*s^3*t^3)/25+x,2-(66*s^2)/25+(16*s^3)/25-(3*t^2)/2+(99*s^2*t^2)/50-(12*s^3*t^2)/25+t^3-(33*s^2*t^3)/25+(8*s^3*t^3)/25+y,-9/10+(27*t)/20-(27*t^2)/40+(3*t^3)/40+z},{2-(66*s^2)/25+(16*s^3)/25-(3*t^2)/2+(99*s^2*t^2)/50-(12*s^3*t^2)/25+t^3-(33*s^2*t^3)/25+(8*s^3*t^3)/25+x,(-84*s)/25+(18*s^2)/25+(16*s^3)/25+(63*s*t^2)/25-(27*s^2*t^2)/50-(12*s^3*t^2)/25-(42*s*t^3)/25+(9*s^2*t^3)/25+(8*s^3*t^3)/25+y,-9/10+(27*t)/20-(27*t^2)/40+(3*t^3)/40+z},{(-84*s)/25+(18*s^2)/25+(16*s^3)/25+(63*s*t^2)/25-(27*s^2*t^2)/50-(12*s^3*t^2)/25-(42*s*t^3)/25+(9*s^2*t^3)/25+(8*s^3*t^3)/25+x,-2+(66*s^2)/25-(16*s^3)/25+(3*t^2)/2-(99*s^2*t^2)/50+(12*s^3*t^2)/25-t^3+(33*s^2*t^3)/25-(8*s^3*t^3)/25+y,-9/10+(27*t)/20-(27*t^2)/40+(3*t^3)/40+z},{8/5-(3*s^2)/10+s^3/5+(21*t)/10+(27*s^2*t)/10-(9*s^3*t)/5-(9*t^2)/10-(9*s^2*t^2)/5+(6*s^3*t^2)/5-t^3/10+(3*s^2*t^3)/10-(s^3*t^3)/5+x,(9*s)/10-(9*s^2)/10+y,-81/40-(27*s^2)/40+(9*s^3)/20+(9*t^3)/40+(27*s^2*t^3)/40-(9*s^3*t^3)/20+z},{3/2+(3*s^2)/10-s^3/5+3*t-(27*s^2*t)/10+(9*s^3*t)/5-(3*t^2)/2+(9*s^2*t^2)/5-(6*s^3*t^2)/5-(3*s^2*t^3)/10+(s^3*t^3)/5+x,(-9*s)/10+(9*s^2)/10+y,-9/4+(27*s^2)/40-(9*s^3)/20+(9*t^3)/20-(27*s^2*t^3)/40+(9*s^3*t^3)/20+z},{27/10+(9*s^2)/10-(3*s^3)/5-(3*t^2)/5-(27*s^2*t^2)/20+(9*s^3*t^2)/10-t^3/10+(3*s^2*t^3)/20-(s^3*t^3)/10+x,(9*s)/10-(9*s^2)/10+y,-9/5+(27*t)/40+(81*s^2*t)/40-(27*s^3*t)/20+(27*t^2)/40-(189*s^2*t^2)/80+(63*s^3*t^2)/40-(9*t^3)/20+(99*s^2*t^3)/80-(33*s^3*t^3)/40+z},{3-(9*s^2)/10+(3*s^3)/5-(21*t^2)/20+(27*s^2*t^2)/20-(9*s^3*t^2)/10-t^3/20-(3*s^2*t^3)/20+(s^3*t^3)/10+x,(-9*s)/10+(9*s^2)/10+y,-9/5+(27*t)/20-(81*s^2*t)/40+(27*s^3*t)/20-(9*t^2)/80+(189*s^2*t^2)/80-(63*s^3*t^2)/40-(3*t^3)/80-(99*s^2*t^3)/80+(33*s^3*t^3)/40+z},{-17/10-(27*t)/10-(9*s^2*t)/2+3*s^3*t+(18*t^2)/5+(81*s^2*t^2)/10-(27*s^3*t^2)/5-(19*t^3)/10-(27*s^2*t^3)/5+(18*s^3*t^3)/5+x,(99*s)/50-(99*s^2)/50-(369*s*t^2)/100+(369*s^2*t^2)/100+(123*s*t^3)/50-(123*s^2*t^3)/50+y,-57/40+(99*s^2)/40-(33*s^3)/20-(81*s^2*t)/40+(27*s^3*t)/20-(81*t^2)/40-(27*s^2*t^2)/10+(9*s^3*t^2)/5+(21*t^3)/20+(9*s^2*t^3)/4-(3*s^3*t^3)/2+z},{-17/10-(21*t)/5+(9*s^2*t)/2-3*s^3*t+(63*t^2)/10-(81*s^2*t^2)/10+(27*s^3*t^2)/5-(37*t^3)/10+(27*s^2*t^3)/5-(18*s^3*t^3)/5+x,(-99*s)/50+(99*s^2)/50+(369*s*t^2)/100-(369*s^2*t^2)/100-(123*s*t^3)/50+(123*s^2*t^3)/50+y,-3/5-(99*s^2)/40+(33*s^3)/20-(27*t)/40+(81*s^2*t)/40-(27*s^3*t)/20-(117*t^2)/40+(27*s^2*t^2)/10-(9*s^3*t^2)/5+(9*t^3)/5-(9*s^2*t^3)/4+(3*s^3*t^3)/2+z},{-27/10-(9*s^2)/5+(6*s^3)/5-(3*t)/10-(9*s^2*t)/8+(3*s^3*t)/4+(27*s^2*t^2)/10-(9*s^3*t^2)/5+t^3/5-(39*s^2*t^3)/40+(13*s^3*t^3)/20+x,(3*s)/4-(3*s^2)/4-(9*s*t^2)/10+(9*s^2*t^2)/10+(3*s*t^3)/5-(3*s^2*t^3)/5+y,-12/5-(9*t)/40-(27*s^2*t)/160+(9*s^3*t)/80+(9*t^2)/40+(27*s^2*t^3)/160-(9*s^3*t^3)/80+z},{-33/10+(9*s^2)/5-(6*s^3)/5-(27*t)/40+(9*s^2*t)/8-(3*s^3*t)/4+(9*t^2)/10-(27*s^2*t^2)/10+(9*s^3*t^2)/5-t^3/8+(39*s^2*t^3)/40-(13*s^3*t^3)/20+x,(-3*s)/4+(3*s^2)/4+(9*s*t^2)/10-(9*s^2*t^2)/10-(3*s*t^3)/5+(3*s^2*t^3)/5+y,-12/5-(9*t)/32+(27*s^2*t)/160-(9*s^3*t)/80+(9*t^2)/40+(9*t^3)/160-(27*s^2*t^3)/160+(9*s^3*t^3)/80+z},{(-12*t)/5+(63*s^2*t)/20-(3*s^3*t)/4+(24*t^2)/5-(63*s^2*t^2)/10+(3*s^3*t^2)/2-(13*t^3)/5+(1707*s^2*t^3)/500-(407*s^3*t^3)/500+x,(81*s*t)/20-(9*s^2*t)/10-(3*s^3*t)/4-(81*s*t^2)/10+(9*s^2*t^2)/5+(3*s^3*t^2)/2+(2193*s*t^3)/500-(243*s^2*t^3)/250-(407*s^3*t^3)/500+y,-63/20+(9*t^2)/10-(9*t^3)/20+z},{(81*s*t)/20-(9*s^2*t)/10-(3*s^3*t)/4-(81*s*t^2)/10+(9*s^2*t^2)/5+(3*s^3*t^2)/2+(2193*s*t^3)/500-(243*s^2*t^3)/250-(407*s^3*t^3)/500+x,(12*t)/5-(63*s^2*t)/20+(3*s^3*t)/4-(24*t^2)/5+(63*s^2*t^2)/10-(3*s^3*t^2)/2+(13*t^3)/5-(1707*s^2*t^3)/500+(407*s^3*t^3)/500+y,-63/20+(9*t^2)/10-(9*t^3)/20+z},{(12*t)/5-(63*s^2*t)/20+(3*s^3*t)/4-(24*t^2)/5+(63*s^2*t^2)/10-(3*s^3*t^2)/2+(13*t^3)/5-(1707*s^2*t^3)/500+(407*s^3*t^3)/500+x,(-81*s*t)/20+(9*s^2*t)/10+(3*s^3*t)/4+(81*s*t^2)/10-(9*s^2*t^2)/5-(3*s^3*t^2)/2-(2193*s*t^3)/500+(243*s^2*t^3)/250+(407*s^3*t^3)/500+y,-63/20+(9*t^2)/10-(9*t^3)/20+z},{(-81*s*t)/20+(9*s^2*t)/10+(3*s^3*t)/4+(81*s*t^2)/10-(9*s^2*t^2)/5-(3*s^3*t^2)/2-(2193*s*t^3)/500+(243*s^2*t^3)/250+(407*s^3*t^3)/500+x,(-12*t)/5+(63*s^2*t)/20-(3*s^3*t)/4+(24*t^2)/5-(63*s^2*t^2)/10+(3*s^3*t^2)/2-(13*t^3)/5+(1707*s^2*t^3)/500-(407*s^3*t^3)/500+y,-63/20+(9*t^2)/10-(9*t^3)/20+z},{-1/5+(33*s^2)/125-(8*s^3)/125-(3*t)/5+(99*s^2*t)/125-(24*s^3*t)/125-(21*t^2)/10+(693*s^2*t^2)/250-(84*s^3*t^2)/125+(8*t^3)/5-(264*s^2*t^3)/125+(64*s^3*t^3)/125+x,(42*s)/125-(9*s^2)/125-(8*s^3)/125+(126*s*t)/125-(27*s^2*t)/125-(24*s^3*t)/125+(441*s*t^2)/125-(189*s^2*t^2)/250-(84*s^3*t^2)/125-(336*s*t^3)/125+(72*s^2*t^3)/125+(64*s^3*t^3)/125+y,-27/10+(9*t)/20-(9*t^2)/20+(3*t^3)/10+z},{(42*s)/125-(9*s^2)/125-(8*s^3)/125+(126*s*t)/125-(27*s^2*t)/125-(24*s^3*t)/125+(441*s*t^2)/125-(189*s^2*t^2)/250-(84*s^3*t^2)/125-(336*s*t^3)/125+(72*s^2*t^3)/125+(64*s^3*t^3)/125+x,1/5-(33*s^2)/125+(8*s^3)/125+(3*t)/5-(99*s^2*t)/125+(24*s^3*t)/125+(21*t^2)/10-(693*s^2*t^2)/250+(84*s^3*t^2)/125-(8*t^3)/5+(264*s^2*t^3)/125-(64*s^3*t^3)/125+y,-27/10+(9*t)/20-(9*t^2)/20+(3*t^3)/10+z},{1/5-(33*s^2)/125+(8*s^3)/125+(3*t)/5-(99*s^2*t)/125+(24*s^3*t)/125+(21*t^2)/10-(693*s^2*t^2)/250+(84*s^3*t^2)/125-(8*t^3)/5+(264*s^2*t^3)/125-(64*s^3*t^3)/125+x,(-42*s)/125+(9*s^2)/125+(8*s^3)/125-(126*s*t)/125+(27*s^2*t)/125+(24*s^3*t)/125-(441*s*t^2)/125+(189*s^2*t^2)/250+(84*s^3*t^2)/125+(336*s*t^3)/125-(72*s^2*t^3)/125-(64*s^3*t^3)/125+y,-27/10+(9*t)/20-(9*t^2)/20+(3*t^3)/10+z},{(-42*s)/125+(9*s^2)/125+(8*s^3)/125-(126*s*t)/125+(27*s^2*t)/125+(24*s^3*t)/125-(441*s*t^2)/125+(189*s^2*t^2)/250+(84*s^3*t^2)/125+(336*s*t^3)/125-(72*s^2*t^3)/125-(64*s^3*t^3)/125+x,-1/5+(33*s^2)/125-(8*s^3)/125-(3*t)/5+(99*s^2*t)/125-(24*s^3*t)/125-(21*t^2)/10+(693*s^2*t^2)/250-(84*s^3*t^2)/125+(8*t^3)/5-(264*s^2*t^3)/125+(64*s^3*t^3)/125+y,-27/10+(9*t)/20-(9*t^2)/20+(3*t^3)/10+z},{(-171*t)/40+(5643*s^2*t)/1000-(171*s^3*t)/125+(81*t^2)/20-(2673*s^2*t^2)/500+(162*s^3*t^2)/125-(51*t^3)/40+(1683*s^2*t^3)/1000-(51*s^3*t^3)/125+x,(-3591*s*t)/500+(1539*s^2*t)/1000+(171*s^3*t)/125+(1701*s*t^2)/250-(729*s^2*t^2)/500-(162*s^3*t^2)/125-(1071*s*t^3)/500+(459*s^2*t^3)/1000+(51*s^3*t^3)/125+y,(-9*t^2)/40+(3*t^3)/40+z},{(3591*s*t)/500-(1539*s^2*t)/1000-(171*s^3*t)/125-(1701*s*t^2)/250+(729*s^2*t^2)/500+(162*s^3*t^2)/125+(1071*s*t^3)/500-(459*s^2*t^3)/1000-(51*s^3*t^3)/125+x,(-171*t)/40+(5643*s^2*t)/1000-(171*s^3*t)/125+(81*t^2)/20-(2673*s^2*t^2)/500+(162*s^3*t^2)/125-(51*t^3)/40+(1683*s^2*t^3)/1000-(51*s^3*t^3)/125+y,(-9*t^2)/40+(3*t^3)/40+z},{(171*t)/40-(5643*s^2*t)/1000+(171*s^3*t)/125-(81*t^2)/20+(2673*s^2*t^2)/500-(162*s^3*t^2)/125+(51*t^3)/40-(1683*s^2*t^3)/1000+(51*s^3*t^3)/125+x,(3591*s*t)/500-(1539*s^2*t)/1000-(171*s^3*t)/125-(1701*s*t^2)/250+(729*s^2*t^2)/500+(162*s^3*t^2)/125+(1071*s*t^3)/500-(459*s^2*t^3)/1000-(51*s^3*t^3)/125+y,(-9*t^2)/40+(3*t^3)/40+z},{(-3591*s*t)/500+(1539*s^2*t)/1000+(171*s^3*t)/125+(1701*s*t^2)/250-(729*s^2*t^2)/500-(162*s^3*t^2)/125-(1071*s*t^3)/500+(459*s^2*t^3)/1000+(51*s^3*t^3)/125+x,(171*t)/40-(5643*s^2*t)/1000+(171*s^3*t)/125-(81*t^2)/20+(2673*s^2*t^2)/500-(162*s^3*t^2)/125+(51*t^3)/40-(1683*s^2*t^3)/1000+(51*s^3*t^3)/125+y,(-9*t^2)/40+(3*t^3)/40+z}};

patches3D[[1]]

patches3D[[17]]

Illustrating the teapot

My thanks to WRI colleague Yu-sung Chang for doing the bulk of work in making this manipulation of the teapot. Warning: This is a simulated teapot. Do not try this on an actual teapot.

potPoints={{1.4,0,2.4},{1.4,-0.784,2.4},{0.784,-1.4,2.4},{0,-1.4,2.4},{1.3375,0,2.53125},{1.3375,-0.749,2.53125},{0.749,-1.3375,2.53125},{0,-1.3375,2.53125},{1.4375,0,2.53125},{1.4375,-0.805,2.53125},{0.805,-1.4375,2.53125},{0,-1.4375,2.53125},{1.5,0,2.4},{1.5,-0.84,2.4},{0.84,-1.5,2.4},{0,-1.5,2.4},{-0.784,-1.4,2.4},{-1.4,-0.784,2.4},{-1.4,0,2.4},{-0.749,-1.3375,2.53125},{-1.3375,-0.749,2.53125},{-1.3375,0,2.53125},{-0.805,-1.4375,2.53125},{-1.4375,-0.805,2.53125},{-1.4375,0,2.53125},{-0.84,-1.5,2.4},{-1.5,-0.84,2.4},{-1.5,0,2.4},{-1.4,0.784,2.4},{-0.784,1.4,2.4},{0,1.4,2.4},{-1.3375,0.749,2.53125},{-0.749,1.3375,2.53125},{0,1.3375,2.53125},{-1.4375,0.805,2.53125},{-0.805,1.4375,2.53125},{0,1.4375,2.53125},{-1.5,0.84,2.4},{-0.84,1.5,2.4},{0,1.5,2.4},{0.784,1.4,2.4},{1.4,0.784,2.4},{0.749,1.3375,2.53125},{1.3375,0.749,2.53125},{0.805,1.4375,2.53125},{1.4375,0.805,2.53125},{0.84,1.5,2.4},{1.5,0.84,2.4},{1.75,0,1.875},{1.75,-0.98,1.875},{0.98,-1.75,1.875},{0,-1.75,1.875},{2.,0,1.35},{2.,-1.12,1.35},{1.12,-2.,1.35},{0,-2.,1.35},{2.,0,0.9},{2.,-1.12,0.9},{1.12,-2.,0.9},{0,-2.,0.9},{-0.98,-1.75,1.875},{-1.75,-0.98,1.875},{-1.75,0,1.875},{-1.12,-2.,1.35},{-2.,-1.12,1.35},{-2.,0,1.35},{-1.12,-2.,0.9},{-2.,-1.12,0.9},{-2.,0,0.9},{-1.75,0.98,1.875},{-0.98,1.75,1.875},{0,1.75,1.875},{-2.,1.12,1.35},{-1.12,2.,1.35},{0,2.,1.35},{-2.,1.12,0.9},{-1.12,2.,0.9},{0,2.,0.9},{0.98,1.75,1.875},{1.75,0.98,1.875},{1.12,2.,1.35},{2.,1.12,1.35},{1.12,2.,0.9},{2.,1.12,0.9},{2.,0,0.45},{2.,-1.12,0.45},{1.12,-2.,0.45},{0,-2.,0.45},{1.5,0,0.225},{1.5,-0.84,0.225},{0.84,-1.5,0.225},{0,-1.5,0.225},{1.5,0,0.15},{1.5,-0.84,0.15},{0.84,-1.5,0.15},{0,-1.5,0.15},{-1.12,-2.,0.45},{-2.,-1.12,0.45},{-2.,0,0.45},{-0.84,-1.5,0.225},{-1.5,-0.84,0.225},{-1.5,0,0.225},{-0.84,-1.5,0.15},{-1.5,-0.84,0.15},{-1.5,0,0.15},{-2.,1.12,0.45},{-1.12,2.,0.45},{0,2.,0.45},{-1.5,0.84,0.225},{-0.84,1.5,0.225},{0,1.5,0.225},{-1.5,0.84,0.15},{-0.84,1.5,0.15},{0,1.5,0.15},{1.12,2.,0.45},{2.,1.12,0.45},{0.84,1.5,0.225},{1.5,0.84,0.225},{0.84,1.5,0.15},{1.5,0.84,0.15},{-1.6,0,2.025},{-1.6,-0.3,2.025},{-1.5,-0.3,2.25},{-1.5,0,2.25},{-2.3,0,2.025},{-2.3,-0.3,2.025},{-2.5,-0.3,2.25},{-2.5,0,2.25},{-2.7,0,2.025},{-2.7,-0.3,2.025},{-3.,-0.3,2.25},{-3.,0,2.25},{-2.7,0,1.8},{-2.7,-0.3,1.8},{-3.,-0.3,1.8},{-3.,0,1.8},{-1.5,0.3,2.25},{-1.6,0.3,2.025},{-2.5,0.3,2.25},{-2.3,0.3,2.025},{-3.,0.3,2.25},{-2.7,0.3,2.025},{-3.,0.3,1.8},{-2.7,0.3,1.8},{-2.7,0,1.575},{-2.7,-0.3,1.575},{-3.,-0.3,1.35},{-3.,0,1.35},{-2.5,0,1.125},{-2.5,-0.3,1.125},{-2.65,-0.3,0.9375},{-2.65,0,0.9375},{-2.,-0.3,0.9},{-1.9,-0.3,0.6},{-1.9,0,0.6},{-3.,0.3,1.35},{-2.7,0.3,1.575},{-2.65,0.3,0.9375},{-2.5,0.3,1.125},{-1.9,0.3,0.6},{-2.,0.3,0.9},{1.7,0,1.425},{1.7,-0.66,1.425},{1.7,-0.66,0.6},{1.7,0,0.6},{2.6,0,1.425},{2.6,-0.66,1.425},{3.1,-0.66,0.825},{3.1,0,0.825},{2.3,0,2.1},{2.3,-0.25,2.1},{2.4,-0.25,2.025},{2.4,0,2.025},{2.7,0,2.4},{2.7,-0.25,2.4},{3.3,-0.25,2.4},{3.3,0,2.4},{1.7,0.66,0.6},{1.7,0.66,1.425},{3.1,0.66,0.825},{2.6,0.66,1.425},{2.4,0.25,2.025},{2.3,0.25,2.1},{3.3,0.25,2.4},{2.7,0.25,2.4},{2.8,0,2.475},{2.8,-0.25,2.475},{3.525,-0.25,2.49375},{3.525,0,2.49375},{2.9,0,2.475},{2.9,-0.15,2.475},{3.45,-0.15,2.5125},{3.45,0,2.5125},{2.8,0,2.4},{2.8,-0.15,2.4},{3.2,-0.15,2.4},{3.2,0,2.4},{3.525,0.25,2.49375},{2.8,0.25,2.475},{3.45,0.15,2.5125},{2.9,0.15,2.475},{3.2,0.15,2.4},{2.8,0.15,2.4},{0,0,3.15},{0,-0.002,3.15},{0.002,0,3.15},{0.8,0,3.15},{0.8,-0.45,3.15},{0.45,-0.8,3.15},{0,-0.8,3.15},{0,0,2.85},{0.2,0,2.7},{0.2,-0.112,2.7},{0.112,-0.2,2.7},{0,-0.2,2.7},{-0.002,0,3.15},{-0.45,-0.8,3.15},{-0.8,-0.45,3.15},{-0.8,0,3.15},{-0.112,-0.2,2.7},{-0.2,-0.112,2.7},{-0.2,0,2.7},{0,0.002,3.15},{-0.8,0.45,3.15},{-0.45,0.8,3.15},{0,0.8,3.15},{-0.2,0.112,2.7},{-0.112,0.2,2.7},{0,0.2,2.7},{0.45,0.8,3.15},{0.8,0.45,3.15},{0.112,0.2,2.7},{0.2,0.112,2.7},{0.4,0,2.55},{0.4,-0.224,2.55},{0.224,-0.4,2.55},{0,-0.4,2.55},{1.3,0,2.55},{1.3,-0.728,2.55},{0.728,-1.3,2.55},{0,-1.3,2.55},{1.3,0,2.4},{1.3,-0.728,2.4},{0.728,-1.3,2.4},{0,-1.3,2.4},{-0.224,-0.4,2.55},{-0.4,-0.224,2.55},{-0.4,0,2.55},{-0.728,-1.3,2.55},{-1.3,-0.728,2.55},{-1.3,0,2.55},{-0.728,-1.3,2.4},{-1.3,-0.728,2.4},{-1.3,0,2.4},{-0.4,0.224,2.55},{-0.224,0.4,2.55},{0,0.4,2.55},{-1.3,0.728,2.55},{-0.728,1.3,2.55},{0,1.3,2.55},{-1.3,0.728,2.4},{-0.728,1.3,2.4},{0,1.3,2.4},{0.224,0.4,2.55},{0.4,0.224,2.55},{0.728,1.3,2.55},{1.3,0.728,2.55},{0.728,1.3,2.4},{1.3,0.728,2.4},{0,0,0},{1.5,0,0.15},{1.5,0.84,0.15},{0.84,1.5,0.15},{0,1.5,0.15},{1.5,0,0.075},{1.5,0.84,0.075},{0.84,1.5,0.075},{0,1.5,0.075},{1.425,0,0},{1.425,0.798,0},{0.798,1.425,0},{0,1.425,0},{-0.84,1.5,0.15},{-1.5,0.84,0.15},{-1.5,0,0.15},{-0.84,1.5,0.075},{-1.5,0.84,0.075},{-1.5,0,0.075},{-0.798,1.425,0},{-1.425,0.798,0},{-1.425,0,0},{-1.5,-0.84,0.15},{-0.84,-1.5,0.15},{0,-1.5,0.15},{-1.5,-0.84,0.075},{-0.84,-1.5,0.075},{0,-1.5,0.075},{-1.425,-0.798,0},{-0.798,-1.425,0},{0,-1.425,0},{0.84,-1.5,0.15},{1.5,-0.84,0.15},{0.84,-1.5,0.075},{1.5,-0.84,0.075},{0.798,-1.425,0},{1.425,-0.798,0}};

pp2=With[{avg=Mean[potPoints]},Map[(#-avg)&,potPoints]];

potPatches={{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16},{4,17,18,19,8,20,21,22,12,23,24,25,16,26,27,28},{19,29,30,31,22,32,33,34,25,35,36,37,28,38,39,40},{31,41,42,1,34,43,44,5,37,45,46,9,40,47,48,13},{13,14,15,16,49,50,51,52,53,54,55,56,57,58,59,60},{16,26,27,28,52,61,62,63,56,64,65,66,60,67,68,69},{28,38,39,40,63,70,71,72,66,73,74,75,69,76,77,78},{40,47,48,13,72,79,80,49,75,81,82,53,78,83,84,57},{57,58,59,60,85,86,87,88,89,90,91,92,93,94,95,96},{60,67,68,69,88,97,98,99,92,100,101,102,96,103,104,105},{69,76,77,78,99,106,107,108,102,109,110,111,105,112,113,114},{78,83,84,57,108,115,116,85,111,117,118,89,114,119,120,93},{121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136},{124,137,138,121,128,139,140,125,132,141,142,129,136,143,144,133},{133,134,135,136,145,146,147,148,149,150,151,152,69,153,154,155},{136,143,144,133,148,156,157,145,152,158,159,149,155,160,161,69},{162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177},{165,178,179,162,169,180,181,166,173,182,183,170,177,184,185,174},{174,175,176,177,186,187,188,189,190,191,192,193,194,195,196,197},{177,184,185,174,189,198,199,186,193,200,201,190,197,202,203,194},{204,204,204,204,207,208,209,210,211,211,211,211,212,213,214,215},{204,204,204,204,210,217,218,219,211,211,211,211,215,220,221,222},{204,204,204,204,219,224,225,226,211,211,211,211,222,227,228,229},{204,204,204,204,226,230,231,207,211,211,211,211,229,232,233,212},{212,213,214,215,234,235,236,237,238,239,240,241,242,243,244,245},{215,220,221,222,237,246,247,248,241,249,250,251,245,252,253,254},{222,227,228,229,248,255,256,257,251,258,259,260,254,261,262,263},{229,232,233,212,257,264,265,234,260,266,267,238,263,268,269,242},{270,270,270,270,279,280,281,282,275,276,277,278,271,272,273,274},{270,270,270,270,282,289,290,291,278,286,287,288,274,283,284,285},{270,270,270,270,291,298,299,300,288,295,296,297,285,292,293,294},{270,270,270,270,300,305,306,279,297,303,304,275,294,301,302,271}};

Bernstein[i_Integer,n_Integer,t_Symbol]:=Binomial[n,i] t^i (1-t)^(n-i)
Bernstein[n_Integer,t_Symbol]:=Map[Bernstein[#,n,t]&,Range[0,n]]

BezierSurface[pts_,u_,v_,m_]:=Module[{dim=Dimensions[pts],avg=Mean[Flatten[pts,1]]},Bernstein[dim[[1]]-1,u].Map[Bernstein[dim[[2]]-1,v].#&,Map[(#)&,pts]]+m*avg]

Manipulate[ParametricPlot3D[Evaluate[BezierSurface[#,s,t,m]&/@(Transpose[Partition[Part[pp2,#],4]]&/@potPatches)],{s,0,1},{t,0,1},Axes->None,Boxed->False,Evaluated->False,PlotPoints->ControlActive[{3,3},{8,8}],PlotStyle->Directive[Specularity[White,50],Yellow],MaxRecursion->0,Mesh->None],{{m,.2,"Explode"},0,1}]

Example 4 (Lichtblau 2006)

This is an algebraic but nonrational parametrization (it uses trigonometric polynomials).

mainvars={x,y};
elimvars={λ,(Subscript[c, a]),(Subscript[c, b]),(Subscript[s, a]),(Subscript[s, b])};polys={-1+(Subscript[c, a]^2)+(Subscript[s, a]^2),-1+(Subscript[c, b]^2)+(Subscript[s, b]^2),x-(Subscript[c, a])-(Subscript[c, b])-(Subscript[c, a]) (Subscript[c, b])-(Subscript[s, a]) (Subscript[s, b]),y-(Subscript[c, a]^5)-(Subscript[c, b]^5)-(Subscript[c, a]^5) (Subscript[c, b]^5)+10 (Subscript[c, a]^3) (Subscript[s, a]^2)+10 (Subscript[c, a]^3) (Subscript[c, b]^5) (Subscript[s, a]^2)-5 (Subscript[c, a]) (Subscript[s, a]^4)-5 (Subscript[c, a]) (Subscript[c, b]^5) (Subscript[s, a]^4)-25 (Subscript[c, a]^4) (Subscript[c, b]^4) (Subscript[s, a]) (Subscript[s, b])+50 (Subscript[c, a]^2) (Subscript[c, b]^4) (Subscript[s, a]^3) (Subscript[s, b])-5 (Subscript[c, b]^4) (Subscript[s, a]^5) (Subscript[s, b])+10 (Subscript[c, b]^3) (Subscript[s, b]^2)+10 (Subscript[c, a]^5) (Subscript[c, b]^3) (Subscript[s, b]^2)-100 (Subscript[c, a]^3) (Subscript[c, b]^3) (Subscript[s, a]^2) (Subscript[s, b]^2)+50 (Subscript[c, a]) (Subscript[c, b]^3) (Subscript[s, a]^4) (Subscript[s, b]^2)+50 (Subscript[c, a]^4) (Subscript[c, b]^2) (Subscript[s, a]) (Subscript[s, b]^3)-100 (Subscript[c, a]^2) (Subscript[c, b]^2) (Subscript[s, a]^3) (Subscript[s, b]^3)+10 (Subscript[c, b]^2) (Subscript[s, a]^5) (Subscript[s, b]^3)-5 (Subscript[c, b]) (Subscript[s, b]^4)-5 (Subscript[c, a]^5) (Subscript[c, b]) (Subscript[s, b]^4)+50 (Subscript[c, a]^3) (Subscript[c, b]) (Subscript[s, a]^2) (Subscript[s, b]^4)-25 (Subscript[c, a]) (Subscript[c, b]) (Subscript[s, a]^4) (Subscript[s, b]^4)-5 (Subscript[c, a]^4) (Subscript[s, a]) (Subscript[s, b]^5)+10 (Subscript[c, a]^2) (Subscript[s, a]^3) (Subscript[s, b]^5)-(Subscript[s, a]^5) (Subscript[s, b]^5),(Subscript[s, a])-25 λ (Subscript[c, a]^4) (Subscript[s, a])+(Subscript[c, b]) (Subscript[s, a])-25 λ (Subscript[c, a]^4) (Subscript[c, b]^5) (Subscript[s, a])+50 λ (Subscript[c, a]^2) (Subscript[s, a]^3)+50 λ (Subscript[c, a]^2) (Subscript[c, b]^5) (Subscript[s, a]^3)-5 λ (Subscript[s, a]^5)-5 λ (Subscript[c, b]^5) (Subscript[s, a]^5)-(Subscript[c, a]) (Subscript[s, b])+25 λ (Subscript[c, a]^5) (Subscript[c, b]^4) (Subscript[s, b])-250 λ (Subscript[c, a]^3) (Subscript[c, b]^4) (Subscript[s, a]^2) (Subscript[s, b])+125 λ (Subscript[c, a]) (Subscript[c, b]^4) (Subscript[s, a]^4) (Subscript[s, b])+250 λ (Subscript[c, a]^4) (Subscript[c, b]^3) (Subscript[s, a]) (Subscript[s, b]^2)-500 λ (Subscript[c, a]^2) (Subscript[c, b]^3) (Subscript[s, a]^3) (Subscript[s, b]^2)+50 λ (Subscript[c, b]^3) (Subscript[s, a]^5) (Subscript[s, b]^2)-50 λ (Subscript[c, a]^5) (Subscript[c, b]^2) (Subscript[s, b]^3)+500 λ (Subscript[c, a]^3) (Subscript[c, b]^2) (Subscript[s, a]^2) (Subscript[s, b]^3)-250 λ (Subscript[c, a]) (Subscript[c, b]^2) (Subscript[s, a]^4) (Subscript[s, b]^3)-125 λ (Subscript[c, a]^4) (Subscript[c, b]) (Subscript[s, a]) (Subscript[s, b]^4)+250 λ (Subscript[c, a]^2) (Subscript[c, b]) (Subscript[s, a]^3) (Subscript[s, b]^4)-25 λ (Subscript[c, b]) (Subscript[s, a]^5) (Subscript[s, b]^4)+5 λ (Subscript[c, a]^5) (Subscript[s, b]^5)-50 λ (Subscript[c, a]^3) (Subscript[s, a]^2) (Subscript[s, b]^5)+25 λ (Subscript[c, a]) (Subscript[s, a]^4) (Subscript[s, b]^5),-(Subscript[c, b]) (Subscript[s, a])+25 λ (Subscript[c, a]^4) (Subscript[c, b]^5) (Subscript[s, a])-50 λ (Subscript[c, a]^2) (Subscript[c, b]^5) (Subscript[s, a]^3)+5 λ (Subscript[c, b]^5) (Subscript[s, a]^5)+(Subscript[s, b])+(Subscript[c, a]) (Subscript[s, b])-25 λ (Subscript[c, b]^4) (Subscript[s, b])-25 λ (Subscript[c, a]^5) (Subscript[c, b]^4) (Subscript[s, b])+250 λ (Subscript[c, a]^3) (Subscript[c, b]^4) (Subscript[s, a]^2) (Subscript[s, b])-125 λ (Subscript[c, a]) (Subscript[c, b]^4) (Subscript[s, a]^4) (Subscript[s, b])-250 λ (Subscript[c, a]^4) (Subscript[c, b]^3) (Subscript[s, a]) (Subscript[s, b]^2)+500 λ (Subscript[c, a]^2) (Subscript[c, b]^3) (Subscript[s, a]^3) (Subscript[s, b]^2)-50 λ (Subscript[c, b]^3) (Subscript[s, a]^5) (Subscript[s, b]^2)+50 λ (Subscript[c, b]^2) (Subscript[s, b]^3)+50 λ (Subscript[c, a]^5) (Subscript[c, b]^2) (Subscript[s, b]^3)-500 λ (Subscript[c, a]^3) (Subscript[c, b]^2) (Subscript[s, a]^2) (Subscript[s, b]^3)+250 λ (Subscript[c, a]) (Subscript[c, b]^2) (Subscript[s, a]^4) (Subscript[s, b]^3)+125 λ (Subscript[c, a]^4) (Subscript[c, b]) (Subscript[s, a]) (Subscript[s, b]^4)-250 λ (Subscript[c, a]^2) (Subscript[c, b]) (Subscript[s, a]^3) (Subscript[s, b]^4)+25 λ (Subscript[c, b]) (Subscript[s, a]^5) (Subscript[s, b]^4)-5 λ (Subscript[s, b]^5)-5 λ (Subscript[c, a]^5) (Subscript[s, b]^5)+50 λ (Subscript[c, a]^3) (Subscript[s, a]^2) (Subscript[s, b]^5)-25 λ (Subscript[c, a]) (Subscript[s, a]^4) (Subscript[s, b]^5)};

It arises in the problem of finding envelope curves for the trignonometric parametrizations depicted below.

trigpoly[n_,a_,b_]=Cos[n a]+Cos[n b]+Cos[n (a-b)];
pair[a_,b_,n_]={trigpoly[1,a,b],trigpoly[n,a,b]};
pointlist[m_]:=Table[2*Pi*RandomReal[{},2],{2^m}]
points8k=pointlist[13];
data[{point__},n_]:=pair[point,n];
datalist8k[n_]:=Map[data[#,n]&,points8k];

Manipulate[ListPlot[datalist8k[n],PlotRange->All],{{n,3,"n"},{2,3,4,5,6,7,8}},SaveDefinitions->True]

Example 5 (Tran 2004)

Here we have a rational parametrization wherein the fourth polynomial is the common denominator. We will implicitize via an old method from (Kalkbrenner 1990).

rational6={-168*t-336*s*t+378*s^2*t+36*t^2-1728*s*t^2+1719*s^2*t^2+32*t^3+2064*s*t^3-2072*s^2*t^3,-100+25*s^2-600*s*t+600*s^2*t+132*t^2-144*s*t^2+111*s^2*t^2-32*t^3+744*s*t^3-736*s^2*t^3,6-6*s+s^2+18*s*t-18*s^2*t+36*s*t^2-36*s^2*t^2-54*s*t^3+54*s^2*t^3,1+6*s*t-6*s^2*t+12*s*t^2-12*s^2*t^2-18*s*t^3+18*s^2*t^3};
denom=Last[rational6];
polys6=Join[vars-Most[rational6]*drecip,{denom*drecip-1}]

Results

Example 1

Timing[gbHoffmann=GroebnerBasis[HoffmannPolys,vars,elims,MonomialOrder->EliminationOrder,Sort->True,Method->{"GroebnerWalk","EarlyEliminate"->True}];]

{21.6854, Null}{21.6854, Null}

The early elimination check is vital for this example; without it the computation takes around 500 seconds.

Example 2, modulo a reasonably large prime

We can do this.

Timing[gbmod=GroebnerBasis[TranPolys,vars,elims,Sort->True,MonomialOrder->EliminationOrder,Modulus->Prime[1000001],Method->"GroebnerWalk"];]

But the early elimination check helps tremendously for this example.

Timing[gbmod=GroebnerBasis[TranPolys,vars,elims,Sort->True,MonomialOrder->EliminationOrder,Modulus->Prime[1000001],Method->{"GroebnerWalk","EarlyEliminate"->True}];]

Of note: the input is already a Gröbner basis with respect to the monomial order given (with variable order s>t>x>y>z) by the weight matrix below.

{{0, 0, 1, 1, 1}, {1, 1, 0, 0, 0}, {0, 0, 0, 0, -1}, {0, 0, 0, -1, 0}, {0, -1, 0, 0, 0}}}

If we use this as initial order for the Gröbner walk then it instead takes longer than I was willing to wait. So we do not use this order.

Example 2 in characteristic zero

How long to do this in characteristic zero?

Timing[gb=GroebnerBasis[TranPolys,vars,elims,Sort->True,MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}];]

Ouch. The gremlin in the machinery is the coefficient swell. We will now see how to keep this at bay.

Example 2 approximated, then rationalized

We can obtain a result in reasonable time using approximation, then rationalization, then verification. For this particular example the total time, including verification, is around 12 minutes. We start by computing a numerical Gröbner basis to high precision.

Timing[gbapprox=GroebnerBasis[TranPolys,vars,elims,MonomialOrder->EliminationOrder,CoefficientDomain->InexactNumbers[600],Method->{"GroebnerWalk","EarlyEliminate"->True}];]

{383.9, Null}{383.9, Null}

implicit=Rationalize[gbapprox[[1]]];
Precision[implicit]

InfinityInfinity

rules=First[Solve[Tranpolys==0,vars]];
Timing[Expand[implicit/.rules]==0]

{340.485, True}{340.485, True}

See (Shirayanagi 1996), (Lichtblau 2000) for discussion regarding computation of numerical Gröbner bases.

Another possibility for example 2

Use a prime modulus, then Hensel lift or else use skeleton of result to solve a linear system.

I tried the latter and got a huge sparse system I was unable to solve without blowing up memory.

Might do better when I get some sparse linear algebra for matrices with arbitrary size integers.

Example 3, modulo a reasonably large prime

We will attempt the implicitization of all patches, using a modulus to keep coefficient swell from occurring.

Table[Timing[implicitmod[j]=GroebnerBasis[patches3D[[j]],vars,elims,MonomialOrder->EliminationOrder,Modulus->Prime[1000001],Method->{"GroebnerWalk","EarlyEliminate"->True}];][[1]],{j,1,Length[patches3D]}]

We do them all, many in "real time". We will show timings for the nonmodular case for all but the four spout parametrizations (patches 17-20), where we will resort to different methods to get the implicitizations.

Example 3 in characteristic zero

Table[Timing[implicit[j]=GroebnerBasis[patches3D[[j]],vars,elims,MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}];][[1]],{j,1,16}]

Table[Timing[implicit[j]=GroebnerBasis[patches3D[[j]],vars,elims,MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}];][[1]],{j,21,Length[patches3D]}]

It is interesting to note that without doing a perturbation of the initial weight vector, several of these take substantially longer (more than a factor of 10). They seem to get caught in large intermediate Gröbner basis computations because ideal initials are not always small.

Plots of a few patches may give some idea of why the implicitization difficulty varies so widely.

patch 1

patch 11

patch 19

One might expect from comparing the modular and characteristic zero timings that the remaining four cases could be handled in reasonable time (an hour each, say). This turns out not to be the case, at least with option settings I have tried. We instead resort to numerical computations followed by rationalization.

Timing[implicit17approx=GroebnerBasis[patches3D[[17]],vars,elims,MonomialOrder->EliminationOrder,CoefficientDomain->InexactNumbers[700],Method->{"GroebnerWalk","EarlyEliminate"->True}];]

{445.356, Null}{445.356, Null}

We rationalize, clear denominators (not necessary, but nice to have result with integer coefficients), plug in original relations and check that it expands to zero.

implicit17attempt=Rationalize[implicit17approx];
imp2=Apply[List,implicit17attempt[[1]]];
imp2=imp2/.{x->1,y->1,z->1};
den=Denominator[imp2];
den=Apply[LCM,den];
implicit17=Expand[den*implicit17attempt[[1]]];
rules17=First[Solve[patches3D[[17]]==0,vars]];
Timing[Expand[implicit17/.rules17]==0]

{280.038, True}{280.038, True}

Example 3 cont'd

We just show timings for computing the remaining three patch implicitizations to 700 digits. That turns out to suffice for rationalizing and recovering the correct exact polynomials.

Timing[implicit18approx=GroebnerBasis[patches3D[[18]],vars,elims,MonomialOrder->EliminationOrder,CoefficientDomain->InexactNumbers[700],Method->{"GroebnerWalk","EarlyEliminate"->True}];]

{465.629, Null}{465.629, Null}

Timing[implicit19approx=GroebnerBasis[patches3D[[19]],vars,elims,MonomialOrder->EliminationOrder,CoefficientDomain->InexactNumbers[700],Method->{"GroebnerWalk","EarlyEliminate"->True}];]

{1216.11, Null}{1216.11, Null}

Timing[implicit20approx=GroebnerBasis[patches3D[[20]],vars,elims,MonomialOrder->EliminationOrder,CoefficientDomain->InexactNumbers[700],Method->{"GroebnerWalk","EarlyEliminate"->True}];]

{1291.09, Null}{1291.09, Null}

Example 4

This next example requires use of a perturbation of the target weight vector. Moreover an a priori computation a la (Tran 2000) is out of the question because the number of variables would cause the exponent vectors to have too many digits to allow the arithmetic to be feasible. We use the method of (Amrheim et al) though the simulated perturbation of (Fukuda et al) would also work here.

Timing[implicitpoly=Factor[First[GroebnerBasis[polys,mainvars,elimvars,Sort->True,MonomialOrder->EliminationOrder,Method->{"GroebnerWalk","EarlyEliminate"->True}]]]]

Note that the result factors nontrivially, that is, it is not prime. So it could not be the implicitization of a rational parametrization.

"The envelope, please."

(I always wanted to say that).

Example 5

Timing[gbrat6=GroebnerBasis[polys6,vars,Prepend[elims,drecip],MonomialOrder->EliminationOrder,Sort->True,Method->{"GroebnerWalk","EarlyEliminate"->True}];]

By using approximate arithmetic of 300 digits one can obtain the same result in about 24 seconds. Not quite "real time" but getting reasonably close...

Summary

What we can do

Implicitization of "difficult" polynomial, rational, and algebraic parametrizations from the literature.

Many other elimination and related problems are also amenable to effective handling via Gröbner walk computations.

Necessary tactics

Perturbation of initial weight vector.

Other than when reaching boundary of final cone, this will help to avoid "large initials" problem.

Always interreduce intermediate bases.

Failure to do so will cause huge blowup due to "false cones" problem.

It will probably also cause crop failures. At least, it has not been proven otherwise.

When indicated, perturbation of final weight vector. Viable possibilities are

Heuristic

Simulated

Early finish when we discover we have attained an elimination ideal

Rarely hurts, never by much, and sometimes helps substantially.

Within variable classes (elimination, not elimination), sort using some sensible heuristic

Usually improves on input order unless you know a "good" variable ordering a priori.

In some cases we must use approximate arithmetic, then rationalize to recover the correct exact result

For some purposes an approximate result, particularly one to high precision, might suffice.

Open bottlenecks

Computation of the initial basis can be costly.

For implicitization one can use initial orders that make the input, in effect, a basis already.

This is not true for nonrational parametrizations.

It generally does not help, and in some cases makes the overall walk much slower.

Provisional conclusion: Use degree reverse lexicographic as initial ordering unless there is a compelling reason not to do so.

Interreduction of intermediate bases is generally the most time consuming step.

Per earlier remarks, it is at present unavoidable.

"Lifting" from initials basis to full basis seems to be second most costly.

Method of (Fukuda et al 2007) might improve on pedestrian method from earlier literature?

The sheer number of cones, large polynomials that complicate the lifting and/or interreduction step, and coefficient swell still pose serious bottlenecks in some problems.

But, as examples above indicate, we are showing progress using the Gröbner walk for implicitization and other computations.

References

B. Amrhein, O. Gloor, and W. Küchlin. On the walk.Theoretical Computer Science 187:179-202. 1997.

S. Collart, M. Kalkbrenner, and D. Mall. Converting bases with the Gröbner walk. Journal of Symbolic Computation 24:465-469. 1997.

D. Cox, J. Little, and D. O'Shea. Ideals, Varieties, and Algorithms: An Introduction to Computational ALgebraic Geometry and Commutative Algebra (third edition). Undergraduate Texts in Mathematics. Springer-Verlag, 2007.

K. Fukuda, A. N. Jensen, N Lauritzen, and R. Thomas. The generic Gröbner walk. Journal of Symbolic Computation 42:198-213. 2007.

M. Kalkbrenner. Implicitization of rational parametric curves and surfaces. AAECC-8, Lecture Notes in Computer Science 508:249-259, S. Sakata, ed. Springer, 1990.

R. Lewis. Private communication, 2007.

D. Lichtblau. Solving finite algebraic systems using numeric Gröbner bases and eigenvalues. Proceedings of SCI 2000 10:555-560. 2000.

D. Lichtblau. Computing curves bounding trigonometric planar maps: symbolic and hybrid methods. Proceedings of the Workshop on Automated and Deductive Geometry (ADG 2004). Lecture Notes in Artificial Intelligence 3763:70-91, H. Hong and D. Wang, eds. Springer-Verlag, 2006.

K. Shirayanagi. Floating point Gröbner bases. Mathematics and Computers in Simulation 42:509-528. 1996.

Q-N Tran. A fast algorithm for Gröbner basis conversion and its applications. Journal of Symbolic Computation 30:451-467. 2000.

Q-N Tran. Efficient Groebner walk conversion for implicitization of geometric objects. Computer Aided Geometric Design 21:837-857. 2004.

Q-N Tran. A new class of term orders for elimination. Journal of Symbolic Computation 42:533-548. 2007.