Another problem with Sympy's solver.

60 views
Skip to first unread message

emanuel.c...@gmail.com

unread,
Dec 21, 2021, 5:25:36 AM12/21/21
to sympy
After :

```
from sympy import symbols, solve
f, q, r, h, i, g = map(symbols, 'fqrhig')
e1 = -4*f**2*q - 5*f**2*r - 5*f*r + 4*q**2 + 3*q*r
e2 = -4*f**2*h*q - 4*f**2*i*q - 8*f*g*q - 10*f*g*r - 8*f*h*q - 20*f*h*r - 10*f*i*r + 8*f*q**2/5 - 5*g*r + 4*h*q**2 - 5*h*r + 4*i*q**2 + 4*q**2/5 + 10*q*r + 5*r**2
```
trying `solve([e1, e2], [f, g])` never returns.

FWIW, in Sage, with equivalent definitions :

```
sage: %time solve([e1, e2], [f, g], solution_dict=True)
CPU times: user 935 ms, sys: 12.1 ms, total: 947 ms
Wall time: 671 ms
[{f: -1/2*(5*r + sqrt(64*q^3 + 128*q^2*r + 5*(12*q + 5)*r^2))/(4*q + 5*r),
  g: -1/5*(5120*h*q^5 - 1024*q^6 + 625*(12*(2*h + i)*q + 10*h + 5*i)*r^4 + 50*(10*(94*h + 38*i - 1)*q^2 - 24*q^3 + 25*(5*h + i)*q)*r^3 + 80*(5*(132*h + 36*i - 1)*q^3 - 44*q^4 + 25*h*q^2)*r^2 + 128*(5*(41*h + 5*i)*q^4 - 26*q^5)*r + (80*(2*h + 2*i + 11)*q^3*r + 64*q^4 + 125*(5*h + 5*i + 18*q)*r^3 + 625*r^4 + 50*(4*(h + i + 12)*q^2 + 5*(h + i)*q)*r^2)*sqrt(64*q^3 + 128*q^2*r + 5*(12*q + 5)*r^2))/(1024*q^5 + 4608*q^4*r + 125*(12*q + 5)*r^4 + 200*(28*q^2 + 5*q)*r^3 + 80*(96*q^3 + 5*q^2)*r^2)},
 {f: -1/2*(5*r - sqrt(64*q^3 + 128*q^2*r + 5*(12*q + 5)*r^2))/(4*q + 5*r),
  g: -1/5*(5120*h*q^5 - 1024*q^6 + 625*(12*(2*h + i)*q + 10*h + 5*i)*r^4 + 50*(10*(94*h + 38*i - 1)*q^2 - 24*q^3 + 25*(5*h + i)*q)*r^3 + 80*(5*(132*h + 36*i - 1)*q^3 - 44*q^4 + 25*h*q^2)*r^2 + 128*(5*(41*h + 5*i)*q^4 - 26*q^5)*r - (80*(2*h + 2*i + 11)*q^3*r + 64*q^4 + 125*(5*h + 5*i + 18*q)*r^3 + 625*r^4 + 50*(4*(h + i + 12)*q^2 + 5*(h + i)*q)*r^2)*sqrt(64*q^3 + 128*q^2*r + 5*(12*q + 5)*r^2))/(1024*q^5 + 4608*q^4*r + 125*(12*q + 5)*r^4 + 200*(28*q^2 + 5*q)*r^3 + 80*(96*q^3 + 5*q^2)*r^2)}]
```

WTF ?

Chris Smith

unread,
Dec 21, 2021, 10:41:06 AM12/21/21
to sympy
Are they using the same versions of SymPy?

e1 is quadratic only if f and e2 is linear in g and quadratic in f...it is a pretty simple system but there are lots of extraneous symbols. Even paring down the symbols to the system `(F1*(F2 + F3*f + f**2), G3*(G4 + G6*f**2 + G7*g + f*(G1 + G5*g - g))` takes a long time to solve and simplify.

/c

Oscar Benjamin

unread,
Dec 21, 2021, 11:33:42 AM12/21/21
to sympy
It solves a lot quicker using check=False:

In [1]: from sympy import symbols, solve
...: f, q, r, h, i, g = map(symbols, 'fqrhig')
...: e1 = -4*f**2*q - 5*f**2*r - 5*f*r + 4*q**2 + 3*q*r
...: e2 = -4*f**2*h*q - 4*f**2*i*q - 8*f*g*q - 10*f*g*r - 8*f*h*q -
20*f*h*r - 10*f*i*r + 8*f*q**2/5 - 5*g*r + 4*h*q**2 - 5*h*r + 4*i*q**2
+
...: 4*q**2/5 + 10*q*r + 5*r**2

In [2]: %time print(solve([e1, e2], [f, g], check=False))
[(-(-1280*h*q**4 - 4960*h*q**3*r - 6800*h*q**2*r**2 - 3000*h*q*r**3 -
500*h*q*r**2 - 625*h*r**3 - 800*i*q**3*r - 2400*i*q**2*r**2 -
1500*i*q*r**3 + 256*q**5 + 512*q**4*r + 240*q**3*r**2 + 80*q**3*r +
1100*q**2*r**2 + 1750*q*r**3 + 625*r**4 + (-sqrt(64*q**3 + 128*q**2*r
+ 60*q*r**2 + 25*r**2)*(160*h*q**3*r + 200*h*q**2*r**2 + 250*h*q*r**2
+ 625*h*r**3 + 160*i*q**3*r + 200*i*q**2*r**2 + 250*i*q*r**2 +
625*i*r**3 + 64*q**4 + 880*q**3*r + 2400*q**2*r**2 + 2250*q*r**3 +
625*r**4)/(5*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 +
5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 +
625*r**4)) + (-80*h*q**2 - 250*h*q*r - 250*h*r**2 - 50*i*q*r -
125*i*r**2 + 16*q**3 + 20*q**2*r)/(5*(4*q + 5*r)**2))*(-1280*q**4 -
4160*q**3*r - 4400*q**2*r**2 - 1500*q*r**3 - 500*q*r**2 -
625*r**3))/(2*(160*h*q**3*r + 200*h*q**2*r**2 + 250*h*q*r**2 +
625*h*r**3 + 160*i*q**3*r + 200*i*q**2*r**2 + 250*i*q*r**2 +
625*i*r**3 + 64*q**4 + 880*q**3*r + 2400*q**2*r**2 + 2250*q*r**3 +
625*r**4)), -sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 +
25*r**2)*(160*h*q**3*r + 200*h*q**2*r**2 + 250*h*q*r**2 + 625*h*r**3 +
160*i*q**3*r + 200*i*q**2*r**2 + 250*i*q*r**2 + 625*i*r**3 + 64*q**4 +
880*q**3*r + 2400*q**2*r**2 + 2250*q*r**3 + 625*r**4)/(5*(1024*q**5 +
4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 +
1500*q*r**4 + 1000*q*r**3 + 625*r**4)) + (-80*h*q**2 - 250*h*q*r -
250*h*r**2 - 50*i*q*r - 125*i*r**2 + 16*q**3 + 20*q**2*r)/(5*(4*q +
5*r)**2)), (-(-1280*h*q**4 - 4960*h*q**3*r - 6800*h*q**2*r**2 -
3000*h*q*r**3 - 500*h*q*r**2 - 625*h*r**3 - 800*i*q**3*r -
2400*i*q**2*r**2 - 1500*i*q*r**3 + 256*q**5 + 512*q**4*r +
240*q**3*r**2 + 80*q**3*r + 1100*q**2*r**2 + 1750*q*r**3 + 625*r**4 +
(sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2)*(160*h*q**3*r +
200*h*q**2*r**2 + 250*h*q*r**2 + 625*h*r**3 + 160*i*q**3*r +
200*i*q**2*r**2 + 250*i*q*r**2 + 625*i*r**3 + 64*q**4 + 880*q**3*r +
2400*q**2*r**2 + 2250*q*r**3 + 625*r**4)/(5*(1024*q**5 + 4608*q**4*r +
7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 +
1000*q*r**3 + 625*r**4)) + (-80*h*q**2 - 250*h*q*r - 250*h*r**2 -
50*i*q*r - 125*i*r**2 + 16*q**3 + 20*q**2*r)/(5*(4*q +
5*r)**2))*(-1280*q**4 - 4160*q**3*r - 4400*q**2*r**2 - 1500*q*r**3 -
500*q*r**2 - 625*r**3))/(2*(160*h*q**3*r + 200*h*q**2*r**2 +
250*h*q*r**2 + 625*h*r**3 + 160*i*q**3*r + 200*i*q**2*r**2 +
250*i*q*r**2 + 625*i*r**3 + 64*q**4 + 880*q**3*r + 2400*q**2*r**2 +
2250*q*r**3 + 625*r**4)), sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 +
25*r**2)*(160*h*q**3*r + 200*h*q**2*r**2 + 250*h*q*r**2 + 625*h*r**3 +
160*i*q**3*r + 200*i*q**2*r**2 + 250*i*q*r**2 + 625*i*r**3 + 64*q**4 +
880*q**3*r + 2400*q**2*r**2 + 2250*q*r**3 + 625*r**4)/(5*(1024*q**5 +
4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 +
1500*q*r**4 + 1000*q*r**3 + 625*r**4)) + (-80*h*q**2 - 250*h*q*r -
250*h*r**2 - 50*i*q*r - 125*i*r**2 + 16*q**3 + 20*q**2*r)/(5*(4*q +
5*r)**2))]
CPU times: user 3.24 s, sys: 29.9 ms, total: 3.27 s
Wall time: 3.3 s
> --
> You received this message because you are subscribed to the Google Groups "sympy" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to sympy+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/2c9fd218-f30f-4e3f-9b7b-4afa3a72e5b4n%40googlegroups.com.

emanuel.c...@gmail.com

unread,
Dec 23, 2021, 3:32:59 AM12/23/21
to sympy
Le mardi 21 décembre 2021 à 16:41:06 UTC+1, smi...@gmail.com a écrit :
Are they using the same versions of SymPy?

Irrelevant : Sage's solution uses Sage's default (i. e. Maxima's) solver. BTW, Sage 9.5.beta8 has sympy 1.9
 

e1 is quadratic only if f and e2 is linear in g and quadratic in f...it is a pretty simple system but there are lots of extraneous symbols. Even paring down the symbols to the system `(F1*(F2 + F3*f + f**2), G3*(G4 + G6*f**2 + G7*g + f*(G1 + G5*g - g))` takes a long time to solve and simplify.

That’s my point… BTW : with check=False, solving with Sympy’s solver is about as fast as with Maxima’s solver ; checking is slower, simplifying is extremely slow :

>>> from sympy import symbols, solve
>>> from time import time as stime
>>> f, q, r, h, i, g = map(symbols, 'fqrhig')
>>> e1 = -4*f**2*q - 5*f**2*r - 5*f*r + 4*q**2 + 3*q*r
>>> e2 = -4*f**2*h*q - 4*f**2*i*q - 8*f*g*q - 10*f*g*r - 8*f*h*q - 20*f*h*r - 10*f*i*r + 8*f*q**2/5 - 5*g*r + 4*h*q**2 - 5*h*r + 4*i*q**2 + 4*q**2/5 + 10*q*r + 5*r**2
>>> t0 = stime() ; Sol = solve([e1, e2], [f, g], check=False, dict=True) ; t1 = stime() ; print(t1 - t0)
0.6349234580993652
>>> t0 = stime() ; Chk = [[e.subs(s).simplify() for e in (e1, e2)] for s in Sol] ; t1 = stime() ; print (t1 - t0)
64.00839900970459
>>> Chk
[[0, 0], [0, 0]]
>>> t0 = stime() ; SSol = [{u:s[u].simplify() for u in s.keys()} for s in Sol] ; t1 = stime() ; print(t1 - t0)
1.9784221649169922
>>> SSol
[{f: -(5*r + sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2))/(8*q + 10*r), g: (-(4*q + 5*r)**2*sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2)*(160*h*q**3*r + 200*h*q**2*r**2 + 250*h*q*r**2 + 625*h*r**3 + 160*i*q**3*r + 200*i*q**2*r**2 + 250*i*q*r**2 + 625*i*r**3 + 64*q**4 + 880*q**3*r + 2400*q**2*r**2 + 2250*q*r**3 + 625*r**4) + (-80*h*q**2 - 250*h*q*r - 250*h*r**2 - 50*i*q*r - 125*i*r**2 + 16*q**3 + 20*q**2*r)*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 + 625*r**4))/(5*(4*q + 5*r)**2*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 + 625*r**4))}, {f: (-5*r + sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2))/(2*(4*q + 5*r)), g: ((4*q + 5*r)**2*sqrt(64*q**3 + 128*q**2*r + 60*q*r**2 + 25*r**2)*(160*h*q**3*r + 200*h*q**2*r**2 + 250*h*q*r**2 + 625*h*r**3 + 160*i*q**3*r + 200*i*q**2*r**2 + 250*i*q*r**2 + 625*i*r**3 + 64*q**4 + 880*q**3*r + 2400*q**2*r**2 + 2250*q*r**3 + 625*r**4) + (-80*h*q**2 - 250*h*q*r - 250*h*r**2 - 50*i*q*r - 125*i*r**2 + 16*q**3 + 20*q**2*r)*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 + 625*r**4))/(5*(4*q + 5*r)**2*(1024*q**5 + 4608*q**4*r + 7680*q**3*r**2 + 5600*q**2*r**3 + 400*q**2*r**2 + 1500*q*r**4 + 1000*q*r**3 + 625*r**4))}]

HTH,

emanuel.c...@gmail.com

unread,
Dec 23, 2021, 4:20:19 AM12/23/21
to sympy

“never returns” is a bit of an over statement. About half an hour is not eternity :

>>> from time import time as stime
>>> t0 = stime() ; Sol2 = solve([e1, e2], [f, g], dict=True) ; t1 = stime() ; print(t1 - t0)
1628.6711568832397

Chris Smith

unread,
Dec 23, 2021, 9:23:17 AM12/23/21
to sympy
The following solves the system in much less time and gives a solution that has 30% operations. It's not that SymPy can't solve the system, it's about knowing when to do something different and what simplification is most apt to be useful.

def test():
    var('a:z')
    e1=(-4*q - 5*r)*(f**2 - 5*f*r/(-4*q - 5*r) + 4*q**2/(
        -4*q - 5*r) + 3*q*r/(-4*q - 5*r))
    e2=10*r*(f**2*(-4*h*q - 4*i*q)/(10*r) + f*(-4*g*q/(5*r
        ) - g - 4*h*q/(5*r) - 2*h - i + 4*q**2/(25*r)
        ) - g/2 + 2*h*q**2/(5*r) - h/2 + 2*i*q**2/(5*r
        ) + 2*q**2/(25*r)+ q + r/2)
    gis=[factor_terms(collect(i, f)) for i in solve(e2,g,
        check=False,simplify=False)]
    fis=solve(e1,f,check=False,simplify=False)
    sol=[]
    for gi in gis:
        for fi in fis:
            fi = factor_terms(fi)
            sol.append({f:fi, g:gi.subs(f, fi)})
    return sol

Issue #15441 has similar problems that have been identified. This case was added there, too.

An improvement to the system-solver was removal of connected components. The next thing that should probably be done is to recursively remove from such subsets equations that uniquely contain a variable of interest or that contain only one variable of interest. e.g. given f(a,b,c), g(a,b), h(a), solve f for c, g for b and h for a.; for a case like f(a,b,c), g(a,b),h(c), solve h for c then solve f and g for a and b.

/c

Oscar Benjamin

unread,
Dec 23, 2021, 3:00:41 PM12/23/21
to sympy
On Thu, 23 Dec 2021 at 08:33, emanuel.c...@gmail.com
<emanuel.c...@gmail.com> wrote:
>
> Le mardi 21 décembre 2021 à 16:41:06 UTC+1, smi...@gmail.com a écrit :
>>
>> e1 is quadratic only if f and e2 is linear in g and quadratic in f...it is a pretty simple system but there are lots of extraneous symbols. Even paring down the symbols to the system `(F1*(F2 + F3*f + f**2), G3*(G4 + G6*f**2 + G7*g + f*(G1 + G5*g - g))` takes a long time to solve and simplify.
>
>
> That’s my point… BTW : with check=False, solving with Sympy’s solver is about as fast as with Maxima’s solver ; checking is slower, simplifying is extremely slow :

I think that the way solve checks solutions is misguided. In a simple
example like this there is no need to check the solutions. There are
cases where it is needed to avoid spurious solutions (e.g. after
unrad) but for a simple polynomial system like this the checking is
unnecessary.

--
Oscar

Omaa Yameend

unread,
Dec 24, 2021, 10:46:59 AM12/24/21
to sympy
yes same issue
Reply all
Reply to author
Forward
0 new messages