Universal Cyclotomic Field implementation using libgap

[vdelecroix]

Sage ships an implementation of the Universal Cyclotomic Field (in sage.rings.universal_cyclotomic_field.*) that is slower and less reliable than what is in gap. We propose in this ticket an implementation based on libgap.

It fixes some issues about the old implementation:

• #14240: Universal cyclotomic field breaks for moderate order
• #16631: Universal Cyclotomic Field doesn't handle int and Integer the same
• #16130: Matrices over Universal Cyclotomic field failing to echelonize/inverse
• #17117: UniversalCyclotomicField elements considered real

And is about 10x faster (see comment:10) except for operations on very sparse elements with very large conductors (see e.g. comment:23, comment:27 and comment:28).

Possible follow up:

• switch to Cython
• implement dense matrices by wrapping libgap matrices (see also #16116)
• slowness in conjugation? #18207 (see also comment:28, this is unavoidable because gap uses a dense representation of elements)
• interface with libgap polynomials (#18266)

[git]

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

​31385f3	Trac 18153: int and infinity from libgap
​23eb30c	Trac 18152: implementation of UCF using libgap
​99c8d80	Trac 18152: remove the old implementation
​7d95734	Trac 18152: Fix coxeter group import

[vdelecroix]

New commits:

​31385f3	Trac 18153: int and infinity from libgap
​23eb30c	Trac 18152: implementation of UCF using libgap
​99c8d80	Trac 18152: remove the old implementation
​7d95734	Trac 18152: Fix coxeter group import

[git]

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

​36fb6e5	Trac 18153: better __int__ and add support for -infinity
​e34e9ca	Trac 18152: fix documentation

[git]

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

​376bcc2

Trac 18152: fix constructor in rings.number_field.number_field

[vdelecroix]

Now all test pass!

[git]

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

​e968468

Trac 18152: better cmp + 32 vs 64 bits in doctest

[vdelecroix]

A test file to compare timings

[vdelecroix]

I did some concrete timings and the conclusion is that in all cases the new implementation is faster by a factor of at least x8 (and it can be up to x20).

[git]

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

​9dd7433

Trac 12158: add a _sub_ method

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​e58f755

Trac 18152: add a _sub_ method

[git]

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

​e4a3524

Trac 18152: speedup (use X._add_(Y) instead of X+Y)

[jipilab]

Hi!

Did you have a look at #16116? Does it improve also the multiplication of generic cyclotomic matrices?

Jean-Philippe

[vdelecroix]

Hello Jean-Philippe,

Replying to jipilab:

Did you have a look at #16116? Does it improve also the multiplication of generic cyclotomic matrices?

I guess so. Should be ten times faster. But I copied-pasted the GAP documentation in the top of the file:

Arithmetical operations are quite expensive, so the use of
internally represented cyclotomics is not recommended for doing
arithmetic over number fields, such as calculations with matrices
of cyclotomics.

Best, Vincent

[stumpc5]

Hi there!

I kind of feel bad that you plan to remove that much code I wrote in weeks of work -- but I guess that's the way it goes if there is something better...

Anyway my first question would be whether it is really desired to remove an algorithm from Sage and instead use the same algorithm from Gap? Don't you think that working on the cython part of my implementation (maybe even moving to plain C of some core functionality) could get it beyond the speed of the black box (from a Sage point of view) Gap implementation?

Secondly, I remember that I did quite some tests back when I implemented the UCF. And that in the end, the Sage implementation was slower on some parts, while faster on others (one I remember where Sage was quicker was the computation of the Galois conjugates of big elements, in case you want to test, one example might be E( 2^11 * 3^4 )). I would really like to redo these tests, but haven't yet been able to get your branch running on my computer. I will get back once that worked out.

I also remember to have found a few bugs in the Gap implementation while implementing UCFs in Sage. That would also not be possible anymore when only relying on the implementation in Gap.

Cheers, Christian

[vdelecroix]

Hello Christian,

Thanks for having a look!

Replying to stumpc5:

I kind of feel bad that you plan to remove that much code I wrote in weeks of work -- but I guess that's the way it goes if there is something better...

Me too. I wondered if something were better with your code. But except the documentation I did not find any case were it was faster.

Anyway my first question would be whether it is really desired to remove an algorithm from Sage and instead use the same algorithm from Gap?

Sage should not reinvent the wheel. Gap is now interfaced as a fairly good level and we should try to maximize its usage at what it is good for.

That being said, it would makes sense to keep your implementation somewhere with more comparisons to the libgap implementation. Especially if you found that the Gap implementation has some bug! But I did not find any relevant example since all tests pass with my branch applied.

Don't you think that working on the cython part of my implementation (maybe even moving to plain C of some core functionality) could get it beyond the speed of the black box (from a Sage point of view) Gap implementation?

No idea. And I will not try. It is much more work than proposing what I did. And if you ask me what I would do if I want even more speed, I would link to Gap at an even lower level.

Secondly, I remember that I did quite some tests back when I implemented the UCF. And that in the end, the Sage implementation was slower on some parts, while faster on others (one I remember were Sage was quicker was the computation of the Galois conjugates of big elements, in case you want to test, one example might be E( 2^11 * 3^4 )). I would really like to redo these tests, but haven't yet been able to get your branch running on my computer. I will get back once that worked out.

Great! I will add this to the ucf_test.py that I wrote. If you have other examples in mind, please propose. I am not sure it will make a difference now since big numbers with pexpect implies a lot of parsing and computation. With libgap it is not the case anymore.

I also remember to have found a few bugs in the Gap implementation while implementing UCFs in Sage. That would also not be possible anymore when only relying on the implementation in Gap.

You do not remember what they are? Did you report them to Gap?

Cheers Vincent

[stumpc5]

Hi Vincent,

Replying to vdelecroix:

Me too. I wondered if something were better with your code. But except the documentation I did not find any case were it was faster.

Did you already redo the Galois conjugate test which I remember was slower in Gap back then?

Great! I will add this to the ucf_test.py that I wrote. If you have other examples in mind, please propose. I am not sure it will make a difference now since big numbers with pexpect implies a lot of parsing and computation. With libgap it is not the case anymore.

I just did

sage: a = E( 2^11 * 3^4 )
sage: %time x = a.galois_conjugates()
CPU times: user 5.52 s, sys: 159 ms, total: 5.68 s
Wall time: 5.66 s
sage: %time x = a.galois_conjugates()
CPU times: user 6.81 s, sys: 170 ms, total: 6.98 s
Wall time: 6.99 s

and compared it to

gap> a := E( 2^11 * 3^4 );;
gap> L := PrimeResidues(165888);;
gap> for i in L do
> x := GaloisCyc(a,i);
> od;

which I had running for 30 sec before I ctrc-c'ed it. (I also believe that in Sage, a big portion of the time is spend of computing the list of coprimes.)

You do not remember what they are? Did you report them to Gap?

I did report them, and they were fixed.

Once I get your branch running, I try to get back other tests I did.

Cheers, Christian

[jipilab]

Hi there!

I will also do some tests on the ticket and let you know.

One thing not to forget is also the number of tickets repaired by this new implementation. If the ticket repairs them and is faster... It's difficult to say no to that.

Grüsse, JP

[vdelecroix]

Hello,

I pushed an experimental branch with both implementations at public/18152. It is certain that most doctest fail but it is not the purpose. With it you can do

sage: UCF_native.<E_native> = UniversalCyclotomicField_native()
sage: UCF_libgap.<E_libgap> = UniversalCyclotomicField_libgap()
sage: UCF_native
Universal Cyclotomic Field (native)
sage: UCF_libgap
Universal Cyclotomic Field (libgap)

and play around with both implementations.

Replying to 20: it is true that GaloisCyc is slow in GAP (and it is not PrimeResidues which takes time). I got with a little bit smaller example

sage: a = E_libgap( 2^9 * 3^4 )
sage: b = E_native( 2^9 * 3^4 )
sage: %time l = a.galois_conjugates()
CPU times: user 6.92 s, sys: 4 ms, total: 6.92 s
Wall time: 6.9 s
sage: %time l = b.galois_conjugates()
CPU times: user 352 ms, sys: 52 ms, total: 404 ms
Wall time: 387 ms

Replying to 21: I do not think that the issues in the four tickets mentioned in the description are that big (and one of them has a fix).

Vincent

[stumpc5]

the branch runs now...

One issue I see is that it would if my code would not break because a few methods disappeared (such as field_order and min_poly. But this problem is minor and easy to fix -- worst case is I update my code.

Redoing the following tests Sage vs. libgap: 1.

sage: a = E(5)
sage: %timeit a.conjugate()
10000 loops, best of 3: 62.4 µs per loop

vs.

sage: a = E(5)
sage: %timeit a.conjugate()
10000 loops, best of 3: 39.9 µs per loop

2.

sage: a = E( 2^10 )
sage: %timeit a.conjugate()
10000 loops, best of 3: 69.1 µs per loop

vs.

sage: sage: a = E( 2^10 )
sage: sage: %timeit a.conjugate()
10000 loops, best of 3: 75.1 µs per loop

3.

sage: a = E( 2^11 * 3^4 )
sage: %timeit a.conjugate()
10000 loops, best of 3: 118 µs per loop

vs.

sage: a = E( 2^11 * 3^4 )
sage: sage: %timeit a.conjugate()
100 loops, best of 3: 9.15 ms per loop

For the Galois conjugation, the libgap implementation doesn't finish in very finite time, see also your example above.

[stumpc5]

conjugate actually uses GaloisCyc as in ​http://www.math.rwth-aachen.de/~Greg.Gamble/gap4r3/doc/htm/ref/CHAP058.htm#SSEC002.4.

[vdelecroix]

I opened #18207 for this slowness.

[vdelecroix]

Replying to vdelecroix:

I opened #18207 for this slowness.

(most time seems to be spent in the function ConvertToBase in cyclotom.c)

[vdelecroix]

On the other hand your code becomes slower when the element is more complicated

sage: a = sum(E_libgap( 2^8 * 3^4, i) for i in range(0, 2**12, 37))
sage: b = sum(E_native( 2^8 * 3^4, i) for i in range(0, 2**12, 37))
sage: %timeit _ = a.conjugate()
1000 loops, best of 3: 262 µs per loop
sage: %timeit _ = b.conjugate()
1000 loops, best of 3: 793 µs per loop

So I would not say that conjugate is slower.

[vdelecroix]

Ok. The whole thing about conjugation just reflects the fact that Gap is using a dense representation where yours is sparse. In a.conjugate() Gap just run through all coefficients, looking for the nonzero ones. For cyclotomic elements with large conductor and few nonzero coefficients, it makes a big difference.

Is there a use case for sparse cyclotomic elements with large conductor?

[jipilab]

Hi,

So, all test pass on my 6.6.rc3. I went through the files and removed some trailing spaces and blank lines.

Replying to 28: I do not have a use case for sparse cyclotomic out of my head. It is also difficult to guess for the use I have in mind (mainly representation of Coxeter groups). Nonetheless, the matrices will not have large conductor. Thus I would say it is not too bad as it is for the purpose of Coxeter groups.

Going through the discussion, the current state of affairs is:

The following tickets are fixed: #14240, #16130, #16631, #17117

The following tickets are follow-ups: #16116,#18207

The ticket looks good to me. It is well doctested with many examples. I would put positive review, nevertheless it would be good if someone else also go over the ticket too, Christian?

Best, JP

---

New commits:

​591409b	Merge branch 'develop' into cyclotomic
​29cb0a7	Corrected some blank lines and trailing spaces

[vdelecroix]

There are merged conflicts on the last beta... Wait a minute, I am doing a rebase. I am doing a rebase.

[jipilab]

Oops, my bad... Should have done it with the 6.6 ...

[vdelecroix]

Salut Jean-Phillippe,

Thanks for your corrections! The branch is a fast forward. So you have to delete your previous branch. Your commit is preserved (it is now 08df5b4). I added three commits for better compatibility with the previous version:

• adding a minpoly and field_order methods as it was requested in comment:23 (see also #18266 for a potential improvement)
• adding a __neg__ (it is now faster to do -E(3))

Vincent

---

New commits:

​0c8af9b	Trac 18152: universal cyclotomic field using libgap
​509acf7	Trac 18152: remove the old implementation
​519dc00	Trac 18152: Fix coxeter group import
​08df5b4	Corrected some blank lines and trailing spaces
​890151f	Trac 18152: add a minpoly method to UCF elements
​70f7895	Trac 18152: field_order as a deprecated alias of conductor
​b8e369b	Trac 18152: __neg__ for UCF elements

[vdelecroix]

Replying to jipilab:

The following tickets are fixed: #14240, #16130, #16631, #17117

And actually, doctests are already there! So these should be closed as soon as this one gets positively reviewed.

Vincent

[vdelecroix]

But there is still the most important question: do we keep both versions (one being sparse the other being dense)?

Vincent

[git]

Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:

​ccfc26e	Trac 18152: universal cyclotomic field using libgap
​1b2d04a	Trac 18152: remove the old implementation
​f3ef31a	Trac 18152: Fix coxeter group import
​bb42614	Corrected some blank lines and trailing spaces
​e3f714e	Trac 18152: add a minpoly method to UCF elements
​8149677	Trac 18152: field_order as a deprecated alias of conductor
​28c23ce	Trac 18152: __neg__ for UCF elements

[vdelecroix]

All right. Updated on the very last 6.7.beta2!

[tscrim]

Replying to vdelecroix:

But there is still the most important question: do we keep both versions (one being sparse the other being dense)?

IMO, we should keep both with the default being the dense.

[vdelecroix]

Replying to tscrim:

Replying to vdelecroix:

But there is still the most important question: do we keep both versions (one being sparse the other being dense)?

IMO, we should keep both with the default being the dense.

I was implictely requiring arguments. We will not keep it just because you think it has to ;-) I do have arguments for both sides, but as I am the ticket author I better not participate.

[nthiery]

The main argument I could see for keeping the code around is as follow: The speed balance may evolve in the future, typically if free modules in Sage get optimized, or if new use cases emerge. Of course we could always revive the code then, but it would probably have rotten in the mean time, and we may have completely forgotten about it. It also can be used for testing purposes for comparing two independent implementations of UCF.

Those are not super strong arguments.

The main argument for not keeping it is that we would not want is wasting time maintaining it.

A sensible approach might be to write a brief comment at the beginning of the file / class stating something like: "at this point this code is not used much (see ...). If maintaining it becomes a bother in the future, it's ok to discard it".

[jipilab]

Perhaps one could add a comment/warning in the beginning of the file also to mention that GAP uses a dense representation and a link to the current ticket.

If there is a use case popping up in the future we could revive the code as Nicolas mentioned.

@Christian: In your actual usage (not tests), is sparse faster than the current dense representation of GAP?

Also: I just updated to 6.7.beta3 on my computer and pulled the ticket. There was a small conflict in the file

src/doc/en/reference/rings_standard/index.rst

Should I push the resolution of the conflict here?

[vdelecroix]

Replying to jipilab:

Also: I just updated to 6.7.beta3 on my computer and pulled the ticket. There was a small conflict in the file

src/doc/en/reference/rings_standard/index.rst

Should I push the resolution of the conflict here?

Yes please.

[jipilab]

Et voilà!

---

New commits:

​dc3de84

Merge branch 'u/vdelecroix/18152' of trac.sagemath.org:sage into SageMath6.7.b3

[jipilab]

Hi,

All test passed on 6.7.beta3!

Could you add a comment/warning at the beginning of the file mentioning the old implementation and a link to the ticket?

Perhaps it would be good to mention this in the documentation as well?

[vdelecroix]

Done. I also simplified some imports.

---

New commits:

​ba5ed57	Trac 18152: add a note about the change of implementation
​b818c83	Trac 18152: remove and simplify some imports

[jipilab]

Hi,

This looks good to go as I see it.

I set it to positive review.