Reintroduce tight enclosure of interval inversion

[videlec]

In #37941 the __invert__ implementation of complex intervals was reverted to the straightforward implementation. This issue stands for reintroducing tight reversion (hopefully correct this time).

@miguelmarco @mezzarobba

[videlec]

I did use the following test function.

def check(re, im):
    precs = [6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 31, 32, 33, 100, 127, 128, 129]
    for prec in precs:
        c = ~ComplexIntervalField(prec)(re, im)
        norm2 = re**2 + im**2

        # check real part
        inv = re / norm2
        left, right = c.real().endpoints()
        left = QQ(left.as_integer_ratio())
        right = QQ(right.as_integer_ratio())
        if left > inv or right < inv:
             print('wrong real at prec={}: left={}~{} right={}~{} value={}~{}'.format(prec, left, left.n(prec=2*prec), right, right.n(prec=2*prec), inv, inv.n(prec=2*prec)))

        # check imaginary part
        inv = - im / norm2
        left, right = c.imag().endpoints()
        left = QQ(left.as_integer_ratio())
        right = QQ(right.as_integer_ratio())
        if left > inv or right < inv:
            print('wrong imag at prec={}: left={}~{} right={}~{} value={}~{}'.format(prec, left, left.n(prec=2*prec), right, right.n(prec=2*prec), inv, inv.n(prec=2*prec)))

It seems to be the case that the tight inversion has

• correct real part always
• wrong imaginary part often

sage: check(1/2, 1/3)
wrong imag at prec=11: left=-473/512~-0.923828 right=-1891/2048~-0.923340 value=-12/13~-0.923077
wrong imag at prec=12: left=-1891/2048~-0.923340 right=-3781/4096~-0.923096 value=-12/13~-0.923077
wrong imag at prec=17: left=-120991/131072~-0.923088074 right=-60495/65536~-0.923080444 value=-12/13~-0.923076923
sage: check(1, 5)
wrong imag at prec=6: left=-25/128~-0.195 right=-25/128~-0.195 value=-5/26~-0.192
wrong imag at prec=7: left=-99/512~-0.193 right=-99/512~-0.193 value=-5/26~-0.192
wrong imag at prec=8: left=-197/1024~-0.1924 right=-197/1024~-0.1924 value=-5/26~-0.1923
wrong imag at prec=9: left=-197/1024~-0.19238 right=-197/1024~-0.19238 value=-5/26~-0.19231
wrong imag at prec=10: left=-197/1024~-0.19238 right=-197/1024~-0.19238 value=-5/26~-0.19231
wrong imag at prec=11: left=-197/1024~-0.192383 right=-197/1024~-0.192383 value=-5/26~-0.192308
wrong imag at prec=12: left=-3151/16384~-0.192322 right=-3151/16384~-0.192322 value=-5/26~-0.192308
wrong imag at prec=13: left=-3151/16384~-0.1923218 right=-3151/16384~-0.1923218 value=-5/26~-0.1923077
wrong imag at prec=14: left=-3151/16384~-0.19232178 right=-3151/16384~-0.19232178 value=-5/26~-0.19230769
...
sage: check(1/4, -1/2)
wrong imag at prec=6: left=13/8~1.62 right=13/8~1.62 value=8/5~1.60
wrong imag at prec=7: left=103/64~1.61 right=103/64~1.61 value=8/5~1.60
wrong imag at prec=8: left=205/128~1.602 right=205/128~1.602 value=8/5~1.600
wrong imag at prec=9: left=205/128~1.6016 right=205/128~1.6016 value=8/5~1.6000
wrong imag at prec=10: left=205/128~1.6016 right=205/128~1.6016 value=8/5~1.6000
...

[unhyperbolic]

I started looking at the implementation of the Rokne-Lancaster method that was just removed in 8d59b12 and found a couple of problems:

1. The following lines should be a rounding down instead of rounding up:
   

   sage/src/sage/rings/complex_interval.pyx

   Line 1251 in d2b1e01

   mpfr_div(imax, c, div1, MPFR_RNDU)

   
   

   sage/src/sage/rings/complex_interval.pyx

   Line 1254 in d2b1e01

   mpfr_div(aux2, d, div2, MPFR_RNDU)

   
   Also note that subtracting from 0 can be avoided by simply switching the sign with mpfr_neg (which can even be done in place).
2. The computation of rmin (the lower bound of the real part of the inverse) does not depend on d (the upper bound of the imaginary part), but it should, consider for example [10,11]+[-3,2], [10,11]+[-3,3], [10,11]+[-3,4]:

   sage/src/sage/rings/complex_interval.pyx

   Lines 1298 to 1305 in d2b1e01

   	if mpfr_cmp(aux, c) >= 0:
   	mpfr_set_str(aux, '-0.5', 10, MPFR_RNDD)
   	mpfr_div(rmin, aux, c, MPFR_RNDD)
   	else:
   	mpfr_mul(a2, a, a, MPFR_RNDD)
   	mpfr_mul(c2, c, c, MPFR_RNDD)
   	mpfr_add(div1, a2, c2, MPFR_RNDD)
   	mpfr_div(rmin, a, div1, MPFR_RNDU)

   
   Also note that the parsing with mpfr_set_str can be avoided by using mpfr_mul_2si or mpfr_div_2si.

I think I understand the Rokne-Lancaster paper well enough that I can give an implementation a shot.

[unhyperbolic]

Feel free to assign this to me.

[videlec]

@unhyperbolic Would be nice if you can provide a (correct) tight enclosure! There is no need for a formal assignment of task. We do not use that much in sagemath.