Opened 3 years ago

# move roots solving of univariate SR polynomial ring

Reported by: Owned by: vdelecroix major sage-8.2 symbolics rws Ralf Stephan N/A u/rws/24942 (Commits) 0c791d057cc5bac4c2d27dd88c385b728445ee18 #25022

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

I'm on it.

### 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:

 ​fcd9773 `24942: 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:

 ​d1c2449 `24942: 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:

 ​8908c9f `24942: 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: ↓ 13 ↓ 31 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: ↓ 16 Changed 3 years ago by vdelecroix

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:

 ​0c791d0 `24942: fixes`

### comment:16 in reply to: ↑ 13 Changed 3 years ago by rws

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 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: ↓ 19 Changed 3 years ago by vdelecroix

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 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 vdelecroix

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 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: ↓ 27 Changed 3 years ago by vdelecroix

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: ↓ 26 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

`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: ↓ 28 Changed 3 years ago by rws

```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 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

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

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: ↓ 33 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

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: ↓ 38 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: ↓ 37 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

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 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 rws

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 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 cdef class SymbolicRing(CommutativeRing): 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)

### 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

... "My" solution consists in removing three lines of code in `SymbolicRing._coerce_map_from_`
Simpler solution is to fix `change_ring` with #24942.