Opened 3 years ago

Last modified 3 years ago

#24942 needs_work defect

move roots solving of univariate SR polynomial ring

Reported by: vdelecroix Owned by:
Priority: major Milestone: sage-8.2
Component: symbolics Keywords:
Cc: rws Merged in:
Authors: Ralf Stephan Reviewers:
Report Upstream: N/A Work issues:
Branch: u/rws/24942 (Commits) Commit: 0c791d057cc5bac4c2d27dd88c385b728445ee18
Dependencies: #25022 Stopgaps:

Description

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).

Change History (43)

comment:1 Changed 3 years ago by rws

I'm on it.

comment:2 Changed 3 years ago by vdelecroix

(I'll review when ready)

comment:3 Changed 3 years ago by rws

  • Branch set to u/rws/move_roots_solving_of_univariate_sr_polynomial_ring

comment:4 Changed 3 years ago by rws

  • Commit set to fcd9773e97ee97dce435316b0983ebd510938f6c

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:

fcd977324942: move roots solving of univariate SR polynomial ring

comment:5 Changed 3 years ago by git

  • Commit changed from fcd9773e97ee97dce435316b0983ebd510938f6c to d1c24497ada050999931d60ecb26e499990b2ee6

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

d1c244924942: fix roots() bug

comment:6 Changed 3 years ago by rws

  • Authors set to Ralf Stephan
  • Status changed from new to needs_review

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

comment:7 Changed 3 years ago by 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)

comment:8 Changed 3 years ago by rws

  • Branch changed from u/rws/move_roots_solving_of_univariate_sr_polynomial_ring to u/rws/24942

comment:9 Changed 3 years ago by rws

  • Commit changed from d1c24497ada050999931d60ecb26e499990b2ee6 to 8908c9f9f67e48b0718b2eb46d666fde19ab2fdc

Still missing fixes for QQbar doctest fails.


New commits:

8908c9f24942: move roots solving of univariate SR polynomial ring

comment:10 Changed 3 years ago by 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.

comment:11 Changed 3 years ago by 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.

comment:12 follow-ups: Changed 3 years ago by 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?

comment:13 in reply to: ↑ 12 ; follow-up: Changed 3 years ago by 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.

Last edited 3 years ago by vdelecroix (previous) (diff)

comment:14 Changed 3 years ago by 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

comment:15 Changed 3 years ago by git

  • Commit changed from 8908c9f9f67e48b0718b2eb46d666fde19ab2fdc to 0c791d057cc5bac4c2d27dd88c385b728445ee18

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

0c791d024942: fixes

comment:16 in reply to: ↑ 13 Changed 3 years ago by 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.

comment:17 follow-up: Changed 3 years ago by 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).

comment:18 in reply to: ↑ 17 ; follow-up: Changed 3 years ago by 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.

comment:19 in reply to: ↑ 18 ; follow-up: Changed 3 years ago by 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.

comment:20 in reply to: ↑ 19 ; follow-up: Changed 3 years ago by 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)

comment:21 in reply to: ↑ 20 ; follow-up: Changed 3 years ago by 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.

comment:22 in reply to: ↑ 21 ; follow-up: Changed 3 years ago by 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

comment:23 Changed 3 years ago by 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.

comment:24 Changed 3 years ago by rws

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

comment:25 follow-up: Changed 3 years ago by 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.

comment:26 in reply to: ↑ 25 Changed 3 years ago by 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.

comment:27 in reply to: ↑ 22 ; follow-up: Changed 3 years ago by 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.

comment:28 in reply to: ↑ 27 ; follow-up: Changed 3 years ago by 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?

comment:29 in reply to: ↑ 28 Changed 3 years ago by 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.

comment:30 Changed 3 years ago by rws

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

comment:31 in reply to: ↑ 12 Changed 3 years ago by 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.

comment:32 follow-up: Changed 3 years ago by rws

  • Dependencies set to #24965

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.

comment:33 in reply to: ↑ 32 Changed 3 years ago by rws

  • Dependencies #24965 deleted

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.

comment:34 Changed 3 years ago by rws

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

comment:35 follow-up: Changed 3 years ago by 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.

comment:36 follow-up: Changed 3 years ago by 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.

comment:37 in reply to: ↑ 36 Changed 3 years ago by 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.

comment:38 in reply to: ↑ 35 ; follow-up: Changed 3 years ago by 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.

comment:39 in reply to: ↑ 38 ; follow-up: Changed 3 years ago by 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.

comment:40 in reply to: ↑ 39 ; follow-up: Changed 3 years ago by 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

    a b cdef class SymbolicRing(CommutativeRing): 
    199199            if ComplexField(mpfr_prec_min()).has_coerce_map_from(R):
    200200                # Almost anything with a coercion into any precision of CC
    201201                return R not in (RLF, CLF)
    202             elif is_PolynomialRing(R) or is_MPolynomialRing(R) or is_FractionField(R) or is_LaurentPolynomialRing(R):
    203                 base = R.base_ring()
    204                 return base is not self and self.has_coerce_map_from(base)
    205202            elif (R is InfinityRing or R is UnsignedInfinityRing
    206203                  or is_RealIntervalField(R) or is_ComplexIntervalField(R)
    207204                  or isinstance(R, RealBallField)

comment:41 Changed 3 years ago by 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.

comment:42 in reply to: ↑ 40 Changed 3 years ago by 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?

comment:43 Changed 3 years ago by vdelecroix

  • Dependencies set to #25022
  • Status changed from needs_review to needs_work

Simpler solution is to fix change_ring with #24942.

Note: See TracTickets for help on using tickets.