Comparisons in quadratic number field

[videlec]

The order of quadratic field with specified embedding is not induced from the order of RR and CC. More precisely we have

sage: K.<sqrt2> = NumberField(x^2 - 2, 'sqrt2', embedding=1.4142)
sage: sqrt2 < 100
False

sage: K.<s> = NumberField(x^3 - 2, 's', embedding=1.26)
sage: s < 100
False

which is not compatible with the order of RR and

sage: K.<i> = QuadraticField(-1)
sage: i > 1
True
sage: 1 > i
True

which is not compatible (!) with the order of CC.

The ticket implements the order for quadratic field (for which comparisons are made using only operations on integers).

Note:

• this patch is partly a duplicate because of comparison with quadratic number field elements #7160
• the modifications for quadratic field field modify the behavior of many commands (especially about output order).
• the patch introduces a new attribute to NumberFieldElement_quadratic called standard_embedding (of type bint)

---

There is no significant lost of speed as shown on the following timings.

sage: K.<sqrt2> = QuadraticField(2, 'sqrt2', embedding=1.4142)
sage: a = (3*sqrt2 + 18)/7
sage: b = (5*sqrt2 + 14)/5
sage: %timeit a < b
1000000 loops, best of 3: 408 ns per loop

sage: K.<s> = QuadraticField(-2)
sage: a = 3*s + 2/4
sage: b = 5/7*s + 1/3
sage: %timeit a < b
1000000 loops, best of 3: 352 ns per loop

Timings without the patch

sage: sage: K.<sqrt2> = QuadraticField(2,'sqrt2',embedding=1.4142)
sage: a = (3*sqrt2 + 18)/7
sage: b = (5*sqrt2 + 14)/5
sage: %timeit a < b
1000000 loops, best of 3: 393 ns per loop

sage: K.<s> = QuadraticField(-2)
sage: a = 3*s+2/4
sage: b = 5/7*s+1/3
sage: %timeit a < b
1000000 loops, best of 3: 344 ns per loop

Apply:

• attachment: trac_13213-quadratic_field_comparison.patch

Component: number fields

Keywords: real number field, comparison

Author: Vincent Delecroix

Reviewer: Burcin Erocal, Volker Braun

Merged: sage-5.10.beta4

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

[a-andre]

comment:1

I had a short look at your patch.
Could you use Python 3 compatible syntax to raise errors. E.g.

raise ValueError("absolute value not implemented for complex embeddings")

[a-andre]

Changed author from vdelecroix to Vincent Delecroix

[videlec]

Description changed:

--- 
+++ 
@@ -9,3 +9,23 @@
 ```
 
 There is a patch for quadratic field. Another one is comming for general number fields.
+
+Timings with the patch
+
+```
+sage: K.<sqrt2> = QuadraticField(2,'sqrt2',embedding=1.4142)
+sage: b = (5*sqrt2 + 14)/5
+sage: a = (3*sqrt2 + 18)/7
+sage: %timeit a < b
+625 loops, best of 3: 6.68 µs per loop
+```
+Timings without the patch
+
+```
+sage: sage: K.<sqrt2> = QuadraticField(2,'sqrt2',embedding=1.4142)
+sage: a = (3*sqrt2 + 18)/7
+sage: b = (5*sqrt2 + 14)/5
+sage: %timeit a < b
+625 loops, best of 3: 478 ns per loop
+```
+}}}

[videlec]

Description changed:

--- 
+++ 
@@ -8,7 +8,8 @@
 True
 ```
 
-There is a patch for quadratic field. Another one is comming for general number fields.
+There is a patch for quadratic field. Another one is comming for general number fields. Note that this patch is partly a duplicate because of #7160.
+
 
 Timings with the patch
 
@@ -28,4 +29,3 @@
 sage: %timeit a < b
 625 loops, best of 3: 478 ns per loop
 ```
-}}}

[videlec]

Description changed:

--- 
+++ 
@@ -6,6 +6,8 @@
 True
 sage: sqrt2 < 2
 True
+sage: sqrt2.floor()
+1
 ```
 
 There is a patch for quadratic field. Another one is comming for general number fields. Note that this patch is partly a duplicate because of #7160.

[burcin]

comment:5

Thanks again for the patch. A few comments:

• reimplementing the comparison of quadratic number field elements is a well defined unit of change itself. I suggest moving floor() and ceil() implementations to a separate patch, even a different ticket.
• Do you really need to use rich comparisons? Working only with __cmp__ should simplify the case comparisons before the return statements and avoid the new python_object.pxi include.
• I guess speed could be further improved by using an approximation (with doubles say) first and falling back to the exact computation if the results are too close.
• left.D is of type Integer. Cython tries to do the right thing for <unsigned int> left.D, but it would be better to use mpz_mul() with left.D.value.
• mpz_sgn() returns 0 if the argument is 0. You don't seem to cover this case, though I couldn't find any example that returns wrong results.
• To avoid code duplication, it might be better change the if statement to:

        if mpz_sgn(i)*mpz_sgn(j) == 1:
            # both i and j are positive or negative
            mpz_mul(i, i, i)
            mpz_mul(j, j, j)
            mpz_mul_ui(j, j, <unsigned int> left.D)

            test = mpz_cmp(i, j)
            if mpz_sgn(i) == -1:
                # both i and j are negative
                test *= -1

[burcin]

Reviewer: Burcin Erocal

[videlec]

comment:6

• reimplementing the comparison of quadratic number field elements is a well defined unit of change itself. I suggest moving floor() and ceil() implementations to a separate patch, even a different ticket.
• left.D is of type Integer. Cython tries to do the right thing for <unsigned int> left.D, but it would be better to use mpz_mul() with left.D.value.

Done and done.

• Do you really need to use rich comparisons? Working only with __cmp__ should simplify the case comparisons before the return statements and avoid the new python_object.pxi include.

We loose speed (~1 micro sec) by doing this. I quickly look at the code of comparison for Integer. The comparison is directly implemented inside __richcmp__ and avoid as much as possible the coercion machinery (ie the type of the input is check with the cython function PY_TYPE_CHECK). Do you think we should proceed that way ?

• I guess speed could be further improved by using an approximation (with doubles say) first and falling back to the exact computation if the results are too close.

What do you mean by too close ? How do I certify that my float comparison is accurate enough ?
Moreover, if we proceed that way we have to perform conversion from mpz_t to float (or double or real_mpfr). A way it may work is to look first if a,b, D and denom have reasonable values (recall that mpz_t have unlimited precision) and if yes, convert it to double and try to compare the arguments.
Last, I'm not sure the gain of speed would be terrific. Comparison requires only 7 integer operations.

• mpz_sgn() returns 0 if the argument is 0. You don't seem to cover this case, though I couldn't find any example that returns wrong results.

Actually there is no problem since we compute the sign of an expression of the form a^2 + b^2 D where a and b are both integers (and one of them is non null) and D is a square free integer. Hence the sign is never 0.

• To avoid code duplication, it might be better change the if statement to:

        if mpz_sgn(i)*mpz_sgn(j) == 1:
            # both i and j are positive or negative
            mpz_mul(i, i, i)
            mpz_mul(j, j, j)
            mpz_mul_ui(j, j, <unsigned int> left.D)

            test = mpz_cmp(i, j)
            if mpz_sgn(i) == -1:
                # both i and j are negative
                test *= -1

Sure but that way we call twice mpz_sgn(i)... I may change it if you prefer.

Actually, the part where a lot of time is lost is the check for embedding

if not RDF.has_coerce_map_from(left._parent):
    ...

By replacing it by

if left.D < 0:
   ...

I pass from 4 micro second to 2 micro seconds. That way comparison without embedding for positive discriminant will be equivalent to comparison with the standard embedding. This is in the last version of the patch. What do you think about that ?

[videlec]

Description changed:

--- 
+++ 
@@ -1,26 +1,27 @@
-Provide the order from RR for number field embedded in RR as in
+Provide the order from RR for real quadratic number field. The current implementation does
 
 ```
 sage: K.<sqrt2> = NumberField(x^2 - 2, 'sqrt2', embedding=1.4142)
-sage: sqrt2 - 1 > 0
-True
-sage: sqrt2 < 2
-True
-sage: sqrt2.floor()
-1
+sage: sqrt2 < 100
+False
 ```
 
 There is a patch for quadratic field. Another one is comming for general number fields. Note that this patch is partly a duplicate because of #7160.
 
-
-Timings with the patch
+The lost of speed is about x5 for positive discriminant and almost nothing for negative ones
 
 ```
 sage: K.<sqrt2> = QuadraticField(2,'sqrt2',embedding=1.4142)
+sage: a = (3*sqrt2 + 18)/7
 sage: b = (5*sqrt2 + 14)/5
-sage: a = (3*sqrt2 + 18)/7
 sage: %timeit a < b
-625 loops, best of 3: 6.68 µs per loop
+625 loops, best of 3: 2.36 µs per loop
+
+sage: K.<s> = QuadraticField(-2)
+sage: a=3*s+2/4
+sage: b=5/7*s+1/3
+sage: %timeit a < b
+625 loops, best of 3: 600 ns per loop
 ```
 Timings without the patch
 
@@ -29,5 +30,11 @@
 sage: a = (3*sqrt2 + 18)/7
 sage: b = (5*sqrt2 + 14)/5
 sage: %timeit a < b
-625 loops, best of 3: 478 ns per loop
+625 loops, best of 3: 491 ns per loop
+
+sage: K.<s> = QuadraticField(-2)
+sage: a=3*s+2/4
+sage: b=5/7*s+1/3
+sage: %timeit a < b
+625 loops, best of 3: 488 ns per loop
 ```