Bug in sqrt in QQbar

[JohnCremona]

See #18836. This bug is holding up that (and also #20028). The following code creates an elemnt d of QQbar and tries to do d.sqrt(). It fails unless you call d.imag().is_zero() first.

K.<i> = QuadraticField(-1)

# define a low-precision embedding from K to CC:                                                                

emb = K.embeddings(CC)[1]

# extend this to the closest embedding into QQbar:                                                              

old_gen = emb(K.gen())
rr = K.defining_polynomial().roots(QQbar, multiplicities=False)
diffs = [(CC(r)-old_gen).abs() for r in rr]
new_gen = rr[diffs.index(min(diffs))]
emb0 = K.hom([new_gen], check=False)

# Take a polynomial with 3 roots in K:                                                                          

e1 = -4+i
e2 = 1+i
e3 = 3-2*i
print("Original ei: %s with parent %s" % ([e1,e2,e3],parent(e1)))
x = polygen(K)
pol = (x-e1)*(x-e2)*(x-e3)

# Find the roots again in QQbar:                                                                                

pol0 = PolynomialRing(QQbar,'x')([emb0(c) for c in list(pol)])
e1, e2, e3 = pol0.roots(QQbar,multiplicities=False)
print("Roots ei: %s with parent %s" % ([e1,e2,e3],parent(e1)))

# Attempt to compute sqrt(e1-e2) from these:                                                                    

d = e1-e2
print("d=%s with parent %s" % (d,d.parent()))
# If the next 2 lines are commented out, an error is raised in the sqrt!                                        
s = d.imag().is_zero()
print("d.imag().is_zero()=%s" % s)
print("d=%s with parent %s" % (d,d.parent()))
d = d.sqrt()
print("d=%s" % d)

Component: numerical

Author: Nils Bruin

Branch/Commit: 48c12ef

Reviewer: John Cremona

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

[JohnCremona]

comment:1

(Extracted from #18836 comment 9)

The fact that _value becomes a real interval element is probably due to this implementation (line 7431 of sage/rings/qqbar.py)

    def _interval_fast(self, prec):
        gen_val = self._generator._interval_fast(prec)
        v = self._value.polynomial()(gen_val)
        if self._exactly_real and is_ComplexIntervalFieldElement(v):
            return v.real()
        return v

It could well be that this is really the intention (and there could be other places where the interval is set to be a real thing!), in which case the bug is indeed in sqrt, which should avoid relying on "argument" being available.

Looking at it a little bit more, I think the error is in fact in __pow__ (sage/rings/qqbar.py:4233). I expect it's not invalid for _value to be a RIF element. That seems to happen quite a bit:

sage: type(QQbar(5)._value)
<type 'sage.rings.real_mpfi.RealIntervalFieldElement'>

It's just that elements where that happens are usually filtered out earlier. So in this case, d isn't recognized as a rational number yet when we enter, but then in the for loop on line 4304, we have that on line 4322 we execute:

               val = self._interval_fast(prec)

So actually, on the first pass val is a CIF element. We apparently get an error on the second pass, when val has been forced through _interval_fast.

So I see two solutions: either ensure that val is put back into CIF or ensure that "argument extraction" is done appropriately in cases where "val" is in RIF.

[JohnCremona]

comment:2

Surely it is simplest to implement arg for RIF elements, returning 0 or Pi depending on sign provided that the interval has constant sign, and raising an error otherwise?

[nbruin]

comment:3

Replying to @JohnCremona:

Surely it is simplest to implement arg for RIF elements, returning 0 or Pi depending on sign provided that the interval has constant sign, and raising an error otherwise?

I would regard that as a "change in design/API", though. I think our real fields currently quite consistently do not have an "argument" method:

sage: CC(-4).argument()
3.14159265358979
sage: RR(-4).argument()
AttributeError: ...
sage: CDF(-4).argument()
3.141592653589793
sage: RDF(-4).argument()
AttributeError: ...

Your proposal is in line with the solution taken in #18337 (to put "real" and "imag" on RIF) so perhaps putting "argument" there is the simplest solution.

From an efficiency point of view: I'm not sure how hard we'd be hitting the coercion system by mixing RIF and CIF elements.

[nbruin]

Branch: u/nbruin/bug_in_sqrt_in_qqbar

[nbruin]

New commits:

48c12ef

trac #20064: introduce "argument" on RealIntervalFieldElement

[nbruin]

Author: Nils Bruin

[nbruin]

Commit: 48c12ef

[JohnCremona]

comment:6

This looks good to me. I have checked that it deals with the original problems I had at #18836, so I am going to set this to positive_review, make that ticket depend on this and set that one to needs_review.

[JohnCremona]

Reviewer: John Cremona

[vbraun]

Changed branch from u/nbruin/bug_in_sqrt_in_qqbar to 48c12ef