New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

factor gives an invalid expression #10161

Closed
asmeurer opened this Issue Nov 19, 2015 · 22 comments

Comments

Projects
None yet
5 participants
@asmeurer
Copy link
Member

asmeurer commented Nov 19, 2015

This was reported on the mailing list a while back

In [29]: h = 2*x*(-2*x + Abs(x))*(x**2 - 1)/Abs(x**2 - 1) + (x/Abs(x) - 2)*Abs(x**2 - 1)

In [30]: print(factor(h))
(x - 1)*(x**4 - 6*x**3*Abs(x) + x**3 - 6*x**2*Abs(x) + x**2 + 2*x*Abs(x) - x + 2*Abs(x))/(Abs(x)*Abs(x**2 - 1))

In [31]: h.subs(x, 2)
Out[31]: -11

In [32]: factor(h).subs(x, 2)
Out[32]: -53/3
@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 19, 2015

Oh, it's very important that x = Symbol('x', real=True).

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 19, 2015

Probably related:

In [24]: together(h).as_numer_denom()[0].expand()
Out[24]:
 5      4          3      2  3      2            2
x  - 6⋅x ⋅│x│ - 2⋅x  + 2⋅x ⋅x  + 8⋅x ⋅│x│ - 2⋅x⋅x  + x - 2⋅│x│

Notice how the x doesn't combine correctly in some of the terms.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 19, 2015

It used to work. I bisected it to

commit 246f52e
Author: Aaron Meurer asmeurer@gmail.com
Date: Mon Jul 23 00:20:26 2012 -0600

Refactor expand() (again)

This is a continuation of the work that was started in
d166080d296249281289ec957ad2ef887dc43ed4.  This commit goes even further than
that one.

expand() now completely handles all the deep logic (i.e., recursion).
_eval_expand_hint() methods now are responsible only for expanding the
top-level expression.  This simplifies the code by a lot, removing duplicate
"if deep: ..." logic and putting it into expand() itself.  Old methods that
still do this should still work, albeit less efficiently, with the exception
that deep is no longer included as a keyword argument to _eval_expand_hint()
methods, so calls like _eval_expand_hint(deep, **hints) (as opposed to
_eval_expand_hint(deep=deep, **hints)) may fail.

Some _eval_expand_* methods had to be refactored a little to make things
still work.

Expr no longer defines no-op _eval_expand_* methods, meaning that only
relevant subclasses of Expr will have _eval_expand_hint defined for any hint.
Therefore, any calls to an _eval_expand_* method outside of an _eval_expand_*
method should be replaced with a call to the public expand() or expand_hint()
functions.

Like its predecessor, this commit does not yet touch expand_complex(), as
there are difficulties there with as_real_imag().

The updated expand() docstring explains more, especially the API section.
@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 19, 2015

I guess this is the real problem

In [1]: x = symbols('x', real=True)

In [2]: x*abs(x)*abs(x)
Out[2]:
   2
x⋅x

which doesn't look like it ever worked correctly.

@certik

This comment has been minimized.

Copy link
Member

certik commented Nov 19, 2015

Nice work figuring out what went wrong!

@smichr

This comment has been minimized.

Copy link
Member

smichr commented Nov 20, 2015

This is one way to fix it:

diff --git a/sympy/core/mul.py b/sympy/core/mul.py
index 2a0e609..6588b7f 100644
--- a/sympy/core/mul.py
+++ b/sympy/core/mul.py
@@ -368,17 +368,27 @@ def flatten(cls, seq):

         # gather exponents of common bases...
         def _gather(c_powers):
-            new_c_powers = []
-            common_b = {}  # b:e
-            for b, e in c_powers:
-                co = e.as_coeff_Mul()
-                common_b.setdefault(b, {}).setdefault(co[1], []).append(co[0])
-            for b, d in common_b.items():
-                for di, li in d.items():
-                    d[di] = Add(*li)
-            for b, e in common_b.items():
-                for t, c in e.items():
-                    new_c_powers.append((b, c*t))
+            while True:
+                redo = False
+                new_c_powers = []
+                common_b = {}  # b:e
+                for b, e in c_powers:
+                    co = e.as_coeff_Mul()
+                    common_b.setdefault(b, {}).setdefault(co[1], []).append(co[0])
+                for b, d in common_b.items():
+                    for di, li in d.items():
+                        d[di] = Add(*li)
+                for b, e in common_b.items():
+                    for t, c in e.items():
+                        bnew, enew = Pow(b, c*t).as_base_exp()
+                        if bnew != b:
+                            redo = True
+                            b, c, t = bnew, enew, 1
+                        new_c_powers.append((b, c*t))
+                if redo:
+                    c_powers = new_c_powers
+                else:
+                    break
             return new_c_powers

This gives

>>> from sympy import *
>>> var('x',real=1)
x
>>> h = 2*x*(-2*x + Abs(x))*(x**2 - 1)/Abs(x**2 - 1) + (x/Abs(x) - 2)*Abs(x**2 - 1)
>>> print(factor(h))
(x - 1)*(x + 1)*(x - 2*Abs(x))*(3*x**2 - 1)/(Abs(x)*Abs(x**2 - 1))
>>> h.subs(x, 2)
-11
>>> factor(h).subs(x, 2)
-11
@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 20, 2015

I'm a little concerned how that would affect the performance of Mul.flatten (which is a very core routine). What is the fundamental cause of the issue?

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 20, 2015

What is the point of bnew, enew = Pow(b, c*t).as_base_exp()? I assume you're intentionally trying to get the Pow constructor to do some simplification, since otherwise bnew, enew would just be b, c*t.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 20, 2015

Oh I get it. It's because abs(x)**2 == x**2. How does this affect the performance? I think a better solution would be to remove that automatic simplification based on assumptions. That's basically already starting to be done in Pow and Abs, but I guess this particular instance happens in Abs._eval_power, and not in the Pow constructor directly.

@Shekharrajak

This comment has been minimized.

Copy link
Member

Shekharrajak commented Nov 20, 2015

It is giving correct answer on latest version ( http://live.sympy.org/ ).

>>> print(factor(h))
-(4*x**4*Abs(x) - 2*x**3*Abs(x)**2 - 4*x**2*Abs(x) + 2*x*Abs(x)**2 - x*Abs(x**2 - 1)**2 + 2*Abs(x)*Abs(x**2 - 1)**2)/(Abs(x)*Abs(x**2 - 1))

and we can check, both are not equal :

>>> e = factor(h)
>>> f = (x - 1)*(x**4 - 6*x**3*Abs(x) + x**3 - 6*x**2*Abs(x) + x**2 + 2*x*Abs(x) - x + 2*Abs(x))/(Abs(x)*Abs(x**2 - 1))
>>> bool(e == f)
False
@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 20, 2015

== will only tell you if they are structurally equal, not mathematically equal. To check that, you need to plug in different values and see if they evaluate to the same result as I have done above, or subtract the two and see if the difference simplifies to 0.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 20, 2015

@smichr did you check if your fix fixes the original problem? I'm not 100% sure that it's the full extent of the issue. The core problem is that the x*x**2, and probably also the Abs(x)**2 == x**2, is confusing the polys, leading it to make invalid assumptions on the independence of the variables (it will pick x and Abs(x) to be the variables of the expression). together and as_numer_denom work because they don't use the polys, but factor and cancel do, and they are the ones giving the wrong answers.

@Shekharrajak

This comment has been minimized.

Copy link
Member

Shekharrajak commented Nov 20, 2015

I checked there is some problem and your guesses look correct.It gives wrong answer when we do

>>> f.subs(x,2)
−53/3

>>> e.subs(x,2)
−11

but gives correct one when we do :

>>> h = 2*x*(-2*x + Abs(x))*(x**2 - 1)/Abs(x**2 - 1) + (x/Abs(x) - 2)*Abs(x**2 - 1)
>>> h.subs(x, 1/2)
−1/4

>>> factor(h).subs(x, 1/2)
−1/4

@smichr

This comment has been minimized.

Copy link
Member

smichr commented Nov 21, 2015

did you check if your fix fixes the original problem

The test expressions from the OP are in my first reply and both give the same answer. I don't know about how performance will be affected -- I cringe to think about iterating over all powers again. There is probably a better way to identify these problems before combining pairs.

@smichr

This comment has been minimized.

Copy link
Member

smichr commented Nov 21, 2015

One can see that Poly will trip over (some) un-gathered powers:

>>> Poly(Mul(2,2,x**2,evaluate=False))
Poly(4*x**2, x, domain='ZZ')
>>> Poly(Mul(x,x**2,evaluate=False))
Poly(x**2, x, domain='ZZ')
>>> Poly(Mul(x,x**2))
Poly(x**3, x, domain='ZZ')
@smichr

This comment has been minimized.

Copy link
Member

smichr commented Nov 21, 2015

The running of the core tests went from 350 (2 run ave) to 360 (2 run ave) with the changes I made. PR forthcoming.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 23, 2015

I'm also worried about the possibility of infinite loops.

@jksuom

This comment has been minimized.

Copy link
Member

jksuom commented Nov 24, 2015

There is another possibility. The critical part in factor seems to be poly_from_expr, in particular this piece:

>>> srepr(x)
"Symbol('x', real=True)"
>>> poly_from_expr(x*Abs(x)*Abs(x))
(Poly(x**2, x, domain='ZZ'), {'polys': False, 'domain': ZZ, 'gens': (x,)})

That can be fixed in polys e.g. in the following way.

diff --git a/sympy/polys/polyutils.py b/sympy/polys/polyutils.py
index 0e74704..9066640 100644
--- a/sympy/polys/polyutils.py
+++ b/sympy/polys/polyutils.py
@@ -262,7 +262,10 @@ def _is_coeff(factor):
                     else:
                         base, exp = decompose_power_rat(factor)

-                    elements[base] = exp
+                    if base in elements:
+                        elements[base] += exp
+                    else:
+                        elements[base] = exp
                     gens.add(base)

             terms.append((coeff, elements))

This might be more appropriate than a change in core, since the issue derives from the polynomial relation Abs(x)**2 - x**2 = 0.

@asmeurer

This comment has been minimized.

Copy link
Member

asmeurer commented Nov 24, 2015

The polys in general are going to choke on generators with algebraic relations on them. That does seem like a simple fix, though. Did you check the tests to see if it doesn't break anything else.

A reminder that a third possible fix would be to remove the automatic simplification of Abs(x)**2. All three of these fixes seem like they could be "right". We shouldn't consider them to be exclusive to one another.

@jksuom

This comment has been minimized.

Copy link
Member

jksuom commented Nov 25, 2015

The fix I proposed is based on the observation that _dict_from_expr currently relies on the monomials being reduced. That is not the case with x*Abs(x)*Abs(x). After Mul gets done wiith it, it is still a product x*x**2 instead of x**3. This goes into the argument of _dict_from_expr. My fix simply collects the exponents 1 + 2 instead of allowing the last one overwrite previous exponent(s).
The tests in polys pass.

(The algebraic relation is not applied here. That is quite independent.)

@smichr

This comment has been minimized.

Copy link
Member

smichr commented Nov 25, 2015

might be more appropriate than a change in core

I am in favor of your change because Poly should correctly deal with unevaluated expressions (I would think).

I believe that the core "flatten" routine should also correctly flatten args and should do the right thing if a power changes form because of autosimplification.

@jksuom

This comment has been minimized.

Copy link
Member

jksuom commented Nov 25, 2015

I agree.

@smichr smichr closed this in #10166 Dec 4, 2015

skirpichev added a commit to skirpichev/diofant that referenced this issue Jul 2, 2016

skirpichev added a commit to skirpichev/diofant that referenced this issue Jul 16, 2016

skirpichev added a commit to skirpichev/diofant that referenced this issue Jul 17, 2016

skirpichev added a commit to skirpichev/diofant that referenced this issue Jul 17, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment