A bug in Sympy's solver

[emanuel.c...@gmail.com]

A bug in Sympy’s solve

Inspired by this ask.sagemath.orgquestion.

This question was about the names of the symbolic variables created by Maxima’s solve to denote arbitrary constants. Exploring the example used in this questions revealed a non-trivial problem with sympy’s solve.

Problem

Solve the system :
$$
\begin{align*}

• a{1}^{3} a{2} + a{1} a{2}^{2} &= 0 \
• 3 a{1}^{2} a{2} b{1} + 2 a{1} a{2} b{2} - a{1} b{2} + a{2}^{2} b{2} &= 0 \
• a{1}^{2} a{2}^{2} + a_{2}^{3} &= 0 \
• 2 a{1}^{2} a{2} b{2} - 2 a{2}^{2} b{1} + 3 a{2}^{2} b_{2} &= 0
  \end{align*}
  $$

# Set up sympy (brutal version)
import sympy
from sympy import *
init_session()
init_printing(pretty_print=False)
# System to solve
a1, a2, b1, b2 = symbols('a1 a2 b1 b2')
Unk = [a1, a2, b1, b2]
eq1 = a1 * a2**2 - a2 * a1**3
eq2 = 2*a1*a2*b2 + b2*a2**2 - 3*a2*a1**2*b1 - a1*b2
eq3 = a2**3 - a2**2*a1**2
eq4 = 3*a2**2*b2 - 2*a2*a1**2*b2 - 2*a2**2*b1
Sys = [eq1, eq2, eq3, eq4]

IPython console for SymPy 1.9 (Python 3.9.9-64-bit) (ground types: gmpy)

These commands were executed:
>>> from __future__ import division
>>> from sympy import *
>>> x, y, z, t = symbols('x y z t')
>>> k, m, n = symbols('k m n', integer=True)
>>> f, g, h = symbols('f g h', cls=Function)
>>> init_printing()

Documentation can be found at https://docs.sympy.org/1.9/

Attempt to use the “automatic” Sympy solver:

Sol = solve(Sys, Unk)
DSol = [dict(zip(Unk, u)) for u in Sol]
DSol

[{a1: 0, a2: 0, b1: b1, b2: b2}, {a1: 0, a2: 0, b1: b1, b2: b2}, {a1: 0, a2: 0, b1: b1, b2: b2}, {a1: 0, a2: 0, b1: b2/2, b2: b2}, {a1: a1, a2: 0, b1: b1, b2: 0}, {a1: -sqrt(a2), a2: a2, b1: 0, b2: 0}, {a1: -2**(5/6)*sqrt(48*2**(1/3) + 624/(1499 + 3*sqrt(303)*I)**(1/3) + 3*2**(2/3)*(1499 + 3*sqrt(303)*I)**(1/3))/6, a2: 16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3, b1: b2/2, b2: b2}, {a1: -sqrt(6)*sqrt(32 - 2**(1/3)*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3) - 416*2**(2/3)/((1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)))/6, a2: 2**(2/3)*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)/((1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)), b1: b2/2, b2: b2}, {a1: -sqrt(6)*sqrt(32 - 416*2**(2/3)/((1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(1/3)*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3))/6, a2: 2**(2/3)*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)/((1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)), b1: b2/2, b2: b2}]

Something went sideways: the six first solutions are okay,but the last three use expressions, some of them being polynomials in b2.

Attempt to check them formally :

Chk=[[u.subs(s) for u in Sys] for s in DSol]
Chk

[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [-2**(5/6)*(16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)**2*sqrt(48*2**(1/3) + 624/(1499 + 3*sqrt(303)*I)**(1/3) + 3*2**(2/3)*(1499 + 3*sqrt(303)*I)**(1/3))/6 + sqrt(2)*(16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)*(48*2**(1/3) + 624/(1499 + 3*sqrt(303)*I)**(1/3) + 3*2**(2/3)*(1499 + 3*sqrt(303)*I)**(1/3))**(3/2)/54, -2**(2/3)*b2*(16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)*(48*2**(1/3) + 624/(1499 + 3*sqrt(303)*I)**(1/3) + 3*2**(2/3)*(1499 + 3*sqrt(303)*I)**(1/3))/12 - 2**(5/6)*b2*(16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)*sqrt(48*2**(1/3) + 624/(1499 + 3*sqrt(303)*I)**(1/3) + 3*2**(2/3)*(1499 + 3*sqrt(303)*I)**(1/3))/3 + b2*(16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)**2 + 2**(5/6)*b2*sqrt(48*2**(1/3) + 624/(1499 + 3*sqrt(303)*I)**(1/3) + 3*2**(2/3)*(1499 + 3*sqrt(303)*I)**(1/3))/6, -2**(2/3)*(16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)**2*(48*2**(1/3) + 624/(1499 + 3*sqrt(303)*I)**(1/3) + 3*2**(2/3)*(1499 + 3*sqrt(303)*I)**(1/3))/18 + (16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)**3, -2**(2/3)*b2*(16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)*(48*2**(1/3) + 624/(1499 + 3*sqrt(303)*I)**(1/3) + 3*2**(2/3)*(1499 + 3*sqrt(303)*I)**(1/3))/9 + 2*b2*(16/3 + 104*2**(2/3)/(3*(1499 + 3*sqrt(303)*I)**(1/3)) + 2**(1/3)*(1499 + 3*sqrt(303)*I)**(1/3)/3)**2], [2**(1/6)*sqrt(3)*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)*(32 - 2**(1/3)*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3) - 416*2**(2/3)/((1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)))**(3/2)/(18*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(5/6)*sqrt(3)*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**2*sqrt(32 - 2**(1/3)*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3) - 416*2**(2/3)/((1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)))/(3*(1 + sqrt(3)*I)**2*(1499 + 3*sqrt(303)*I)**(2/3)), -2*2**(1/6)*sqrt(3)*b2*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)*sqrt(32 - 2**(1/3)*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3) - 416*2**(2/3)/((1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)))/(3*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 3*2**(2/3)*b2*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)*(16/3 - 2**(1/3)*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)/6 - 208*2**(2/3)/(3*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)))/(2*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) + sqrt(6)*b2*sqrt(32 - 2**(1/3)*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3) - 416*2**(2/3)/((1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)))/6 + 2*2**(1/3)*b2*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**2/((1 + sqrt(3)*I)**2*(1499 + 3*sqrt(303)*I)**(2/3)), -2*2**(1/3)*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**2*(16/3 - 2**(1/3)*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)/6 - 208*2**(2/3)/(3*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)))/((1 + sqrt(3)*I)**2*(1499 + 3*sqrt(303)*I)**(2/3)) + 4*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**3/((1 + sqrt(3)*I)**3*(1499 + 3*sqrt(303)*I)), -2*2**(2/3)*b2*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)*(16/3 - 2**(1/3)*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)/6 - 208*2**(2/3)/(3*(1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)))/((1 + sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) + 4*2**(1/3)*b2*(-208/3 + (1 + sqrt(3)*I)*(32 + (-1 - sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**2/((1 + sqrt(3)*I)**2*(1499 + 3*sqrt(303)*I)**(2/3))], [-2**(5/6)*sqrt(3)*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**2*sqrt(32 - 416*2**(2/3)/((1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(1/3)*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3))/(3*(1 - sqrt(3)*I)**2*(1499 + 3*sqrt(303)*I)**(2/3)) + 2**(1/6)*sqrt(3)*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)*(32 - 416*2**(2/3)/((1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(1/3)*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3))**(3/2)/(18*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)), sqrt(6)*b2*sqrt(32 - 416*2**(2/3)/((1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(1/3)*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3))/6 + 2*2**(1/3)*b2*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**2/((1 - sqrt(3)*I)**2*(1499 + 3*sqrt(303)*I)**(2/3)) - 3*2**(2/3)*b2*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)*(16/3 - 208*2**(2/3)/(3*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(1/3)*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)/6)/(2*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2*2**(1/6)*sqrt(3)*b2*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)*sqrt(32 - 416*2**(2/3)/((1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(1/3)*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3))/(3*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)), 4*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**3/((1 - sqrt(3)*I)**3*(1499 + 3*sqrt(303)*I)) - 2*2**(1/3)*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**2*(16/3 - 208*2**(2/3)/(3*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(1/3)*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)/6)/((1 - sqrt(3)*I)**2*(1499 + 3*sqrt(303)*I)**(2/3)), -2*2**(2/3)*b2*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)*(16/3 - 208*2**(2/3)/(3*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) - 2**(1/3)*(1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)/6)/((1 - sqrt(3)*I)*(1499 + 3*sqrt(303)*I)**(1/3)) + 4*2**(1/3)*b2*(-208/3 + (1 - sqrt(3)*I)*(32 + (-1 + sqrt(3)*I)*(2998 + 6*sqrt(303)*I)**(1/3))*(2998 + 6*sqrt(303)*I)**(1/3)/12)**2/((1 - sqrt(3)*I)**2*(1499 + 3*sqrt(303)*I)**(2/3))]]

Some of these expressions are polynomial in b2 whose coefficients of not-null dregree are nort obviously null.
A bit of hit-and-miss trials leads to this attempt at numerical check :

def chz(x):
    r = x.factor().is_zero
    if r is None: return x.coeff(b2).n()
    return r
[[chz(u) for u in v] for v in Chk]

[[True, True, True, True],
 [True, True, True, True],
 [True, True, True, True],
 [True, True, True, True],
 [True, True, True, True],
 [True, True, True, True],
 [True, -223.427427555427 + 0.e-25*I, True, True],
 [True, -0.388012232966408 - 5.25038392589403e-29*I, True, True],
 [True, -0.e-137 + 1.12445821450677e-139*I, True, True]]

Solution #8 may be exactr, but #6 and #7 cannot.

Manual solution.

eq1 and eq3 give us solutions for a1 and a2 :

S13=solve([eq3, eq1], [a2, a1], dict=True) ; S13

[{a1: -sqrt(a2)}, {a1: sqrt(a2)}, {a2: 0}]

which we prefer to rewrite as :

S13 = [{a2:a1**2}, {a2:0}]; S13

[{a2: a1**2}, {a2: 0}]

Substituting these values in eq4 gives us :

E4 = [eq4.subs(s) for s in S13] ; E4

[-2*a1**4*b1 + a1**4*b2, 0]

The second solution tells us that S13[1] is also,a solution to [eq1, eq3, eq4].

S134=S13[1:] ; S134

[{a2: 0}]

Substituting S13[0] ineq4and solving for the variables of the resulting polynomial augments the setS134of solutions of[eq1, eq3`, eq4]’ :

V4=E4[0].free_symbols
S4=[solve(E4[0], v, dict=True) for v in V4] ; S4
for S in S4:
    for s in S:
        S0={u:S13[0][u].subs(s) for u in S13[0].keys()}
        S134 += [S0.copy()|s]
S134

[{a2: 0}, {a2: a1**2, b1: b2/2}, {a1: 0, a2: 0}, {a2: a1**2, b2: 2*b1}]

Again, we prefer to rewrite it in a simpler (and shorter) fashion :

S134=[{a2:0}, {a2:a1**2, b2:2*b1}] ; S134

[{a2: 0}, {a2: a1**2, b2: 2*b1}]

Substituting in eq3 gives :

[eq2.subs(s) for s in S134]

[-a1*b2, -a1**4*b1 + 4*a1**3*b1 - 2*a1*b1]

Solving these equations for their free symbols a,d merging with the previous partial solutions gives us the solutions of the full system :

S1234=[]
for S in S134:
    # print("S=",S)
    E=eq2.subs(S)
    # print("E=",E)
    S1=[solve(E, v, dict=True) for v in E.free_symbols]
    # print("S1=",S1)
    for s in flatten(S1):
        # print("s=",s)
        S0={u:S[u].subs(s) if "subs" in dir(S[u]) else S[u] for u in S.keys()}
        S1234+=[S0.copy()|s]
S1234

[{a1: 0, a2: 0}, {a2: 0, b2: 0}, {a2: a1**2, b1: 0, b2: 0}, {a1: 0, a2: 0, b2: 2*b1}, {a1: 4/3 + (-1/2 - sqrt(3)*I/2)*(37/27 + sqrt(303)*I/9)**(1/3) + 16/(9*(-1/2 - sqrt(3)*I/2)*(37/27 + sqrt(303)*I/9)**(1/3)), a2: (4/3 + (-1/2 - sqrt(3)*I/2)*(37/27 + sqrt(303)*I/9)**(1/3) + 16/(9*(-1/2 - sqrt(3)*I/2)*(37/27 + sqrt(303)*I/9)**(1/3)))**2, b2: 2*b1}, {a1: 4/3 + 16/(9*(-1/2 + sqrt(3)*I/2)*(37/27 + sqrt(303)*I/9)**(1/3)) + (-1/2 + sqrt(3)*I/2)*(37/27 + sqrt(303)*I/9)**(1/3), a2: (4/3 + 16/(9*(-1/2 + sqrt(3)*I/2)*(37/27 + sqrt(303)*I/9)**(1/3)) + (-1/2 + sqrt(3)*I/2)*(37/27 + sqrt(303)*I/9)**(1/3))**2, b2: 2*b1}, {a1: 4/3 + 16/(9*(37/27 + sqrt(303)*I/9)**(1/3)) + (37/27 + sqrt(303)*I/9)**(1/3), a2: (4/3 + 16/(9*(37/27 + sqrt(303)*I/9)**(1/3)) + (37/27 + sqrt(303)*I/9)**(1/3))**2, b2: 2*b1}]

Again, some solutions, substituted in eq2, give first-degree monomials in b1 whose oefficient cannot be shownt to be null byis_zero :

[[e.subs(s).is_zero for e in Sys] for s in S1234]

[[True, True, True, True],
 [True, True, True, True],
 [True, True, True, True],
 [True, True, True, True],
 [True, None, True, True],
 [True, None, True, True],
 [True, None, True, True]]

But, this time, the numerical check points to a probably null result :

[Sys[1].subs(S1234[u]).coeff(b1).n() for u in range(3,6)]

[0, 0.e-125 + 0.e-127*I, 0.e-125 - 0.e-127*I]

​

[Chris Smith]

`nonlinsolve(Sys, Unk)` gives

```

{(0, 0, b1, b2), (0, 0, b2/2, b2), (a1, 0, b1, 0), (-sqrt(a2), a2, 0, 0), (sqrt(a2), a2, 0, 0)}

```

[Chris Smith]

If you factor the equations before passing them to solve then the solution is

```

>>> fsol

[{a2: 0, a1: 0}, {a2: 0, a1: 0}, {a2: 0, b2: 0}, {a2: 0, b1: b2/2, a1:

0}, {a2: (64 + (-8 + (1 - sqrt(3)*I)*(37 + 3*sqrt(303)*I)**(1/3))*(1 -

sqrt(3)*I)*(37 + 3*sqrt(303)*I)**(1/3))**2/(36*(1 - sqrt(3)*I)**2*(37

+ 3*sqrt(303)*I)**(2/3)), b1: b2/2, a1: (-64 + (1 - sqrt(3)*I)*(8 +

(-1 + sqrt(3)*I)*(37 + 3*sqrt(303)*I)**(1/3))*(37 +

3*sqrt(303)*I)**(1/3))/(6*(1 - sqrt(3)*I)*(37 +

3*sqrt(303)*I)**(1/3))}, {a2: (64 + (-8 + (1 + sqrt(3)*I)*(37 +

3*sqrt(303)*I)**(1/3))*(1 + sqrt(3)*I)*(37 +

3*sqrt(303)*I)**(1/3))**2/(36*(1 + sqrt(3)*I)**2*(37 +

3*sqrt(303)*I)**(2/3)), b1: b2/2, a1: (-64 + (1 + sqrt(3)*I)*(8 - (1 +

sqrt(3)*I)*(37 + 3*sqrt(303)*I)**(1/3))*(37 +

3*sqrt(303)*I)**(1/3))/(6*(1 + sqrt(3)*I)*(37 +

3*sqrt(303)*I)**(1/3))}, {a2: (16 + (4 + (37 +

3*sqrt(303)*I)**(1/3))*(37 + 3*sqrt(303)*I)**(1/3))**2/(9*(37 +

3*sqrt(303)*I)**(2/3)), b1: b2/2, a1: 4/3 + 16/(3*(37 +

3*sqrt(303)*I)**(1/3)) + (37 + 3*sqrt(303)*I)**(1/3)/3}, {a2: a1**2,

b1: 0, b2: 0}]

```

and all of them are true solutions:

```

>>> [[e.subs(s).simplify() for e in Sys] for s in fsol]

[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]

```

[Chris Smith]

It just looks like the OP issue has two extra solutions (when the initial equations are not factored) which couldn't be eliminated (easily) and the solver sends them back to be safe.

NOTE: `fsol` in the previous was obtained with `solve([factor(i) for i in Sys], Unk, simplify=False)`.

If checking is turned off then two invalid solutions are included, but the solution involving `sqrt(303)` is missing:

```python

>>> s=solve([factor(i) for i in Sys], Unk, check=0,dict=1); s

[{a2: 0, a1: 0}, {a2: 0, a1: 0}, {a1: 0, a2: 0}, {a1: 0, b2: 0}, {a2: 0, b2: 0}, {a2: a1**2, b1: b2*(a1**3 + 2*a1**2 - 1)/(3*a1**3)}]

>>> [[e.subs(si).simplify() for e in Sys] for si in s]

[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, a2**3, -2*a2**2*b1], [0, 0, 0, 0], [0, 0, 0, a1*b2*(a1**3 - 4*a1**2 + 2)/3]]

```

/c

[emanuel.c...@gmail.com]

It’s a bit more involved than that; see the parallel discussion of the same problem in this post on Sagemath support list.

​

[emanuel.c...@gmail.com]

It turns ut that Chris' sugestion is effective in sympy 1.9 (current), but not in sympy 1.8 (still used by Sage 9.5.beta7). I still think that this equation preprocessing shouldn't be necessary...

Sorry for the noise.