From: Daniel Lichtblau <danl@wolfram.com>
Newsgroups: sci.math.symbolic
Subject: Re: Q:implicitation of curves
Date: Mon, 21 Jul 1997 12:27:21 -0500

Andreas Unterkircher wrote:
> 
> Say I have plane curve c=(f_1(t),f_2(t)) and I want to compute the
> implicit form p(x,y) of this curve. This could be done with Groebner
> Basis techniques. Unfortunately this approach often is not useful due to
> the complexity of this algorithm. My question is: are there any other
> techniques to compute the implicit form of c (maybe some special cases).
> I am especially interested in the case where c is a Bezier curve or some
> spline. Any hints are welcome !
> 
> Andreas
> 
> --
> Andreas Unterkircher             Institut fuer Umformtechnik
> ETH Zentrum / CLA, Tannenstr. 3,CH-8092 Zuerich, Switzerland
> Phone: ++43 1 632 6621                       www.ifu.ethz.ch


The reply below is based on e-mail back and forth between AU and myself.

A first attempt at an answer is that it depends on the degree of the
data and type of coefficients. For low degree, Groebner techniques
should suffice. Alternatively one can form a resultant. Below I show a
simple example with Mathematica, running on a Pentium-Pro under Linux.

In[1]:= f[t] = Array[a,4,0] . t^Range[0,3];

In[2]:= g[t] = Array[b,4,0] . t^Range[0,3];

In[3]:= Timing[imp1 = GroebnerBasis[{x-f[t], y-g[t]}, {x,y}, t,
        MonomialOrder->EliminationOrder][[1]];]
Out[3]= {0.16 Second, Null}

In[4]:= Timing[imp2 = Resultant[x-f[t], y-g[t], t];] (* gives -imp1 *)
Out[4]= {0.14 Second, Null}


For much higher degree input, it is not likely that these will be
viable using symbolic coefficients. In our development version of
Mathematica one can find a Groebner basis with specified bignum
precision coefficients, so for inexact data this might be a reasonable
approach.

In[6]:= SeedRandom[1111];

In[7]:= f[t] = Table[Random[Real, 10, 150], {10}] . t^Range[0,9];

In[8]:= g[t] = Table[Random[Real, 10, 150], {10}] . t^Range[0,9];

In[9]:= Timing[imp = GroebnerBasis[{x-f[t], y-g[t]}, {x,y}, t,
        MonomialOrder->EliminationOrder,
        CoefficientDomain->InexactNumbers[150]][[1]];]
Out[9]= {3.75 Second, Null}


Alternative methods for implicitization might include interpolating the 
polynomial using some a priori degree bounds. One can check the text
"Effective Polynomial Computation" by Zippel for more information on
this. I think this reference may also provide information about the
resultant-based approach.


Based on this preliminary response, and an invitation to send an example
of interest, AU sent me the following with a note to the effect that
Mathematica 3.0 could not handle it.

polys = {x - (110*t^2 - 495*t^3 + 1320*t^4 - 2772*t^5 + 5082*t^6 -
7590*t^7 +
	  8085*t^8 - 5555*t^9 + 2189*t^10 - 374*t^11),
	y - (22*t - 110*t^2 + 330*t^3 - 1848*t^5 + 3696*t^6 - 3300*t^7 +
1650*t^8 -
	  550*t^9 + 88*t^10 + 22*t^11)};

The development version of Mathematica will handle this provided one
starts at 250 or so digits of precision (200 did not suffice). While I
have not tried an exact computation, AU says it runs over 3 hours with
no result. So clearly the inexact method is of interest here. Results
below are abbreviated for purposes of, well, brevity.


In[1]:= InputForm[Timing[poly250 = First[GroebnerBasis[
	polys, {x,y}, t, MonomialOrder->EliminationOrder,
	  CoefficientDomain->InexactNumbers[250]]]]]
    
Out[1]//InputForm= 
  {20.689999999999998*Second,
1.97552413410743382899217474906499999999999\
99999999999580685`46.7148*^30*x^3 - 
     6.7273185232570481684919210047500000000000000000000000198397`47.46\
78*^29*x^4 +
1.753245020195242250921029593500000000000000000000000003732\
38`47.7702*^29*x^5 -
1.954692034344520501241951172400000000000000000000000\
000215795`48.9375*^28*x^6 +
1.13993498927601163681082731799999999999999999\
99999999999999887325`52.0201*^27*x^7 - ....

The result contains decimals that appear to end generally in long
sequences of 0's or 9's. After playing a bit with Rationalize, it seemed
that the common denominator is 1024. In fact we do get a polynomial with
integer coefficients below.

In[7]:= newpoly250 = Rationalize[Expand[poly250*1024]] // InputForm

Out[7]//InputForm= 
  2022936713326012240887986943042560*x^3 - 
   688877416781521732453572710886400*x^4 + 
   179532290067992806494313430374400*x^5 - 
   20016046431687889932717580005376*x^6 + ...

This is in fact the correct exact implicit polynomial in x and y for the
curve defined parametrically by polys. First evidence that this is so
comes when one repeats the computation at 500 digits of precision,
rationalizes, and gets the same polynomial. More importantly, one can
prove that this is the implicit polynomial using the following facts.
(i) The polynomial is irreducible (Factor[] will show this)
(ii) Polynomial reduction with respect to a Groebner basis will reduce
it to zero. As it happens, the input is already a Groebner basis with
respect to the lexicographic term order in the variables {x,y,t} (this
is the case for ideals defined by parametric equations). A simple
reduction with respect to polys, in this term order, gives zero.


In[14]:= PolynomialReduce[newpoly250,
        {x - (110*t^2 - 495*t^3 + 1320*t^4 - 2772*t^5 + 5082*t^6 -
7590*t^7 +
          8085*t^8 - 5555*t^9 + 2189*t^10 - 374*t^11),
        y - (22*t - 110*t^2 + 330*t^3 - 1848*t^5 + 3696*t^6 - 3300*t^7 +
1650*t^8 -
          550*t^9 + 88*t^10 + 22*t^11)}, {x,y,t}] [[2]]

Out[14]= 0


If the coefficients of the inexact polynomial did not lend themselves to
easy rationalization, one might still use its "skeleton," that is, the
power-products in x and y that show up with non-zero coefficients. A bit
of linear algebra, described next, will allow one to determine the exact
coefficients. This is pretty much the method of basis conversion, and
the idea goes back to Bruno Buchberger's '85 work "Groebner bases: an
algorithmic method in polynomial ideal theory." One multiplies each
power-product by a distinct new indeterminate to form a new polynomial
poly. One then calls PolynomialReduce on the whole thing (with respect
to any Groebner bases; polys will suffice), and gets a "reduced"
polynomial in the original variables, with coefficients in terms of the
new indeterminates. One sets this to zero, which means each coefficient
must be zero; this gives a system of equations in the indeterminates.
Solving this system and back-substituting the solutions into poly gives
the exact polynomial in x and y.

These last paragraphs were a bit of a digression on how one can use
inexact methods to arrive at an exact solution. In practice the
approximate solutions might suffice for one's work, and indeed the input
data might well be approximate in the setting described (Bezier or
spline curves).


Daniel Lichtblau
Wolfram Research
danl@wolfram.com