Repository: sympy/sympy

Issue:
Matrices should not use is_zero for a zero matrix
Integrating a zero matrix returns the scalar zero. I think this is because matrices define is_zero as True for a zero matrix:
```julia
In [1]: M = Matrix([[0, 0]])                                                                                                      

In [2]: M                                                                                                                         
Out[2]: [0  0]

In [3]: integrate(M, t)                                                                                                           
Out[3]: 0

In [4]: M.is_zero                                                                                                                 
Out[4]: True
```
I think the error here is in Matrix.is_zero (which should always be False) rather than integrate. In the assumptions system is_zero means the scalar zero so is_zero can never be True for a matrix.

Discussion:
sympy-bot: :white_check_mark:

Hi, I am the [SymPy bot](https://github.com/sympy/sympy-bot) (v149). I'm here to help you write a release notes entry. Please read the [guide on how to write release notes](https://github.com/sympy/sympy/wiki/Writing-Release-Notes).



Your release notes are in good order.

Here is what the release notes will look like:
* matrices
  - `Matrix.is_zero` will always return `False` and the original `is_zero` is renamed to `is_zero_matrix` for assumptions consistency. ([#18464](https://github.com/sympy/sympy/pull/18464) by [@sylee957](https://github.com/sylee957))

This will be added to https://github.com/sympy/sympy/wiki/Release-Notes-for-1.6.

Note: This comment will be updated with the latest check if you edit the pull request. You need to reload the page to see it. <details><summary>Click here to see the pull request description that was parsed.</summary>

    <!-- Your title above should be a short description of what
    was changed. Do not include the issue number in the title. -->

    #### References to other Issues or PRs
    <!-- If this pull request fixes an issue, write "Fixes #NNNN" in that exact
    format, e.g. "Fixes #1234" (see
    https://tinyurl.com/auto-closing for more information). Also, please
    write a comment on that issue linking back to this pull request once it is
    open. -->

    Fixes #17131

    #### Brief description of what is fixed or changed


    #### Other comments


    #### Release Notes

    <!-- Write the release notes for this release below. See
    https://github.com/sympy/sympy/wiki/Writing-Release-Notes for more information
    on how to write release notes. The bot will check your release notes
    automatically to see if they are formatted correctly. -->

    <!-- BEGIN RELEASE NOTES -->
    - matrices
      - `Matrix.is_zero` will always return `False` and the original `is_zero` is renamed to `is_zero_matrix` for assumptions consistency.
    <!-- END RELEASE NOTES -->

</details><p>


Patch:
diff --git a/sympy/calculus/util.py b/sympy/calculus/util.py
index 7174d92c28..c81b77fe5e 100644
--- a/sympy/calculus/util.py
+++ b/sympy/calculus/util.py
@@ -1338,10 +1338,9 @@ def __pow__(self, other):
                                 return AccumBounds(0, real_root(self.max, den))
                     return AccumBounds(real_root(self.min, den),
                                        real_root(self.max, den))
-                if den!=1:
-                    num_pow = self**num
-                    return num_pow**(1 / den)
-            return AccumBounds(-oo, oo)
+                num_pow = self**num
+                return num_pow**(1 / den)
+            return Pow(self, other, evaluate=False)
 
         return NotImplemented
 
diff --git a/sympy/geometry/plane.py b/sympy/geometry/plane.py
index 1807ab7045..33a7a8871c 100644
--- a/sympy/geometry/plane.py
+++ b/sympy/geometry/plane.py
@@ -494,7 +494,7 @@ def is_parallel(self, l):
         elif isinstance(l, Plane):
             a = Matrix(l.normal_vector)
             b = Matrix(self.normal_vector)
-            if a.cross(b).is_zero:
+            if a.cross(b).is_zero_matrix:
                 return True
             else:
                 return False
@@ -527,7 +527,7 @@ def is_perpendicular(self, l):
         if isinstance(l, LinearEntity3D):
             a = Matrix(l.direction_ratio)
             b = Matrix(self.normal_vector)
-            if a.cross(b).is_zero:
+            if a.cross(b).is_zero_matrix:
                 return True
             else:
                 return False
diff --git a/sympy/holonomic/holonomic.py b/sympy/holonomic/holonomic.py
index 7594ef5fda..9ced290153 100644
--- a/sympy/holonomic/holonomic.py
+++ b/sympy/holonomic/holonomic.py
@@ -598,7 +598,7 @@ def __add__(self, other):
 
         # if a solution is not obtained then increasing the order by 1 in each
         # iteration
-        while sol.is_zero:
+        while sol.is_zero_matrix:
             dim += 1
 
             diff1 = (gen * rowsself[-1])
@@ -986,7 +986,7 @@ def __mul__(self, other):
         sol = (NewMatrix(lin_sys).transpose()).gauss_jordan_solve(homo_sys)
 
         # until a non trivial solution is found
-        while sol[0].is_zero:
+        while sol[0].is_zero_matrix:
 
             # updating the coefficients Dx^i(f).Dx^j(g) for next degree
             for i in range(a - 1, -1, -1):
@@ -1211,19 +1211,19 @@ def composition(self, expr, *args, **kwargs):
         homogeneous = Matrix([[S.Zero for i in range(a)]]).transpose()
         sol = S.Zero
 
-        while sol.is_zero:
+        while True:
             coeffs_next = [p.diff(self.x) for p in coeffs]
             for i in range(a - 1):
                 coeffs_next[i + 1] += (coeffs[i] * diff)
-
             for i in range(a):
                 coeffs_next[i] += (coeffs[-1] * subs[i] * diff)
             coeffs = coeffs_next
-
             # check for linear relations
             system.append(coeffs)
             sol, taus = (Matrix(system).transpose()
                 ).gauss_jordan_solve(homogeneous)
+            if sol.is_zero_matrix is not True:
+                break
 
         tau = list(taus)[0]
         sol = sol.subs(tau, 1)
diff --git a/sympy/holonomic/linearsolver.py b/sympy/holonomic/linearsolver.py
index 5e37e6b35e..6eb2e0c97d 100644
--- a/sympy/holonomic/linearsolver.py
+++ b/sympy/holonomic/linearsolver.py
@@ -66,7 +66,7 @@ def gauss_jordan_solve(self, b, freevar=False):
 
         # check for existence of solutions
         # rank of aug Matrix should be equal to rank of coefficient matrix
-        if not v[rank:, 0].is_zero:
+        if not v[rank:, 0].is_zero_matrix:
             raise ValueError("Linear system has no solution")
 
         # Get index of free symbols (free parameters)
diff --git a/sympy/matrices/common.py b/sympy/matrices/common.py
index 4617511626..f5753115ca 100644
--- a/sympy/matrices/common.py
+++ b/sympy/matrices/common.py
@@ -1114,7 +1114,7 @@ def _eval_is_diagonal(self):
     # sure the names don't clash by adding `_matrix_` in name.
     def _eval_is_matrix_hermitian(self, simpfunc):
         mat = self._new(self.rows, self.cols, lambda i, j: simpfunc(self[i, j] - self[j, i].conjugate()))
-        return mat.is_zero
+        return mat.is_zero_matrix
 
     def _eval_is_Identity(self) -> FuzzyBool:
         def dirac(i, j):
@@ -1141,9 +1141,9 @@ def _eval_is_symbolic(self):
 
     def _eval_is_symmetric(self, simpfunc):
         mat = self._new(self.rows, self.cols, lambda i, j: simpfunc(self[i, j] - self[j, i]))
-        return mat.is_zero
+        return mat.is_zero_matrix
 
-    def _eval_is_zero(self):
+    def _eval_is_zero_matrix(self):
         if any(i.is_zero == False for i in self):
             return False
         if any(i.is_zero is None for i in self):
@@ -1623,7 +1623,7 @@ def is_upper(self):
                    for j in range(min(i, self.cols)))
 
     @property
-    def is_zero(self):
+    def is_zero_matrix(self):
         """Checks if a matrix is a zero matrix.
 
         A matrix is zero if every element is zero.  A matrix need not be square
@@ -1641,17 +1641,17 @@ def is_zero(self):
         >>> c = Matrix([[0, 1], [0, 0]])
         >>> d = Matrix([])
         >>> e = Matrix([[x, 0], [0, 0]])
-        >>> a.is_zero
+        >>> a.is_zero_matrix
         True
-        >>> b.is_zero
+        >>> b.is_zero_matrix
         True
-        >>> c.is_zero
+        >>> c.is_zero_matrix
         False
-        >>> d.is_zero
+        >>> d.is_zero_matrix
         True
-        >>> e.is_zero
+        >>> e.is_zero_matrix
         """
-        return self._eval_is_zero()
+        return self._eval_is_zero_matrix()
 
     def values(self):
         """Return non-zero values of self."""
diff --git a/sympy/matrices/dense.py b/sympy/matrices/dense.py
index 681d6636d2..8ce9042fce 100644
--- a/sympy/matrices/dense.py
+++ b/sympy/matrices/dense.py
@@ -787,7 +787,8 @@ def zip_row_op(self, i, k, f):
 
         self._mat[i0: i0 + self.cols] = [f(x, y) for x, y in zip(ri, rk)]
 
-    # Utility functions
+    is_zero = False
+
 
 MutableMatrix = Matrix = MutableDenseMatrix
 
diff --git a/sympy/matrices/immutable.py b/sympy/matrices/immutable.py
index 841bdb86c6..6ab214f68c 100644
--- a/sympy/matrices/immutable.py
+++ b/sympy/matrices/immutable.py
@@ -2,7 +2,7 @@
 
 from typing import Callable
 
-from sympy.core import Basic, Dict, Integer, S, Tuple, sympify
+from sympy.core import Basic, Dict, Integer, S, Tuple
 from sympy.core.cache import cacheit
 from sympy.core.sympify import converter as sympify_converter
 from sympy.matrices.dense import DenseMatrix
@@ -91,8 +91,11 @@ def _eval_Eq(self, other):
         if isinstance(other, MatrixExpr) and not isinstance(
                 other, ImmutableDenseMatrix):
             return None
-        diff = self - other
-        return sympify(diff.is_zero)
+        diff = (self - other).is_zero_matrix
+        if diff is True:
+            return S.true
+        elif diff is False:
+            return S.false
 
     def _eval_extract(self, rowsList, colsList):
         # self._mat is a Tuple.  It is slightly faster to index a
@@ -122,17 +125,6 @@ def is_diagonalizable(self, reals_only=False, **kwargs):
     is_diagonalizable = cacheit(is_diagonalizable)
 
 
-# XXX: The assumptions system correctly infers that is_zero is always False
-# because commutative is False and zero is commutative. This assignment of
-# the property has to take place after the class definition so that the
-# ManagedProperties metaclass that handles the assumptions has finished
-# post-processing the class first. This should be removed because a Matrix
-# should always have is_zero=False.
-# https://github.com/sympy/sympy/issues/7213
-# https://github.com/sympy/sympy/issues/17131
-# mypy (correctly) complains about this line:
-ImmutableDenseMatrix.is_zero = DenseMatrix.is_zero  # type: ignore
-
 # make sure ImmutableDenseMatrix is aliased as ImmutableMatrix
 ImmutableMatrix = ImmutableDenseMatrix
 
diff --git a/sympy/matrices/matrices.py b/sympy/matrices/matrices.py
index dad7eebc6c..7c6c6054f1 100644
--- a/sympy/matrices/matrices.py
+++ b/sympy/matrices/matrices.py
@@ -1970,7 +1970,7 @@ def gauss_jordan_solve(self, B, freevar=False, dotprodsimp=None):
 
         # check for existence of solutions
         # rank of aug Matrix should be equal to rank of coefficient matrix
-        if not v[rank:, :].is_zero:
+        if not v[rank:, :].is_zero_matrix:
             raise ValueError("Linear system has no solution")
 
         # Get index of free symbols (free parameters)
@@ -3259,7 +3259,7 @@ def _eval_pinv_full_rank(self):
         This routine can apply for both cases by checking the shape
         and have small decision.
         """
-        if self.is_zero:
+        if self.is_zero_matrix:
             return self.H
 
         if self.rows >= self.cols:
@@ -3274,7 +3274,7 @@ def _eval_pinv_rank_decomposition(self):
         rank matrices, and each matrix can take pseudoinverse
         individually.
         """
-        if self.is_zero:
+        if self.is_zero_matrix:
             return self.H
 
         B, C = self.rank_decomposition()
@@ -3290,7 +3290,7 @@ def _eval_pinv_diagonalization(self):
         This routine can sometimes fail if SymPy's eigenvalue
         computation is not reliable.
         """
-        if self.is_zero:
+        if self.is_zero_matrix:
             return self.H
 
         A = self
@@ -3364,7 +3364,7 @@ def pinv(self, method='RD'):
 
         """
         # Trivial case: pseudoinverse of all-zero matrix is its transpose.
-        if self.is_zero:
+        if self.is_zero_matrix:
             return self.H
 
         if method == 'RD':
diff --git a/sympy/matrices/sparse.py b/sympy/matrices/sparse.py
index ee856a4047..b99ac0c278 100644
--- a/sympy/matrices/sparse.py
+++ b/sympy/matrices/sparse.py
@@ -1431,3 +1431,5 @@ def zip_row_op(self, i, k, f):
 
         """
         self.row_op(i, lambda v, j: f(v, self[k, j]))
+
+    is_zero = False
diff --git a/sympy/matrices/subspaces.py b/sympy/matrices/subspaces.py
index 204d107f9c..e78cac9d94 100644
--- a/sympy/matrices/subspaces.py
+++ b/sympy/matrices/subspaces.py
@@ -202,7 +202,7 @@ def perp_to_subspace(vec, basis):
     ret  = []
     vecs = list(vecs) # make sure we start with a non-zero vector
 
-    while len(vecs) > 0 and vecs[0].is_zero:
+    while len(vecs) > 0 and vecs[0].is_zero_matrix:
         if rankcheck is False:
             del vecs[0]
         else:
@@ -211,7 +211,7 @@ def perp_to_subspace(vec, basis):
     for vec in vecs:
         perp = perp_to_subspace(vec, ret)
 
-        if not perp.is_zero:
+        if not perp.is_zero_matrix:
             ret.append(cls(perp))
         elif rankcheck is True:
             raise ValueError("GramSchmidt: vector set not linearly independent")
diff --git a/sympy/polys/multivariate_resultants.py b/sympy/polys/multivariate_resultants.py
index c5ea3f4f65..24755aec45 100644
--- a/sympy/polys/multivariate_resultants.py
+++ b/sympy/polys/multivariate_resultants.py
@@ -194,7 +194,7 @@ def KSY_precondition(self, matrix):
         the last column. For the precondition to hold the last non-zero
         row of the rref matrix should be of the form [0, 0, ..., 1].
         """
-        if matrix.is_zero:
+        if matrix.is_zero_matrix:
             return False
 
         m, n = matrix.shape
@@ -213,8 +213,10 @@ def KSY_precondition(self, matrix):
 
     def delete_zero_rows_and_columns(self, matrix):
         """Remove the zero rows and columns of the matrix."""
-        rows = [i for i in range(matrix.rows) if not matrix.row(i).is_zero]
-        cols = [j for j in range(matrix.cols) if not matrix.col(j).is_zero]
+        rows = [
+            i for i in range(matrix.rows) if not matrix.row(i).is_zero_matrix]
+        cols = [
+            j for j in range(matrix.cols) if not matrix.col(j).is_zero_matrix]
 
         return matrix[rows, cols]
 
diff --git a/sympy/polys/polytools.py b/sympy/polys/polytools.py
index af1db9115e..1f15d8d98c 100644
--- a/sympy/polys/polytools.py
+++ b/sympy/polys/polytools.py
@@ -6593,7 +6593,7 @@ def cancel(f, *gens, **args):
     Examples
     ========
 
-    >>> from sympy import cancel, sqrt, Symbol, together
+    >>> from sympy import cancel, sqrt, Symbol
     >>> from sympy.abc import x
     >>> A = Symbol('A', commutative=False)
 
@@ -6601,14 +6601,6 @@ def cancel(f, *gens, **args):
     (2*x + 2)/(x - 1)
     >>> cancel((sqrt(3) + sqrt(15)*A)/(sqrt(2) + sqrt(10)*A))
     sqrt(6)/2
-
-    Note: due to automatic distribution of Rationals, a sum divided by an integer
-    will appear as a sum. To recover a rational form use `together` on the result:
-
-    >>> cancel(x/2 + 1)
-    x/2 + 1
-    >>> together(_)
-    (x + 2)/2
     """
     from sympy.core.exprtools import factor_terms
     from sympy.functions.elementary.piecewise import Piecewise


Test patch:
diff --git a/sympy/calculus/tests/test_util.py b/sympy/calculus/tests/test_util.py
index 5035cceb4e..a3d34a6169 100644
--- a/sympy/calculus/tests/test_util.py
+++ b/sympy/calculus/tests/test_util.py
@@ -452,14 +452,14 @@ def test_AccumBounds_pow():
     assert AccumBounds(-1, 2)**oo == AccumBounds(-oo, oo)
     assert AccumBounds(-2, S.Half)**oo == AccumBounds(-oo, oo)
 
-    assert AccumBounds(1, 2)**x == Pow(AccumBounds(1, 2), x)
+    assert AccumBounds(1, 2)**x == Pow(AccumBounds(1, 2), x, evaluate=False)
 
     assert AccumBounds(2, 3)**(-oo) is S.Zero
     assert AccumBounds(0, 2)**(-oo) == AccumBounds(0, oo)
     assert AccumBounds(-1, 2)**(-oo) == AccumBounds(-oo, oo)
 
     assert (tan(x)**sin(2*x)).subs(x, AccumBounds(0, pi/2)) == \
-        Pow(AccumBounds(-oo, oo), AccumBounds(0, 1))
+        Pow(AccumBounds(-oo, oo), AccumBounds(0, 1), evaluate=False)
 
 
 def test_comparison_AccumBounds():
diff --git a/sympy/functions/combinatorial/tests/test_comb_factorials.py b/sympy/functions/combinatorial/tests/test_comb_factorials.py
index 54f3e89a2b..09e62d05dd 100644
--- a/sympy/functions/combinatorial/tests/test_comb_factorials.py
+++ b/sympy/functions/combinatorial/tests/test_comb_factorials.py
@@ -147,9 +147,7 @@ def test_ff_eval_apply():
     for i in range(100):
         x = -500 + 500 * random.random()
         k = -500 + 500 * random.random()
-        a = mpmath_ff(x, k)
-        b = ff(x, k)
-        assert (abs(a - b) < abs(a) * 10**(-15))
+        assert (abs(mpmath_ff(x, k) - ff(x, k)) < 10**(-15))
 
 
 def test_rf_ff_eval_hiprec():
diff --git a/sympy/functions/elementary/tests/test_piecewise.py b/sympy/functions/elementary/tests/test_piecewise.py
index fd609c9f71..17b0bacf32 100644
--- a/sympy/functions/elementary/tests/test_piecewise.py
+++ b/sympy/functions/elementary/tests/test_piecewise.py
@@ -2,7 +2,7 @@
     adjoint, And, Basic, conjugate, diff, expand, Eq, Function, I, ITE,
     Integral, integrate, Interval, KroneckerDelta, lambdify, log, Max, Min,
     oo, Or, pi, Piecewise, piecewise_fold, Rational, solve, symbols, transpose,
-    cos, sin, exp, Abs, Ne, Not, Symbol, S, sqrt, Sum, Tuple, zoo, Float,
+    cos, sin, exp, Abs, Ne, Not, Symbol, S, sqrt, Sum, Tuple, zoo,
     DiracDelta, Heaviside, Add, Mul, factorial, Ge, Contains)
 from sympy.core.expr import unchanged
 from sympy.functions.elementary.piecewise import Undefined, ExprCondPair
@@ -1345,9 +1345,3 @@ def test_identical_conds_issue():
     # Result is quite big, so not really important here (and should ideally be
     # simpler). Should not give an exception though.
     density(u1 + u2)
-
-
-def test_issue_7370():
-    f = Piecewise((1, x <= 2400))
-    v = integrate(f, (x, 0, Float("252.4", 30)))
-    assert str(v) == '252.400000000000000000000000000'
diff --git a/sympy/matrices/tests/test_commonmatrix.py b/sympy/matrices/tests/test_commonmatrix.py
index 5ae0cfaa39..d29af83226 100644
--- a/sympy/matrices/tests/test_commonmatrix.py
+++ b/sympy/matrices/tests/test_commonmatrix.py
@@ -440,14 +440,14 @@ def test_is_hessenberg():
 
 
 def test_is_zero():
-    assert PropertiesOnlyMatrix(0, 0, []).is_zero
-    assert PropertiesOnlyMatrix([[0, 0], [0, 0]]).is_zero
-    assert PropertiesOnlyMatrix(zeros(3, 4)).is_zero
-    assert not PropertiesOnlyMatrix(eye(3)).is_zero
-    assert PropertiesOnlyMatrix([[x, 0], [0, 0]]).is_zero == None
-    assert PropertiesOnlyMatrix([[x, 1], [0, 0]]).is_zero == False
+    assert PropertiesOnlyMatrix(0, 0, []).is_zero_matrix
+    assert PropertiesOnlyMatrix([[0, 0], [0, 0]]).is_zero_matrix
+    assert PropertiesOnlyMatrix(zeros(3, 4)).is_zero_matrix
+    assert not PropertiesOnlyMatrix(eye(3)).is_zero_matrix
+    assert PropertiesOnlyMatrix([[x, 0], [0, 0]]).is_zero_matrix == None
+    assert PropertiesOnlyMatrix([[x, 1], [0, 0]]).is_zero_matrix == False
     a = Symbol('a', nonzero=True)
-    assert PropertiesOnlyMatrix([[a, 0], [0, 0]]).is_zero == False
+    assert PropertiesOnlyMatrix([[a, 0], [0, 0]]).is_zero_matrix == False
 
 
 def test_values():
diff --git a/sympy/matrices/tests/test_matrices.py b/sympy/matrices/tests/test_matrices.py
index a0800e5cef..cce489d755 100644
--- a/sympy/matrices/tests/test_matrices.py
+++ b/sympy/matrices/tests/test_matrices.py
@@ -3083,17 +3083,17 @@ def test_zeros_eye():
 
 
 def test_is_zero():
-    assert Matrix().is_zero
-    assert Matrix([[0, 0], [0, 0]]).is_zero
-    assert zeros(3, 4).is_zero
-    assert not eye(3).is_zero
-    assert Matrix([[x, 0], [0, 0]]).is_zero == None
-    assert SparseMatrix([[x, 0], [0, 0]]).is_zero == None
-    assert ImmutableMatrix([[x, 0], [0, 0]]).is_zero == None
-    assert ImmutableSparseMatrix([[x, 0], [0, 0]]).is_zero == None
-    assert Matrix([[x, 1], [0, 0]]).is_zero == False
+    assert Matrix().is_zero_matrix
+    assert Matrix([[0, 0], [0, 0]]).is_zero_matrix
+    assert zeros(3, 4).is_zero_matrix
+    assert not eye(3).is_zero_matrix
+    assert Matrix([[x, 0], [0, 0]]).is_zero_matrix == None
+    assert SparseMatrix([[x, 0], [0, 0]]).is_zero_matrix == None
+    assert ImmutableMatrix([[x, 0], [0, 0]]).is_zero_matrix == None
+    assert ImmutableSparseMatrix([[x, 0], [0, 0]]).is_zero_matrix == None
+    assert Matrix([[x, 1], [0, 0]]).is_zero_matrix == False
     a = Symbol('a', nonzero=True)
-    assert Matrix([[a, 0], [0, 0]]).is_zero == False
+    assert Matrix([[a, 0], [0, 0]]).is_zero_matrix == False
 
 
 def test_rotation_matrices():
diff --git a/sympy/series/tests/test_limits.py b/sympy/series/tests/test_limits.py
index a001378e8a..93eb22b3d8 100644
--- a/sympy/series/tests/test_limits.py
+++ b/sympy/series/tests/test_limits.py
@@ -608,6 +608,3 @@ def test_issue_17671():
 
 def test_issue_18306():
     assert limit(sin(sqrt(x))/sqrt(sin(x)), x, 0, '+') == 1
-
-def test_issue_18442():
-    assert limit(tan(x)**(2**(sqrt(pi))), x, oo, dir='-') == AccumBounds(-oo, oo)