Opened 5 years ago

Closed 5 years ago

#24378 closed defect (fixed)

complex_root_of uses inexact index

Reported by: Jeroen Demeyer Owned by:
Priority: blocker Milestone: sage-8.2
Component: symbolics Keywords:
Cc: Ralf Stephan, Emmanuel Charpentier Merged in:
Authors: Ralf Stephan Reviewers: Jeroen Demeyer
Report Upstream: N/A Work issues:
Branch: 277a403 (Commits, GitHub, GitLab) Commit: 277a403ef00692247f9853d89322f9a464abd471
Dependencies: Stopgaps:

Status badges

Description (last modified by Jeroen Demeyer)

When numerically evaluating complex_root_of(a, k), the index k is passed as floating-point number:

sage: complex_root_of(x^8-1, 7).n(2)
---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-1-4eba84b7c14c> in <module>()
----> 1 complex_root_of(x**Integer(8)-Integer(1), Integer(7)).n(Integer(2))

/home/patchbot/sage-patchbot/src/sage/structure/element.pyx in sage.structure.element.Element.n (build/cythonized/sage/structure/element.c:8131)()
    863             0.666666666666667
    864         """
--> 865         return self.numerical_approx(prec, digits, algorithm)
    866 
    867     N = deprecated_function_alias(13055, n)

/home/patchbot/sage-patchbot/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression.numerical_approx (build/cythonized/sage/symbolic/expression.cpp:36209)()
   5818         kwds = {'parent': R, 'algorithm': algorithm}
   5819         try:
-> 5820             x = x._convert(kwds)
   5821         except TypeError: # numerical approximation for real number failed
   5822             pass          # try again with complex

/home/patchbot/sage-patchbot/src/sage/symbolic/expression.pyx in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:10663)()
   1257         sig_on()
   1258         try:
-> 1259             res = self._gobj.evalf(0, kwds)
   1260         finally:
   1261             sig_off()

/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/functions/other.pyc in _evalf_(self, poly, index, parent, algorithm)
   2918         except AttributeError:
   2919             prec = 53
-> 2920         sobj = CRootOf(Poly(poly._sympy_()), int(index))
   2921         return sobj.n(ceil(prec*3/10))._sage_()
   2922 

/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sympy/polys/rootoftools.pyc in __new__(cls, f, x, index, radicals, expand)
    130         if index < -degree or index >= degree:
    131             raise IndexError("root index out of [%d, %d] range, got %d" %
--> 132                              (-degree, degree - 1, index))
    133         elif index < 0:
    134             index += degree

IndexError: root index out of [-8, 7] range, got 8

When gmpy2 is installed, the index is even passed as complex number:

sage -t --long src/sage/functions/other.py
**********************************************************************
File "src/sage/functions/other.py", line 2862, in sage.functions.other.Function_crootof
Failed example:
    c.n()
Exception raised:
    Traceback (most recent call last):
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 517, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 920, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.functions.other.Function_crootof[1]>", line 1, in <module>
        c.n()
      File "sage/structure/element.pyx", line 865, in sage.structure.element.Element.n (build/cythonized/sage/structure/element.c:8131)
        return self.numerical_approx(prec, digits, algorithm)
      File "sage/symbolic/expression.pyx", line 5824, in sage.symbolic.expression.Expression.numerical_approx (build/cythonized/sage/symbolic/expression.cpp:36286)
        x = x._convert(kwds)
      File "sage/symbolic/expression.pyx", line 1259, in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:10663)
        res = self._gobj.evalf(0, kwds)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/functions/other.py", line 2920, in _evalf_
        sobj = CRootOf(Poly(poly._sympy_()), int(index))
      File "sage/rings/complex_number.pyx", line 1058, in sage.rings.complex_number.ComplexNumber.__int__ (build/cythonized/sage/rings/complex_number.c:10918)
        raise TypeError("can't convert complex to int; use int(abs(z))")
    TypeError: can't convert complex to int; use int(abs(z))
**********************************************************************
File "src/sage/functions/other.py", line 2864, in sage.functions.other.Function_crootof
Failed example:
    c.n(100)
Exception raised:
    Traceback (most recent call last):
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 517, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 920, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.functions.other.Function_crootof[2]>", line 1, in <module>
        c.n(Integer(100))
      File "sage/structure/element.pyx", line 865, in sage.structure.element.Element.n (build/cythonized/sage/structure/element.c:8131)
        return self.numerical_approx(prec, digits, algorithm)
      File "sage/symbolic/expression.pyx", line 5824, in sage.symbolic.expression.Expression.numerical_approx (build/cythonized/sage/symbolic/expression.cpp:36286)
        x = x._convert(kwds)
      File "sage/symbolic/expression.pyx", line 1259, in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:10663)
        res = self._gobj.evalf(0, kwds)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/functions/other.py", line 2920, in _evalf_
        sobj = CRootOf(Poly(poly._sympy_()), int(index))
      File "sage/rings/complex_number.pyx", line 1058, in sage.rings.complex_number.ComplexNumber.__int__ (build/cythonized/sage/rings/complex_number.c:10918)
        raise TypeError("can't convert complex to int; use int(abs(z))")
    TypeError: can't convert complex to int; use int(abs(z))
**********************************************************************
File "src/sage/functions/other.py", line 2866, in sage.functions.other.Function_crootof
Failed example:
    (c^6 + c + 1).n(100) < 1e-25
Exception raised:
    Traceback (most recent call last):
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 517, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 920, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.functions.other.Function_crootof[3]>", line 1, in <module>
        (c**Integer(6) + c + Integer(1)).n(Integer(100)) < RealNumber('1e-25')
      File "sage/structure/element.pyx", line 865, in sage.structure.element.Element.n (build/cythonized/sage/structure/element.c:8131)
        return self.numerical_approx(prec, digits, algorithm)
      File "sage/symbolic/expression.pyx", line 5824, in sage.symbolic.expression.Expression.numerical_approx (build/cythonized/sage/symbolic/expression.cpp:36286)
        x = x._convert(kwds)
      File "sage/symbolic/expression.pyx", line 1259, in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:10663)
        res = self._gobj.evalf(0, kwds)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/functions/other.py", line 2920, in _evalf_
        sobj = CRootOf(Poly(poly._sympy_()), int(index))
      File "sage/rings/complex_number.pyx", line 1058, in sage.rings.complex_number.ComplexNumber.__int__ (build/cythonized/sage/rings/complex_number.c:10918)
        raise TypeError("can't convert complex to int; use int(abs(z))")
    TypeError: can't convert complex to int; use int(abs(z))
**********************************************************************
File "src/sage/functions/other.py", line 2908, in sage.functions.other.Function_crootof._evalf_
Failed example:
    complex_root_of(x^2-2, 1).n()
Exception raised:
    Traceback (most recent call last):
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 517, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 920, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.functions.other.Function_crootof._evalf_[0]>", line 1, in <module>
        complex_root_of(x**Integer(2)-Integer(2), Integer(1)).n()
      File "sage/structure/element.pyx", line 865, in sage.structure.element.Element.n (build/cythonized/sage/structure/element.c:8131)
        return self.numerical_approx(prec, digits, algorithm)
      File "sage/symbolic/expression.pyx", line 5824, in sage.symbolic.expression.Expression.numerical_approx (build/cythonized/sage/symbolic/expression.cpp:36286)
        x = x._convert(kwds)
      File "sage/symbolic/expression.pyx", line 1259, in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:10663)
        res = self._gobj.evalf(0, kwds)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/functions/other.py", line 2920, in _evalf_
        sobj = CRootOf(Poly(poly._sympy_()), int(index))
      File "sage/rings/complex_number.pyx", line 1058, in sage.rings.complex_number.ComplexNumber.__int__ (build/cythonized/sage/rings/complex_number.c:10918)
        raise TypeError("can't convert complex to int; use int(abs(z))")
    TypeError: can't convert complex to int; use int(abs(z))
**********************************************************************
2 items had failures:
   3 of   5 in sage.functions.other.Function_crootof
   1 of   3 in sage.functions.other.Function_crootof._evalf_
    [622 tests, 4 failures, 17.27 s]  
sage -t --long src/sage/interfaces/sympy.py
**********************************************************************
File "src/sage/interfaces/sympy.py", line 617, in sage.interfaces.sympy._sympysage_crootof
Failed example:
    (sols[0]+1)._sage_().n()
Exception raised:
    Traceback (most recent call last):
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 517, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/doctest/forker.py", line 920, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.interfaces.sympy._sympysage_crootof[6]>", line 1, in <module>
        (sols[Integer(0)]+Integer(1))._sage_().n()
      File "sage/structure/element.pyx", line 865, in sage.structure.element.Element.n (build/cythonized/sage/structure/element.c:8131)
        return self.numerical_approx(prec, digits, algorithm)
      File "sage/symbolic/expression.pyx", line 5824, in sage.symbolic.expression.Expression.numerical_approx (build/cythonized/sage/symbolic/expression.cpp:36286)
        x = x._convert(kwds)
      File "sage/symbolic/expression.pyx", line 1259, in sage.symbolic.expression.Expression._convert (build/cythonized/sage/symbolic/expression.cpp:10663)
        res = self._gobj.evalf(0, kwds)
      File "/home/patchbot/sage-patchbot/local/lib/python2.7/site-packages/sage/functions/other.py", line 2920, in _evalf_
        sobj = CRootOf(Poly(poly._sympy_()), int(index))
      File "sage/rings/complex_number.pyx", line 1058, in sage.rings.complex_number.ComplexNumber.__int__ (build/cythonized/sage/rings/complex_number.c:10918)
        raise TypeError("can't convert complex to int; use int(abs(z))")
    TypeError: can't convert complex to int; use int(abs(z))
**********************************************************************
1 item had failures:
   1 of   8 in sage.interfaces.sympy._sympysage_crootof
    [218 tests, 1 failure, 3.27 s]

Change History (54)

comment:1 Changed 5 years ago by Jeroen Demeyer

Cc: Emmanuel Charpentier added

comment:2 Changed 5 years ago by Jeroen Demeyer

Description: modified (diff)
Summary: Doctest failures related to sympyDoctest failures related to sympy CRootOf

comment:3 Changed 5 years ago by Ralf Stephan

Branch: u/rws/doctest_failures_related_to_sympy_crootof

comment:4 Changed 5 years ago by Ralf Stephan

Authors: Ralf Stephan
Commit: a93662faed523eb117f4ce0197b3178855849fd2
Status: newneeds_review

That exact arguments are made inexact is a bug in itself. But that they are made complex (EDIT: in some machines) is beyond me.


New commits:

a93662f24378: Doctest failures related to sympy
Last edited 5 years ago by Ralf Stephan (previous) (diff)

comment:5 Changed 5 years ago by Jeroen Demeyer

Do you have an idea why this depends on the system?

comment:6 in reply to:  4 Changed 5 years ago by Jeroen Demeyer

Replying to rws:

That exact arguments are made inexact is a bug in itself. But that they are made complex (EDIT: in some machines) is beyond me.

Do you know where this happens? It seems to me that it is Pynac doing this.

comment:7 Changed 5 years ago by Jeroen Demeyer

Description: modified (diff)

comment:8 Changed 5 years ago by Vincent Delecroix

Isn't it related to presence/absence of gmpy2 (that has been upgraded in 8.2.beta0)? Some sympy behavior depends on gmpy2.

Version 0, edited 5 years ago by Vincent Delecroix (next)

comment:9 Changed 5 years ago by Jeroen Demeyer

Description: modified (diff)
Summary: Doctest failures related to sympy CRootOfcomplex_root_of uses inexact index

comment:10 in reply to:  4 ; Changed 5 years ago by Jeroen Demeyer

Replying to rws:

That exact arguments are made inexact is a bug in itself.

Right, this ticket should just fix that.

comment:11 Changed 5 years ago by Jeroen Demeyer

Status: needs_reviewneeds_work

comment:12 Changed 5 years ago by Jeroen Demeyer

Description: modified (diff)

comment:13 Changed 5 years ago by Jeroen Demeyer

Description: modified (diff)

comment:14 Changed 5 years ago by Jeroen Demeyer

Description: modified (diff)

comment:15 in reply to:  10 Changed 5 years ago by Vincent Delecroix

Replying to jdemeyer:

Replying to rws:

That exact arguments are made inexact is a bug in itself.

Right, this ticket should just fix that.

Indeed numerical_approximation of symbolic expression is lazily converting all arguments to the same parent. Which is very wrong when a function is defined on R x N.

Perhaps it would be good to add in _evalf_ the check

if not isinstance(index, numbers.Integral):
    raise ValueError

Though I am curious about this sympy misbehavior.

Last edited 5 years ago by Vincent Delecroix (previous) (diff)

comment:16 Changed 5 years ago by Vincent Delecroix

Indeed there is a sympy bug that is now tracked at #24380.

comment:17 Changed 5 years ago by Vincent Delecroix

numerical_approx is very wrong when evaluating. It converts x^2 - 2 to x^2 - 2.0 and hence loosing precision for no good reason.

comment:18 Changed 5 years ago by Vincent Delecroix

There was indeed a bug in sympy: #24380 needs review!

comment:19 Changed 5 years ago by Ralf Stephan

I suggest accepting my quick fix to resolve the blocker. Meanwhile I'm looking into Pynac.

comment:20 Changed 5 years ago by Jeroen Demeyer

If your only goal is to resolve the blocker doctest failures, then I think it is better to use # known bug. The fact that the index is made floating-point is the real bug, whether it is real or complex is just a small variation.

comment:21 in reply to:  17 ; Changed 5 years ago by Ralf Stephan

Replying to vdelecroix:

numerical_approx is very wrong when evaluating. It converts x^2 - 2 to x^2 - 2.0 and hence loosing precision for no good reason.

I disagree in this case. The result is exactly what was asked for.

comment:22 in reply to:  21 ; Changed 5 years ago by Vincent Delecroix

Replying to rws:

Replying to vdelecroix:

numerical_approx is very wrong when evaluating. It converts x^2 - 2 to x^2 - 2.0 and hence loosing precision for no good reason.

I disagree in this case. The result is exactly what was asked for.

The problem under discussion is that you have a number specified by an equation like x^2 - 2 = 0. Now, you want an numerical approximation of it. I claim that there is not much reason to transform the equation into 1.0 * x^2 - 2.0 for doing that. In particular, this is not what the user wants. This kind of conversion is a can of worms. If you want 100 bits of your number, and start by converting the coefficients to floating points with 100 bits of precision you might loose your roots.

Note that in sympy they solve the original equation with appropriate root finding algorithms on rational numbers.

comment:23 in reply to:  22 Changed 5 years ago by Ralf Stephan

Replying to vdelecroix:

Replying to rws:

Replying to vdelecroix:

numerical_approx is very wrong when evaluating. It converts x^2 - 2 to x^2 - 2.0 and hence loosing precision for no good reason.

I disagree in this case. The result is exactly what was asked for.

The problem under discussion is that you have a number specified by an equation like x^2 - 2 = 0. Now, you want an numerical approximation of it.

No, that is solving. N() just converts symbolic expressions---if they happen to contain symbols you won't get a number, you get an error. _convert() helps with that.

Note that in sympy they solve the original equation with appropriate root finding algorithms on rational numbers.

This has nothing to do with this ticket. This ticket is about losing precision in expressions that can be FP evaluated because they contain no symbols but may contain functions. While the result precision from the function at hand complex_root_of is in principle not affected by the exactness of the index there may be existing/future functions where it is. This can be fixed by not FP converting exact function arguments and letting the function itself (i.e. its evalf member) decide what to do with exact arguments. I may well be wrong, and the best solution instead is a type database for function arguments...

comment:24 in reply to:  22 Changed 5 years ago by Ralf Stephan

Replying to vdelecroix:

Note that in sympy they solve the original equation with appropriate root finding algorithms on rational numbers.

Yes but we don't do the solving. We let SymPy do it.

comment:25 Changed 5 years ago by Ralf Stephan

Branch: u/rws/doctest_failures_related_to_sympy_crootofu/rws/24378

comment:26 Changed 5 years ago by Ralf Stephan

Branch: u/rws/24378u/rws/doctest_failures_related_to_sympy_crootof
Status: needs_workneeds_review

comment:27 Changed 5 years ago by Ralf Stephan

Branch: u/rws/doctest_failures_related_to_sympy_crootofu/rws/24378
Commit: a93662faed523eb117f4ce0197b3178855849fd2ac7681458cdeb8b6118080f81f3d9601db485021

New commits:

ac7681424378: avoid doctest fails

comment:28 Changed 5 years ago by Ralf Stephan

See #24397 for the fix.

comment:29 Changed 5 years ago by Jeroen Demeyer

Status: needs_reviewneeds_work

You should add the ticket number: # known bug (Trac #24397).

comment:30 Changed 5 years ago by git

Commit: ac7681458cdeb8b6118080f81f3d9601db48502149c6d614b4293059911f19ccedb87a3c5e96df35

Branch pushed to git repo; I updated commit sha1. New commits:

49c6d6124378: cosmetics

comment:31 Changed 5 years ago by Ralf Stephan

Status: needs_workneeds_review

comment:32 Changed 5 years ago by Ralf Stephan

Branch: u/rws/24378u/rws/24378-1

comment:33 Changed 5 years ago by Ralf Stephan

Commit: 49c6d614b4293059911f19ccedb87a3c5e96df35ae3a1e0be21a0aecbb46791fc2cacf84495eee5c

Fortunately this can be resolved easily.


New commits:

ae3a1e024387: complex_root_of uses inexact index

comment:34 Changed 5 years ago by Ralf Stephan

I have added the information to https://github.com/pynac/pynac/wiki/%7C-floating-point-evaluation, but Sage Function initialization options need better docs too.

comment:35 Changed 5 years ago by Jeroen Demeyer

Status: needs_reviewneeds_work

If this fixes the actual underlying issue, could you add

sage: complex_root_of(x^8-1, 7).n(2)

as a doctest?

comment:36 Changed 5 years ago by Ralf Stephan

It fixes the issue but there are other problems at super-low precision.

comment:37 Changed 5 years ago by git

Commit: ae3a1e0be21a0aecbb46791fc2cacf84495eee5c8c3f951c7b5f726260ccadb74368713d90740d93

Branch pushed to git repo; I updated commit sha1. New commits:

8c3f95124378: set minimum digits; convert to given prec

comment:38 Changed 5 years ago by Ralf Stephan

Status: needs_workneeds_review

comment:39 Changed 5 years ago by Jeroen Demeyer

Reviewers: Jeroen Demeyer
Status: needs_reviewneeds_work

The expression prec*3/10 will behave differently on Python 2 and Python 3. You should either use // or from __future__ import division such that it works the same on Python 2 and Python 3.

comment:40 Changed 5 years ago by Jeroen Demeyer

What is the meaning of the factor 3/10 anyway?

comment:41 in reply to:  40 ; Changed 5 years ago by Ralf Stephan

Replying to jdemeyer:

What is the meaning of the factor 3/10 anyway?

An approximation to log(2,10) (EDITED) which is good enough until about 1000 bits precision.

Last edited 5 years ago by Ralf Stephan (previous) (diff)

comment:42 in reply to:  41 Changed 5 years ago by Jeroen Demeyer

Replying to rws:

Replying to jdemeyer:

What is the meaning of the factor 3/10 anyway?

An approximation to log(2,10) (EDITED) which is good enough until about 1000 bits precision.

Then why don't use the actual value 0.3010299956639812? That avoids the problem with division also.

comment:43 Changed 5 years ago by git

Commit: 8c3f951c7b5f726260ccadb74368713d90740d938e65608d35c2527bcf6a7ebba8eadc7eb0377c38

Branch pushed to git repo; I updated commit sha1. New commits:

8e6560824378: don't prematurely optimize

comment:44 Changed 5 years ago by Ralf Stephan

Status: needs_workneeds_review

comment:45 Changed 5 years ago by Ralf Stephan

Ping?

comment:46 Changed 5 years ago by Vincent Delecroix

ceil(prec*log(2,10))) is terrible. Why not using 0.3010299956639812 as suggested in comment:42? A possible alternative would be

from math import log10
int(prec * log10(2))

You definitely do not want to use two symbolic functions here.

comment:47 Changed 5 years ago by git

Commit: 8e65608d35c2527bcf6a7ebba8eadc7eb0377c3814ed7c6f210f4e809c3458dc54315542bbf3547b

Branch pushed to git repo; I updated commit sha1. New commits:

14ed7c624378: fix previous commit

comment:48 Changed 5 years ago by Jeroen Demeyer

Or better yet, use the sympy function prec_to_dps which is meant to convert bits to digits:

from sympy.core.evalf import prec_to_dps

comment:49 Changed 5 years ago by Ralf Stephan

Thanks both. I'll make a fresh branch.

comment:50 Changed 5 years ago by Ralf Stephan

Branch: u/rws/24378-1u/rws/24378-2

comment:51 Changed 5 years ago by Ralf Stephan

Commit: 14ed7c6f210f4e809c3458dc54315542bbf3547b277a403ef00692247f9853d89322f9a464abd471

A minor note is that the results from prec_dps(n) are not identical to the previously proposed formulae at small values, even with a suitable shift.


New commits:

277a40324378: complex_root_of uses inexact index

comment:52 Changed 5 years ago by Jeroen Demeyer

Good for me if this passes patchbot testing.

comment:53 Changed 5 years ago by Ralf Stephan

Status: needs_reviewpositive_review

comment:54 Changed 5 years ago by Volker Braun

Branch: u/rws/24378-2277a403ef00692247f9853d89322f9a464abd471
Resolution: fixed
Status: positive_reviewclosed
Note: See TracTickets for help on using tickets.