move roots solving of univariate SR polynomial ring

[vdelecroix]

Move the root finding code of polynomial ring over SR from the generic method roots in rings/polynomial/polynomial_element.pyx to the specialized _roots_univariate_polynomial that has to be put directly in the base ring SR (code in symbolic/ring.pyx).

[rws]

I'm on it.

[vdelecroix]

(I'll review when ready)

[rws]

Although this is a straight move, with addition of two examples, now 3 doctests fail:

File "src/sage/symbolic/ring.pyx", line 227, in sage.symbolic.ring.SymbolicRing._element_constructor_
Failed example:
    a + sin(x)
Expected:
    I*sqrt(3) + sin(x)
Got:
    sin(x) + 1.732050807568878?*I

File "src/sage/rings/polynomial/polynomial_element.pyx", line 7038, in sage.rings.polynomial.polynomial_element.Polynomial.roots
Failed example:
    f.roots(SR)
Expected:
    [(-I*sqrt(2), 1), (I*sqrt(2), 1)]
Got:
    []

File "src/sage/rings/polynomial/polynomial_element.pyx", line 7040, in sage.rings.polynomial.polynomial_element.Polynomial.roots
Failed example:
    f.roots(SR, multiplicities=False)
Expected:
    [-I*sqrt(2), I*sqrt(2)]
Got:
    []

There must be some fuzziness on if SR is the base ring, or the ring argument of roots().

---

New commits:

​fcd9773

24942: move roots solving of univariate SR polynomial ring

[git]

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

​d1c2449

24942: fix roots() bug

[rws]

Found and fixed a bug in roots(). Please review.

[vdelecroix]

The following is wrong

-        if hasattr(K, '_roots_univariate_polynomial'):
+        if hasattr(L, '_roots_univariate_polynomial'):

You do want to delegate to the base ring of the polynomial (not to the target for the roots). You can check that at the end of the roots code there is

return self.change_ring(L).roots(multiplicities=multiplicities, algorithm=algorithm)

[rws]

Still missing fixes for QQbar doctest fails.

---

New commits:

​8908c9f

24942: move roots solving of univariate SR polynomial ring

[vdelecroix]

I am not sure we want to support the ring argument when different from None/SR. Do you? If so, you should add proper examples.

BTW, in the doctest (x^2 + x + 1).roots(SR) the SR should not be needed.

[rws]

There may be some serious inefficiencies involved up to now because the polys that get transferred in the call to SR._roots... have coefficients in SR BUT they are actually embedded algebraic numbers. I need to have a further look.

[rws]

Agree with ignoring the ring argument. Here is a problem:

sage: QQbar(sqrt(2)).minpoly().change_ring(SR)
x^2 - 2
sage: _.coefficients()
[x^2 - 2]

That is, the result from change_ring(SR) has degree 0. Is this expected?

[vdelecroix]

Replying to rws:

Agree with ignoring the ring argument. Here is a problem:

sage: QQbar(sqrt(2)).minpoly().change_ring(SR)
x^2 - 2
sage: _.coefficients()
[x^2 - 2]

That is, the result from change_ring(SR) has degree 0. Is this expected?

Another bug (mostly disjoint from this ticket)

sage: PSR = SR['x']
sage: PSR('x^2 - 1').degree()   # fine
2
sage: PSR([-1, 0, 1]).degree()  # fine
2
sage: PSR(polygen(ZZ)^2 - 1).degree()  # weird!
0

I guess that the _element_constructor of PolynomialRing has to be blame here. The point is that SR is able to merge x^2 - 1 in its base ring so the the constructor get confused. We should special case when the input is a univariate polynomial. But please in another ticket.

[rws]

I think I can work around it, and the following:

sage: R.<x> = SR[]
sage: SR(x^2 + x + 1)
...
/home/ralf/sage/local/lib/python2.7/site-packages/sage/structure/parent.pyx in sage.structure.parent.Parent.__call__ (build/cythonized/sage/structure/parent.c:9671)()
    918         if mor is not None:
    919             if no_extra_args:
--> 920                 return mor._call_(x)
    921             else:
    922                 return mor._call_with_args(x, args, kwds)

/home/ralf/sage/local/lib/python2.7/site-packages/sage/rings/polynomial/polynomial_element.pyx in sage.rings.polynomial.polynomial_element.ConstantPolynomialSection._call_ (build/cythonized/sage/rings/polynomial/polynomial_element.c:102674)()
  11117                 return <Element>((<Polynomial>x).constant_coefficient())
  11118         else:
> 11119             raise TypeError("not a constant polynomial")
  11120 
  11121 cdef class PolynomialBaseringInjection(Morphism):

TypeError: not a constant polynomial

[git]

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

​0c791d0

24942: fixes

[rws]

Replying to vdelecroix:

Another bug (mostly disjoint from this ticket)

Not disjoint because I need a special case to work with these polynomials that get thrown at the _roots...() function. Please have a look. Tests pass so far.

[rws]

Note there is a performance hit, QQbar(I).radical_expression() now takes 4-5x the time (mostly because code in polynomial_elements.pyx:roots() is executed that was cut short before (if L is SR is now missing).

[vdelecroix]

Replying to rws:

Note there is a performance hit, QQbar(I).radical_expression() now takes 4-5x the time (mostly because code in polynomial_elements.pyx:roots() is executed that was cut short before (if L is SR is now missing).

The roots code first calls _roots_univariate_polynomial. So there must be something wrong.

[rws]

Replying to vdelecroix:

Replying to rws:

Note there is a performance hit, QQbar(I).radical_expression() now takes 4-5x the time (mostly because code in polynomial_elements.pyx:roots() is executed that was cut short before (if L is SR is now missing).

The roots code first calls _roots_univariate_polynomial. So there must be something wrong.

No, the result from QQbar(I).minpoly() does not have a _roots_univariate_polynomial member.

[vdelecroix]

Replying to rws:

Replying to vdelecroix:

Replying to rws:

Note there is a performance hit, QQbar(I).radical_expression() now takes 4-5x the time (mostly because code in polynomial_elements.pyx:roots() is executed that was cut short before (if L is SR is now missing).

The roots code first calls _roots_univariate_polynomial. So there must be something wrong.

No, the result from QQbar(I).minpoly() does not have a _roots_univariate_polynomial member.

_roots_univariate_polynomial is a method of the base ring not of the polynomial. That being said, the code in radical_expression would better be changed to

- roots = poly.roots(SR, multiplicities=False)
+ roots = poly.change_ring(SR).roots(multiplicities=False)

[rws]

Replying to vdelecroix: That being said, the code in radical_expression would better be changed to

- roots = poly.roots(SR, multiplicities=False)
+ roots = poly.change_ring(SR).roots(multiplicities=False)

Is nearly as slow, despite being more correct. I tried converting to SR, or calling _symbolic_ on the poly but it's even slower. If these are Flint polys, I could put code in Pynac for fast conversion.

[vdelecroix]

Replying to rws:

Replying to vdelecroix: That being said, the code in radical_expression would better be changed to

- roots = poly.roots(SR, multiplicities=False)
+ roots = poly.change_ring(SR).roots(multiplicities=False)

Is nearly as slow, despite being more correct. I tried converting to SR, or calling _symbolic_ on the poly but it's even slower. If these are Flint polys, I could put code in Pynac for fast conversion.

Conversion from ZZ can already be improved

sage: a = ZZ(23)
sage: %timeit SR(a)
100000 loops, best of 3: 3.51 µs per loop
sage: %timeit RR(a)
1000000 loops, best of 3: 218 ns per loop

which is basically the same x15 ratio as in

sage: p = polygen(ZZ)^2 + 1
sage: %timeit p.change_ring(SR)
1000 loops, best of 3: 322 µs per loop
sage: %timeit p.change_ring(RR)
100000 loops, best of 3: 12.9 µs per loop

[rws]

I'm able to reduce SR(a) to 300ns by simply rearranging the code in SR._element_constructor(), see #24952. I'm fighting for my life with Cython to get this implemented however.

[rws]

Note that speeding up SR(a) with #24952 does not affect change_ring(SR) timing.

[rws]

p.change_ring(SR) uses the generic polynomial_element.pyx:_call_ calling subs() while p.change_ring(RR) uses the generic Parent.__call__. It seems from profiling that thousands of dict queries are done for a single call to p.change_ring(SR). Either way, implementing a dedicated polynomial_integer_dense_flint.pyx:_symbolic_() seems to me the optimal way to deal with this.

[vdelecroix]

Replying to rws:

p.change_ring(SR) uses the generic polynomial_element.pyx:_call_ calling subs() while p.change_ring(RR) uses the generic Parent.__call__. It seems from profiling that thousands of dict queries are done for a single call to p.change_ring(SR). Either way, implementing a dedicated polynomial_integer_dense_flint.pyx:_symbolic_() seems to me the optimal way to deal with this.

These are different operations. You seem to confuse polynomials with coefficients in SR like SR['x'].gen()^2 + 1 and the corresponding element of SR SR.var('x')^2 + 1.

If you have a polynomial p, the operations SR(p) and p.change_ring(SR) are different. In the second case, you change the ring of the coefficients. Not the fact that it is a Sage polynomial.

[rws]

Replying to vdelecroix:

sage: p = polygen(ZZ)^2 + 1
sage: %timeit p.change_ring(SR)
1000 loops, best of 3: 322 µs per loop
sage: %timeit p.change_ring(RR)
100000 loops, best of 3: 12.9 µs per loop

The biggest part (230µs) of these 300µs is spent in Expression.__nonzero__ disproving x^2 + 1 == 0. So again, it's calling the wrong function when comparing an expression with zero, that leads to performance issues. You can easily confirm this by inserting return False at the start of Expression.__nonzero__ which brings the time down to 70µs.

[vdelecroix]

Replying to rws:

Replying to vdelecroix:

sage: p = polygen(ZZ)^2 + 1
sage: %timeit p.change_ring(SR)
1000 loops, best of 3: 322 µs per loop
sage: %timeit p.change_ring(RR)
100000 loops, best of 3: 12.9 µs per loop

The biggest part (230µs) of these 300µs is spent in Expression.__nonzero__ disproving x^2 + 1 == 0. So again, it's calling the wrong function when comparing an expression with zero, that leads to performance issues. You can easily confirm this by inserting return False at the start of Expression.__nonzero__ which brings the time down to 70µs.

Why on earth the code would you check that x^2 + 1 == 0!? You just need to consider the coefficients, that is 1, 0 and 1. Have you read 25?

[rws]

Replying to vdelecroix:

Have you read 25?

I have, and I have traced what happens. I will document this clearly and confirmable, but this takes time.

[rws]

Please see #24965, which we might depend on here.

[rws]

Replying to rws:

That is, the result from change_ring(SR) has degree 0. Is this expected?

Note that this is the reason x^2 + 1 is tested as documented in #24965.

[rws]

It seems I mentally did put your comment:13 in a different slot than comment:12. Anyway, we should depend on the coming fix in #24965 and revert the commit 0c791d0 here.

[rws]

Replying to rws:

Anyway, we should depend on the coming fix in #24965 and revert the commit 0c791d0 here.

This was a cul-de-sac. So the alternative, doing something about the comparison of objects with zero (where the comparison explicitly needs no simplification) looks like an alternative. The problem is outlined in ​https://trac.sagemath.org/wiki/symbolics/nonzero and #21201 looks like the best start at it.

[rws]

With #21201 and patching polynomial_element.pyx to use special comparison I can get the ring change of comment:22 down to 90us.

[vdelecroix]

Actually, moving the code to SR._roots_univariate_polynomial is nearly impossible for the following reason. As you can see the code that is now in roots uses a silent change of variables with name vname = 'do_not_use_this_name_in_a_polynomial_coefficient'. The reason for this is that you want to support

sage: x = polygen(SR)
sage: (x^3 - SR.var('x')^3).roots()
[(1/2*x*(I*sqrt(3) - 1), 1), (-1/2*x*(I*sqrt(3) + 1), 1), (x, 1)]

Now if you convert first to SR as a symbolic expression you get 0. Which is certainly not what you want.

My conclusion (for now): SR is terribly intrusive and does not play nicely with anything else from Sage.

[rws]

The only intrusion I see is in roots, and well, if you request symbolic output then you need a special case. QQbar could simply convert to SR, and call solve on the resulting expression.

[vdelecroix]

Replying to rws:

The only intrusion I see is in roots, and well, if you request symbolic output then you need a special case. QQbar could simply convert to SR, and call solve on the resulting expression.

It is in roots and you can not remove it. It is terribly intrusive. All of the following works correctly

sage: x = polygen(ZZ)
sage: (x^3 - 2).roots(QQ)
[]
sage: (x^3 - 2).roots(RR)
[(1.25992104989487, 1)]
sage: (x^3 - 2).roots(CC)
[(1.25992104989487, 1),
 (-0.629960524947437 - 1.09112363597172*I, 1),
 (-0.629960524947437 + 1.09112363597172*I, 1)]
sage: (x^3 - 2).roots(AA)
[(1.259921049894873?, 1)]
sage: (x^3 - 2).roots(QQbar)
[(1.259921049894873?, 1),
 (-0.6299605249474365? - 1.091123635971722?*I, 1),
 (-0.6299605249474365? + 1.091123635971722?*I, 1)]

and there is no special casing for them.

[vdelecroix]

Replying to vdelecroix:

My conclusion (for now): SR is terribly intrusive and does not play nicely with anything else from Sage.

I should have added: unless we go for my solution one in #24965.

[rws]

Replying to vdelecroix:

Replying to vdelecroix:

My conclusion (for now): SR is terribly intrusive and does not play nicely with anything else from Sage.

I should have added: unless we go for my solution one in #24965.

I have no problems with that, but can't contribute much.

[vdelecroix]

Replying to rws:

Replying to vdelecroix:

Replying to vdelecroix:

My conclusion (for now): SR is terribly intrusive and does not play nicely with anything else from Sage.

I should have added: unless we go for my solution one in #24965.

I have no problems with that, but can't contribute much.

You can at least give your opinion. "My" solution consists in removing three lines of code in SymbolicRing._coerce_map_from_

src/sage/symbolic/ring.pyx

@@ -199 +199 @@
             if ComplexField(mpfr_prec_min()).has_coerce_map_from(R):
                 # Almost anything with a coercion into any precision of CC
                 return R not in (RLF, CLF)
-            elif is_PolynomialRing(R) or is_MPolynomialRing(R) or is_FractionField(R) or is_LaurentPolynomialRing(R):
-                base = R.base_ring()
-                return base is not self and self.has_coerce_map_from(base)
             elif (R is InfinityRing or R is UnsignedInfinityRing
                   or is_RealIntervalField(R) or is_ComplexIntervalField(R)
                   or isinstance(R, RealBallField)

[rws]

I see. The main consequence is that some symbolic functions that traditionally supported polynomials as arguments will need a separation into an interface taking all sorts of arguments, and the symbolic function class. Accidentally, for the orthogonal poly functions there is already a ticket that does that (#24554 needing review). For binomial there is already an interface (in arith) and a symbolic function, the only problem being that the import of the symbolic version overwrites the arith version on startup. There is #24178 for that but Jeroen opposed it. The split into global interface / symbolic function is inherently necessary because the symbolic function creation/call code restricts what is allowed as argument (even which kwds may be given), and I get no support to change that.

The other doctest failures I get are minor. All in all I support this change.

[rws]

Replying to vdelecroix:

... "My" solution consists in removing three lines of code in SymbolicRing._coerce_map_from_

Are you going to upload a branch with this? What's the plan?

[vdelecroix]

Simpler solution is to fix change_ring with #24942.