Improve refine_root()

[jdemeyer]

In particular, allow both real and complex input. Also implement bisection if Newton-Raphson cannot be used.

Component: algebra

Author: Jeroen Demeyer

Branch/Commit: u/jdemeyer/ticket/19362 @ 9d5f99c

Issue created by migration from https://trac.sagemath.org/ticket/19362

[jdemeyer]

Branch: u/jdemeyer/ticket/19362

[jdemeyer]

Commit: 3b34e49

[jdemeyer]

New commits:

2049f5a	Move refine_root() to refine_root.pyx
3b34e49	Improve refine_root()

[sagetrac-git]

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

950eb84

Implement bisection

[sagetrac-git]

Changed commit from 3b34e49 to 950eb84

[jdemeyer]

comment:5

sage: from sage.rings.polynomial.refine_root import refine_root
sage: RIF = RealIntervalField(106)
sage: R.<x> = RIF[]
sage: pol = 10*x^6 - 10*x^2 + 1
sage: refine_root(pol, pol.derivative(), RIF(3/4,9/8), RIF)
Traceback (most recent call last):
...
ArithmeticError: cannot refine polynomial root (not enough steps?)

[sagetrac-git]

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

972f8e0

Once converging, keep converging

[sagetrac-git]

Changed commit from 950eb84 to 972f8e0

[sagetrac-git]

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

eb61557

When smashing, stop converging

[sagetrac-git]

Changed commit from 972f8e0 to eb61557

[sagetrac-git]

Changed commit from eb61557 to 7c689b6

[sagetrac-git]

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

7c689b6

When converging, take the intersection of old and new interval

[jdemeyer]

comment:10

I am going to continue working on this a bit more. If anybody feels like reviewing this code, let me know and I'll put up the current version for review.

[jdemeyer]

Changed dependencies from #19361 to #19361, #19364

[sagetrac-git]

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

388d495	Add edges() and endpoints() methods to intervals
8c10c00	Merge branch 't/19364/add_edges___method_to_rif_and_cif_elements' into t/19362/ticket/19362
8ab38c4	Trim instead of bisect

[sagetrac-git]

Changed commit from 7c689b6 to 8ab38c4

[videlec]

comment:13

Did you notice that some code in QQbar actually duplicates this: method _real_refine_interval/_complex_refine_interval of the ANRoot class.

[jdemeyer]

comment:14

Replying to @videlec:

Did you notice that some code in QQbar actually duplicates this: method _real_refine_interval/_complex_refine_interval of the ANRoot class.

Yes, there are way too many re-implementations of this in Sage. My eventual goal is to replace all "real/complex root refining" code in Sage by this new refine_root() function. But I'm not there yet...

[jdemeyer]

comment:15

Making it support all use cases of the various existing implementations is also what makes it tricky.

[sagetrac-git]

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

9d5f99c

Improve convergence check

[sagetrac-git]

Changed commit from 8ab38c4 to 9d5f99c

[JohnCremona]

comment:17

What is the work which needs doing?

[jdemeyer]

comment:18

I don't really remember myself, I do remember that it wasn't ready. It was trickier than I initially estimated. I will need to look at it again.

[jdemeyer]

Changed dependencies from #19361, #19364 to none

[mezzarobba]

comment:19

You probably know that, but just in case: it is not strictly true that you cannot use the interval Newton method and must switch to bisection when the slope interval contains zero. Another option is to interpret 1/[-a,b] as a kind of ”projective interval” containing ∞ (which gives rise to two disjoint complex intervals after intersecting with the previous estimate, so you'll still need to deal with several pieces).

[videlec]

comment:20

IMHO, refine_root is a bad name. I would expect such function to refine an interval that is already known to contain a (possibly unique) root... What about isolating_interval_from_approximate_root or something similar?

On the other hand, there are some possible alternative in the real case that does not involve convergence of Newton algorithm (and might already be implemented elsewhere). Namely p has a unique root in an interval [a, b] if the following holds:

1. the sign of sgn(p(a)) * sgn(p(b)) = -1
2. the interval p'([a,b]) does not contain zero
   It is also possible to replace 2 with Descartes rules of sign (which is a more expensive).