Unexpected result after polynomial zero-vector multiplication

[0xsamgreen]

Hi,

When I have two or more zero polynomials in a vector, I am not able to simplify the answer. This is causing me problems downstream in my logic and it makes printing difficult. Here is an example:

import sympy as sym

x = sym.Symbol('x')

b1 = sym.Poly(3*x**2 + x, x)
b2 = sym.Poly(2*x**3 + x**2 + 4, x)
b = sym.Matrix([b1,b2])

c1 = sym.Poly(5*x, x)
c2 = sym.Poly(x**7+x**5+3, x)
c = sym.Matrix([c1,c2])

# Calculate b dot c. Works as expected.
d = b.T*c
print('d orig',d)
# Does not simplify automatically. Simplify it manually.
d = d.applyfunc(lambda p: p.collect(x))
print('d simp',d)

# Create zero vector
e1 = sym.Poly(0, x)
e2 = sym.Poly(0, x)
e = sym.Matrix([e1,e2])

# Does not work as expected
f = (b.T*e) + sym.Matrix([b1])
print('\nf orig',f)
f = f.applyfunc(lambda p: p.collect(x))
print('f simp',f)

# Expected
g = sym.Matrix([b[0]*e[0] + b[1]*e[1] + b1])
print('\ng expe',g)

This prints:

d orig Matrix([[Poly(15*x**3 + 5*x**2, x, domain='ZZ') + Poly(2*x**10 + x**9 + 2*x**8 + 5*x**7 + 4*x**5 + 6*x**3 + 3*x**2 + 12, x, domain='ZZ')]])
d simp Matrix([[Poly(2*x**10 + x**9 + 2*x**8 + 5*x**7 + 4*x**5 + 21*x**3 + 8*x**2 + 12, x, domain='ZZ')]])

f orig Matrix([[Poly(3*x**2 + x, x, domain='ZZ') + 2*Poly(0, x, domain='ZZ')]])
f simp Matrix([[Poly(3*x**2 + x, x, domain='ZZ') + 2*Poly(0, x, domain='ZZ')]])

g expe Matrix([[Poly(3*x**2 + x, x, domain='ZZ')]])

Is this a bug? I would like the 2*Poly(0, x, domain='ZZ') to be removed from the result.

Thanks!

[sylee957]

Now it looks like there is an another regression for printing polynomial addition Add(Poly(2*x**10 + x**9 + 2*x**8 + 5*x**7 + 4*x**5 + 6*x**3 + 3*x**2 + 12, x, domain='ZZ'), Poly(15*x**3 + 5*x**2, x, domain='ZZ'))

[oscarbenjamin]

I don't think a poly should be in an Add. The arguments to Add should be Expr and Poly is not Expr (any more). https://github.com/sympy/sympy/blob/6de346061821d2bc5b7e11be0382baaf2b22083d/sympy/core/operations.py#L34-L37

I also don't think that a Poly should be in a Matrix for similar reasons. We can have a class for matrices of polynomials. In fact there is one. Using that in the example above and removing collect:

import sympy as sym
from sympy.polys.polymatrix import PolyMatrix

x = sym.Symbol('x')

b1 = sym.Poly(3*x**2 + x, x)
b2 = sym.Poly(2*x**3 + x**2 + 4, x)
b = PolyMatrix([b1,b2])

c1 = sym.Poly(5*x, x)
c2 = sym.Poly(x**7+x**5+3, x)
c = PolyMatrix([c1,c2])

# Calculate b dot c. Works as expected.
d = b.T*c
print('d orig',d)
# Does not simplify automatically. Simplify it manually.
#d = d.applyfunc(lambda p: p.collect(x))
print('d simp',d)

# Create zero vector
e1 = sym.Poly(0, x)
e2 = sym.Poly(0, x)
e = PolyMatrix([e1,e2])

# Does not work as expected
f = (b.T*e) + PolyMatrix([b1])
print('\nf orig',f)
#f = f.applyfunc(lambda p: p.collect(x))
print('f simp',f)

# Expected
g = PolyMatrix([b[0]*e[0] + b[1]*e[1] + b1])
print('\ng expe',g)

This outputs:

d orig Matrix([[Poly(2*x**10 + x**9 + 2*x**8 + 5*x**7 + 4*x**5 + 21*x**3 + 8*x**2 + 12, x, domain='ZZ')]])
d simp Matrix([[Poly(2*x**10 + x**9 + 2*x**8 + 5*x**7 + 4*x**5 + 21*x**3 + 8*x**2 + 12, x, domain='ZZ')]])

f orig Matrix([[Poly(3*x**2 + x, x, domain='ZZ')]])
f simp Matrix([[Poly(3*x**2 + x, x, domain='ZZ')]])

g expe Matrix([[Poly(3*x**2 + x, x, domain='ZZ')]])

[sylee957]

I don't think any error should be raised on printing side. But what's the reason to disallow Poly in Add or Mul, even if matrix objects that cannot be handled easily, are supported there with constructor postprocessor.

[oscarbenjamin]

The printing error is this:

>>> from sympy import *
>>> x = Symbol('x')
>>> Add(Poly(2*x**10 + x**9 + 2*x**8 + 5*x**7 + 4*x**5 + 6*x**3 + 3*x**2 + 12, x, domain='ZZ'), Poly(15*x**3 + 5*x**2, x, domain='ZZ'))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/enojb/current/sympy/sympy/sympy/core/basic.py", line 422, in __repr__
    return sstr(self, order=None)
  File "/Users/enojb/current/sympy/sympy/sympy/printing/str.py", line 909, in sstr
    s = p.doprint(expr)
  File "/Users/enojb/current/sympy/sympy/sympy/printing/printer.py", line 251, in doprint
    return self._str(self._print(expr))
  File "/Users/enojb/current/sympy/sympy/sympy/printing/printer.py", line 289, in _print
    return getattr(self, printmethod)(expr, **kwargs)
  File "/Users/enojb/current/sympy/sympy/sympy/printing/str.py", line 51, in _print_Add
    terms = self._as_ordered_terms(expr, order=order)
  File "/Users/enojb/current/sympy/sympy/sympy/printing/printer.py", line 308, in _as_ordered_terms
    return expr.as_ordered_terms(order=order)
  File "/Users/enojb/current/sympy/sympy/sympy/core/expr.py", line 1147, in as_ordered_terms
    terms, gens = self.as_terms()
  File "/Users/enojb/current/sympy/sympy/sympy/core/expr.py", line 1177, in as_terms
    coeff, _term = term.as_coeff_Mul()
AttributeError: 'Poly' object has no attribute 'as_coeff_Mul'

The cause of the error is precisely the fact that the args of the Add are not of type Expr and therefore don't have Expr methods.

If you want to add Polys you can add them directly with +. If you want to put them in an Add then use the as_expr methods:

>>> Add(Poly(2*x**10 + x**9 + 2*x**8 + 5*x**7 + 4*x**5 + 6*x**3 + 3*x**2 + 12, x, domain='ZZ').as_expr(), Poly(15*x**3 + 5*x**2, x, domain='ZZ').as_expr())
2*x**10 + x**9 + 2*x**8 + 5*x**7 + 4*x**5 + 21*x**3 + 8*x**2 + 12

It is not necessary for every object that can be added with + to be allowed to be in the args of an Add. It is also not sustainable to write code for something like Add without having any constraints on what the args should be.

SymPy is full of bugs like this where the wrong type of object is used in the wrong place. The combination of Python not having any inbuilt type checking and the flexible architecture of Basic and the printing code means that we end up with all kinds of nonsensical objects that look reasonable but don't actually work.

Until we start actually limiting the scope of things like Add we will be constantly plagued by bugs where the wrong kind of object is used in the wrong place:

>>> obj = Add((2, 3), {2}, S.true, Poly(x), 1.0, {1:2})
>>> repr(obj)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/enojb/current/sympy/sympy/sympy/core/basic.py", line 422, in __repr__
    return sstr(self, order=None)
  File "/Users/enojb/current/sympy/sympy/sympy/printing/str.py", line 909, in sstr
    s = p.doprint(expr)
  File "/Users/enojb/current/sympy/sympy/sympy/printing/printer.py", line 251, in doprint
    return self._str(self._print(expr))
  File "/Users/enojb/current/sympy/sympy/sympy/printing/printer.py", line 289, in _print
    return getattr(self, printmethod)(expr, **kwargs)
  File "/Users/enojb/current/sympy/sympy/sympy/printing/str.py", line 51, in _print_Add
    terms = self._as_ordered_terms(expr, order=order)
  File "/Users/enojb/current/sympy/sympy/sympy/printing/printer.py", line 308, in _as_ordered_terms
    return expr.as_ordered_terms(order=order)
  File "/Users/enojb/current/sympy/sympy/sympy/core/expr.py", line 1147, in as_ordered_terms
    terms, gens = self.as_terms()
  File "/Users/enojb/current/sympy/sympy/sympy/core/expr.py", line 1177, in as_terms
    coeff, _term = term.as_coeff_Mul()
AttributeError: 'BooleanTrue' object has no attribute 'as_coeff_Mul'

[0xsamgreen]

Thank you @oscarbenjamin and @sylee957. Using PolyMatrix solves the collect issue, but it presents another issue. My polynomials are over a field. When I use sym.Matrix, the product is also over the field. However, when I use PolyMatrix, the modulus is dropped and the product is over the integers; this causes future calculations to also be over the integers. Manual calculation of the dot product works as expected.

import sympy as sym
from sympy.polys.polymatrix import PolyMatrix

x = sym.Symbol('x')

field_size = 3

# Create B and C vectors with elements in GF(3)
b1 = sym.Poly(3*x**2 + x, x, domain=sym.GF(field_size, symmetric=False))
b2 = sym.Poly(2*x**3 + x**2 + 4, x, domain=sym.GF(field_size, symmetric=False))
B_M = sym.Matrix([b1,b2])
B_PM = PolyMatrix([b1,b2])

c1 = sym.Poly(5*x, x, domain=sym.GF(field_size, symmetric=False))
c2 = sym.Poly(x**7+x**5+3, x, domain=sym.GF(field_size, symmetric=False))
C_M = sym.Matrix([c1,c2])
C_PM = PolyMatrix([c1,c2])

# Calculate B_M dot C_M. Modulus works as expected, but no simplification.
d1 = B_M.T*C_M
print('d1:',d1)
# Add the dot product to itself. Result still in GF(3)? Not sure because
#   coefficients outside the polynomials
print('d1+d1:',d1+d1)

# Calculate B_PM dot C_PM. Result simplified, but we leave GF(3).
d2 = B_PM.T*C_PM
print('\nd2:',d2)
# Add the dot product to itself. Result no longer in GF(3).
print('d2+d2',d2+d2)

# Manually calculate dot product. Modulus works as expected and simplified.
d3 = b1*c1 + b2*c2
print('\nd3:',d3)
# Add the dot product to itself. Result still in GF(3).
print('d3+d3:',d3+d3)

This outputs:

d1: Matrix([[Poly(2*x**2, x, modulus=3) + Poly(2*x**10 + x**9 + 2*x**8 + 2*x**7 + x**5, x, modulus=3)]])
d1+d1: Matrix([[2*Poly(2*x**2, x, modulus=3) + 2*Poly(2*x**10 + x**9 + 2*x**8 + 2*x**7 + x**5, x, modulus=3)]])

d2: Matrix([[Poly(2*x**10 + x**9 + 2*x**8 + 2*x**7 + x**5 + 2*x**2, x, domain='ZZ')]])
d2+d2 Matrix([[Poly(4*x**10 + 2*x**9 + 4*x**8 + 4*x**7 + 2*x**5 + 4*x**2, x, domain='ZZ')]])

d3: Poly(2*x**10 + x**9 + 2*x**8 + 2*x**7 + x**5 + 2*x**2, x, modulus=3)
d3+d3: Poly(x**10 + 2*x**9 + x**8 + x**7 + 2*x**5 + x**2, x, modulus=3)

[oscarbenjamin]

This is caused by the use of sum and can be fixed in one way with this:

diff --git a/sympy/polys/polymatrix.py b/sympy/polys/polymatrix.py
index 3cc71358d2..10f817b128 100644
--- a/sympy/polys/polymatrix.py
+++ b/sympy/polys/polymatrix.py
@@ -1,5 +1,7 @@
 from __future__ import print_function

+from functools import reduce
+
 from sympy.matrices.dense import MutableDenseMatrix
 from sympy.polys.polytools import Poly

@@ -73,7 +75,7 @@ def _eval_matrix_mul(self, other):
                 col_indices = range(col, other_len, other_cols)
                 vec = (mat[a]*other_mat[b] for a,b in zip(row_indices, col_indices))
                 # 'Add' shouldn't be used here
-                new_mat[i] = sum(vec)
+                new_mat[i] = reduce(lambda a, b: a + b, vec)

         return self.__class__(new_mat_rows, new_mat_cols, new_mat, copy=False)

The problem is that sum adds with zero and:

In [1]: Poly(x, x, modulus=3)                                                                                                     
Out[1]: Poly(x, x, modulus=3)

In [2]: Poly(x, x, modulus=3) + 0                                                                                                 
Out[2]: Poly(x, x, domain='ZZ')

Alternatively maybe adding zero to the poly should not change the domain.

[0xsamgreen]

Thanks @oscarbenjamin, that fixes everything. Out of curiosity, will you be doing a pull request on this change?

[oscarbenjamin]

I'm reopening this. The point of the issue is to stay open until everything is fixed in sympy. I left the code there so that it would be easy for someone else to open a PR.

The outstanding question is whether to fix this in PolyMatrix with the diff shown or to fix it in Poly by changing the behaviour of Poly(x, x, modulus=3) + 0.