# Solving polynomial equations using Grobner basis technique

Given a list of 24 polynomials over $\mathbb{Z}[x,y,z]$ namely Lf, computed by sage using the code

R.<x,y,z>=QQ['x, y, z']
I=(x*y^3-z)*R
Q= R.quo(I)
N=15
p=3
q=5
L=[(0, 0, 0),(0, 0, 1),(0, 0, 2),(0, 0, 3),(0, 0, 4),(0, 0, 5),(0, 0, 6),(0, 0, 7),(0, 1, 0),(0, 1, 1),(0, 1, 2),(0, 2, 0),(0,2,1),
(0, 2, 2),(1, 0, 0),(1, 0, 1),(1, 0, 2),(1, 0, 3),(1, 0, 4),(1, 1, 0),(1, 1, 1),(1, 1, 2),(2, 0, 0),(2, 0, 1)]
X=floor(N^2)
Y=floor(N^2)
Z=X*Y^3
phi=((p^4-1)*(q^4-1))/((p-1)*(q-1))
e=79
d=79
f=x*(y^3+(N+1)*y^2+(N^2-2*N+1)*y+N^3-N^2-N+1)+1
M=Matrix(0,len(L))
for (k,i,j) in L:
g=Q(x^i*y^j*f^k*e^(2-k)).lift()
h=g(x*X,y*Y,z*Z)
L1=vector(h[x^i*y^j*z^k] for (k,i,j) in L)
M = M.stack(L1)
N=M.LLL()

L2=[x^i*y^j*z^k for (k,i,j) in L]
r=0
Lf=[]
Li=[]
for i in srange(24):
D=N[i,:]*transpose(matrix([L2]))
r=D[0,0]
Li=[r]
Lf.extend(Li)
print(Lf)


In the theory, I know that there are three polynomials from the list Lf say $p(x,y,z)$, $q(x,y,z)$, and $r(x,y,z)$ having the same root (x,y,z). To extract such root, it suffices to compute the grobner basis of the three polynomials which yields three polynomilas and solving them to obtain the common root. I have computed the grobner basis by sage, here it is denoted by $E$.

for i in srange(24):
for j in srange(i+1,24):
for k in srange(j+1,24):
J = R.ideal([Lf[i],Lf[j],Lf[k]])
E=J.groebner_basis()
if len(E)==3:
print(E)


So far, E consists of the three polynomials we want. My probelm is what is the command to solve them.?

edit retag close merge delete

What do you want to do if J.dimension() > 0 (i.e. there are not only common points but there is a common curve or surface)? And where do you seek the point? Over the complex numbers CC?

( 2024-06-19 12:03:18 +0200 )edit

Sort by » oldest newest most voted

The set of solutions, when it's finite, can be computed with .variety() method. In your problem, there are many cases when the ideal does not depend on one of the variables, say z, essentially making it free. To take care of those cases, variety can be computed over the subring depending only on the rest of the variables (say, x and y). The corresponding code is given below.

For example, when it prints {x: 0, y: 0} it means that solutions are given by x=0, y=0 and any z.

In cases when variety still fails, it reports an error like

ERROR: Dimension: 1 Basis: [x - 3796875/8*z - 1/5400, y]
ERROR: Dimension: 1 Basis: [x^2 - 21515625/143*x*z + 14416259765625/143*z^2 - 62/7425*x + 11250/143*z + 1/65154375, y]


In some cases, one variable can be expressed in terms of another (like x in terms of z in the first error above), while others (like the second error above) may lead to quadratic Pell-type, elliptic, or higher-degree equations. One may try to address them with an appropriate functionality in Sage if needed.

from sage.misc.verbose import set_verbose
set_verbose(-1)

for fff in Combinations(Lf,3):
B = R.ideal(fff).groebner_basis()
V = set.union(*(set(f.variables()) for f in B))
R1 = PolynomialRing(QQ,len(V),', '.join(map(str,V)))
J = R1.ideal(B)
try:
S = J.variety()
except:
print('ERROR: Dimension:',J.dimension(),'Basis:',B)
continue
if S:
print(S)


ADDED. Here is a reduced list of solutions (without repeats) resulted from cases not involving quadratic or higher-degree equations in each of the variables:

{y: 0}
{y: 0, x: 0}
{y: 0, x: -16/2475}
{y: 0, x: -1/525}
{y: 0, x: -23/225}
{y: 0, x: -19/4050}
{y: 0, x: 3796875/8*z + 1/5400}
{z: -21499/568107421875, y: 0, x: 0}
{z: -79/7688671875, y: 0, x: 0}
{z: -1/2562890625, y: 0, x: 0}
{z: 24964/1606932421875, y: 0, x: -16/2475}
{z: 6241/414333984375, y: 0, x: 2401/327375}
{z: (-6733125*x^2 + 472050*x + 32)/(20182763671875*x - 82012500000), y: 0}
{z: 118579/2145139453125, y: 0, x: -23/225}
{z: 68651/2493692578125, y: 0, x: 967/72975}
{z: -79/5638359375, y: 0, x: -16/2475}
{z: 755161/57249850781250, y: 0, x: -19/4050}
{z: -79/7688671875, y: 0, x: -19/4050}
{z: -553/2562890625, y: 0, x: -23/225}
{z: 6241/843191015625, y: 0, x: -1/525}
{z: -79/17940234375, y: 0, x: -1/525}
{z: -1/2562890625, y: -19/3600, x: 0}
{z: 493/48694921875, y: 19/3600, x: 0}
{z: -2209/1355769140625, y: -23/2700, x: 0}
{z: -79/7688671875, y: -4/1425, x: 0}
{z: 1381/2050312500, y: 2/675, x: 3/25}
{z: -12637/5207793750000, y: -4/675, x: 3/3175}

more

A lazy alternative to @Max Alekseyev's excellent answer is to use one of Sage's symbolic solvers, which already take care of non-zero-dimension solutions, and to filter out the undeterminates returned by these solvers.

After running :

SymUnk=list(map(SR, (x, y, z)))
def foo():
Sols=[]
for fff in Combinations(Lf, 3):
if len((B:=RR.ideal(fff).groebner_basis()))==3:
SSol=solve([SR(u)==0 for u in B],
SymUnk,
solution_dict=True)
Sol=[]
for S in SSol:
SRR={}
for u in S.keys():
try:
SRR|={RR(u):RR(S[u])}
except:
continue
Sol+=[SRR]
Sols+=[Sol]
return Sols


Running :

sage: %time Sols=foo()
CPU times: user 2min 1s, sys: 105 ms, total: 2min 1s
Wall time: 1min 48s
sage: len(Sols)
245


returns 245 solutions. However, many of them are duplicates, as can be shown by indexing them on the groebner_basis() solved :

def bar():
Sols={}
for fff in Combinations(Lf, 3):
if len((B:=RR.ideal(fff).groebner_basis()))==3:
SSol=solve([SR(u)==0 for u in B],
SymUnk,
solution_dict=True)
Sol=[]
for S in SSol:
SRR={}
for u in S.keys():
try:
SRR|={RR(u):RR(S[u])}
except:
continue
Sol+=[SRR]
Sols|={B:Sol}
return Sols

sage: %time DSols=bar()
CPU times: user 2min 2s, sys: 84.4 ms, total: 2min 2s
Wall time: 1min 49s
sage: len(DSols)
156


These solutions to 156 distinct Groebner bases give 260 solutions :

sage: Set(flatten(list(DSols.values()))).cardinality()
260


In the present case, the laziness isn't too costly (2 minutes is way less than the time necessary to rig up a general solution with polynomial rings tools...).

UPDATE :

It turns out that Sympy's solver does a better job than Sage's for dimension>0 ideals. An alternative way to write the solutions is :

 # List of pairs (triplets, base) with Groebnar base of length 3
Triplets = [[u, b]
for u in Combinations(Lf, 3)
if len(b:=Rr.ideal(u).groebner_basis())==3]
# List of *distinct* bases
Bases = list(uniq([u[1] for u in Triplets]))
# Solve them *IN SR* using sympy :
# Precompute the unknowns
Unk = list(map(lambda u:SR(u), Rr.gens()))
# Solve
SSols = [solve([SR(u) for u in v], Unk, algorithm="sympy") for v in Bases]


FWIW :

sage: %time SSols = [solve([SR(u) for u in v], Unk, algorithm="sympy") for v in Bases]
CPU times: user 3.27 s, sys: 0 ns, total: 3.27 s
Wall time: 3.27 s


Sympy's solver seems also faster than Sage's...

Since neither the triplets nor the bases can be used as dictionary keys ("unhashable"), we need a bit of record-keeping :

# Record the base-solutions link
BS=list(map(list, zip(Bases, SSols)))
# Link triplets to their solutions
TBS=[[u[0]]+BS[Bases.index(u[1])] for u in Triplets]


Expressing the SR solutions in the original ring is lazily left to the reader...

HTH

more