On Tuesday, March 1, 2016 at 5:53:11 PM UTC+2, Oscar wrote: > > On 26 February 2016 at 17:15, Aaron Meurer <[email protected] > <javascript:>> 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? > >
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. -- 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 [email protected]. To post to this group, send email to [email protected]. 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/d2760dfb-323d-4e5a-be20-2f42b88df01a%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
