Opened 7 years ago

Closed 3 years ago

# Computing all roots is faster than computing a single one

Reported by: Owned by: gagern major sage-8.9 number fields qqbar polynomial_root roots performance sbrandhorst Vincent Delecroix Simon Brandhorst N/A d5105b6 d5105b625060d4b2a9cfe9cee4316d8b7c96ea46 #28191

The following example comes from comment:27:ticket:16964 via comment:2:ticket:17886.

```sage: x,y = polygens(QQ,"x,y")
sage: p1 = x^5 + 6*x^4 - 42*x^3 - 142*x^2 + 467*x + 422
sage: p2 = p1(x=(x-1)^2)
sage: p3 = p2(x=x*y).resultant(p2,x).univariate_polynomial()
sage: p4, = [f[0] for f in p3.factor() if f[0].degree() == 80]
sage: ival = CIF((0.77, 0.78), (-0.08, -0.07))
sage: %time z1, = [r for r in p4.roots(QQbar, False) if r in ival]
CPU times: user 1.43 s, sys: 195 ms, total: 1.62 s
Wall time: 1.47 s
sage: %time z2 = QQbar.polynomial_root(p4, ival)
CPU times: user 1min 5s, sys: 212 ms, total: 1min 5s
Wall time: 1min 5s
```

The computation for `z1` works reasonably well and completes in under 2 seconds, but the one for `z2` takes over a minute. Which is definitely wrong, since computing all roots and then choosing the right one should be more work, not less than just computing a single one!

The reason for this is the time spent in the square-free decomposition (called from `sage/rings/polynomial/complex_roots.py`) which behave differently depending whether the polynomial is defined with coefficients in `QQ` or `AA`

```sage: %time _ = p4.squarefree_decomposition()
CPU times: user 807 µs, sys: 0 ns, total: 807 µs
Wall time: 883 µs
sage: %time _ = p4.change_ring(AA).squarefree_decomposition()
CPU times: user 40.1 s, sys: 3.21 ms, total: 40.1 s
Wall time: 40.1 s
```

(this somehow independent problem about gcd is tracked at #28199).

With the branch attached the timings become

```sage: sage: %time z1, = [r for r in p4.roots(QQbar, False) if r in ival]
CPU times: user 411 ms, sys: 152 µs, total: 411 ms
Wall time: 410 ms
sage: sage: %time z2 = QQbar.polynomial_root(p4, ival)
CPU times: user 294 ms, sys: 1 µs, total: 294 ms
Wall time: 296 ms
```

### comment:1 Changed 7 years ago by pernici

If the interval is close enough to the root, QQbar.polynomial_root` is fast

```sage: z1
0.7745442035157198? - 0.07841573936877289?*I
sage: ival = CIF((0.77454420351571, 0.77454420351572), (-0.078415739368773, -0.078415739368772))
sage: %time z2 = QQbar.polynomial_root(p4, ival); z2._more_precision()
Wall time: 26.7 ms
sage: z2
0.7745442035157198? - 0.07841573936877289?*I
```
Version 1, edited 7 years ago by pernici (previous) (next) (diff)

### comment:2 Changed 3 years ago by vdelecroix

For the record, here is what profiling shows

```   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
6320   23.769    0.004   23.769    0.004 {built-in method _operator.sub}
6481   17.266    0.003   17.266    0.003 {built-in method _operator.mul}
2    0.225    0.112    0.231    0.116 {method 'roots' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
16711    0.188    0.000    0.188    0.000 sets_cat.py:974(_element_constructor_from_element_class)
89    0.069    0.001    0.070    0.001 {sage.rings.polynomial.refine_root.refine_root}
15053    0.037    0.000    0.272    0.000 qqbar.py:3074(__init__)
15287    0.031    0.000    0.227    0.000 qqbar.py:5723(_interval_fast)
1    0.028    0.028   41.429   41.429 fields.py:218(_gcd_univariate_polynomial)
14889    0.024    0.000   41.096    0.003 qqbar.py:7736(an_binop_rational)
15052    0.024    0.000    0.036    0.000 qqbar.py:5653(__init__)
47328    0.022    0.000    0.022    0.000 {built-in method builtins.isinstance}
15052    0.017    0.000    0.297    0.000 qqbar.py:4712(__init__)
6481    0.014    0.000   17.451    0.003 qqbar.py:3240(_mul_)
6320    0.012    0.000   23.946    0.004 qqbar.py:3284(_sub_)
15449    0.009    0.000    0.009    0.000 {built-in method sage.rings.real_mpfi.RealIntervalField}
15287    0.006    0.000    0.008    0.000 qqbar.py:4730(_ensure_real)
2088    0.003    0.000    0.022    0.000 qqbar.py:3271(_add_)
15450    0.003    0.000    0.003    0.000 {sage.rings.complex_interval.is_ComplexIntervalFieldElement}
448/160    0.002    0.000    0.004    0.000 complex_field.py:389(_element_constructor_)
4    0.002    0.000    0.024    0.006 {method 'derivative' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
2088    0.002    0.000    0.002    0.000 {built-in method _operator.add}
484    0.001    0.000    0.005    0.000 qqbar.py:4023(interval)
1322/1000    0.001    0.000    0.005    0.000 complex_interval_field.py:427(__call__)
484    0.001    0.000    0.003    0.000 qqbar.py:3999(interval_diameter)
216    0.001    0.000    0.003    0.000 {method 'n' of 'sage.structure.element.Element' objects}
1    0.001    0.001   41.766   41.766 qqbar.py:6487(_complex_refine_interval)
4    0.001    0.000    0.005    0.001 {built-in method builtins.sorted}
81    0.001    0.000    0.001    0.000 qqbar.py:5809(invert)
628/324    0.001    0.000    0.005    0.000 complex_field.py:352(__call__)
561    0.001    0.000    0.001    0.000 {built-in method builtins.hasattr}
240    0.001    0.000    0.004    0.000 misc.py:5558(_key_complex_for_display)
144    0.001    0.000    0.001    0.000 {method '__copy__' of 'sage.categories.map.Map' objects}
719    0.000    0.000    0.000    0.000 {method 'diameter' of 'sage.rings.real_mpfi.RealIntervalFieldElement' objects}
15    0.000    0.000    0.007    0.000 polynomial_element_generic.py:1052(__init__)
161/155    0.000    0.000    0.001    0.000 homset.py:84(Hom)
312    0.000    0.000    0.001    0.000 copy.py:66(copy)
2    0.000    0.000    0.086    0.043 complex_roots.py:47(interval_roots)
235    0.000    0.000    0.001    0.000 qqbar.py:3905(_more_precision)
1    0.000    0.000   41.757   41.757 complex_roots.py:147(complex_roots)
162    0.000    0.000    0.001    0.000 qqbar.py:5244(real_number)
161    0.000    0.000    0.001    0.000 qqbar.py:3970(interval_fast)
1067    0.000    0.000    0.000    0.000 complex_field.py:289(_real_field)
11    0.000    0.000    0.000    0.000 homset.py:582(__init__)
81    0.000    0.000    0.004    0.000 qqbar.py:3262(__invert__)
235    0.000    0.000    0.001    0.000 qqbar.py:4751(_more_precision)
2    0.000    0.000    0.000    0.000 polynomial_ring.py:234(__init__)
469    0.000    0.000    0.001    0.000 <frozen importlib._bootstrap>:1009(_handle_fromlist)
85    0.000    0.000    0.000    0.000 qqbar.py:3386(__bool__)
2    0.000    0.000    0.000    0.000 complex_field.py:193(__init__)
82    0.000    0.000    0.001    0.000 qqbar.py:720(_element_constructor_)
671    0.000    0.000    0.000    0.000 complex_interval_field.py:266(prec)
160    0.000    0.000    0.000    0.000 {method 'real' of 'sage.rings.complex_number.ComplexNumber' objects}
1    0.000    0.000   41.766   41.766 qqbar.py:6320(refine_interval)
218    0.000    0.000    0.001    0.000 complex_field.py:91(ComplexField)
80    0.000    0.000    0.003    0.000 fields.py:767(inverse_of_unit)
1    0.000    0.000    0.001    0.001 {method 'change_ring' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
216    0.000    0.000    0.000    0.000 {method 'abs' of 'sage.structure.element.RingElement' objects}
6    0.000    0.000    0.001    0.000 polynomial_ring_constructor.py:52(PolynomialRing)
1    0.000    0.000   41.436   41.436 fields.py:441(_squarefree_decomposition_univariate_polynomial)
160    0.000    0.000    0.000    0.000 {method 'imag' of 'sage.rings.complex_number.ComplexNumber' objects}
1    0.000    0.000    0.001    0.001 complex_roots.py:85(intervals_disjoint)
17    0.000    0.000    0.000    0.000 dynamic_class.py:127(dynamic_class)
160    0.000    0.000    0.000    0.000 {method 'real' of 'sage.rings.complex_interval.ComplexIntervalFieldElement' objects}
1    0.000    0.000    0.000    0.000 {method 'monic' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
9    0.000    0.000    0.008    0.001 polynomial_ring.py:320(_element_constructor_)
80    0.000    0.000    0.000    0.000 complex_roots.py:118(row_disjoint)
6    0.000    0.000    0.000    0.000 polynomial_ring_constructor.py:677(_single_variate)
1    0.000    0.000    0.000    0.000 dynamic_class.py:346(dynamic_class_internal)
6/5    0.000    0.000    0.000    0.000 {method 'has_coerce_map_from' of 'sage.structure.parent.Parent' objects}
496    0.000    0.000    0.000    0.000 complex_interval_field.py:255(is_exact)
2    0.000    0.000    0.000    0.000 category.py:499(__repr_object_names)
91    0.000    0.000    0.000    0.000 weakref.py:341(__init__)
6    0.000    0.000    0.000    0.000 homset.py:40(RingHomset)
468    0.000    0.000    0.000    0.000 {method 'precision' of 'sage.rings.real_mpfi.RealIntervalFieldElement' objects}
8    0.000    0.000    0.000    0.000 rings.py:322(_Hom_)
1    0.000    0.000   41.757   41.757 qqbar.py:6060(complex_roots)
160    0.000    0.000    0.000    0.000 {method 'imag' of 'sage.rings.complex_interval.ComplexIntervalFieldElement' objects}
36    0.000    0.000    0.000    0.000 complex_roots.py:112(column_disjoint)
85    0.000    0.000    0.000    0.000 {method 'contains_zero' of 'sage.rings.real_mpfi.RealIntervalFieldElement' objects}
3    0.000    0.000    0.000    0.000 {method '_coerce_map_via' of 'sage.structure.parent.Parent' objects}
97    0.000    0.000    0.000    0.000 weakref.py:336(__new__)
5    0.000    0.000    0.001    0.000 rings.py:845(__getitem__)
80    0.000    0.000    0.001    0.000 misc.py:5634(<lambda>)
152    0.000    0.000    0.000    0.000 {built-in method builtins.issubclass}
8    0.000    0.000    0.000    0.000 pushout.py:817(__init__)
97    0.000    0.000    0.000    0.000 {built-in method __new__ of type object at 0x7fc5f0f21a00}
6    0.000    0.000    0.000    0.000 {built-in method sage.structure.category_object.normalize_names}
333    0.000    0.000    0.000    0.000 {method 'append' of 'list' objects}
155    0.000    0.000    0.000    0.000 homset.py:1174(domain)
318    0.000    0.000    0.000    0.000 {method 'get' of 'dict' objects}
5    0.000    0.000    0.000    0.000 polynomial_ring.py:661(_coerce_map_from_)
306    0.000    0.000    0.000    0.000 complex_field.py:43(late_import)
1    0.000    0.000   41.768   41.768 {built-in method builtins.exec}
1    0.000    0.000    0.000    0.000 qqbar.py:6545(<listcomp>)
144    0.000    0.000    0.000    0.000 {method 'precision' of 'sage.rings.real_mpfr.RealNumber' objects}
5    0.000    0.000    0.000    0.000 {method 'coerce_map_from' of 'sage.structure.parent.Parent' objects}
2    0.000    0.000    0.000    0.000 polynomial_singular_interface.py:355(can_convert_to_singular)
4    0.000    0.000    0.000    0.000 {method '_populate_coercion_lists_' of 'sage.structure.parent.Parent' objects}
162    0.000    0.000    0.000    0.000 {method 'precision' of 'sage.rings.real_mpfr.RealField_class' objects}
3    0.000    0.000    0.000    0.000 complex_field.py:431(_coerce_map_from_)
1    0.000    0.000    0.000    0.000 sequence.py:80(Sequence)
8    0.000    0.000    0.000    0.000 polynomial_ring.py:628(construction)
2    0.000    0.000    0.000    0.000 polynomial_ring.py:1650(__init__)
3    0.000    0.000    0.000    0.000 rational_field.py:349(_coerce_map_from_)
16    0.000    0.000    0.000    0.000 <frozen importlib._bootstrap>:416(parent)
1    0.000    0.000   41.768   41.768 <string>:1(<module>)
1    0.000    0.000    0.000    0.000 qqbar.py:6533(<listcomp>)
146    0.000    0.000    0.000    0.000 {built-in method builtins.getattr}
200    0.000    0.000    0.000    0.000 {built-in method builtins.len}
2    0.000    0.000    0.000    0.000 polynomial_ring.py:1826(__init__)
4    0.000    0.000    0.000    0.000 polynomial_ring_constructor.py:672(_save_in_cache)
155    0.000    0.000    0.000    0.000 homset.py:1189(codomain)
149    0.000    0.000    0.000    0.000 copy.py:111(_copy_immutable)
2    0.000    0.000    0.000    0.000 category.py:525(_repr_object_names)
1    0.000    0.000   41.757   41.757 qqbar.py:6060(complex_roots)
1    0.000    0.000   41.768   41.768 qqbar.py:6171(__init__)
1    0.000    0.000   41.436   41.436 {method 'squarefree_decomposition' of 'sage.rings.polynomial.polynomial_element.Polynomial' objects}
1    0.000    0.000   41.768   41.768 qqbar.py:1361(polynomial_root)
1    0.000    0.000   41.757   41.757 qqbar.py:6637(_complex_isolate_interval)
```

### comment:3 Changed 3 years ago by vdelecroix

• Description modified (diff)
Last edited 3 years ago by vdelecroix (previous) (diff)

### comment:4 Changed 3 years ago by vdelecroix

• Authors set to Vincent Delecroix
• Branch set to u/vdelecroix/17895
• Commit set to 2984ecb978aca6e1f3d3926f3e4e9afbe954e3f7
• Milestone changed from sage-6.6 to sage-8.9
• Status changed from new to needs_review

With the change in the branch, timings go back to normal! Though, I have some strange doctest failures in `qqbar.py`. Before trying to understand them I would like to see if I did break anything else. So I am just waiting for a more complete patchbot report.

New commits:

 ​3893ba3 `rewrite root refinement in qqbar` ​2984ecb `allow rational polys in AlgebraicPolynomialTracker`

### comment:5 Changed 3 years ago by vdelecroix

• Dependencies set to #28191

### comment:6 Changed 3 years ago by git

• Commit changed from 2984ecb978aca6e1f3d3926f3e4e9afbe954e3f7 to 659b7a6ce931a998df29fca979772f604a78f677

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

 ​242a41a `allow rational polys in AlgebraicPolynomialTracker` ​659b7a6 `fix doctests in qqbar.py`

### comment:7 Changed 3 years ago by vdelecroix

• Description modified (diff)

### comment:9 Changed 3 years ago by vdelecroix

• Description modified (diff)

### comment:10 Changed 3 years ago by git

• Commit changed from 659b7a6ce931a998df29fca979772f604a78f677 to 49f917a3b74b3784e5ced901404b3778c94ecf92

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

 ​49f917a `also handle ZZ`

### comment:12 follow-up: ↓ 14 Changed 3 years ago by sbrandhorst

seems to look good. One question:

Is there a reason to hardcode `QQy` in several places?

### comment:13 Changed 3 years ago by sbrandhorst

• Reviewers set to Simon Brandhorst
• Status changed from needs_review to needs_info

### comment:14 in reply to: ↑ 12 Changed 3 years ago by vdelecroix

seems to look good. One question:

Is there a reason to hardcode `QQy` in several places?

Mainly for consistency and the fact that you don't want to multiply parents. This is an internal class, so the name of the variable does not matter from the user perspective. Before, it used to be canoncalized with a `y` (after the call to complex roots). I kept this convention.

### comment:15 Changed 3 years ago by sbrandhorst

You could add a doctest to show the speedup. Then I am happy.

### comment:16 Changed 3 years ago by git

• Commit changed from 49f917a3b74b3784e5ced901404b3778c94ecf92 to d5105b625060d4b2a9cfe9cee4316d8b7c96ea46

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

 ​d5105b6 `add a doctest for speedup`

### comment:17 Changed 3 years ago by vdelecroix

• Status changed from needs_info to needs_review

### comment:18 Changed 3 years ago by sbrandhorst

• Status changed from needs_review to positive_review

Thank you.

### comment:20 Changed 3 years ago by vbraun

• Branch changed from u/vdelecroix/17895 to d5105b625060d4b2a9cfe9cee4316d8b7c96ea46
• Resolution set to fixed
• Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.