Repository: sympy/sympy
Pull request: core: fixed _eval_is_zero() and added test
<!-- 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
Fixes #17130
<!-- 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. -->


#### Brief description of what is fixed or changed
_eval_is_zero() would give wrong result with imaginary numbers,
resolved as directed by @oscarbenjamin, and added
appropriate test.

#### Other comments
Please suggest any extra tests I should add to it.

#### 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 -->
* core
  * Fixed `_eval_is_zero()` functionality when imaginary numbers are involved.
<!-- END RELEASE NOTES -->

Problem statement:
IndexError from jordan_form
I get an error computing the Jordan form of this Matrix:
```julia
In [6]: M = Matrix([ 
   ...: [exp(t**3/3 - t**2/2)/2 + exp(t**3/3 + t**2/2)/2, exp(t**3/3 - t**2/2)/2 - exp(t**3/3 + t**2/2)/2], 
   ...: [exp(t**3/3 - t**2/2)/2 - exp(t**3/3 + t**2/2)/2, exp(t**3/3 - t**2/2)/2 + exp(t**3/3 + t**2/2)/2]])                      

In [7]: M                                                                                                                         
Out[7]: 
⎡  3    2     3    2    3    2     3    2⎤
⎢ t    t     t    t    t    t     t    t ⎥
⎢ ── - ──    ── + ──   ── - ──    ── + ──⎥
⎢ 3    2     3    2    3    2     3    2 ⎥
⎢ℯ          ℯ         ℯ          ℯ       ⎥
⎢──────── + ────────  ──────── - ────────⎥
⎢   2          2         2          2    ⎥
⎢                                        ⎥
⎢  3    2     3    2    3    2     3    2⎥
⎢ t    t     t    t    t    t     t    t ⎥
⎢ ── - ──    ── + ──   ── - ──    ── + ──⎥
⎢ 3    2     3    2    3    2     3    2 ⎥
⎢ℯ          ℯ         ℯ          ℯ       ⎥
⎢──────── - ────────  ──────── + ────────⎥
⎣   2          2         2          2    ⎦

In [8]: M.jordan_form()                                                                                                           
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-8-192ccf6137fc> in <module>
----> 1 M.jordan_form()

~/current/sympy/sympy/sympy/matrices/matrices.py in jordan_form(self, calc_transform, **kwargs)
   1806             if not calc_transform:
   1807                 return restore_floats(jordan_mat)
-> 1808             jordan_basis = [eig_mat(eig, 1).nullspace()[0] for eig in blocks]
   1809             basis_mat = mat.hstack(*jordan_basis)
   1810             return restore_floats(basis_mat, jordan_mat)

~/current/sympy/sympy/sympy/matrices/matrices.py in <listcomp>(.0)
   1806             if not calc_transform:
   1807                 return restore_floats(jordan_mat)
-> 1808             jordan_basis = [eig_mat(eig, 1).nullspace()[0] for eig in blocks]
   1809             basis_mat = mat.hstack(*jordan_basis)
   1810             return restore_floats(basis_mat, jordan_mat)

IndexError: list index out of range
```
Looks like nullspace returns an empty list:
```julia
In [17]: M = Matrix([ 
    ...: [sqrt((exp(2*t**2) - 2*exp(t**2) + 1)*exp(t**2*(2*t/3 - 1)))/2 - (exp(t**2) + 1)*exp(t**2*(2*t - 3)/6)/2 + exp(t**3/3 - t
    ...: **2/2)/2 + exp(t**3/3 + t**2/2)/2,                                                                                       
    ...:                     exp(t**3/3 - t**2/2)/2 - exp(t**3/3 + t**2/2)/2], 
    ...: [                                                                                                          exp(t**3/3 - t
    ...: **2/2)/2 - exp(t**3/3 + t**2/2)/2, sqrt((exp(2*t**2) - 2*exp(t**2) + 1)*exp(t**2*(2*t/3 - 1)))/2 - (exp(t**2) + 1)*exp(t*
    ...: *2*(2*t - 3)/6)/2 + exp(t**3/3 - t**2/2)/2 + exp(t**3/3 + t**2/2)/2]])                                                   

In [18]: M                                                                                                                        
Out[18]: 
⎡     _____________________________________                 2               3    2     3    2                                   
⎢    ╱                         2 ⎛2⋅t    ⎞                 t ⋅(2⋅t - 3)    t    t     t    t                                    
⎢   ╱  ⎛    2      ⎛ 2⎞    ⎞  t ⋅⎜─── - 1⎟    ⎛ ⎛ 2⎞    ⎞  ────────────    ── - ──    ── + ──                                   
⎢  ╱   ⎜ 2⋅t       ⎝t ⎠    ⎟     ⎝ 3     ⎠    ⎜ ⎝t ⎠    ⎟       6          3    2     3    2                                    
⎢╲╱    ⎝ℯ     - 2⋅ℯ     + 1⎠⋅ℯ                ⎝ℯ     + 1⎠⋅ℯ               ℯ          ℯ                                          
⎢────────────────────────────────────────── - ───────────────────────── + ──────── + ────────                                   
⎢                    2                                    2                  2          2                                       
⎢                                                                                                                               
⎢                                      3    2     3    2                                            ____________________________
⎢                                     t    t     t    t                                            ╱                         2 ⎛
⎢                                     ── - ──    ── + ──                                          ╱  ⎛    2      ⎛ 2⎞    ⎞  t ⋅⎜
⎢                                     3    2     3    2                                          ╱   ⎜ 2⋅t       ⎝t ⎠    ⎟     ⎝
⎢                                    ℯ          ℯ                                              ╲╱    ⎝ℯ     - 2⋅ℯ     + 1⎠⋅ℯ    
⎢                                    ──────── - ────────                                       ─────────────────────────────────
⎣                                       2          2                                                               2            

     3    2     3    2                                     ⎤
    t    t     t    t                                      ⎥
    ── - ──    ── + ──                                     ⎥
    3    2     3    2                                      ⎥
   ℯ          ℯ                                            ⎥
   ──────── - ────────                                     ⎥
      2          2                                         ⎥
                                                           ⎥
_________                 2               3    2     3    2⎥
2⋅t    ⎞                 t ⋅(2⋅t - 3)    t    t     t    t ⎥
─── - 1⎟    ⎛ ⎛ 2⎞    ⎞  ────────────    ── - ──    ── + ──⎥
 3     ⎠    ⎜ ⎝t ⎠    ⎟       6          3    2     3    2 ⎥
            ⎝ℯ     + 1⎠⋅ℯ               ℯ          ℯ       ⎥
───────── - ───────────────────────── + ──────── + ────────⎥
                        2                  2          2    ⎦

In [19]: M.nullspace()                                                                                                            
Out[19]: []
```
After simplifying the nullspace is found:
```julia
In [20]: M.simplify()                                                                                                             

In [21]: M                                                                                                                        
Out[21]: 
⎡      __________________________________                                          ⎤
⎢     ╱                           3                                3    2          ⎥
⎢    ╱                         2⋅t     2                          t    t           ⎥
⎢   ╱   ⎛    2      ⎛ 2⎞    ⎞  ──── - t              ⎛     ⎛ 2⎞⎞  ── - ──          ⎥
⎢  ╱    ⎜ 2⋅t       ⎝t ⎠    ⎟   3                    ⎜     ⎝t ⎠⎟  3    2           ⎥
⎢╲╱     ⎝ℯ     - 2⋅ℯ     + 1⎠⋅ℯ                      ⎝1 - ℯ    ⎠⋅ℯ                 ⎥
⎢────────────────────────────────────────            ────────────────────          ⎥
⎢                   2                                         2                    ⎥
⎢                                                                                  ⎥
⎢                                                __________________________________⎥
⎢                        3    2                 ╱                           3      ⎥
⎢                       t    t                 ╱                         2⋅t     2 ⎥
⎢          ⎛     ⎛ 2⎞⎞  ── - ──               ╱   ⎛    2      ⎛ 2⎞    ⎞  ──── - t  ⎥
⎢          ⎜     ⎝t ⎠⎟  3    2               ╱    ⎜ 2⋅t       ⎝t ⎠    ⎟   3        ⎥
⎢          ⎝1 - ℯ    ⎠⋅ℯ                   ╲╱     ⎝ℯ     - 2⋅ℯ     + 1⎠⋅ℯ          ⎥
⎢          ────────────────────            ────────────────────────────────────────⎥
⎣                   2                                         2                    ⎦

In [22]: M.nullspace()                                                                                                            
Out[22]: 
⎡⎡                        3    2          ⎤⎤
⎢⎢                       t    t           ⎥⎥
⎢⎢          ⎛     ⎛ 2⎞⎞  ── - ──          ⎥⎥
⎢⎢          ⎜     ⎝t ⎠⎟  3    2           ⎥⎥
⎢⎢         -⎝1 - ℯ    ⎠⋅ℯ                 ⎥⎥
⎢⎢────────────────────────────────────────⎥⎥
⎢⎢      __________________________________⎥⎥
⎢⎢     ╱                           3      ⎥⎥
⎢⎢    ╱                         2⋅t     2 ⎥⎥
⎢⎢   ╱   ⎛    2      ⎛ 2⎞    ⎞  ──── - t  ⎥⎥
⎢⎢  ╱    ⎜ 2⋅t       ⎝t ⎠    ⎟   3        ⎥⎥
⎢⎢╲╱     ⎝ℯ     - 2⋅ℯ     + 1⎠⋅ℯ          ⎥⎥
⎢⎢                                        ⎥⎥
⎣⎣                   1                    ⎦⎦
```
I guess this is a problem in nullspace then.

Hints from issue discussion:
Same here, with much simplier matrix. Consider
```
from numpy import array
from sympy import Matrix, pprint

S0 = array([ [-0.7, 0.0, 0.0, 0.5],
             [0.5, -0.2, 0.0, 0.25],
             [0.0, 0.2, -0.2, 0.0],
             [0.2, 0.0, 0.2, -0.75]
           ])

S1 = array([ [-1.25, 0.5, 0., 0.],
             [1., -0.5, 0., 0.5],
             [0.25, 0., -1., 0.5],
             [0., 0., 1., -1.]
           ])

m_s = Matrix(S0)
pprint(m_s.nullspace())
```
For other matrices with det=0 (such as S1) it prints column of nullspace elements. For matrix S0 it prints '[]'. Broken nullspace() breaks jordan_form(), which in turn breaks dsolve() for the equivalent system of ODEs.
I'm using sympy 1.6.2.
EDIT: _iszero() problem? Pull #15169 maybe related, custom iszerofunc=my_iszero from that pull can't help. At least with iszerofunc like:
```
def my_iszero():
    return True
```
nullspace() returns entire identity matrix  ;)
UPDATE: For those who want to use dsolve() || Matrix.jordan_form() || Matrix.nullspace() with constant quotitent matrices I suggest to use this workaround:
```
from sympy import Matrix
from scipy.linalg import null_space
from numpy import array, cdouble

def nullspace_sci(matrix):
    ''' SciPy version of nullspace() '''
    return [Matrix(null_space(array(matrix).astype(cdouble)))]

Matrix.nullspace = nullspace_sci
```
Very fragile, but works with S0.
Using `nsimplify` to convert the floats to rational works for your case:
```
In [4]: m_s.applyfunc(nsimplify).nullspace()
Out[4]: 
⎡⎡5/7⎤⎤
⎢⎢   ⎥⎥
⎢⎢85 ⎥⎥
⎢⎢── ⎥⎥
⎢⎢28 ⎥⎥
⎢⎢   ⎥⎥
⎢⎢85 ⎥⎥
⎢⎢── ⎥⎥
⎢⎢28 ⎥⎥
⎢⎢   ⎥⎥
⎣⎣ 1 ⎦⎦
```
It's ok with nullspace() but I wonder how can I propagate this behaviour up to dsolve().
AFAIK, dsolve() deals with list of Eq's. Currently I'm preparing ICs and Eq's by hand with given matrix. And dsolve() decomposes it back to matrix...
Also this Issue #17218 definitely related with my problem, because
```
pprint(Matrix(S0).rref()) 
pprint(Matrix(S0).rref(iszerofunc=lambda x: abs(x) < 10**-10))     # Different, more plausible output 
```
So to fix dsolve() I need to just set iszerofunc in Matrix.nullspace() (which will passed to .rref()) to deal with it w/o SciPy.
You can use `nsimplify` on your original equations. It is generally better to avoid floats especially for something like `nullspace`:
https://en.wikipedia.org/wiki/Kernel_(linear_algebra)#Floating_point_computation
Even if I apply nsimplify to every equation in system (with rational=True, checking results after that with pprint()), aforementioned IndexError in dsolve()->jordan_form()->nullspace()[0] still persists.
Also, replacing iszerofunc default value not a solution since dsolve() calls it from different contexts and subtle errors appear. 
Do you have an example where this happens with rational numbers? I think that the original matrix in the OP is also something that came up from `dsolve`.
The OP examples seem to work fine now for `jordan_form` and `nullspace`.

I think that this issue is fixed and the new example you've introduced is really another case of #17218 which is about matrices with floats.
This is somewhat truncated example (S0 and S1 same as above):
```
from numpy import array
from sympy import Eq, Derivative, symbols, Function, sympify, Add, dsolve, nsimplify
from sympy import pprint

S0 = array([ [-0.7, 0.0, 0.0, 0.5],
             [0.5, -0.2, 0.0, 0.25],
             [0.0, 0.2, -0.2, 0.0],
             [0.2, 0.0, 0.2, -0.75]
           ])

S1 = array([ [-1.25, 0.5, 0., 0.],
             [1., -0.5, 0., 0.5],
             [0.25, 0., -1., 0.5],
             [0., 0., 1., -1.]
           ])
S = S0  # S1 OK, S0 FAIL
ICS = sympify('{P1(0):0, P2(0):0, P3(0):0, P4(0):1}')

def sym_sol():
    dim = S.shape[0]
    t_s = symbols('t')
    p_i_syms = symbols('P1, P2, P3, P4', cls=Function)
    eq_sys = [None] * S.shape[0]
    for i, p_i in enumerate(p_i_syms):
        rhs_i = 0
        for j in range(dim):
            rhs_i = Add(rhs_i, S[i][j]*p_i_syms[j](t_s))
        eq_sys[i] = nsimplify(Eq(Derivative(p_i(t_s)), rhs_i), rational=True)
    eq_sys = tuple(eq_sys)
    #pprint(eq_sys)
    # Hint obtained by benchmarking
    sol = dsolve(eq_sys, ics=ICS, hint='1st_power_series')
    print('Solution:')
    for equ in sol:
        pprint(equ)

sym_sol()
```
Okay so the matrix there is
```python
M = Matrix([
[-7/10,    0,    0,  1/2],
[  1/2, -1/5,    0,  1/4],
[    0,  1/5, -1/5,    0],
[  1/5,    0,  1/5, -3/4]]).applyfunc(nsimplify)
```
That does not have floats and gives:
```python
In [20]: M
Out[20]: 
⎡-7/10   0     0    1/2 ⎤
⎢                       ⎥
⎢ 1/2   -1/5   0    1/4 ⎥
⎢                       ⎥
⎢  0    1/5   -1/5   0  ⎥
⎢                       ⎥
⎣ 1/5    0    1/5   -3/4⎦

In [21]: M.jordan_form()
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-21-192ccf6137fc> in <module>
----> 1 M.jordan_form()

~/current/sympy/sympy/sympy/matrices/matrices.py in jordan_form(self, calc_transform, **kwargs)
    410 
    411     def jordan_form(self, calc_transform=True, **kwargs):
--> 412         return _jordan_form(self, calc_transform=calc_transform, **kwargs)
    413 
    414     def left_eigenvects(self, **flags):

~/current/sympy/sympy/sympy/matrices/eigen.py in _jordan_form(M, calc_transform, chop)
   1185             return restore_floats(jordan_mat)
   1186 
-> 1187         jordan_basis = [eig_mat(eig, 1).nullspace()[0]
   1188                 for eig in blocks]
   1189         basis_mat    = mat.hstack(*jordan_basis)

~/current/sympy/sympy/sympy/matrices/eigen.py in <listcomp>(.0)
   1185             return restore_floats(jordan_mat)
   1186 
-> 1187         jordan_basis = [eig_mat(eig, 1).nullspace()[0]
   1188                 for eig in blocks]
   1189         basis_mat    = mat.hstack(*jordan_basis)

IndexError: list index out of range
> /Users/enojb/current/sympy/sympy/sympy/matrices/eigen.py(1187)<listcomp>()
   1185             return restore_floats(jordan_mat)
   1186 
-> 1187         jordan_basis = [eig_mat(eig, 1).nullspace()[0]
   1188                 for eig in blocks]
   1189         basis_mat    = mat.hstack(*jordan_basis)
```
Internally this computes nullspaces e.g.:
```python
In [29]: [len((M - e*eye(4)).nullspace()) for e in M.eigenvals()]
Out[29]: [1, 0, 0, 1]
```
So it looks like it fails for the 2nd and 3rd eigenvalues. Note these have complicated expressions from the cubic formula:
```python
In [30]: e1, e2, e3, e4 = M.eigenvals()

In [31]: e2
Out[31]: 
                        ______________                                      
       ⎛  1   √3⋅ⅈ⎞    ╱ 69⋅√165   23                                       
       ⎜- ─ + ────⎟⋅3 ╱  ─────── + ──                                       
  37   ⎝  2    2  ⎠ ╲╱     1600    40                     23                
- ── - ─────────────────────────────── - ───────────────────────────────────
  60                  3                                       ______________
                                             ⎛  1   √3⋅ⅈ⎞    ╱ 69⋅√165   23 
                                         240⋅⎜- ─ + ────⎟⋅3 ╱  ─────── + ── 
                                             ⎝  2    2  ⎠ ╲╱     1600    40 

In [32]: e3
Out[32]: 
                                                              ______________
                                             ⎛  1   √3⋅ⅈ⎞    ╱ 69⋅√165   23 
                                             ⎜- ─ - ────⎟⋅3 ╱  ─────── + ── 
  37                    23                   ⎝  2    2  ⎠ ╲╱     1600    40 
- ── - ─────────────────────────────────── - ───────────────────────────────
  60                        ______________                  3               
           ⎛  1   √3⋅ⅈ⎞    ╱ 69⋅√165   23                                   
       240⋅⎜- ─ - ────⎟⋅3 ╱  ─────── + ──                                   
           ⎝  2    2  ⎠ ╲╱     1600    40 
```
There is I guess a bug somewhere in zero detection during nullspace.

Note that there is a new implementation of eigenvectors on the master branch (not yet released) that can calculate the eigenvectors much faster and correctly:
```python
In [58]: %time M.eigenvects()
CPU times: user 29.9 ms, sys: 274 µs, total: 30.2 ms
Wall time: 30.5 ms
```
That implementation should also be extended to `jordan_form` although it is still a work in progress.

The bug in the old implementation still warrants investigation.
The bug is that zero testing fails on this pivot:
```python
expr = (
    (
    -512109718324422794624654112703630525511*2645**(S(1)/3)*sqrt(3) +
    512109718324422794624654112703630525511*23**(S(2)/3)*sqrt(3)*5**(S(1)/3) -
    512109718324422794624654112703630525511*23**(S(2)/3)*5**(S(1)/3)*I +
    512109718324422794624654112703630525511*2645**(S(1)/3)*I)
    /(11232853364312579514316204343197764426240000*sqrt(3)*5**(S(2)/3)*(3*sqrt(165)
        + 40)**(S(1)/3) +
        13117146726057831427621087078765368441868800*sqrt(11)*5**(1/6)*(3*sqrt(165)
            + 40)**(S(1)/3) +
        4372382242019277142540362359588456147289600*sqrt(33)*5**(1/6)*I*(3*sqrt(165)
            + 40)**(S(1)/3) +
        11232853364312579514316204343197764426240000*5**(S(2)/3)*I*(3*sqrt(165) +
            40)**(S(1)/3))
    )

print('expr.n():', expr.n())
print('expr.is_zero:', expr.is_zero)
```
That outputs:
```python
expr.n(): 0.e-144 - 0.e-146*I
expr.is_zero: False
```
That in turn boils down to this:
```python
In [8]: e = parse_expr('-2645**(1/3)*sqrt(3) + 23**(2/3)*sqrt(3)*5**(1/3) - 23**(2/3)*5**(1/3)*I + 2645**(1/3)*I')

In [9]: e
Out[9]: 
  3 ______        2/3    3 ___     2/3 3 ___     3 ______  
- ╲╱ 2645 ⋅√3 + 23   ⋅√3⋅╲╱ 5  - 23   ⋅╲╱ 5 ⋅ⅈ + ╲╱ 2645 ⋅ⅈ

In [10]: e.n()
Out[10]: 0.e-124 + 0.e-124⋅ⅈ

In [11]: e.is_zero
Out[11]: False
```
The faulty logic is here:
https://github.com/sympy/sympy/blob/c9dbb014993f6248be0a4a3d514545cbf7d4c1c1/sympy/core/add.py#L610-L615

Looking at the code suggests a simpler way to show the bug:
```python
In [1]: x = Symbol('x', positive=True)

In [2]: e = Add(x, -x, I, -I, evaluate=False)

In [3]: e.is_zero
Out[3]: False

In [4]: e
Out[4]: -x + x - ⅈ + ⅈ
```

The fix is
```diff
diff --git a/sympy/core/add.py b/sympy/core/add.py
index 772def1604..aab9e3f5b8 100644
--- a/sympy/core/add.py
+++ b/sympy/core/add.py
@@ -588,7 +588,7 @@ def _eval_is_zero(self):
         nz = []
         z = 0
         im_or_z = False
-        im = False
+        im = 0
         for a in self.args:
             if a.is_extended_real:
                 if a.is_zero:
@@ -598,7 +598,7 @@ def _eval_is_zero(self):
                 else:
                     return
             elif a.is_imaginary:
-                im = True
+                im += 1
             elif (S.ImaginaryUnit*a).is_extended_real:
                 im_or_z = True
             else:
@@ -609,10 +609,11 @@ def _eval_is_zero(self):
             return None
         b = self.func(*nz)
         if b.is_zero:
-            if not im_or_z and not im:
-                return True
-            if im and not im_or_z:
-                return False
+            if not im_or_z:
+                if im == 0:
+                    return True
+                elif im == 1:
+                    return False
         if b.is_zero is False:
             return False
```
With that diff it seems that `jordan_form` takes a two minutes to complete but comes up with a presumably correct answer (it's slow because it's trying really hard to simplfy the expressions to figure out whether or not they are zero).

The newer implementation of `jordan_form` would presumably finish in milliseconds like the newer eigenvects implementation so that's a fixable problem.

Commit message: core: fixed _eval_is_zero() and added test, fixing issue #17130

_eval_is_zero() would give wrong result with imaginary numbers,
resolved as directed by @oscarbenjamin, and added
appropriate test.

Patch to review:
diff --git a/sympy/core/add.py b/sympy/core/add.py
index 6a3d10ac0e66..301d0132c5ea 100644
--- a/sympy/core/add.py
+++ b/sympy/core/add.py
@@ -604,7 +604,7 @@ def _eval_is_zero(self):
         nz = []
         z = 0
         im_or_z = False
-        im = False
+        im = 0
         for a in self.args:
             if a.is_extended_real:
                 if a.is_zero:
@@ -614,7 +614,7 @@ def _eval_is_zero(self):
                 else:
                     return
             elif a.is_imaginary:
-                im = True
+                im += 1
             elif (S.ImaginaryUnit*a).is_extended_real:
                 im_or_z = True
             else:
@@ -625,10 +625,12 @@ def _eval_is_zero(self):
             return None
         b = self.func(*nz)
         if b.is_zero:
-            if not im_or_z and not im:
-                return True
-            if im and not im_or_z:
-                return False
+            if not im_or_z:
+                if im == 0:
+                    return True
+                elif im == 1:
+                    return False
+            return True
         if b.is_zero is False:
             return False
 
diff --git a/sympy/core/tests/test_arit.py b/sympy/core/tests/test_arit.py
index 2cd5e484bb98..27d40ffcf432 100644
--- a/sympy/core/tests/test_arit.py
+++ b/sympy/core/tests/test_arit.py
@@ -2307,3 +2307,8 @@ def test__neg__():
 
 def test_issue_18507():
     assert Mul(zoo, zoo, 0) is nan
+
+
+def test_issue_17130():
+    e = Add(b, -b, I, -I, evaluate=False)
+    assert e._eval_is_zero() == True


Code review comments:
diff_hunk: @@ -625,10 +625,12 @@ def _eval_is_zero(self):
             return None
         b = self.func(*nz)
         if b.is_zero:
-            if not im_or_z and not im:
-                return True
-            if im and not im_or_z:
-                return False
+            if not im_or_z:
+                if im == 0:
+                    return True
+                elif im == 1:
+                    return False
+            return True
path: sympy/core/add.py
text: @user1:
This return True was not there in my suggestion in https://github.com/sympy/sympy/issues/17130#issuecomment-749117670
Why is it needed?

@author:
I tried doing it without this `return True` line, but then the test I wrote as suggested by you in #17130 was failing.

Since this function was not returning anything when `b.is_zero` is true and `im_or_z` is true, I added this line so as to return true when `b.is_zero` is true and `im_or_z` is true.


@user1:
What test was failing?

The return True here will give erroneous results:
```python
In [2]: Add(1, -1, I, -I, I, evaluate=False).is_zero
Out[2]: True
```

@author:
It was failing the test you suggested in the [comment](https://github.com/sympy/sympy/issues/17130#issuecomment-749117670). 
It is:
```python
def test_issue_17130():
    e = Add(b, -b, I, -I, evaluate=False)
    assert e._eval_is_zero() == True
```

@user1:
It's fine if that case gives `None`. In the assumptions system we just need to make sure that we don't incorrectly return True or False. There will always be cases that are not resolvable because an assumptions check like `is_zero` needs to be efficient.

The test could just be for None but with a comment indicating why:
```
    assert e._eval_is_zero() is None # ideally this would be True
```

@author:
I see, thanks for the help. I will change the code accordingly.