Re: commutative polynomial ideals in Mathematica and Maple



Roman Pearce wrote:
> Maple 10 added some new functionality for working with commutative
> polynomial ideals. This is an ongoing project and obviously we have a
> long ways to go, but for now I was hoping to have a friendly comparison
> with Mathematica and perhaps some other systems. Since it wouldn't be
> fair to compare top level commands, light coding is allowed. Below are
> the examples I posted to comp.soft-sys.math.maple, along with times for
> a 1.4 GHz Athlon. I have some other machines I can use as well. If
> you have any of your own example problems feel free to add them to this
> list. Everything is appreciated :)
>
> Roman
>
>
> 1) compute a lexicographic Groebner basis for Katsura-5 with x > y > z
> > t > u > v :
> {x*y+y*z+2*z*t+2*t*u+2*u*v-u, 2*x^2+2*y^2+2*z^2+2*t^2+2*u^2+ v^2-v,
> 2*x*z+2*y*t+2*z*u+u^2+2*v*t-t, 2*x*t+2*u*y+2*t*u+2*
> z*v-z, t^2+2*x*v+2*y*v+2*z*v-y, 2*x+2*y+2*z+2*t+2*u+v-1}
> Maple: 6.5 sec
>
> 2) compute a lexicographic Groebner basis for the Butcher system with w
> > v > u > t > z > y > x :
> {2*z*u+2*y*v+2*t*w-2*w^2-w-1,
> 3*z*u^2+3*y*v^2-3*t*w^2+3*w^3+ 3*w^2-t+4*w,
> 6*x*z*v-6*t*w^2+6*w^3-3*t*w+6* w^2-t+4*w, 4*z*
> u^3+4*y*v^3+4*t*w^3-4*w^4-6*w^ 3+4*t*w-10*w^2-w-1,
> 8*x*z*u*v+8*t*w^3-8*w^4+4*t*w^ 2-12*w^3+4*t*w-14*w^2-3*w-1, 12*x*z*v^
> 2+12*t*w^3-12*w^4+12*t*w^2-18* w^3+8*t*w-14*w^2-w-1,
> -24*t*w^3+24*w^4-24*t*w^2+36*w ^3-8*t*w+26*w^2+7*w+1}
> Maple: 9.5 sec
>
> 3) test if this ideal (weispfenning94) is maximal in Q[x,y,z]. ie: is
> the quotient ring Q[x,y,z]/I a field ?
> {x*y^4+y*z^4-2*x^2*y-3, y^4+x*y^2*z+x^2-2*x*y+y^2+z^2,
> -x^3*y^2+x*y*z^3+y^4+x*y^2*z-2 *x*y}
> Maple: 4.4 sec
>
> 4) Test if this ideal is radical (kinema):
> {z1^2+z2^2+z3^2-12*z1-68, z4^2+z5^2+z6^2-12*z5-68,
> z7^2+z8^2+z9^2-24*z8-12*z9+100 , z1*z4+z2*z5+z3*z6-6*z1-6*z5-52 , z1*
> z7+z2*z8+z3*z9-6*z1-12*z8-6*z9 +64, z1+z3-2*z4+z5-z6+2*z7-2*z8+8,
> 2*z2+2*z3-z4-z5-2*z6-z7-z9+18, z1+z2+2*z3+2*z4+2*z6-2
> *z7+z8-z9-38, z4*z7+z5*z8+z6*z9-6*z5-12*z8-6 *z9+32}
> Maple: 55 sec
>
> 5) Compute a primary decomposition for schwarz7:
> {x4^2+2*x3*x5+2*x6*x2+2*x1*x7+ x7, 2*x3*x4+2*x2*x5+2*x1*x6+x7^2+x 6,
> x3^2+2*x2*x4+2*x1*x5+2*x6*x7+x 5, 2*x2*x3+2*x1*x4+x6
> ^2+2*x7*x5+x4, x2^2+2*x1*x3+2*x5*x6+2*x7*x4+x 3,
> 2*x1*x2+x5^2+2*x6*x4+2*x7*x3+x 2, x1^2+2*x4*x5+2*x6*x3+2*x7*x2+x 1}
> Maple: 14 sec
>
> 6) Compute a primary decomposition of schwarz9:
> {2*x4*x5+2*x6*x3+2*x7*x2+2*x8* x1+x9^2+x8,
> x5^2+2*x6*x4+2*x7*x3+2*x8*x2+2 *x1*x9+x9,
> x4^2+2*x3*x5+2*x6*x2+2*x1*x7+2 *x8*
> x9+x7, 2*x3*x4+2*x2*x5+2*x1*x6+x8^2+2 *x7*x9+x6,
> x3^2+2*x2*x4+2*x1*x5+2*x7*x8+2 *x6*x9+x5, 2*x2*x3+2*x1*x4+x7^2+2*x8*x6+
>
> 2*x5*x9+x4, x2^2+2*x1*x3+2*x6*x7+2*x8*x5+2 *x4*x9+x3,
> 2*x1*x2+x6^2+2*x7*x5+2*x8*x4+2 *x3*x9+x2, x1^2+2*x5*x6+2*x7*x4+2*
> x8*x3+2*x2*x9+x1}
> Maple: 6 minutes
>
> 7) test whether katsura6 is radical
> {2*z*t+2*u*y+2*x*v+2*w*y-v, x^2+2*y^2+2*z^2+2*t^2+2*u^2+2* v^2+2*w^2-x,
>
> x+2*y+2*z+2*t+2*u+2*v+2*w-1, 2*y*z+2*x*t+2*u*y+
> 2*z*v+2*t*w-t, 2*x*y+2*y*z+2*z*t+2*t*u+2*u*v+ 2*v*w-y,
> y^2+2*x*z+2*y*t+2*z*u+2*v*t+2* w*u-z, z^2+2*y*t+2*u*x+2*y*v+2*w*z-u}
> Maple: 2 minutes, 6 seconds.


First I would like to thank Roman Pearce for reposting these to
sci.math.symbolic, where people such as myself can show how to do such
problems in Mathematica or other programs.

Second, I would like to point out that his presentation of results
shows substantial good work. I'm not sure if this was obvious to
everyone, but it is clear enough to me.

Third in the order of business is to make some mention of what one
might hope to accomplish with these examples in Mathematica. Timings
are on a 3 GHz pentium 4 processor running under Linux. I'm using a
development version of Mathematica but I think all these will work
pretty much the same with the current street version.

(1)
polys1 = {x*y+y*z+2*z*t+2*t*u+2*u*v-u,
2*x^2+2*y^2+2*z^2+2*t^2+2*u^2+v^2-v,
2*x*z+2*y*t+2*z*u+u^2+2*v*t-t, 2*x*t+2*u*y+2*t*u+2*z*v-z,
t^2+2*x*v+2*y*v+2*z*v-y, 2*x+2*y+2*z+2*t+2*u+v-1};
vars1 = {x,y,z,t,u,v};

I think this will take a very long time if done via Buchberger's
algorithm, so I'll specify to use the Groebner walk as method.

In[4]:= Timing[gb1 = GroebnerBasis[polys1, vars1,
Method->"GroebnerWalk"];]
Out[4]= {1.03484 Second, Null}

For purposes of checking, I'll show the first element. As coefficients
are several hundred digits I'll numericize for further brevity.

In[5]:= InputForm[N[First[gb1]]]
Out[5]//InputForm=
2.5222979497716460426446941626391`15.954589770191005*^562 +
1.908352507544669496226853414803113`15.954589770191005*^553*u -
1.384149770805177959851895757420161`15.954589770191005*^564*v +
3.6019533202604346464555471`15.954589770191005*^565*v^2 -
5.92788627635073525155352762`15.954589770191005*^566*v^3 +
6.940018931690148135058458423`15.954589770191005*^567*v^4 -
6.1624458473675889986719332145`15.954589770191005*^568*v^5 +
4.32003940602465869825805929657`15.954589770191005*^569*v^6 -
2.456768255554960314331979076083`15.954589770191005*^570*v^7 +
1.1555431360706851654193403495708`15.954589770191005*^571*v^8 -
4.5595277857633470697226226917442`15.954589770191005*^571*v^9 +
1.52529830660595045486123698678602`15.954589770191005*^572*v^10 -
4.36015516836721746128269543542333`15.954589770191005*^572*v^11 +
1.071127646747865863266676290053613`15.954589770191005*^573*v^12 -
2.270252897150901624111358338042538`15.954589770191005*^573*v^13 +
4.161137953154269822621625243844157`15.954589770191005*^573*v^14 -
6.60104777330259696382011134936249`15.954589770191005*^573*v^15 +
9.056756973469580669152325655887127`15.954589770191005*^573*v^16 -
1.072177208506528387903001672519329`15.954589770191005*^574*v^17 +
1.0904896896057841944056401271675306`15.954589770191005*^574*v^18 -
9.464254506482286169862515115784789`15.954589770191005*^573*v^19 +
6.93772884459763884920911232719225`15.954589770191005*^573*v^20 -
4.229577137103989977820623991692447`15.954589770191005*^573*v^21 +
2.092700064229363597075341734062981`15.954589770191005*^573*v^22 -
8.0526201027071615750471460084194`15.954589770191005*^572*v^23 +
2.20151869815276075799705176996865`15.954589770191005*^572*v^24 -
3.1492790235782082491287672760722`15.954589770191005*^571*v^25 -
3.644114656332327607048066919726`15.954589770191005*^570*v^26 +
2.9557727650010426719940027889`15.954589770191005*^570*v^27 -
5.59653280653038756139324776093`15.954589770191005*^569*v^28 -
1.646205595091591492065754689`15.954589770191005*^567*v^29 +
1.5893617212514762591095667423`15.954589770191005*^568*v^30 -
1.761144344042372288854511546`15.954589770191005*^567*v^31


(2)
polys2 = {2*z*u+2*y*v+2*t*w-2*w^2-w-1,
3*z*u^2+3*y*v^2-3*t*w^2+3*w^3+3*w^2-t+4*w,
6*x*z*v-6*t*w^2+6*w^3-3*t*w+6*w^2-t+4*w,
4*z*u^3+4*y*v^3+4*t*w^3-4*w^4-6*w^3+4*t*w-10*w^2-w-1,
8*x*z*u*v+8*t*w^3-8*w^4+4*t*w^2-12*w^3+4*t*w-14*w^2-3*w-1,
12*x*z*v^2+12*t*w^3-12*w^4+12*t*w^2-18*w^3+8*t*w-14*w^2-w-1,
-24*t*w^3+24*w^4-24*t*w^2+36*w^3-8*t*w+26*w^2+7*w+1};
vars2 = {w,v,u,t,z,y,x};

Same idea: use the Groebner walk as method.

In[10]:= InputForm[Timing[
gb2 = GroebnerBasis[polys2, vars2, Method->"GroebnerWalk"]]]

Out[10]//InputForm=
{0.5379189999999999*Second, {-1200636*x*z - 1200636*t*x*z +
206829*x^2*z +
206829*t*x^2*z + 1855828*x^3*z + 1855828*t*x^3*z + 3962844*x*y*z +
3962844*t*x*y*z, 81*x*z + 81*t*x*z - 243*x^2*z - 243*t*x^2*z +
129*x^3*z +
129*t*x^3*z + 356*x^4*z + 356*t*x^4*z, -202284*x*z - 202284*t*x*z -
542115*x^2*z - 542115*t*x^2*z - 614812*x^3*z - 614812*t*x^3*z +
1320948*x*z^2 + 1320948*t*x*z^2, 4077 + 8154*t + 4077*t^2 -
51378*x*z -
51378*t*x*z + 66128*x^2*z + 66128*t*x^2*z + 119616*x^3*z +
119616*t*x^3*z,
587088*y + 587088*t*y - 880632*x^2*y - 880632*t*x^2*y - 783081*x*z -
783081*t*x*z + 2015748*x^2*z + 2015748*t*x^2*z + 2634772*x^3*z +
2634772*t*x^3*z - 880632*u*x*z^2 + 1320948*u*x^3*z^2,
-213993576*y - 213993576*t*y + 213993576*x*y + 213993576*t*x*y -
855974304*y^2 - 855974304*t*y^2 + 512876502*x*z + 512876502*t*x*z -
1483609725*x^2*z - 1483609725*t*x^2*z - 1583657812*x^3*z -
1583657812*t*x^3*z + 320990364*u*x*z^2 - 320990364*u*x^2*z^2 +
1283961456*u*x*y*z^2, -71331192*y - 71331192*t*y + 71331192*x*y +
71331192*t*x*y + 678530817*x*z + 678530817*t*x*z - 1022191362*x^2*z
-
1022191362*t*x^2*z - 1273988008*x^3*z - 1273988008*t*x^3*z -
570649536*y*z - 570649536*t*y*z + 106996788*u*x*z^2 -
106996788*u*x^2*z^2 + 855974304*u*x*z^3, 4077*u + 4077*t*u -
47118*x*z -
47118*t*x*z + 36604*x^2*z + 36604*t*x^2*z + 75472*x^3*z +
75472*t*x^3*z,
293544*x + 293544*t*x - 2348352*y - 2348352*t*y - 3522528*x*y -
3522528*t*x*y + 8064009*x*z + 8064009*t*x*z + 220158*u^2*x*z -
13940271*x^2*z - 13940271*t*x^2*z - 20078044*x^3*z -
20078044*t*x^3*z +
3522528*u*x*z^2 + 5283792*u*x^2*z^2, 1467720*y + 1467720*t*y -
3376296*x*z - 3376296*t*x*z + 4815081*x^2*z + 4815081*t*x^2*z +
6093652*x^3*z + 6093652*t*x^3*z + 440316*u^2*y*z + 440316*u^2*z^2 -
1320948*u*x*z^2, 110079 + 110079*t - 1174176*y - 1174176*t*y +
36693*v*y -
1761264*x*y - 1761264*t*x*y + 36693*u*z + 3539322*x*z +
3539322*t*x*z -
6161292*x^2*z - 6161292*t*x^2*z - 8263472*x^3*z - 8263472*t*x^3*z +
1761264*u*x*z^2 + 2641896*u*x^2*z^2, 8154 + 8154*t - 97848*y -
97848*t*y -
97848*x*y - 97848*t*x*y + 249462*x*z + 249462*t*x*z + 12231*v*x*z -
471626*x^2*z - 471626*t*x^2*z - 591672*x^3*z - 591672*t*x^3*z +
146772*u*x*z^2 + 146772*u*x^2*z^2, 36693*v + 36693*t*v - 374850*x*z
-
374850*t*x*z + 542088*x^2*z + 542088*t*x^2*z + 712000*x^3*z +
712000*t*x^3*z, -48924 - 48924*t + 293544*y + 293544*t*y +
293544*x*y +
293544*t*x*y - 36693*u^2*z + 36693*u*v*z - 623142*x*z - 623142*t*x*z
+
1315704*x^2*z + 1315704*t*x^2*z + 1530800*x^3*z + 1530800*t*x^3*z -
440316*u*x*z^2 - 440316*u*x^2*z^2, -85617 - 97848*t + 12231*w +
1174176*y + 1174176*t*y + 1174176*x*y + 1174176*t*x*y - 3132756*x*z
-
3132756*t*x*z + 5666556*x^2*z + 5666556*t*x^2*z + 7260976*x^3*z +
7260976*t*x^3*z - 1761264*u*x*z^2 - 1761264*u*x^2*z^2}}


(3)
For maximality testing I chose to proceed as follows. First get a
lexicographic Groebner basis. Next show the ideal is finite (a
necessary condition), then show the univariate polynomial does not
factor (another necessary condition) and all remaining polynomials are
linear. This last is a sufficient condition. If the ideal is radical
but not in general position we would actually need to do a random (that
is, "generic") linear change in the "smallest" variable. I do not do
this here, but it is simple to code.

isMaximal[polys_] := Module[
{gb, n, vars=Variables[polys]},
n = Length[vars];
gb = GroebnerBasis[GroebnerBasis[polys, vars,
Method->"GroebnerWalk"], vars];
If [Length[gb]!=n ||
Variables[gb[[1]]]=!={Last[vars]} ||
Length[FactorList[gb[[1]]]]=!=2,
Return[False]];
For [j=2, j<=n, j++,
If [Complement[Variables[gb[[j]]],{vars[[n-j+1]],Last[vars]}]=!={} ||
Exponent[gb[[j]],vars[[n-j+1]]]=!=1, Return[False]]];
True
]

polys3 = {x*y^4+y*z^4-2*x^2*y-3, y^4+x*y^2*z+x^2-2*x*y+y^2+z^2,
-x^3*y^2+x*y*z^3+y^4+x*y^2*z-2*x*y};
vars3 = Variables[polys3];

In[14]:= Timing[isMaximal[polys3]]
Out[14]= {0.6609 Second, True}


(5), (6)
I do not have code to do a primary decomposition. Offhand I do not know
how easy or difficult that might be to implement, but I suspect it is
not particularly easy. Needless to say, if anyone has code to do this
and wants to share it, feel free to let me know (or post in this
thread).


(4)
Showing in general that an ideal is radical is not trivial. In the
special case where the variety is zero dimensional it can be done as
follows. One gets minimal univariate polynomials in each variable. (One
exists for each variable by assumption of a finite solution set.) The
ideal is radical iff there is no nontrivial factorization of any of
these univariates.

I chose to implement this in a way that is relatively easy to code but
is not so good in regard to speed. First I produce a Groebner basis by
degree reverse lexicographic order. Check the finiteness condition
(easy). Then for each variable we convert to lex order with that
variable last (hence produce a univariate in that polynomial). If it
factors nontrivially, return False immediately. If we get through all
variables, return True.

The potential improvement is to use an FGLM type of conversion to find
the univariate polynomials via linear algebra. This actually predates
both the "FL" and "GM" work, and goes back to Buchberger, I think in
the early-to-mid 80's. If I get time I may code this variant for
comparison.

The asture reader will notice a GroebnerBasis[GroebnerBasis[...]] sort
of construction in the part that actually gets the univariates. This is
to work around a silly deficiency wherein the Groebner walk method
fails to sort its result by term order. I intend to repair this. In any
case the impact on timing is not large. There are other potential
savings e.g. by not recomputing graded reverse lex basis, but again,
the savings will be small at least in the examples we discuss below.

zeroDimensionalRadicalQ[polys_] := Catch[Module[
{gb, vars=Variables[polys], len, leadvecs, pureterms, dtl},
len = Length[vars];
gb = GroebnerBasis[polys, vars,
MonomialOrder->DegreeReverseLexicographic];
dtl = First[Internal`DistributedTermsList[gb, vars,
MonomialOrder->DegreeReverseLexicographic]];
leadvecs = Map[#[[1,1]]&, dtl];
pureterms = Select[leadvecs, Count[#,a_ /;a=!=0]===1&];
If [Length[pureterms]<len, Throw[False]];
checkUnivariatePolynomials[gb,vars]
]]

checkUnivariatePolynomials[basis_,vars_] := Module[
{unipolys, len=Length[vars], newgb, poly},
Do [
var = vars[[j]];
newvars = RotateRight[vars,j];
newgb = GroebnerBasis[basis, newvars,
Method->{"GroebnerWalk"}];
poly = First[GroebnerBasis[newgb, newvars]];
fax = Rest[FactorList[poly]];
If [Length[fax]>1 || fax[[1,2]]!=1, Throw[False]];
, {j,len}];
True
]

polys4 = {z1^2+z2^2+z3^2-12*z1-68, z4^2+z5^2+z6^2-12*z5-68,
z7^2+z8^2+z9^2-24*z8-12*z9+100, z1*z4+z2*z5+z3*z6-6*z1-6*z5-52,
z1*z7+z2*z8+z3*z9-6*z1-12*z8-6*z9+64, z1+z3-2*z4+z5-z6+2*z7-2*z8+8,
2*z2+2*z3-z4-z5-2*z6-z7-z9+18, z1+z2+2*z3+2*z4+2*z6-2*z7+z8-z9-38,
z4*z7+z5*z8+z6*z9-6*z5-12*z8-6*z9+32};

In[10]:= Timing[zeroDimensionalRadicalQ[polys4]]
Out[10]= {18.28 Second, False}


(7)
polys7 = {2*z*t+2*u*y+2*x*v+2*w*y-v,
x^2+2*y^2+2*z^2+2*t^2+2*u^2+2*v^2+2*w^2-x,
x+2*y+2*z+2*t+2*u+2*v+2*w-1, 2*y*z+2*x*t+2*u*y+2*z*v+2*t*w-t,
2*x*y+2*y*z+2*z*t+2*t*u+2*u*v+2*v*w-y,
y^2+2*x*z+2*y*t+2*z*u+2*v*t+2*w*u-z, z^2+2*y*t+2*u*x+2*y*v+2*w*z-u};

In[12]:= Timing[zeroDimensionalRadicalQ[polys7]]
Out[12]= {20.19 Second, False}


Daniel Lichtblau
Wolfram Research

.

Relevant Pages

  • Re: d|a[n]
    ... and easier to use than those complicated representing polynomials. ... ideals I_0,I_1,... ... a regular basis that all ideals _m belong to the main class. ... If one forms the product of all prime ideals of the same degree, ...
    (sci.math)
  • Groebner bases over the integers
    ... I'm currently working with computations in the ring Zof ... Algorithm for Finitely Generated Ideals in Rings", Buchberger, 1983. ... by adding all polynomials that can be obtained by ... a Groebner basis of the ideal I', where I' is the smallest ideal ...
    (sci.math.symbolic)
  • Re: Prime ideals of C[x]?
    ... The polynomials of degree at most m-1 ... | But what are ideals of this ring? ... |are prime ideals in this ring? ...
    (sci.math)
  • Re: hilberts nullstellensatz
    ... The argument I used is much simpler than the Nullstellensatz. ... Let Vthe variety of zeros of all polynomials in I ... Since radis the intersection of all prime ideals in R containing I, ... Let M be a maximal ideal in the finitely generated k-algebra ...
    (sci.math)
  • Re: High Precision Approximate GCD
    ... EpsilonGCD and QuasiGCD are probably the first packages available. ... Maple 9) is much more ... UvGCD is reported in "The approximate GCD of inexact polynomials. ... > The package has commands like EpsilonGCD and QuasiGCD and others. ...
    (sci.math.symbolic)