Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

(1/x).evalf(subs={x: 0}) produces 0 #8242

Closed
immerrr opened this issue Oct 14, 2014 · 8 comments · Fixed by #22211
Closed

(1/x).evalf(subs={x: 0}) produces 0 #8242

immerrr opened this issue Oct 14, 2014 · 8 comments · Fixed by #22211
Labels
core.evalf Wrong Result The output produced by SymPy is mathematically incorrect.

Comments

@immerrr
Copy link
Contributor

immerrr commented Oct 14, 2014

Here's how it looks like:

In [1]: sp.__version__
Out[1]: '0.7.5-git'

In [2]: x = sp.abc.x

In [3]: (1/x).evalf(subs={x: 0})
Out[3]: 0

In [4]: (1/x).subs({x: 0})
Out[4]: zoo

In [5]: (1/x).xreplace({x: 0})
Out[5]: zoo

I did investigate it some time ago and, AFAIR, it happens in Pow evaluation that short-circuits to 0 whenever the base is 0 while it should return infs if the exponent is negative.

@pelegm
Copy link
Contributor

pelegm commented Oct 15, 2014

Another interesting and related issue: while

In [3]: x = Symbol('x', real=True)

In [4]: (1/x).subs({x: 0.})
Out[4]: +inf

In [5]: (1/x).xreplace({x: 0.})
Out[5]: +inf

the following

(1/x).evalf(subs={x: 0.})

raises

Traceback (most recent call last):
  File "<ipython-input-6-226388149257>", line 1, in <module>
    (1/x).evalf(subs={x: 0.})
  File "sympy/core/evalf.py", line 1326, in evalf
    result = evalf(self, prec + 4, options)
  File "sympy/core/evalf.py", line 1220, in evalf
    r = rf(x, prec, options)
  File "sympy/core/evalf.py", line 622, in evalf_pow
    return mpf_pow_int(re, p, target_prec), None, target_prec, None
  File "sympy/mpmath/libmp/libmpf.py", line 1044, in mpf_pow_int
    if n == -1: return mpf_div(fone, s, prec, rnd)
  File "sympy/mpmath/libmp/libmpf.py", line 934, in mpf_div
    raise ZeroDivisionError
ZeroDivisionError

@smichr
Copy link
Member

smichr commented Jul 7, 2017

Valid for sympy-1.0 and sympy-1.1-58-gd250ff0, too, as noted in #12904. There are many related issues:

#10662
#9807
#7867
#7192
#6974

@asmeurer asmeurer added the Wrong Result The output produced by SymPy is mathematically incorrect. label Jul 7, 2017
@smichr
Copy link
Member

smichr commented Jul 7, 2017

I get these results

>>> (x*y).n(subs={y:0})
0
>>> log(x).evalf(subs={x: 0})
zoo
>>> (x-1).evalf(subs={x:1})
0.e-125
>>> (x**I).evalf(subs={x:0})
nan

with these mods (which are suggestive as to where the problem lies)

diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index c794656..233bbe5 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -821,6 +821,8 @@ def evalf_log(expr, prec, options):
 
     imaginary_term = (mpf_cmp(xre, fzero) < 0)
 
+    if xre == (0,)*4:
+        return S.ComplexInfinity
     re = mpf_log(mpf_abs(xre), prec, rnd)
     size = fastlog(re)
     if prec - size > workprec and re != fzero:
@@ -1284,13 +1286,25 @@ def evalf(x, prec, options):
     try:
         rf = evalf_table[x.func]
         r = rf(x, prec, options)
-    except KeyError:
+        if r is S.ComplexInfinity:
+            return r
+        if x and r[0] == r[1] == None:
+            assert not 'subs' in options
+    except (KeyError, AssertionError):
         try:
             # Fall back to ordinary evalf if possible
             if 'subs' in options:
                 x = x.subs(evalf_subs(prec, options['subs']))
-            xe = x._eval_evalf(prec)
-            re, im = xe.as_real_imag()
+            if x is S.ComplexInfinity:
+                return x
+            if x:
+                if not x.is_number:
+                    return (None,)*4
+                x = evalf(x, prec, options)
+                xe = x._eval_evalf(prec)
+                re, im = xe.as_real_imag()
+            else:
+                return (0, 0, 0, 0), None, None, None
             if re.has(re_) or im.has(im_):
                 raise NotImplementedError
             if re == 0:
@@ -1306,6 +1320,8 @@ def evalf(x, prec, options):
                 im = im._to_mpmath(prec, allow_ints=False)._mpf_
                 imprec = prec
             r = re, im, reprec, imprec
+            if re == im == None:
+                return (0, 0, 0, 0), None, None, None
         except AttributeError:
             raise NotImplementedError
     if options.get("verbose"):
@@ -1392,6 +1408,8 @@ def evalf(self, n=15, subs=None, maxn=100, chop=False, strict=False, quad=None,
             options['quad'] = quad
         try:
             result = evalf(self, prec + 4, options)
+            if result is S.ComplexInfinity:
+                return result
         except NotImplementedError:
             # Fall back to the ordinary evalf
             v = self._eval_evalf(prec)

Whoever tackles this beware: there is the evalf method and the evalf function (which references the evalf_table (of special evalf functions) and the _eval_evalf method. This is not a beginner's issue...unless you are a special beginner!

@asmeurer asmeurer marked this as a duplicate of #12904 Jul 25, 2017
@asmeurer
Copy link
Member

I believe the problem is here. evalf_pow assumes that 0 to an integer power is 0 (Nones represent exact zeros). But this is only true if the exponent is > 0. The question is, what should it return. I guess we want it to give zoo, but what value to we return from that function to give that (there is no floating point version of zoo). I guess your patch is in that direction, although I don't like that it makes the evalf functions return just zoo instead of a tuple.

We could also make it return inf quite easily:

diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index c79465651..9437c9c17 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -670,7 +670,12 @@ def evalf_pow(v, prec, options):
                 return None, mpf_neg(z), None, target_prec
         # Zero raised to an integer power
         if not re:
-            return None, None, None, None
+            if exp > 0:
+                return None, None, None, None
+            elif exp < 0:
+                return mpmath_inf._mpf_, None, target_prec, None
+            else:
+                return mpf(1)._mpf_, None, target_prec, None
         # General complex number to arbitrary integer power
         re, im = libmp.mpc_pow_int((re, im), p, prec)
         # Assumes full accuracy in input

Maybe we should apply the above as a simpler fix for now (inf is not as correct as zoo, but it's definitely less wrong than 0).

@smichr
Copy link
Member

smichr commented Aug 12, 2017

SO reports this issue, too.

@eagleoflqj
Copy link
Member

I want to work on this. Ideally mpmath should provide a tuple that represents zoo, but it doesn't. And there is no benefit that we assign a tuple as zoo then check this tuple everywhere, so I'd like to follow @smichr 's idea that directly return the zoo when appropriate.

eagleoflqj added a commit to eagleoflqj/sympy that referenced this issue Oct 3, 2021
eagleoflqj added a commit to eagleoflqj/sympy that referenced this issue Oct 4, 2021
eagleoflqj added a commit to eagleoflqj/sympy that referenced this issue Oct 5, 2021
@bessavagner
Copy link

bessavagner commented Sep 4, 2022

I've found the same issue when using numeric evaluation on parser with evaluate=False:

import sympy as sp
parser = sp.parsing.sympy_parser.parse_expr
expr = '1/0'
sp.N(parser(expr, evaluate=False))

>>> 0

sympy==1.11.1
python==3.9.13

@oscarbenjamin
Copy link
Collaborator

sympy==1.11.1

Are you sure about this? I tried with 1.11.1 and also with current master:

In [1]: import sympy as sp
   ...: parser = sp.parsing.sympy_parser.parse_expr
   ...: expr = '1/0'
   ...: sp.N(parser(expr, evaluate=False))
Out[1]: zoo

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
core.evalf Wrong Result The output produced by SymPy is mathematically incorrect.
Projects
None yet
7 participants