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
comment:2 Changed 3 years ago by
(I'll review when ready)
comment:3 Changed 3 years ago by
- Branch set to u/rws/move_roots_solving_of_univariate_sr_polynomial_ring
comment:4 Changed 3 years ago by
- 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:
fcd9773 | 24942: move roots solving of univariate SR polynomial ring
|
comment:5 Changed 3 years ago by
- Commit changed from fcd9773e97ee97dce435316b0983ebd510938f6c to d1c24497ada050999931d60ecb26e499990b2ee6
Branch pushed to git repo; I updated commit sha1. New commits:
d1c2449 | 24942: fix roots() bug
|
comment:6 Changed 3 years ago by
- Status changed from new to needs_review
Found and fixed a bug in roots(). Please review.
comment:7 Changed 3 years ago by
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
- Branch changed from u/rws/move_roots_solving_of_univariate_sr_polynomial_ring to u/rws/24942
comment:9 Changed 3 years ago by
- Commit changed from d1c24497ada050999931d60ecb26e499990b2ee6 to 8908c9f9f67e48b0718b2eb46d666fde19ab2fdc
Still missing fixes for QQbar doctest fails.
New commits:
8908c9f | 24942: move roots solving of univariate SR polynomial ring
|
comment:10 Changed 3 years ago by
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
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: ↓ 13 ↓ 31 Changed 3 years ago by
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: ↓ 16 Changed 3 years ago by
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.
comment:14 Changed 3 years ago by
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
- Commit changed from 8908c9f9f67e48b0718b2eb46d666fde19ab2fdc to 0c791d057cc5bac4c2d27dd88c385b728445ee18
Branch pushed to git repo; I updated commit sha1. New commits:
0c791d0 | 24942: fixes
|
comment:16 in reply to: ↑ 13 Changed 3 years ago by
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: ↓ 18 Changed 3 years ago by
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: ↓ 19 Changed 3 years ago by
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: ↓ 20 Changed 3 years ago by
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: ↓ 21 Changed 3 years ago by
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: ↓ 22 Changed 3 years ago by
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: ↓ 27 Changed 3 years ago by
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
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
Note that speeding up SR(a)
with #24952 does not affect change_ring(SR)
timing.
comment:25 follow-up: ↓ 26 Changed 3 years ago by
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
Replying to rws:
p.change_ring(SR)
uses the genericpolynomial_element.pyx:_call_
callingsubs()
whilep.change_ring(RR)
uses the genericParent.__call__
. It seems from profiling that thousands of dict queries are done for a single call top.change_ring(SR)
. Either way, implementing a dedicatedpolynomial_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: ↓ 28 Changed 3 years ago by
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: ↓ 29 Changed 3 years ago by
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 loopThe biggest part (230µs) of these 300µs is spent in
Expression.__nonzero__
disprovingx^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 insertingreturn False
at the start ofExpression.__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
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
Please see #24965, which we might depend on here.
comment:31 in reply to: ↑ 12 Changed 3 years ago by
comment:32 follow-up: ↓ 33 Changed 3 years ago by
- 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
- 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
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: ↓ 38 Changed 3 years ago by
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: ↓ 37 Changed 3 years ago by
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
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: ↓ 39 Changed 3 years ago by
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: ↓ 40 Changed 3 years ago by
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: ↓ 42 Changed 3 years ago by
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): 199 199 if ComplexField(mpfr_prec_min()).has_coerce_map_from(R): 200 200 # Almost anything with a coercion into any precision of CC 201 201 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)205 202 elif (R is InfinityRing or R is UnsignedInfinityRing 206 203 or is_RealIntervalField(R) or is_ComplexIntervalField(R) 207 204 or isinstance(R, RealBallField)
comment:41 Changed 3 years ago by
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
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
- Dependencies set to #25022
- Status changed from needs_review to needs_work
Simpler solution is to fix change_ring
with #24942.
I'm on it.