Enhanced checks for sign() with qqbar elements with addition and subtraction

[tscrim]

The following are painfully slow:

sage: x1 = AA(2^(1/100))
sage: x2 = AA(2^(1/100))
sage: y = x1 - x2
sage: z = x1 - x2
sage: y == y
sage: y == 0
sage: y == z

Note that for testing, you need to recreate y before each comparison as otherwise it has been made exact.

We avoid computing the minimum polynomial (as it calls exactify()) in the rich comparison until we need to. We also avoid calling exactify() for a binary operation in sign() until we absolutely need to. Additionally, the identity test is not done before reaching _richcmp_.

[tscrim]

I have not run any timings on this, but things to check for speed regressions would be low degree minimal polynomials. Suggestions for good timing tests would be useful.

I don't completely understand why I need to update some of the test for the eigenvalues, but the real precision seems to have significantly improved. So there might be some other doctests that need fixing; waiting for the patchbot.

I also took the change to mark a test as # long time that was taking nearly 2 seconds on my computer.

---

New commits:

​fc64ba2

Improvements to sign() and comparisons for qqbar elements.

[tscrim]

See also #31768 for a followup subclassing ANUnaryExpr for some finer control and trivial simplifications.

[vdelecroix]

This proposal looks good to me but

+            # Rationals
+            if type(sd._left._descr) is ANRational and type(sd._right._descr) is ANRational:
+                ret = sd._op(sd._left._descr._value, sd._right._descr._value)
+                if ret == 0:
+                    self._set_descr(ANRational(QQ.zero()))
+                    return 0
+                return ret.sign()

Do you have an instance where the above code path is used? When both operands are rationals the descriptor should be automatically updated. We don't want to repeat this block in all possible places cmp, sign, _add_, ...

Replying to tscrim:

I have not run any timings on this, but things to check for speed regressions would be low degree minimal polynomials. Suggestions for good timing tests would be useful.

I don't completely understand why I need to update some of the test for the eigenvalues, but the real precision seems to have significantly improved. So there might be some other doctests that need fixing; waiting for the patchbot.

When precision increases it is typically that more enclosure refinement occurred. With the new code it looks to me that we should have less refinement. If my diagnostic is correct, this is bad. Would be nice to check what are the diameters of the intervals before and after the change.

[tscrim]

Replying to vdelecroix:

This proposal looks good to me but

+            # Rationals
+            if type(sd._left._descr) is ANRational and type(sd._right._descr) is ANRational:
+                ret = sd._op(sd._left._descr._value, sd._right._descr._value)
+                if ret == 0:
+                    self._set_descr(ANRational(QQ.zero()))
+                    return 0
+                return ret.sign()

Do you have an instance where the above code path is used? When both operands are rationals the descriptor should be automatically updated. We don't want to repeat this block in all possible places cmp, sign, _add_, ...

This is occurring after a call to exactify(). Thus, we could have something that ends up being a rational that didn't know it before. For example, b = a - (a + 1) for some irrational a, then taking b + 1 == 0). Adding a print("rat check") just after the if statement gives:

sage: a = AA(sqrt(2))
sage: b = a - (a + 1)
sage: (b + 1).sign()
rat check
0

I will add this as a test and push. I also need to figure out why the test in src/sage/geometry/polyhedron/backend_field.py is failing too...

Replying to tscrim:

I have not run any timings on this, but things to check for speed regressions would be low degree minimal polynomials. Suggestions for good timing tests would be useful.

I don't completely understand why I need to update some of the test for the eigenvalues, but the real precision seems to have significantly improved. So there might be some other doctests that need fixing; waiting for the patchbot.

When precision increases it is typically that more enclosure refinement occurred. With the new code it looks to me that we should have less refinement. If my diagnostic is correct, this is bad. Would be nice to check what are the diameters of the intervals before and after the change.

Here is what is happening. Because we are sorting the (complex) roots, we compare things lex on the (real, imag) pair. Eventually it has to check that the real parts of the two numbers are not equal (see case 3:). This calls the real _richcmp_ method. In here, before it computed the minpoly() (which called exactify() of the treal parts), whereas now only gets the minpoly if they are known, so it first tries the refinement before ultimately ended up calling exactify() (as both are supposed to actually be 0). This gets called when comparing these roots with a root that is explicitly known with a real part of 0.

I don't think this is a bad thing as I believe the refinement is a relatively fast operation, especially if we ultimately compute exactify() on most elements. Also, getting rid of calls to exactify() in cases when it is not needed is a larger benefit IMO. However, some more actual use-cases could be useful here for testing timings.

[git]

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

​b069b81

Fixing bug with addition x + (-x), adding doctests, removing redundancy.

[git]

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

​3b8e217

Removing one other redundancy.

[git]

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

​5da8b27

Removing same redundant checks for the QQbar elements.

[tscrim]

I have fixed a bug with addition x + (-x) giving the wrong sign. This was the issue with backend_field.py. I added in a few other doctests and removed some redundant checks.

I just don't know how much this will affect timings refining the precision instead of calling the exactify() via the minpoly(). Something else that might be useful is adding a way to compute the minpoly without computing exactify() for elements that have an idea about what their minpoly is. Mainly here I am thinking elements using ANRoot. We could also fashion something that passes around the minpoly when we do operations that preserve it, such as conjugation. This can be done on a separate ticket. What do you think?

[vdelecroix]

Replying to tscrim:

I have fixed a bug with addition x + (-x) giving the wrong sign. This was the issue with backend_field.py. I added in a few other doctests and removed some redundant checks.

I just don't know how much this will affect timings refining the precision instead of calling the exactify() via the minpoly(). Something else that might be useful is adding a way to compute the minpoly without computing exactify() for elements that have an idea about what their minpoly is. Mainly here I am thinking elements using ANRoot. We could also fashion something that passes around the minpoly when we do operations that preserve it, such as conjugation. This can be done on a separate ticket. What do you think?

Computing minpoly is doable and desirable for all descriptors without exactification. This is the reason why was introduced the function composed_op on polynomials

sage: x = polygen(QQ)
sage: p = x^3 + 2*x - 1
sage: q = x^3 - x + 3
sage: a = p.roots(QQbar, multiplicities=False)[0]
sage: b = q.roots(QQbar, multiplicities=False)[0]
sage: r = p.composed_op(q, operator.add)
sage: r(a + b)  # not necessary minpoly, but at least a multiple
0.?e-15

More generally, I think that QQbar elements beyond degree 8 (let say) should not try to find a reasonable representation of the number field they represent (what is exactify currently doing via polred which is costly). They should just be considered as the root of some irreducible polynomial. Computing minpoly of a + b, abs(a), etc is doable.

Though, all this logic is now in ​Fredrik's calcium and I am not sure how much it is desirable to improve the version in Sage compared to having a calcium version of AA/QQbar.

[gh-Ilia-Smilga]

Replying to tscrim:

New commits:

​fc64ba2

Improvements to sign() and comparisons for qqbar elements.

Why not use (good enough refinements of) approximate values to compute the sign of x-y when x and y are known to have different minimal polynomials? (cf. ​https://groups.google.com/g/sage-devel/c/OaUlf5VHBbc/m/doZC1KpMBAAJ )

[tscrim]

comment:9: I agree that it is better to not invest a lot of time in improving AA/QQbar since it will almost certainly be replaced by Calcium. However, I think we can do a few improvements around the current code, and some of these changes might still be applicable once we switch the backend over (unlike, say, #31768). Perhaps we can include the current patch, maybe I will do a followup to try and lessen the impact of the changes.

Fredrik, do you have any comments about this?

comment:10: Once you have different minimal polynomials, then I believe you essentially have everything. Perhaps I have misunderstood something?

Green patchbot (the pyflakes warning existed before this ticket).

[vdelecroix]

Actually, to construct a better refinement of a + b knowing the minimal polynomials m_a(x) and m_b(x) there are at least two strategies

• compute m_{a+b}(x) (eg via composed_op) and do refinement with it
• computer better refinement of a and b and add them

Depending on deg(m_a), deg(m_b), deg(m_{a+b}) it might be better to use the second option.

[fredrik.johansson]

To compare (the real parts of) two qqbars in Calcium, I just compare enclosures at a couple of different levels of precision, otherwise do an exact subtraction and check the sign of the difference. Most of the work will then be in the subtraction (which factors the resultant polynomial and then refines the precision of the inputs if needed to get the right root). This is not really optimized, but so far it has not been a bottleneck that I needed to fix. In fact, this issue reminds me that I did not even optimize for both inputs being real, in which case it will suffice to increase the precision until convergence...

More generally, it's an interesting problem to optimize computing sign(a1+a2+...+an) or proving a1+a2+...+an == 0 without computing the entire sum explicitly.

[gh-Ilia-Smilga]

Replying to tscrim:

comment:10: Once you have different minimal polynomials, then I believe you essentially have everything. Perhaps I have misunderstood something?

Consider the following computation:

x = AA(1+100*2^(-1000))^(1/100)
y = AA(1+101*2^(-1000))^(1/101)
x < y

(which calls sign(y-x) internally). With your patch (if I read it correctly; I did not actually run it), Sage will notice that the minimal polynomials are different, but will not take advantage of that information. It will still trigger exact computation of y-x (which takes forever).

I suggest adding for this case (different minimal polynomials) the code

while self._value.contains_zero():
    self._more_precision()
return self._value.unique_sign()

Is there any reason not to do so? (For the record: in this particular example, this loop takes 82ms on my machine.)

[tscrim]

Replying to gh-Ilia-Smilga:

I suggest adding for this case (different minimal polynomials) the code

while self._value.contains_zero():
    self._more_precision()
return self._value.unique_sign()

Is there any reason not to do so? (For the record: in this particular example, this loop takes 82ms on my machine.)

I see what you're saying, but that is a more invasive change. It probably is better to just do that here. However, I don't quite understand Dima's suggestion on the sage-devel thread. It probably is something simple, and I am just too tired to see it right now.

[vdelecroix]

In this part of the new code

        if sd._left.minpoly() == sd._right.minpoly():
            # Negating the element does not change the minpoly
            right = sd._right if sd._op is operator.sub else -sd._right
            c = cmp_elements_with_same_minpoly(sd._left, right, sd._left.minpoly())
            if c == 0:
                self._set_descr(ANRational(QQ.zero()))
                return 0
            elif c is not None:
                return c

you seemed to ignore that the operator could be mul or truediv.

[soehms]

I wonder if we can improve exactification instead of avoiding it. At the moment any two instances of AlgebraicGenerator are recognized to be different, even if they describe the same root of the same polynomial. Is there a special reason for this? What is the idea behind the attribute _index?

Howsoever, I did the following attempt modifying the method _richcmp_ of this class and changed self is other in union to self == other and accordingly super is self to super == self in super_poly (diff shown with respect to stable 9.3):

src/sage/rings/qqbar.py

@@ -2911 +2911 @@
             sage: gen > qq_generator
             True
         """
+        sf = self._field
+        of = other._field
+        if type(sf) != type(of):
+           return rich_to_bool(op, 1)
+        if sf != of:
+           return richcmp_not_equal(sf, of, op)
+        sa = self.root_as_algebraic()
+        oa = other.root_as_algebraic()
+        if sa in QQ and oa in QQ:
+            return richcmp(QQ(sa), QQ(oa), op)
+        c = cmp_elements_with_same_minpoly(sa, oa, sf.polynomial())
+        if c is not None:
+            return rich_to_bool(op, c)
         return richcmp(self._index, other._index, op)
 
     def is_complex(self):
@@ -3108 +3121 @@
             return other
         if other._trivial:
             return self
-        if self is other:
+        if self == other:
             return self
         if other in self._unions:
             return self._unions[other].parent
@@ -3207 +3220 @@
         if checked is None:
             checked = {}
         checked[self] = True
-        if super is self:
+        if super == self:
             return self._field.gen()
         for u in self._unions.values():
             if u.parent in checked:

This change passes through all doctests (including --long) and makes exactification possible for the testcase of the ticket description:

sage: x1 = AA(2^(1/100))
....: x2 = AA(2^(1/100))
....: y = x1 - x2
....: z = x1 - x2
sage: %time y._exact_field()
CPU times: user 160 ms, sys: 16 ms, total: 176 ms
Wall time: 174 ms
Trivial generator
sage: %time y == z
CPU times: user 260 ms, sys: 0 ns, total: 260 ms
Wall time: 259 ms
True

Of course this does not solve the testcase of ​comment:14. Concerning this it seems, that the function composed_op could do a good job here:

sage: x = AA(1+2*2^(-1000))^(1/2); y = AA(1+3*2^(-1000))^(1/3)
sage: p = x.minpoly()
sage: q = y.minpoly()
sage: %time r = p.composed_op(q, operator.sub)
CPU times: user 4 ms, sys: 0 ns, total: 4 ms
Wall time: 5.53 ms
sage: %time F = r.factor()
CPU times: user 0 ns, sys: 0 ns, total: 0 ns
Wall time: 2.65 ms

I would prefer such a solution since it would not just fix this issue, but also improve exactification in general. On the other hand I agree with Vincent that this only makes sense (in view of the upcoming Calcium), if it will be easy to realize.

But note, that this issue is different from the former one in not just hanging in Pari's nffactor but already one step before in the construction of the Pari number field (of degree two!):

sage: AG = x._exact_field()
sage: AG.pari_field()
^C---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-8-2b02365498e1> in <module>
----> 1 AG.pari_field()

~/devel/sage/local/lib/python3.9/site-packages/sage/rings/qqbar.py in pari_field(self)
   3025         if self._pari_field is None:
   3026             pari_pol = self._field.pari_polynomial("y")
-> 3027             self._pari_field = pari_pol.nfinit(1)
   3028         return self._pari_field
   3029

cypari2/auto_gen.pxi in cypari2.gen.Gen_base.nfinit()

KeyboardInterrupt:

Thus, I would take this for bug in Pari, isn't it?

[tscrim]

Replying to vdelecroix:

In this part of the new code

        if sd._left.minpoly() == sd._right.minpoly():
            # Negating the element does not change the minpoly
            right = sd._right if sd._op is operator.sub else -sd._right
            c = cmp_elements_with_same_minpoly(sd._left, right, sd._left.minpoly())
            if c == 0:
                self._set_descr(ANRational(QQ.zero()))
                return 0
            elif c is not None:
                return c

you seemed to ignore that the operator could be mul or truediv.

Those cases already would have been picked up earlier in this block:

         elif type(sd) is ANBinaryExpr:
             ls = sd._left.sign()
             rs = sd._right.sign()
             if sd._op is operator.mul or sd._op is operator.truediv:
-                return sd._left.sign() * sd._right.sign()
+                return ls * rs

[tscrim]

@soehms That seems like a good idea. However, we now have to come up with a better hash. My guess behind the _index idea is that they wanted very fast comparison checks and hashing that is consistent with equality. This line actually dates back to the initial implementation of algebraic reals (commit 07704df513). So there is room for improvement here.

[vdelecroix]

What about postponing further improvements to other tickets? This one is ready and I think it is worth setting to positive review.

[vdelecroix]

Replying to gh-Ilia-Smilga:

Replying to tscrim:

comment:10: Once you have different minimal polynomials, then I believe you essentially have everything. Perhaps I have misunderstood something?

Consider the following computation:

x = AA(1+100*2^(-1000))^(1/100)
y = AA(1+101*2^(-1000))^(1/101)
x < y

I agree that for a single + operation on two elements for which you know the minimal polynomial, infinite refinement is a way to go.

However, you need someting to compute minimal polynomials in order to handle

sage: a = AA(1 + 2^(-1000))
sage: x = a^(1/100) + a^(1/104) - a^(1/102) - a^(1/103)

I opened #31889 for that purpose.

[soehms]

Replying to vdelecroix:

What about postponing further improvements to other tickets? This one is ready and I think it is worth setting to positive review.

Agreed!

[git]

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

​526bb7d	Merge branch 'public/performance/qqbar_sign-31767' of git://trac.sagemath.org/sage into public/performance/qqbar_sign-31767
​642b7c2	Fixing pyflakes warning.

[tscrim]

The bot was green modulo one pyflakes before I pushed. The last change just removes an unused import. Tests still pass for me.

[tscrim]

Thank you.