Opened 6 years ago

Closed 6 years ago

#20064 closed defect (fixed)

Bug in sqrt in QQbar

Reported by: cremona Owned by:
Priority: major Milestone: sage-7.1
Component: numerical Keywords:
Cc: Merged in:
Authors: Nils Bruin Reviewers: John Cremona
Report Upstream: N/A Work issues:
Branch: 48c12ef (Commits, GitHub, GitLab) Commit: 48c12efbf46cbd4143b1c605004bd6224461b3d8
Dependencies: Stopgaps:

Status badges

Description

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)

Change History (7)

comment:1 Changed 6 years ago by cremona

(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.

comment:2 follow-up: Changed 6 years ago by cremona

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?

comment:3 in reply to: ↑ 2 Changed 6 years ago by nbruin

Replying to cremona:

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.

comment:4 Changed 6 years ago by nbruin

  • Branch set to u/nbruin/bug_in_sqrt_in_qqbar

comment:5 Changed 6 years ago by nbruin

  • Authors set to Nils Bruin
  • Commit set to 48c12efbf46cbd4143b1c605004bd6224461b3d8
  • Component changed from number fields to numerical
  • Status changed from new to needs_review

New commits:

48c12eftrac #20064: introduce "argument" on RealIntervalFieldElement

comment:6 Changed 6 years ago by cremona

  • Reviewers set to John Cremona
  • Status changed from needs_review to positive_review

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.

comment:7 Changed 6 years ago by vbraun

  • Branch changed from u/nbruin/bug_in_sqrt_in_qqbar to 48c12efbf46cbd4143b1c605004bd6224461b3d8
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.