The PolynomialIdeals Package
This worksheet demonstrates some features of the PolynomialIdeals package and their use in solving large polynomial systems.
Warning: Some examples may require a powerful computer.
restart;withPolynomialIdeals:
Introduction
We begin by examining a cyclic-5 system. A good first step for working with any system is to compute the Hilbert dimension. If the dimension is zero, then the system has a finite number of solutions.
cyclic5:=x⁢y+y⁢z+z⁢t+t⁢u+u⁢x,x⁢y⁢z+y⁢z⁢t+z⁢t⁢u+t⁢u⁢x+u⁢x⁢y,x+y+z+t+u,x⁢y⁢z⁢t+y⁢z⁢t⁢u+z⁢t⁢u⁢x+t⁢u⁢x⁢y+u⁢x⁢y⁢z,x⁢y⁢z⁢t⁢u−1
cyclic5:=x⁢y+y⁢z+z⁢t+t⁢u+u⁢x,x⁢y⁢z⁢t⁢u−1,x+y+z+t+u,x⁢y⁢z⁢t+y⁢z⁢t⁢u+z⁢t⁢u⁢x+t⁢u⁢x⁢y+u⁢x⁢y⁢z,x⁢y⁢z+y⁢z⁢t+z⁢t⁢u+t⁢u⁢x+u⁢x⁢y
HilbertDimension⁡cyclic5
0
An ideal is radical if it has no repeated solutions. It is primary if it does not split into two or more components. An ideal is maximal if it is not contained in any larger (proper) ideal. We can test all of the following using the commands below.
IsRadical⁡cyclic5
true
IsPrimary⁡cyclic5
false
IsMaximal⁡cyclic5
To do these computations, PolynomialIdeals uses Groebner bases. We can see which ones were computed using the IdealInfo[KnownGroebnerBases] command. Once computed, Groebner bases are remembered by the system and used again whenever possible.
IdealInfoKnownGroebnerBasescyclic5
tdeg⁡u,y,t,x,z,plex⁡u,y,t,x,z
All of the major algorithms make use of lexicographic Groebner bases. Below we will factor the polynomial system into its prime components. Because the ideal is zero-dimensional, each component corresponds to a finite set of points in affine space.
P:=PrimeDecomposition⁡cyclic5:
nops⁡P
20
A lexicographic basis is known for every component in the decomposition. The Simplify command replaces the ideal generators by the Groebner basis.
Simplify⁡P1
−z2+x,−z3+t,1+z+z2+z3+y,z4+z3+z2+z+1,−1+u
The form below is particularly useful.
map⁡lprint,sort⁡map⁡GroebnerBasis,P,plexz,y,x,t,u,length
[1+3*u+u^2, 3+u+t, -1+x, -1+y, z-1] [-1+u, 1+3*t+t^2, -1+x, -1+y, t+3+z] [-1+u, -1+t, 1+3*x+x^2, 3+x+y, z-1] [-1+u, -1+t, -1+x, 1+3*y+y^2, 3+z+y] [1+3*u+u^2, -1+t, 3+x+u, -1+y, z-1] [1+u+u^2+u^3+u^4, -u+t, -1+u-u^2+x, 1+2*u+u^2+y, z-u] [1+u+u^2+u^3+u^4, -1+u-u^2+t, -u+x, -u+y, 1+2*u+u^2+z] [1+u+u^2+u^3+u^4, -u+t, -u+x, 1+2*u+u^2+y, -1+u-u^2+z] [1+u+u^2+u^3+u^4, -u+t, 1+2*u+u^2+x, -1+u-u^2+y, z-u] [1+u+u^2+u^3+u^4, 1+2*u+u^2+t, -u+x, -u+y, -1+u-u^2+z] [1+u+u^2+u^3+u^4, -u+t, -u+x, -1+u-u^2+y, 1+2*u+u^2+z] [-1+u, 1+t^2+t+t^3+t^4, 1+t^2+t+t^3+x, -t^3+y, -t^2+z] [1+u+u^2+u^3+u^4, -1+t, -u^2+x, -u^3+y, 1+u+u^2+u^3+z] [1+u+u^2+u^3+u^4, -u^3+t, 1+u+u^2+u^3+x, -u^2+y, z-1] [1+u+u^2+u^3+u^4, -u^2+t, -1+x, u^3+1+u+u^2+y, -u^3+z] [1+u+u^2+u^3+u^4, 1+u+u^2+u^3+t, -u^3+x, -1+y, -u^2+z] [1-4*u+6*u^2+u^3+u^4, -3-64*u-15*u^2-9*u^3+11*t, 1+25*u+5*u^2+3*u^3+11*x, 1+25*u+5*u^2+3*u^3+11*y, 1+25*u+5*u^2+3*u^3+11*z] [1-4*u+6*u^2+u^3+u^4, 1+25*u+5*u^2+3*u^3+11*t, -3-64*u-15*u^2-9*u^3+11*x, 1+25*u+5*u^2+3*u^3+11*y, 1+25*u+5*u^2+3*u^3+11*z] [1+u+6*u^2-4*u^3+u^4, 6-67*u+39*u^2-9*u^3+11*t, -2+26*u-13*u^2+3*u^3+11*x, -2+26*u-13*u^2+3*u^3+11*y, -2+26*u-13*u^2+3*u^3+11*z] [1+u+6*u^2-4*u^3+u^4, -2+26*u-13*u^2+3*u^3+11*t, 6-67*u+39*u^2-9*u^3+11*x, -2+26*u-13*u^2+3*u^3+11*y, -2+26*u-13*u^2+3*u^3+11*z]
Separating Solutions
In the example above, several solutions - corresponding to distinct points on the affine variety - share the same value of u. For example, there are 4 different solutions with u=1. For some applications, it may be desirable to separate all of the solutions, so that they differ in at least one single coordinate.
One way to do this is to introduce a new variable equal to a "random" linear combination of the original variables. In almost all cases, the solutions will have distinct values in the new variable.
f:=W−x+y−2⁢z+2⁢t−3⁢u
cyc5W:=Add⁡cyclic5,f:
IdealInfoKnownGroebnerBasescyc5W
One difficulty with this method is that as soon as the new polynomial is added to the ideal all of our Groebner bases are lost. However, note that if G is a Groebner basis for the original ideal, then [G, f] is a Groebner basis for the new ideal under an appropriate product order.
T≔IdealInfoDefaultMonomialOrdercyclic5,'tdeg';G≔GroebnerBasiscyclic5,T:
T:=tdeg⁡u,y,t,x,z
Gnew:=op⁡G,f:Tnew:=prod⁡plex⁡W,T
Tnew:=prod⁡plex⁡W,tdeg⁡u,y,t,x,z
cyc5W:=PolynomialIdeal⁡Add⁡cyclic5,f,known_groebner_bases=Tnew=Gnew:
Now that a Groebner basis is known, we can compute the univariate polynomial with respect to the new variable using only linear algebra.
Wpoly:=factor⁡UnivariatePolynomial⁡W,cyc5W
Wpoly:=W2−6⁢W−11⁢W4+13⁢W3+69⁢W2+77⁢W+121⁢W4−7⁢W3+69⁢W2−143⁢W+121⁢W2−W−11⁢W4+8⁢W3+34⁢W2+77⁢W+121⁢W4−7⁢W3+34⁢W2−88⁢W+121⁢W2+14⁢W+44⁢W4−12⁢W3+64⁢W2−88⁢W+1936⁢W4−2⁢W3+64⁢W2−528⁢W+1936⁢W2−6⁢W+4⁢W4+8⁢W3+24⁢W2−8⁢W+16⁢W4−2⁢W3+24⁢W2+32⁢W+16⁢W2−W−31⁢W4−12⁢W3+94⁢W2−403⁢W+961⁢W4+13⁢W3+94⁢W2+372⁢W+961⁢W4+13⁢W3+54⁢W2+92⁢W+181⁢W4−12⁢W3+74⁢W2−183⁢W+181⁢W4−7⁢W3−6⁢W2+112⁢W+181⁢W4+8⁢W3+54⁢W2+97⁢W+181⁢W4−2⁢W3+39⁢W2−118⁢W+181
nops⁡Wpoly
degree⁡Wpoly
70
Each factor above corresponds to a prime component of the original ideal. There are 70 distinct solutions of this system in C[x,y,z,t,u]. To verify that the solutions are completely separated, we check that a Groebner basis with W << {x, y, z, t, u} is linear in {x, y, z, t, u}. Again, since a Groebner basis is known this will be computed with only linear algebra.
map⁡degree,GroebnerBasiscyc5W,lexdeg⁡x,y,z,t,u,W,x,y,z,t,u
0,1,1,1,1,1
Infolevel[GroebnerBasis]
The Groebner basis command can output various amounts of information about its progress, as specified by infolevel[GroebnerBasis]. This is especially useful on large systems.
trinks:=−9⁢w+15⁢p⁢t+20⁢z⁢s,99⁢w−11⁢s⁢b+3⁢b2,35⁢p+40⁢z+25⁢t−27⁢s,45⁢p+35⁢s−165⁢b−36,w⁢p+2⁢z⁢t−11⁢b3,15⁢w+25⁢p⁢s+30⁢z−18⁢t−165⁢b2
trinks:=35⁢p+40⁢z+25⁢t−27⁢s,15⁢w+25⁢p⁢s+30⁢z−18⁢t−165⁢b2,99⁢w−11⁢s⁢b+3⁢b2,−9⁢w+15⁢p⁢t+20⁢z⁢s,w⁢p+2⁢z⁢t−11⁢b3,45⁢p+35⁢s−165⁢b−36
Infolevel 2 times the various algorithms and indicates their execution.
infolevelGroebnerBasis≔2:
GroebnerBasistrinks,lexdeg⁡z,t,w,p,s,b:
-> FGb total time: 0.031 sec -> FGLM algorithm total monomials : 16 total time (sec) 0.15e-1
Infolevel 3 outputs more comprehensive statistics about algorithm performance.
infolevelGroebnerBasis:=3:
GroebnerBasistrinks,plex⁡z,t,w,p,s,b:
-> FGLM algorithm rational coefficients total monomials : 16 normalforms (sec) 0.15e-1 linearsolve (sec) 0. total time (sec) 0.15e-1
-----------------------------
Infolevel 4 prints intermediate information in a compact form. This is useful for monitoring large computations.
infolevelGroebnerBasis:=4:
GroebnerBasistrinks,plex⁡b,p,s,t,w,z:
-> FGLM algorithm rational coefficients current monomial 1 current monomial z current monomial z^2 current monomial z^3 current monomial z^4 current monomial z^5 current monomial z^6 current monomial z^7 current monomial z^8
current monomial z^9 current monomial z^10 current monomial w current monomial t current monomial s current monomial p current monomial b total monomials : 16 normalforms (sec) 0.94e-1 linearsolve (sec) 0.16e-1 total time (sec) .110 -----------------------------
Infolevel 5 prints detailed information about individual steps of the algorithms. This is useful for checking for the point where a computation fails.
infolevelGroebnerBasis:=5:
GroebnerBasistrinks,tdeg⁡p,s,b,z,t,w:
-> FGb domain: rat_int_cof Lin Bk ignored [NEW lib]/S:0 -> 8/ [1](2x6)100/ [2](10x18)100/ [3](14x25)100/ [3](20x30)100/ [4](41x43)100/ [4](30x34)50.0/100/ Mingbasis2 (13x23)100/restore Z1 Copy 18.25 for 23/126 exposants {92.31 per completed} SWAP Z1/2 Memory usage (estimate): 0.000 /S:1 -> 16/ Mingbasis2 restore Z1 Compute_upper bound 5 {done} SWAP Z1/2 Memory usage (estimate): 0.000
13 polynomials, 126 terms total time: 0.156 sec ------------------------------
You can also assign the infolevel of individual algorithms separately using infolevel[buchberger], infolevel[fglm], and infolevel[walk]. A higher value of infolevel[GroebnerBasis] will take precedence over these values.
infolevelGroebnerBasis:=1:
Positive-Dimensional Systems
It is sometimes possible to simplify non zero-dimensional ideals by removing solutions. In the example below, we would like to compute the smallest positive value of m. If we remove the solutions corresponding to m=0, we get an ideal which is zero-dimensional.
circles:=x2+y2−m2,2−2⁢z+z2+w2+2⁢w⁢m−2⁢w−2⁢m,z2+1−2⁢y−2⁢m+y2+2⁢y⁢m−3⁢m2,5−6⁢m−4⁢z+m2+4⁢m⁢z+z2+4⁢w2−4⁢w−4⁢y⁢w−4⁢w⁢m+2⁢y+y2+2⁢y⁢m,1−2⁢x−2⁢m+x2+2⁢x⁢m+w2
circles:=5−6⁢m−4⁢z+m2+4⁢m⁢z+z2+4⁢w2−4⁢w−4⁢y⁢w−4⁢w⁢m+2⁢y+y2+2⁢y⁢m,x2+y2−m2,2−2⁢z+z2+w2+2⁢w⁢m−2⁢w−2⁢m,1−2⁢x−2⁢m+x2+2⁢x⁢m+w2,z2+1−2⁢y−2⁢m+y2+2⁢y⁢m−3⁢m2
HilbertDimension⁡circles
1
circles1:=Saturate⁡circles,m−0
circles1:=−178798791+513052018⁢m+205926593⁢y−25508929⁢w+134703⁢z+2552037⁢x+120046030⁢m2−240278726⁢y⁢m+95369357⁢w⁢m−279141279⁢m⁢z−94004920⁢y2−16810987⁢y⁢z−5620181⁢x⁢y+65711447⁢w⁢z−17337367⁢x⁢w+73039735⁢m3+46037710⁢y⁢m2−191591206⁢w⁢m2+43664773⁢w⁢m3,x⁢m−w⁢m+y⁢m−m2−x+z+w−y−m,x2+y2−m2,17826⁢z+22615⁢y⁢m+142⁢x⁢y+3890⁢y2−5015⁢w⁢z+269⁢x−523⁢x⁢w−12060⁢y+444⁢y2⁢x+7442⁢w−539⁢m−13760⁢w⁢m−19646⁢m2−5973+8992⁢w⁢m2−10446⁢y⁢m2+9549⁢m3−8217⁢m⁢z,w2+1−y2+2⁢w⁢m−2⁢y⁢m+3⁢m2−2⁢z−2⁢w+2⁢y,255⁢z+41⁢y⁢m−3⁢x⁢y−16⁢y2−10⁢w⁢z−230+5⁢x−21⁢x⁢w+38⁢y+95⁢w+378⁢m−74⁢y⁢z−109⁢w⁢m−166⁢m2−55⁢w⁢m2+74⁢z⁢y⁢m−80⁢y⁢m2+168⁢m3−316⁢m⁢z,−12846953871⁢z−21056954801⁢y⁢m−310564157⁢x⁢y−5940128812⁢y2+6666214216⁢w⁢z−378236976+103112603⁢x−839297383⁢x⁢w+13531707252⁢y−7640771911⁢w+13704648040⁢m−4342320⁢y⁢z+392982957⁢m4+15704929507⁢w⁢m+19839428668⁢m2−14106629405⁢w⁢m2+6745180152⁢y⁢m2−4883346564⁢m3−2949654150⁢m⁢z,1158+1136⁢m+2880⁢y−2309⁢w−3861⁢z+7⁢x+4178⁢m2−5167⁢y⁢m+4391⁢w⁢m+342⁢m⁢z−1088⁢y2−19⁢x⁢y+1466⁢w⁢z−59⁢x⁢w−3037⁢w⁢m2+666⁢z⁢m2+2256⁢y⁢m2−1452⁢m3,−5+x⁢z+12⁢z+10⁢y⁢m+3⁢y2−3⁢w⁢z−2⁢x−8⁢y+8⁢w−y⁢z−8⁢w⁢m−15⁢m2−m⁢z,−14382⁢z−9853⁢y⁢m+110⁢x⁢y−4322⁢y2+6077⁢w⁢z−23⁢x+289⁢x⁢w+6423+5760⁢y+1332⁢y3−6542⁢w−1687⁢m+2664⁢y⁢z+6932⁢w⁢m+15386⁢m2−2152⁢w⁢m2−2370⁢y⁢m2+1203⁢m3+1683⁢m⁢z,−27324⁢z−33467⁢y⁢m−32⁢x⁢y−8878⁢y2+1332⁢y2⁢z+8067+6097⁢w⁢z−181⁢x+479⁢x⁢w+20484⁢y−12652⁢w+5401⁢m+23134⁢w⁢m+33922⁢m2−11810⁢w⁢m2+12738⁢y⁢m2−10455⁢m3+6903⁢m⁢z,−531⁢z−1837⁢y⁢m−19⁢x⁢y−422⁢y2+134⁢w⁢z+7⁢x−59⁢x⁢w+882⁢y−311⁢w+803⁢m+1061⁢w⁢m+333⁢y2⁢m+848⁢m2−174−1039⁢w⁢m2+924⁢y⁢m2−1119⁢m3+342⁢m⁢z,−802080060+5346880081⁢m+4152985503⁢y−2042981563⁢w−3182078763⁢z+40573910⁢x+5603587108⁢m2−6076124063⁢y⁢m+4362338179⁢w⁢m−1848236604⁢m⁢z−1874980288⁢y2−40125501⁢y⁢z−107818307⁢x⁢y+1973431840⁢w⁢z−305585482⁢x⁢w−934748403⁢m3+1582183248⁢y⁢m2−4255227986⁢w⁢m2+130994319⁢y⁢m3,y⁢w−y2−m⁢z+3⁢w⁢m−2⁢y⁢m+2⁢m2−z−w+y+m,−93⁢z−203⁢y⁢m−5⁢x⁢y−76⁢y2−29⁢w⁢z−4⁢x+2⁢x⁢w+162⁢y+35⁢w+334⁢m−83⁢w⁢m+118⁢m2+7⁢w⁢m2+27⁢y⁢m2−75+21⁢m3−132⁢m⁢z+111⁢z⁢w⁢m,z2+1−2⁢y−2⁢m+y2+2⁢y⁢m−3⁢m2
HilbertDimension⁡circles1
UnivariatePolynomial⁡m,circles1
819200−14172160⁢m+114038784⁢m2−563649536⁢m3+7918461504⁢m6+98015844⁢m16−1398966480⁢m13−4564076288⁢m5+1899131648⁢m4−9722063488⁢m7+1180129⁢m18+1145811528⁢m14−11436428⁢m17−462103584⁢m15−1038261808⁢m10−2960321792⁢m9+1526909568⁢m11+7803109440⁢m8+227573920⁢m12
For more general systems we can use the ZeroDimensionalDecomposition command. This command factors the system into a sequence of zero-dimensional ideals. In the process, variables are placed into the coefficient ring to make the resulting ideals zero-dimensional. To recover the original system, we must contract each ideal back to the original ring before computing its intersection. This technique allows us to run zero-dimensional algorithms on ideals of positive dimension, provided they interact nicely with the extension/contraction process.
cyclohexane:=−15−34⁢z−23⁢z2−34⁢x−20⁢x⁢z−23⁢x2+9⁢x2⁢z2+6⁢z⁢x2+6⁢x⁢z2,−15−23⁢y2−20⁢y⁢z−34⁢y−23⁢z2−34⁢z+9⁢y2⁢z2+6⁢y⁢z2+6⁢y2⁢z,−9⁢z3⁢x−9⁢z3⁢y−9⁢x⁢z2−9⁢y⁢z2−6⁢z3+20⁢x⁢y+41⁢x⁢z+41⁢y⁢z+18⁢z2+21⁢x+21⁢y+54⁢z+18,27⁢z4⁢x+108⁢z3⁢x+18⁢x⁢z2−212⁢x⁢z−69⁢x+27⁢z4⁢y+108⁢z3⁢y+18⁢y⁢z2−212⁢y⁢z+18⁢z4−69⁢y−284⁢z2−400⁢z−102
cyclohexane:=−15−34⁢z−23⁢z2−34⁢x−20⁢x⁢z−23⁢x2+9⁢x2⁢z2+6⁢x2⁢z+6⁢x⁢z2,−9⁢x⁢z3−9⁢y⁢z3−9⁢x⁢z2−9⁢y⁢z2−6⁢z3+20⁢x⁢y+41⁢x⁢z+41⁢y⁢z+18⁢z2+21⁢x+21⁢y+54⁢z+18,−15−23⁢y2−20⁢y⁢z−34⁢y−23⁢z2−34⁢z+9⁢y2⁢z2+6⁢y⁢z2+6⁢y2⁢z,27⁢x⁢z4+108⁢x⁢z3+18⁢x⁢z2−212⁢x⁢z−69⁢x+27⁢y⁢z4+108⁢y⁢z3+18⁢y⁢z2−212⁢y⁢z+18⁢z4−69⁢y−284⁢z2−400⁢z−102
IdealInfoVariablescyclohexane
z,x,y
HilbertDimension⁡cyclohexane
S:=Simplify⁡ZeroDimensionalDecomposition⁡cyclohexane
S:=−34+6⁢z2−20⁢z+−23+9⁢z2+6⁢z⁢x+−23+9⁢z2+6⁢z⁢y,−15−34⁢z−23⁢z2+−23+9⁢z2+6⁢z⁢x2+−34+6⁢z2−20⁢z⁢x,3⁢y2+4⁢x⁢y+3⁢x2−5⁢y⁢z−5⁢x⁢z+3⁢y+3⁢x−6⁢z,3⁢x⁢y⁢z+8⁢x⁢y+2⁢y⁢z+2⁢x⁢z−2⁢z2+3⁢y+3⁢x−3⁢z,9⁢y⁢x2−3⁢x⁢z2+12⁢x⁢y−9⁢x2+18⁢y⁢z−z2+9⁢y−9⁢x+6⁢z,9⁢x3+12+6⁢x⁢z2+3⁢x⁢y+37⁢x2−3⁢y⁢z+25⁢x⁢z+41⁢x+14⁢z,9⁢x2⁢z+9+3⁢x⁢z2−3⁢x⁢y+15⁢x2+3⁢y⁢z+27⁢x⁢z+7⁢z2+21⁢x+21⁢z,3⁢y⁢z2+3+3⁢x⁢z2+x⁢y+9⁢y⁢z+9⁢x⁢z+4⁢z2+3⁢y+3⁢x+10⁢z,3⁢z3+3⁢x⁢y−3⁢y⁢z−3⁢x⁢z+13⁢z2+3⁢z
We know immediately that the ideal is not prime because it has factored.
seq⁡IdealInfoVariablesi,i=S
x,y,z,x,y
seq⁡HilbertDimension⁡i,i=S
0,0
S1,S2:=selectremove⁡has,S,variables=x,y,zS1:=op⁡S1;S2:=op⁡S2
S1,S2:=3⁢y2+4⁢x⁢y+3⁢x2−5⁢y⁢z−5⁢x⁢z+3⁢y+3⁢x−6⁢z,3⁢x⁢y⁢z+8⁢x⁢y+2⁢y⁢z+2⁢x⁢z−2⁢z2+3⁢y+3⁢x−3⁢z,9⁢y⁢x2−3⁢x⁢z2+12⁢x⁢y−9⁢x2+18⁢y⁢z−z2+9⁢y−9⁢x+6⁢z,9⁢x3+12+6⁢x⁢z2+3⁢x⁢y+37⁢x2−3⁢y⁢z+25⁢x⁢z+41⁢x+14⁢z,9⁢x2⁢z+9+3⁢x⁢z2−3⁢x⁢y+15⁢x2+3⁢y⁢z+27⁢x⁢z+7⁢z2+21⁢x+21⁢z,3⁢y⁢z2+3+3⁢x⁢z2+x⁢y+9⁢y⁢z+9⁢x⁢z+4⁢z2+3⁢y+3⁢x+10⁢z,3⁢z3+3⁢x⁢y−3⁢y⁢z−3⁢x⁢z+13⁢z2+3⁢z,−34+6⁢z2−20⁢z+−23+9⁢z2+6⁢z⁢x+−23+9⁢z2+6⁢z⁢y,−15−34⁢z−23⁢z2+−23+9⁢z2+6⁢z⁢x2+−34+6⁢z2−20⁢z⁢x
S1:=3⁢y2+4⁢x⁢y+3⁢x2−5⁢y⁢z−5⁢x⁢z+3⁢y+3⁢x−6⁢z,3⁢x⁢y⁢z+8⁢x⁢y+2⁢y⁢z+2⁢x⁢z−2⁢z2+3⁢y+3⁢x−3⁢z,9⁢y⁢x2−3⁢x⁢z2+12⁢x⁢y−9⁢x2+18⁢y⁢z−z2+9⁢y−9⁢x+6⁢z,9⁢x3+12+6⁢x⁢z2+3⁢x⁢y+37⁢x2−3⁢y⁢z+25⁢x⁢z+41⁢x+14⁢z,9⁢x2⁢z+9+3⁢x⁢z2−3⁢x⁢y+15⁢x2+3⁢y⁢z+27⁢x⁢z+7⁢z2+21⁢x+21⁢z,3⁢y⁢z2+3+3⁢x⁢z2+x⁢y+9⁢y⁢z+9⁢x⁢z+4⁢z2+3⁢y+3⁢x+10⁢z,3⁢z3+3⁢x⁢y−3⁢y⁢z−3⁢x⁢z+13⁢z2+3⁢z
S2:=−34+6⁢z2−20⁢z+−23+9⁢z2+6⁢z⁢x+−23+9⁢z2+6⁢z⁢y,−15−34⁢z−23⁢z2+−23+9⁢z2+6⁢z⁢x2+−34+6⁢z2−20⁢z⁢x
The first ideal lies in the original polynomial ring, so we can factor it directly.
P1:=Simplify⁡PrimeDecomposition⁡S1
P1:=−1+y,3+x,z+3,9⁢x+7,1+3⁢y,3⁢z+1,1+3⁢y,3⁢z+1,1+3⁢x,3+y,z+3,−1+x,3⁢z+1,9⁢y+7,1+3⁢x,3+y,3+x,z+3
In the second ideal, the variable z has been made into a parameter. We verify that this ideal in Q(z)[x, y] is prime.
tord≔IdealInfoDefaultMonomialOrderS2,'plex';G≔GroebnerBasisS2,tord
tord:=plex⁡y,x
G:=−15−34⁢z−23⁢z2+−23+9⁢z2+6⁢z⁢x2+−34+6⁢z2−20⁢z⁢x,−34+6⁢z2−20⁢z+−23+9⁢z2+6⁢z⁢x+−23+9⁢z2+6⁢z⁢y
factor⁡G1
−15−34⁢z−23⁢z2−34⁢x−20⁢x⁢z−23⁢x2+9⁢x2⁢z2+6⁢x2⁢z+6⁢x⁢z2
map⁡degree,G,convert⁡tord,'list'
2,1
G[1] is irreducible, and G[2] is linear so it is irreducible mod G[1]. Therefore, this ideal is prime and if we contract, we will recover a prime component of the cyclohexane system.
IsPrime⁡S2
P2:=Contract⁡S2,x,y,z
P2:=−15−34⁢z−23⁢z2−34⁢x−20⁢x⁢z−23⁢x2+9⁢x2⁢z2+6⁢x2⁢z+6⁢x⁢z2,−34−23⁢x−23⁢y−20⁢z+6⁢x⁢z+6⁢y⁢z+6⁢z2+9⁢x⁢z2+9⁢y⁢z2,1+2⁢z+2⁢x+3⁢x⁢z+2⁢y+3⁢y⁢z+3⁢x⁢y
We check using the PrimeDecomposition command:
Simplify⁡PrimeDecomposition⁡cyclohexane
−15−34⁢z−23⁢z2−34⁢x−20⁢x⁢z−23⁢x2+9⁢x2⁢z2+6⁢x2⁢z+6⁢x⁢z2,−34−23⁢x−23⁢y−20⁢z+6⁢x⁢z+6⁢y⁢z+6⁢z2+9⁢x⁢z2+9⁢y⁢z2,1+2⁢z+2⁢x+3⁢x⁢z+2⁢y+3⁢y⁢z+3⁢x⁢y,1+3⁢y,3⁢z+1,1+3⁢x,3+y,z+3,−1+x,9⁢x+7,1+3⁢y,3⁢z+1,−1+y,3+x,z+3,3⁢z+1,9⁢y+7,1+3⁢x,3+y,3+x,z+3
Advanced Example: Computing the radical
In this section we will write our own command to compute the radical of an ideal. We begin by writing a program for zero-dimensional ideals. For the zero-dimensional case, for each variable we can compute a univariate polynomial in our ideal. We make all those univariate polynomials square-free and add them as generators. This is demonstrated below.
J:=x2−y2,x⁢y−1
f:=UnivariatePolynomial⁡x,J
f:=1−2⁢x3+x6
convert⁡f,'sqrfree'
−1+x32
Our first program, a subroutine, will take polynomials and make them square-free.
squarefreePart:=procf,xnormal⁡f/gcd⁡f,diff⁡f,xend proc:
squarefreePart⁡f,x
−1+x3
Now our radical program just has to do this for the univariate polynomial in each variable.
zradical:=procJlocalx,P;P:=seq⁡squarefreePart⁡UnivariatePolynomial⁡x,J,x,x=IdealInfoVariables⁡J;PolynomialIdeal⁡J,Pend proc:
zradical⁡J
−1+x3,−1+y3,x2−y2,x⁢y−1
Simplify⁡,tdegx,y
y2−x,x⁢y−1,x2−y
Simplify⁡Radical⁡J,tdeg⁡x,y
Now we need to make it work on non zero-dimensional ideals. As in the previous section, we will use the ZeroDimensionalDecomposition command to reduce the problem to the zero-dimensional case. The ideal below is not radical. We will demonstrate the technique before writing the program.
K:=x⁢y−z,y⁢x2−z3
K:=y⁢x2−z3,x⁢y−z
First we decompose the system into a sequence of zero-dimensional ideals, some of which live in an extension ring where some of the the variables are considered to be part of the coefficient field.
zdd:=ZeroDimensionalDecomposition⁡K
zdd:=x⁢y−z,−z+y⁢z2,−y+y2⁢z,z3−x⁢z,−x⁢y+z,x2
seq⁡IdealInfoVariablesi,i=zdd
z,y,z,x
seq⁡HilbertDimension⁡i,i=zdd
Now we can use our radical program above to make each component radical.
rdd:=seq⁡zradical⁡i,i=zdd
rdd:=x⁢y−z,−y+x⁢y3,−z+y⁢z2,−y+y2⁢z,z3−x⁢z,z,−x⁢y+z,x2,x
Simplify⁡
−x⁢y+z,−y+x⁢y3,z,x
Finally we contract these components back to the original ring and intersect them. The result is the radical of the original ideal.
seq⁡Contract⁡i,x,y,z,i=rdd
x⁢y−z,−z+y⁢z2,−y+y2⁢z,z3−x⁢z,z,x
R:=Intersect⁡
R:=x⁢y−z,−z+y⁢z2,z3−x⁢z
IdealContainment⁡R,Radical⁡K,R
Here is the program. Notice that we have stored the zero-dimensional decomposition as a list, in case there is only one ideal.
MyRadical:=procJlocalvars,zdd;ifHilbertDimension⁡J=0thenzradical⁡Jelsezdd:=ZeroDimensionalDecomposition⁡J;zdd:=map⁡zradical,zdd;zdd:=map⁡Contract,zdd,IdealInfoVariables⁡J;Intersect⁡op⁡zddend ifend proc:
MyRadical⁡K
x⁢y−z,−z+y⁢z2,z3−x⁢z
We will try it on a more difficult example. This system arises from trying to compute the singular locus of the discriminant of a generic monic polynomial of degree four.
f:=discrim⁡x4+a⁢x3+b⁢x2+c⁢x+d,x
f:=18⁢a3⁢c⁢b⁢d+256⁢d3−27⁢c4−6⁢a2⁢c2⁢d−192⁢a⁢c⁢d2+18⁢a⁢c3⁢b+144⁢b⁢a2⁢d2+b2⁢a2⁢c2−4⁢b3⁢a2⁢d+144⁢d⁢c2⁢b−4⁢a3⁢c3−128⁢b2⁢d2+16⁢b4⁢d−4⁢b3⁢c2−27⁢a4⁢d2−80⁢a⁢c⁢b2⁢d
The singular locus of f is the set of points where f and its partial derivatives vanish.
discriminant4:=f,seq⁡∂∂i⁢f,i=a,b,c,d
discriminant4:=27⁢a2⁢c⁢b⁢d−6⁢a⁢c2⁢d−96⁢c⁢d2+9⁢c3⁢b+144⁢b⁢a⁢d2+b2⁢a⁢c2−4⁢b3⁢a⁢d−6⁢a2⁢c3−54⁢a3⁢d2−40⁢c⁢b2⁢d,9⁢a3⁢c⁢d+9⁢a⁢c3+72⁢a2⁢d2+b⁢a2⁢c2−6⁢b2⁢a2⁢d+72⁢d⁢c2−128⁢b⁢d2+32⁢b3⁢d−6⁢b2⁢c2−80⁢a⁢c⁢b⁢d,9⁢a3⁢b⁢d−54⁢c3−6⁢a2⁢c⁢d−96⁢a⁢d2+27⁢a⁢c2⁢b+b2⁢a2⁢c+144⁢d⁢c⁢b−6⁢a3⁢c2−4⁢b3⁢c−40⁢a⁢b2⁢d,9⁢a3⁢c⁢b+384⁢d2−3⁢a2⁢c2−192⁢a⁢c⁢d+144⁢b⁢a2⁢d−2⁢b3⁢a2+72⁢c2⁢b−128⁢b2⁢d+8⁢b4−27⁢a4⁢d−40⁢a⁢c⁢b2,18⁢a3⁢c⁢b⁢d+256⁢d3−27⁢c4−6⁢a2⁢c2⁢d−192⁢a⁢c⁢d2+18⁢a⁢c3⁢b+144⁢b⁢a2⁢d2+b2⁢a2⁢c2−4⁢b3⁢a2⁢d+144⁢d⁢c2⁢b−4⁢a3⁢c3−128⁢b2⁢d2+16⁢b4⁢d−4⁢b3⁢c2−27⁢a4⁢d2−80⁢a⁢c⁢b2⁢d
The system is not zero-dimensional. The singular locus is two-dimensional surface in the four-dimensional parameter space.
HilbertDimension⁡discriminant4
2
Our radical program is powerful enough to remove the multiple solutions from this system. The resulting system will be easier to work with and solve.
R:=MyRadical⁡discriminant4
R:=b2⁢c2−4⁢b3⁢d−3⁢a⁢c3+14⁢a⁢c⁢b⁢d−18⁢a2⁢d2−6⁢d⁢c2+16⁢b⁢d2,−9⁢a3⁢d+c⁢b⁢a2+32⁢a⁢b⁢d+3⁢a⁢c2−48⁢d⁢c−4⁢c⁢b2,a2⁢b2−3⁢c⁢a3−4⁢b3+14⁢c⁢b⁢a−6⁢d⁢a2−18⁢c2+16⁢b⁢d,a⁢c2⁢b−4⁢a⁢b2⁢d+3⁢a2⁢c⁢d−9⁢c3+32⁢d⁢c⁢b−48⁢a⁢d2,a2⁢c2−3⁢b⁢a2⁢d−3⁢c2⁢b+8⁢b2⁢d+4⁢a⁢c⁢d−32⁢d2
solve⁡Generators⁡R
b=0,c=0,a=a,d=0,d=164⁢a4−18⁢b⁢a2+14⁢b2,a=a,b=b,c=12⁢a⁢b−18⁢a3,c=RootOf⁡108⁢_Z2+−108⁢a⁢b+27⁢a3⁢_Z+32⁢b3−9⁢a2⁢b2,d=14⁢RootOf⁡108⁢_Z2+−108⁢a⁢b+27⁢a3⁢_Z+32⁢b3−9⁢a2⁢b2⁢a−112⁢b2,a=a,b=b
It looks like there are two distinct pieces. We can verify this by factoring the system.
PrimeDecomposition⁡R
3⁢a⁢c−12⁢d−b2,27⁢c2−108⁢b⁢d−b3+27⁢d⁢a2,−108⁢d⁢c⁢b+108⁢a⁢d2+27⁢c3−b3⁢c+9⁢a⁢b2⁢d,−432⁢d3−72⁢b2⁢d2+108⁢d⁢c2⁢b−27⁢c4−3⁢b4⁢d+b3⁢c2,c⁢a2−4⁢c⁢b+8⁢a⁢d,c2⁢b−4⁢b2⁢d+2⁢a⁢c⁢d+16⁢d2,c3−4⁢d⁢c⁢b+8⁢a⁢d2,b⁢a2−4⁢b2+2⁢a⁢c+16⁢d,a⁢c2−4⁢a⁢b⁢d+8⁢d⁢c,d⁢a2−c2,8⁢c−4⁢a⁢b+a3
map⁡HilbertDimension,
2,2
So the locus splits into two irreducible surfaces.
One Last Example
PolynomialIdeals can factor and solve large systems. The system below has 96 prime components. It will take about two and a half minutes to run on a 700 MHz Pentium III.
schwarz8:=2⁢x4⁢x5+2⁢x6⁢x3+2⁢x7⁢x2+2⁢x8⁢x1+x8,x42+2⁢x3⁢x5+2⁢x6⁢x2+2⁢x1⁢x7+x82+x7,2⁢x3⁢x4+2⁢x2⁢x5+2⁢x1⁢x6+2⁢x7⁢x8+x6,x32+2⁢x2⁢x4+2⁢x1⁢x5+x72+2⁢x8⁢x6+x5,2⁢x2⁢x3+2⁢x1⁢x4+2⁢x6⁢x7+2⁢x8⁢x5+x4,x22+2⁢x1⁢x3+x62+2⁢x7⁢x5+2⁢x8⁢x4+x3,2⁢x1⁢x2+2⁢x5⁢x6+2⁢x7⁢x4+2⁢x8⁢x3+x2,x12+x52+2⁢x6⁢x4+2⁢x7⁢x3+2⁢x8⁢x2+x1
schwarz8:=2⁢x4⁢x5+2⁢x6⁢x3+2⁢x7⁢x2+2⁢x8⁢x1+x8,x42+2⁢x3⁢x5+2⁢x6⁢x2+2⁢x1⁢x7+x82+x7,2⁢x3⁢x4+2⁢x2⁢x5+2⁢x1⁢x6+2⁢x7⁢x8+x6,2⁢x1⁢x2+2⁢x5⁢x6+2⁢x7⁢x4+2⁢x8⁢x3+x2,2⁢x2⁢x3+2⁢x1⁢x4+2⁢x6⁢x7+2⁢x8⁢x5+x4,x22+2⁢x1⁢x3+x62+2⁢x7⁢x5+2⁢x8⁢x4+x3,x32+2⁢x2⁢x4+2⁢x1⁢x5+x72+2⁢x8⁢x6+x5,x12+x52+2⁢x6⁢x4+2⁢x7⁢x3+2⁢x8⁢x2+x1
tord≔IdealInfoDefaultMonomialOrderschwarz8,'plex';time⁡GroebnerBasisschwarz8,tord
tord:=plex⁡x1,x2,x7,x5,x8,x4,x6,x3
-> FGb domain: rat_int_cof Lin Bk ignored [NEW lib]/S:0 -> 8/ [3](68x159)25.0/100/ [4](247x377)100/ [5](485x606)1.6/53.2/100/ [6](729x819)3.3/38.5/73.6/100/ [7](723x800)8.2/40.8/73.5/100/ [8](593x709)8.6/54.3/100/ [9](391x580)0.0/100/ [10](146x302)100/ Mingbasis2 (126x380)6.3/57.1/100/restore Z1 Copy 21.62 for 380/1758 exposants {done} SWAP Z1/2 Memory usage (estimate): 0.000 126 polynomials, 1758 terms total time: 0.281 sec ------------------------------ -> FGLM algorithm rational coefficients current monomial 1 current monomial x3 current monomial x3^2 current monomial x3^3 current monomial x3^4 current monomial x3^5 current monomial x3^6 current monomial x3^7
current monomial x3^8 current monomial x3^9 current monomial x3^10 current monomial x3^11 current monomial x3^12 current monomial x3^13 current monomial x3^14 current monomial x3^15 current monomial x3^16 current monomial x3^17 current monomial x3^18 current monomial x3^19 current monomial x3^20 current monomial x3^21 current monomial x3^22 current monomial x3^23 current monomial x3^24 current monomial x3^25 current monomial x6 current monomial x6*x3 current monomial x6*x3^2 current monomial x6*x3^3 current monomial x6*x3^4 current monomial x6*x3^5 current monomial x6*x3^6 current monomial x6*x3^7
current monomial x6*x3^8 current monomial x6*x3^9 current monomial x6*x3^10 current monomial x6*x3^11 current monomial x6*x3^12 current monomial x6*x3^13 current monomial x6*x3^14 current monomial x6*x3^15 current monomial x6*x3^16 current monomial x6*x3^17 current monomial x6*x3^18 current monomial x6*x3^19
current monomial x6*x3^20 current monomial x6*x3^21 current monomial x6^2 current monomial x6^2*x3 current monomial x6^2*x3^2 current monomial x6^2*x3^3 current monomial x6^2*x3^4 current monomial x6^2*x3^5 current monomial x6^2*x3^6 current monomial x6^2*x3^7 current monomial x6^2*x3^8 current monomial x6^2*x3^9 current monomial x6^2*x3^10 current monomial x6^2*x3^11 current monomial x6^2*x3^12 current monomial x6^2*x3^13 current monomial x6^2*x3^14 current monomial x6^2*x3^15 current monomial x6^2*x3^16 current monomial x6^2*x3^17 current monomial x6^2*x3^18 current monomial x6^2*x3^19 current monomial x6^2*x3^20 current monomial x6^3 current monomial x6^3*x3 current monomial x6^3*x3^2 current monomial x6^3*x3^3
current monomial x6^3*x3^4 current monomial x6^3*x3^5 current monomial x6^3*x3^6 current monomial x6^3*x3^7 current monomial x6^3*x3^8 current monomial x6^3*x3^9 current monomial x6^3*x3^10 current monomial x6^3*x3^11 current monomial x6^3*x3^12 current monomial x6^3*x3^13 current monomial x6^3*x3^14 current monomial x6^3*x3^15 current monomial x6^3*x3^16 current monomial x6^3*x3^17 current monomial x6^3*x3^18 current monomial x6^3*x3^19 current monomial x6^3*x3^20 current monomial x6^4 current monomial x6^4*x3 current monomial x6^4*x3^2