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

bad ceiling(sqrt(big integer)) #10323

Closed
swn1 opened this Issue Dec 26, 2015 · 9 comments

Comments

Projects
None yet
3 participants
@swn1
Copy link

swn1 commented Dec 26, 2015

import sympy
ll = 1206577996382235787095214 # lower limit for x**2
x = sympy.ceiling(sympy.sqrt(ll))
assert (x-1)**2 < ll and x**2 >= ll

fails in a python3 notebook on jupyter.org.

@swn1

This comment has been minimized.

Copy link

swn1 commented Dec 26, 2015

In hindsight, "ll" is a poor choice for inclusion in this repro case. instead of "lower limit" it shouts "eleven". Sigh.

@swn1

This comment has been minimized.

Copy link

swn1 commented Dec 26, 2015

a more thorough set of asserts.

import sympy
ll = 1206577996382235787095214
x = sympy.ceiling(sympy.sqrt(ll))
y = sympy.floor(sympy.sqrt(ll))
assert (x == y) == (x**2 == ll)
assert (x-1)**2 < ll and x**2 >= ll
assert (y+1)**2 > ll and y**2 <= ll
@bjodah

This comment has been minimized.

Copy link
Member

bjodah commented Dec 27, 2015

I can confirm this bug on Python 2.7, sympy 0.7.6.1, and even for smaller numbers (inside the range of ordinary 32-bit integers):

>>> def check(x):
...     c = sympy.ceiling(sympy.sqrt(x))
...     f = sympy.floor(sympy.sqrt(x))
...     assert (c - 1)**2 < x and c**2 >= x
...     assert (f + 1)**2 > x and f**2 <= x
... 
>>> check(2**30 + 1)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "<stdin>", line 5, in check
AssertionError
@smichr

This comment has been minimized.

Copy link
Member

smichr commented Dec 28, 2015

I get

>>> ll=1206577996382235787095214
>>> x=ceiling(sqrt(ll))
>>> assert x**2 >= ll
>>> 

with the following change (correcting what I believe to be a misuse of and with negated args):

diff --git a/sympy/core/evalf.py b/sympy/core/evalf.py
index 2fc8d98..cc25240 100644
--- a/sympy/core/evalf.py
+++ b/sympy/core/evalf.py
@@ -321,7 +321,8 @@ def calc_part(expr, nexpr):
         from sympy import Add
         nint = int(to_int(nexpr, rnd))
         n, c, p, b = nexpr
-        if (c != 1 and p != 0) or p < 0:
+        is_one = c == 1 and p == 0
+        if not is_one or p < 0:
             expr = Add(expr, -nint, evaluate=False)
             x, _, x_acc, _ = evalf(expr, 10, options)
             try:

Would you like to make the change and add a test?

@swn1

This comment has been minimized.

Copy link

swn1 commented Dec 28, 2015

Captain DeMorgan strikes again! I don't know the logic flow intimately but I think your sense of the coding error is likely correct. If so, it would seem more unit tests are needed closer to this level.

Adding the test case seems like good practice whether or not the fix is the right one.

But I'm a newbie here.

@bjodah

This comment has been minimized.

Copy link
Member

bjodah commented Dec 31, 2015

@swn1 yes, we should definitely add a unit test for this fix to avoid regressions. You may add these changes yourself and open a PR (you can find instructions here and you can always ask us for help), or let us know if you want someone else to add @smichr's fix.

@swn1

This comment has been minimized.

Copy link

swn1 commented Dec 31, 2015

I have some doubts about the suggested fix.
is_one = c == 1 and p == 0
is equivalent to not(c != 1 or p != 0) by De Morgan's law.
then not is_zero or p < 0 is just not is_zero by associativity and p<0 ==> p!=0.
So the p<0 clause is no longer doing anything.

Unless I'm smoking crack, of course. But I get a new failure running bin/test evalf with the patch in.

So I'm going to have to understand this code a bit more deeply.

@swn1

This comment has been minimized.

Copy link

swn1 commented Dec 31, 2015

adding a print statement to the calc_part function we see a discontinuity in the structure of the call at the point where the incorrect results kick in. The problem is not (just) in calc_part.

In [1]: for o in range(10,0,-1):
   ...:     print(ceiling(sqrt(2**30+o))**2)
   ...:     
re(sqrt(1073741834)) (0, 268435457, -13, 29)
1073807361
re(sqrt(1073741833)) (0, 268435457, -13, 29)
1073807361
re(18*sqrt(3314018)) (0, 268435457, -13, 29)
1073807361
re(sqrt(1073741831)) (0, 536870913, -14, 30)
1073807361
re(sqrt(1073741830)) (0, 536870913, -14, 30)
1073807361
re(sqrt(1073741829)) (0, 536870913, -14, 30)
1073807361
re(2*sqrt(268435457)) (0, 536870913, -14, 30)
1073807361
re(sqrt(1073741827)) (0, 1, 15, 1)
1073741824
re(sqrt(1073741826)) (0, 1, 15, 1)
1073741824
re(5*sqrt(42949673)) (0, 1, 15, 1)
1073741824

skirpichev added a commit to diofant/diofant that referenced this issue Jan 1, 2016

Use mpmath's floor/ceil to calculate round/ceiling
Drop get_integer_part() helper function.

Fixes sympy/sympy#10323

skirpichev added a commit to diofant/diofant that referenced this issue Jan 1, 2016

Use mpmath's floor/ceil to calculate round/ceiling
Drop get_integer_part() helper function.

Fixes sympy/sympy#10323
@smichr

This comment has been minimized.

Copy link
Member

smichr commented Jan 2, 2016

There is an alternative approach but I am not sure of the relative merits. The author of mpmath is the one that added the get_integer_parts in the first place, so I assume there is some efficiency in this function.

@smichr smichr closed this in #10346 Jan 7, 2016

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment