Bug in factor() and simplify()

[Yingchi Li]

I have found a bug in factor() and simplify()
please look at the image below

the correct answer is q2 + 0.1

[Oscar Benjamin]

Looks like a recent regression since I don't see it with 0.7.1 but I
do with master:

$ isympy
IPython console for SymPy 0.7.7.dev (Python 2.7.3-64-bit) (ground types: python)
...
In [1]: q1, q2 = symbols('q1 q2')

In [2]: b = q2*sin(q1)**2 + q2*cos(q1)**2 + 0.1*sin(q1)**2 + 0.1*cos(q1)**2

In [3]: simplify(b)
Out[3]: 0.1⋅q₂ + 0.01

It doesn't happen if you split the sum in two:

In [1]: q1, q2 = symbols('q1 q2')

In [2]: a = q2*sin(q1)**2 + q2*cos(q1)**2

In [3]: b = 0.1*sin(q1)**2 + 0.1*cos(q1)**2

In [4]: simplify(a)
Out[4]: q₂

In [5]: simplify(b)
Out[5]: 0.100000000000000

In [6]: simplify(a + b)
Out[6]: 0.1⋅q₂ + 0.01

In [7]: simplify(a + b) == simplify(a) + simplify(b)
Out[7]: False

--
Oscar

[Yingchi Li]

I found that use nsimplify() can bypass this problem

c = b = q2*sin(q1)**2 + q2*cos(q1)**2 + nsimplify(0.1)*sin(q1)**2 + nsimplify(0.1)*cos(q1)**2
simplify(c) will give result q2 + 1/10

[Oscar Benjamin]

On 25 February 2016 at 15:30, Oscar Benjamin <oscar.j....@gmail.com> wrote:
> On 25 February 2016 at 10:44, Yingchi Li <liyi...@gmail.com> wrote:
>>
>> I have found a bug in factor() and simplify()
>> please look at the image below
>>
>> the correct answer is q2 + 0.1
>
> Looks like a recent regression since I don't see it with 0.7.1 but I
> do with master:

I traced this to here:

In [1]: expr = x*sin(y)**2 + x*cos(y)**2 + 0.1*sin(y)**2 + 0.1*cos(y)**2

In [2]: p = Poly(expr, x, cos(y), sin(y), domain='RR')

In [3]: p
Out[3]: Poly(1.0*x*cos(y)**2 + 1.0*x*sin(y)**2 + 0.1*cos(y)**2 +
0.1*sin(y)**2, x, cos(y), sin(y), domain='RR')

In [4]: p.rep
Out[4]: DMP([[[1.0], [], [1.0, 0.0, 0.0]], [[0.1], [], [0.1, 0.0,
0.0]]], RR, None)

In [5]: p.rep.factor_list()
Out[5]:
(1.0,
[(DMP([[[0.1], [], [0.1, 0.0, 0.0]]], RR, None), 1),
(DMP([[[1.0]], [[0.1]]], RR, None), 1)])

I think the bug is in there somehow but I have no idea what it means.
What is p.rep and how do you interpret it?

--
Oscar

[Aaron Meurer]

On Fri, Feb 26, 2016 at 11:28 AM, Oscar Benjamin

It's the internal representation of the polynomial. It's a list of
lists of lists of the coefficients (one level deep for each variable
of the polynomial). factor_list() has a bug with floating point
coefficients, it seems. Can you open an issue, if you haven't already?

Aaron Meurer

>
> --
> Oscar
>
> --
> 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 post to this group, send email to sy...@googlegroups.com.
> Visit this group at https://groups.google.com/group/sympy.
> To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/CAHVvXxTc6KYEwEBP1GqN5eh-W0JgmsGxLD8DpR1HCAzD29%3DTEw%40mail.gmail.com.
> For more options, visit https://groups.google.com/d/optout.

[Oscar Benjamin]

On 26 February 2016 at 17:15, Aaron Meurer <asme...@gmail.com> wrote:
>> I traced this to here:
>>
>> In [1]: expr = x*sin(y)**2 + x*cos(y)**2 + 0.1*sin(y)**2 + 0.1*cos(y)**2
>>
>> In [2]: p = Poly(expr, x, cos(y), sin(y), domain='RR')
>>
>> In [3]: p
>> Out[3]: Poly(1.0*x*cos(y)**2 + 1.0*x*sin(y)**2 + 0.1*cos(y)**2 +
>> 0.1*sin(y)**2, x, cos(y), sin(y), domain='RR')
>>
>> In [4]: p.rep
>> Out[4]: DMP([[[1.0], [], [1.0, 0.0, 0.0]], [[0.1], [], [0.1, 0.0,
>> 0.0]]], RR, None)
>>
>> In [5]: p.rep.factor_list()
>> Out[5]:
>> (1.0,
>> [(DMP([[[0.1], [], [0.1, 0.0, 0.0]]], RR, None), 1),
>> (DMP([[[1.0]], [[0.1]]], RR, None), 1)])
>>
>> I think the bug is in there somehow but I have no idea what it means.
>> What is p.rep and how do you interpret it?
>
> It's the internal representation of the polynomial. It's a list of
> lists of lists of the coefficients (one level deep for each variable
> of the polynomial). factor_list() has a bug with floating point
> coefficients, it seems. Can you open an issue, if you haven't already?

I will do soon as I've almost fixed it. Actually I've found the bug in
dmp_factor_list in sympy/polys/factortools.py.

A simple example is

In [1]: x**2 - 0.01*y**2
Out[1]:
2 2
x - 0.01⋅y

In [2]: factor(x**2 - 0.01*y**2)
Out[2]: 1.0⋅(0.1⋅x - 0.01⋅y)⋅(0.1⋅x + 0.01⋅y)

In [3]: expand(factor(x**2 - 0.01*y**2))
Out[3]:
2 2
0.01⋅x - 0.0001⋅y

The bug is due to this:
https://github.com/sympy/sympy/blob/master/sympy/polys/factortools.py#L1272

The first pulls out a factor of 1/100 from the polynomial so that we get

(1/100) and (100*x**2 - y**2)

the second part is factored to

(10*x + y)*(10*x - y)

The factor of 100 is brought back in here:
https://github.com/sympy/sympy/blob/master/sympy/polys/factortools.py#L1304

However it is divided from BOTH factors to give

(10*x + y)/100 * (10*x - y)/100

So the result is 100x smaller than it should be. This will happen any
time you have non-integer valued floating point coefficients and there
are at least 2 polynomial factors in the expression.

I have an easy fix in dmp_factor_list:

$ git diff master
diff --git a/sympy/polys/factortools.py b/sympy/polys/factortools.py
index ebbabb8..ccecda6 100644
--- a/sympy/polys/factortools.py
+++ b/sympy/polys/factortools.py
@@ -1302,7 +1302,8 @@ def dmp_factor_list(f, u, K0):
coeff = coeff/denom
else:
for i, (f, k) in enumerate(factors):
- f = dmp_quo_ground(f, denom, u, K0)
+ if i == 0:
+ f = dmp_quo_ground(f, denom, u, K0)
f = dmp_convert(f, u, K0, K0_inexact)
factors[i] = (f, k)

I think I can see the same bug in dup_factor_list:

https://github.com/sympy/sympy/blob/master/sympy/polys/factortools.py#L1197
https://github.com/sympy/sympy/blob/master/sympy/polys/factortools.py#L1225

But I don't know how to trigger a code path that would go through
dup_factor_list. Under what circumstances would I hit that path?

--
Oscar

[Kalevi Suominen]

dup_factor_list is called automatically from dmp_factor_list for univariate polynomials on this line:

https://github.com/sympy/sympy/blob/master/sympy/polys/factortools.py#L1253

The parameter 'u' is the number of variables minus one.