search

sympy / sympy

A computer algebra system written in pure Python
https://sympy.org/
Other
12.85k stars 4.41k forks source link

radsimp() makes expressions more complex. #3557

Closed certik closed 12 years ago

certik commented 16 years ago

bc.. a,b,c,d,e,f,g,h,i,x,y,z = symbols('abcdefghixyz')

eq1 = x_a + y_b + z_c == 1 eq2 = x_d + y_e + z_f == 1 eq3 = x_g + y_h + z_i == 1 solve([eq1,eq2,eq3],[x,y,z]) Traceback (most recent call last): File "", line 1, in File "/usr/lib/python2.5/site-packages/sympy/solvers/solvers.py", line 112, in solve return solve_linear_system(matrix, syms, simplified) File "/usr/lib/python2.5/site-packages/sympy/solvers/solvers.py", line 205, in solve_linearsystem solutions[syms[k]] = simplify(content) File "/usr/lib/python2.5/site-packages/sympy/simplify/simplify.py", line 1007, in simplify return q + factor(r) / factor(d) File "/usr/lib/python2.5/site-packages/sympy/polynomials/ wrapper.py", line 38, in factor return Mul([ff.sympyexpr for ff in factor.factor(f, var, order)]) File "/usr/lib/python2.5/site-packages/sympy/polynomials/ factor_.py", line 194, in factor factors = kroneckermv(f) File "/usr/lib/python2.5/site-packages/sympy/polynomials/ factor.py", line 336, in kronecker_mv gfactors = factor(g)[1:] # Don't use constant factor. File "/usr/lib/python2.5/site-packages/sympy/polynomials/ factor.py", line 189, in factor factors = fast.intpoly.factor(Polynomial2IntPoly(f)) File "/usr/lib/python2.5/site-packages/sympy/polynomials/fast/ intpoly.py", line 347, in factor sqf_part = squarefree_part(pp) File "/usr/lib/python2.5/site-packages/sympy/polynomials/fast/ intpoly.py", line 334, in squarefree_part g = gcd(f, f.diff()) File "/usr/lib/python2.5/site-packages/sympy/polynomials/fast/ intpoly.py", line 55, in gcd_small_primes B = int(math.ceil(2*_n_A_b_math.sqrt(n+1))) OverflowError: long int too large to convert to float

Maxima:

eq1 : x_a + y_b + z_c = 1; eq2 : x_d + y_e + z_f = 1; eq3 : x_g + y_h + z*i = 1;

(%i18) solve([eq1,eq2,eq3],[x,y,z]);

              b (i - f) - e i + f h + c (e - h)

(%o18) [[x = ---------------------------------------------, a (f h - e i) + b (d i - f g) + c (e g - d h) a (i - f) - d i + f g + c (d - g) y = - ---------------------------------------------, a (f h - e i) + b (d i - f g) + c (e g - d h) a (h - e) - d h + e g + b (d - g) z = ---------------------------------------------]] a (f h - e i) + b (d i - f g) + c (e g - d h)

p. Original issue for "#3557":https://github.com/sympy/sympy/issues/3557: "http://code.google.com/p/sympy/issues/detail?id=458":http://code.google.com/p/sympy/issues/detail?id=458

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

sympy-issue-migrator commented 16 years ago

bc.. The point of failure here ist the attempt of multivariate factorization.

Multivariate factorization is not at all robust, because it just tries to solve the problem by factorization of a new univariate polynomial (of high degree, depending on the number of variables). This is probably, were the overflow is generated.

Factorization is called from simplify in the solver. I guess that the solving wouldn't pose any great problems though. Ideally, the solver would recognize this as linear equations and call the appropriate matrix algorithm, right?

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c1":http://code.google.com/p/sympy/issues/detail?id=458#c1

p. Original author: "https://code.google.com/u/116277859734487035258/":https://code.google.com/u/116277859734487035258/

certik commented 16 years ago

bc.. Yes, in this case it should be converted to matrix algorithm automatically by the solver. So it should be actually pretty easy to fix this bug.

The only problem is how to easily automatically recognize that we deal with a linear system.

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c2":http://code.google.com/p/sympy/issues/detail?id=458#c2

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

certik commented 16 years ago

bc.. Labels: -Priority-Medium Priority-High

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c3":http://code.google.com/p/sympy/issues/detail?id=458#c3

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

certik commented 16 years ago

bc.. You can apply the following patch:

$ hg di diff -r 48db9d4aaf29 sympy/solvers/solvers.py --- a/sympy/solvers/solvers.py Fri Dec 14 09:14:00 2007 +0300 +++ b/sympy/solvers/solvers.py Fri Dec 14 22:19:10 2007 +0100 @@ -202,7 +202,7 @@ def solve_linear_system(system, symbols, content -= matrix[k, j]*solutions[syms[j]]

         if simplified == True:

That doesn't use "simplify", but just expand() and then it works, even though the result is complicated:

$ bin/isympy Python 2.4.4 console for SymPy 0.5.8-hg. These commands were executed:

from future import division from sympy import * x, y, z = symbols('xyz') k, m, n = symbols('kmn', integer=True) f = Function("f") Basic.set_repr_level(2) # pretty print output; Use "1" for python output pprint_try_use_unicode() # use unicode pretty print when available

In [1]: a,b,c,d,e,f,g,h,i,x,y,z = symbols('abcdefghixyz')

In [2]: eq1 = x_a + y_b + z*c == 1

In [3]: eq2 = x_d + y_e + z*f == 1

In [4]: eq3 = x_g + y_h + z*i == 1

In [5]: print solve([eq1,eq2,eq3],[x,y,z]) {y:
1 d f
─────── - ─────────── - ────────────────────────────────────────────────────── b_d ⎛ b_d⎞ ⎛ b_d⎞ ⎛ c_g f_h b_f_g c_dh e - ─── a⎜e - ───⎟ ⎜e - ───⎟_⎜i - ─── - ─────── + ─────────── + ───────── a ⎝ a ⎠ ⎝ a ⎠ ⎜ a b_d ⎛ bd⎞ ⎛ b ⎜ e - ─── a⎜e - ───⎟ a⎜e - ── ⎝ a ⎝ a ⎠ ⎝ a

[...]

(and many more lines like that)

Status: Started

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c4":http://code.google.com/p/sympy/issues/detail?id=458#c4

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

certik commented 16 years ago

bc.. It's a bad idea to call simplify() in the solver. Imho. But it's clear, that it needs to try to simplify things, compare the result from sympy and from maxima.

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c5":http://code.google.com/p/sympy/issues/detail?id=458#c5

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

mattpap commented 16 years ago

bc.. With some effort sympy can compute the system in finite time:

In [1]: var('a b c d e f g h i') Out[1]: (a, b, c, d, e, f, g, h, i)

In [2]: f_1 = x_a + y_b + z*c - 1

In [3]: f_2 = x_d + y_e + z*f - 1

In [4]: f_3 = x_g + y_h + z*i - 1

In [5]: s = solve([f_1,f_2,f_3], x,y,z, simplified=False)

In [6]: s_x = Poly.cancel(s[x])

In [7]: s_y = Poly.cancel(s[y])

In [8]: s_z = Poly.cancel(s[z])

In [9]: s_x Out[9]: b_f + c_h + e_i - b_i - c_e - f_h
───────────────────────────────────────────── a_e_i + b_f_g + c_d_h - a_f_h - b_d_i - c_e_g

In [10]: s_y Out[10]: a_f + c_g + d_i - a_i - c_d - f_g
───────────────────────────────────────────── a_f_h + b_d_i + c_e_g - a_e_i - b_f_g - c_d_h

In [11]: s_z Out[11]: a_e + b_g + d_h - a_h - b_d - e_g
───────────────────────────────────────────── a_e_i + b_f_g + c_d_h - a_f_h - b_d_i - c_e_g

It's a bad idea to call simplify() in the solver.

It's bad to use radsimp() (if fact match()) anywhere:

In [12]: s_x.count_ops(False) Out[12]: 36

In [13]: s[x].count_ops(False) Out[13]: 137

In [14]: powsimp(s[x]).count_ops(False) Out[14]: 137

In [15]: radsimp(s[x]).count_ops(False) Out[15]: 1152

For now a solution is to remove radsimp() from simplify(). If fact functionality of radsimp() was reimplemented in Poly.cancel() for Symbol-less case (where radsimp() gives any reasonable results), so we won't loose anything with this change.

Summary: radsimp() makes expressions more complex.

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c6":http://code.google.com/p/sympy/issues/detail?id=458#c6

p. Original author: "https://code.google.com/u/101069955704897915480/":https://code.google.com/u/101069955704897915480/

certik commented 16 years ago

bc.. So it's radsimp redundant? Let's remove it then, or call Poly.cancel() from it.

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c7":http://code.google.com/p/sympy/issues/detail?id=458#c7

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

certik commented 16 years ago

bc.. The original problem was fixed and a test written.

In [1]: var('a b c d e f g h i') Out[1]: (a, b, c, d, e, f, g, h, i)

In [2]: f_1 = x_a + y_b + z*c - 1

In [3]: f_2 = x_d + y_e + z*f - 1

In [4]: f_3 = x_g + y_h + z*i - 1

In [5]: solve([f_1,f_2,f_3], x,y,z) Out[5]: ⎧ a_e + b_g + d_h - a_h - b_d - e_g a_f + c_g + d_i - ⎨z: ─────────────────────────────────────────────, y: ──────────────────────── ⎩ a_e_i + b_f_g + c_d_h - a_f_h - b_d_i - c_e_g a_f_h + b_d_i + c_e*g -

a_i - c_d - f_g b_f + c_h + e_i - b_i - c_e - f_h ⎫ ─────────────────────, x: ─────────────────────────────────────────────⎬ a_e_i - b_f_g - c_d_h a_e_i + b_f_g + c_d_h - a_f_h - b_d_i - c_e*g⎭

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c8":http://code.google.com/p/sympy/issues/detail?id=458#c8

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

certik commented 16 years ago

bc.. So I think this is not a high priority anymore.

Labels: -Priority-High Priority-Medium

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c9":http://code.google.com/p/sympy/issues/detail?id=458#c9

p. Original author: "https://code.google.com/u/104039945248245758823/":https://code.google.com/u/104039945248245758823/

mattpap commented 16 years ago

bc.. > So it's radsimp redundant?

radsimp() is needed, but:

[1] Basic.match must be fixed, there are many other issues concerning it [2] radsimp() must know about sqrtdenest(), to handle trivial cases fast and also non-trivial input [3] it can't be applied blindly as it was in simplify()

p. Original comment: "http://code.google.com/p/sympy/issues/detail?id=458#c10":http://code.google.com/p/sympy/issues/detail?id=458#c10

p. Original author: "https://code.google.com/u/101069955704897915480/":https://code.google.com/u/101069955704897915480/

asmeurer commented 14 years ago 
Chris, do your radsimp() changes in t (commit from issue 4924 ) fix any of this?  Or is that issue a duplicate of this one?

**Cc:** smichr

Referenced issues: #4924 Original comment: http://code.google.com/p/sympy/issues/detail?id=458#c12 Original author: https://code.google.com/u/asmeurer@gmail.com/

smichr commented 12 years ago 
c8 pointed out the the OP issue was fixed. 1825 fixes the explosion-in-size issue and a pull is forthcoming.

    >>> ok = var('a:z')
    >>> eq1 = x*a + y*b + z*c - 1
    >>> eq2 = x*d + y*e + z*f - 1
    >>> eq3 = x*g + y*h + z*i - 1
    >>> s = solve([eq1,eq2,eq3],[x,y,z])
    >>> for k,v in s.items():
    ...  print k
    ...  pprint(v)
    ...  print
    ...
    x
          b*f - b*i - c*e + c*h + e*i - f*h
    ---------------------------------------------
    a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g

    y
          -a*f + a*i + c*d - c*g - d*i + f*g
    ---------------------------------------------
    a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g

    z
          a*e - a*h - b*d + b*g + d*h - e*g
    ---------------------------------------------
    a*e*i - a*f*h - b*d*i + b*f*g + c*d*h - c*e*g

    >>> s[x].count_ops()
    29
    >>> radsimp(s[x]).count_ops()
    29

Original comment: http://code.google.com/p/sympy/issues/detail?id=458#c13 Original author: https://code.google.com/u/117933771799683895267/

smichr commented 12 years ago 
https://github.com/sympy/sympy/pull/660

**Labels:** NeedsReview smichr

Original comment: http://code.google.com/p/sympy/issues/detail?id=458#c14 Original author: https://code.google.com/u/117933771799683895267/

smichr commented 12 years ago 
Actually, the explosion of size was already fixed in master, but the other issues (which I will enumerate in the commit message) are addressed by the pull request.

Original comment: http://code.google.com/p/sympy/issues/detail?id=458#c15 Original author: https://code.google.com/u/117933771799683895267/

smichr commented 12 years ago 
**Status:** Fixed

Original comment: http://code.google.com/p/sympy/issues/detail?id=458#c16 Original author: https://code.google.com/u/117933771799683895267/