Slightly fasten is_cyclotomic()

[bgrenet]

In src/sage/rings/polynomial/polynomial_element.pyx, the method is_cyclotomic() begins by checking the irreducibility of the input polynomial, and then checks the leading and constant coefficients. I simply propose to invert those tests to fasten a bit the test for polynomials that are not monic or whose constant coefficient is not 1.

CC: @fchapoton

Component: number theory

Keywords: cyclotomic polynomials

Author: Bruno Grenet

Branch/Commit: a2a6f02

Reviewer: Frédéric Chapoton

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

[bgrenet]

Branch: u/bruno/slightly_fasten_is_cyclotomic__

[bgrenet]

Commit: f79eeaa

[bgrenet]

Changed keywords from none to cyclotomic polynomials

[videlec]

comment:3

Hello Bruno,

Actually, I discussed this with Karim Belabas in Bordeaux few weeks ago. We should either call Pari or stick to their implementation. It is much faster:

sage: p = cyclotomic_polynomial(150)
sage: pp = p._gp_()
sage: timeit("pp.poliscyclo()")
625 loops, best of 3: 697 µs per loop
sage: timeit("p.is_cyclotomic()")
125 loops, best of 3: 2.78 ms per loop

Moreover they do have the function poliscycloprod that checks if all roots are roots of unity (which is even faster).

Vincent

[bgrenet]

comment:5

Replying to @videlec:

Hello Bruno,

Actually, I discussed this with Karim Belabas in Bordeaux few weeks ago. We should either call Pari or stick to their implementation. It is much faster:

sage: p = cyclotomic_polynomial(150)
sage: pp = p._gp_()
sage: timeit("pp.poliscyclo()")
625 loops, best of 3: 697 µs per loop
sage: timeit("p.is_cyclotomic()")
125 loops, best of 3: 2.78 ms per loop

Moreover they do have the function poliscycloprod that checks if all roots are roots of unity (which is even faster).

Vincent

It is not surprising indeed that Pari is faster for this. I don't understand what you mean by

either call Pari or stick to their implementation.

Do you mean there are two ways to use Pari's implementation?

[videlec]

comment:6

Replying to @bgrenet:

Replying to @videlec:
It is not surprising indeed that Pari is faster for this. I don't understand what you mean by

either call Pari or stick to their implementation.

Do you mean there are two ways to use Pari's implementation?

Nope. Either we should copy their implementation within Sage (I do not think this is justified) or just call Pari.

And by the way there are two ways of using pari. Either with gp (this is a pexpect interface with a lot of latency in communication) and with libpari (which is very efficient in C). But sadly there is not yet poliscyclo and poliscycloprod in the libpari interface (i.e. in the file libs/pari/gen.pyx). I guess that we should wait for #17631 and not modify gen.pyx directly.

Vincent

[fchapoton]

comment:7

A side remark:

I would like to have a certificate (something like : this is the cyclotomic polynomial Phi_20)

Is this allowed by the pari function ?

[videlec]

comment:8

Replying to @fchapoton:

A side remark:

I would like to have a certificate (something like : this is the cyclotomic polynomial Phi_20)

Is this allowed by the pari function ?

Of course ;-) From pari documentation

poliscyclo(f): returns 0 if f is not a cyclotomic polynomial,
               and n > 0 if f = Phi_n, the n-th cyclotomic polynomial.

and you can already use it

sage: R.<x> = ZZ[]
sage: (x^2+x+1)._gp_().poliscyclo()
3
sage: (x^2-1)._gp_().poliscyclo()
0

Vincent

[bgrenet]

comment:9

Replying to @videlec:

Replying to @bgrenet:

Replying to @videlec:
It is not surprising indeed that Pari is faster for this. I don't understand what you mean by

either call Pari or stick to their implementation.

Do you mean there are two ways to use Pari's implementation?

Nope. Either we should copy their implementation within Sage (I do not think this is justified) or just call Pari.

And by the way there are two ways of using pari. Either with gp (this is a pexpect interface with a lot of latency in communication) and with libpari (which is very efficient in C). But sadly there is not yet poliscyclo and poliscycloprod in the libpari interface (i.e. in the file libs/pari/gen.pyx). I guess that we should wait for #17631 and not modify gen.pyx directly.

Vincent

For concreteness, what do you propose to do now? If I understand correctly, you propose to wait for #17631 and then make is_cyclotomic() directly call the pari functions, and not introduce the tiny change I made? I am fine with that, just wanted to be sure!

[videlec]

comment:10

Replying to @bgrenet:

For concreteness, what do you propose to do now? If I understand correctly, you propose to wait for #17631 and then make is_cyclotomic() directly call the pari functions, and not introduce the tiny change I made? I am fine with that, just wanted to be sure!

For now, we can already do the thing I proposed with the gp interface (and return the certificate instead of True/False). And we can also add a TODO section in the documentation saying that we should move to the pari interface as soon as #17631 is ready.

Vincent

[bgrenet]

comment:11

Replying to @videlec:

Replying to @bgrenet:

For concreteness, what do you propose to do now? If I understand correctly, you propose to wait for #17631 and then make is_cyclotomic() directly call the pari functions, and not introduce the tiny change I made? I am fine with that, just wanted to be sure!

For now, we can already do the thing I proposed with the gp interface (and return the certificate instead of True/False). And we can also add a TODO section in the documentation saying that we should move to the pari interface as soon as #17631 is ready.

I am not sure to agree since it actually makes the function slower:

sage: p = cyclotomic_polynomial(150)
sage: %timeit p.is_cyclotomic()
100 loops, best of 3: 2.75 ms per loop
sage: %timeit p._gp_().poliscyclo()
100 loops, best of 3: 3.96 ms per loop

As an aside, isn't it a problem for a method named is_cyclotomic() to return an integer rather than a boolean? Of course, the certificate is interesting to have, and one can still test with is_cyclotomic() as before. It is just that I expect a method is_...() to return a boolean. Maybe I am wrong...

[videlec]

comment:12

Replying to @bgrenet:

Replying to @videlec:

Replying to @bgrenet:

For concreteness, what do you propose to do now? If I understand correctly, you propose to wait for #17631 and then make is_cyclotomic() directly call the pari functions, and not introduce the tiny change I made? I am fine with that, just wanted to be sure!

For now, we can already do the thing I proposed with the gp interface (and return the certificate instead of True/False). And we can also add a TODO section in the documentation saying that we should move to the pari interface as soon as #17631 is ready.

I am not sure to agree since it actually makes the function slower:

sage: p = cyclotomic_polynomial(150)
sage: %timeit p.is_cyclotomic()
100 loops, best of 3: 2.75 ms per loop
sage: %timeit p._gp_().poliscyclo()
100 loops, best of 3: 3.96 ms per loop

Right. Most of the time seems to be lost in the conversion

sage: l = [cyclotomic_polynomial(i) for i in range(2,50)]
sage: timeit("for p in l: _ = p.is_cyclotomic()")
5 loops, best of 3: 51.7 ms per loop
sage: timeit("for p in l: _ = p._gp_().poliscyclo()")
5 loops, best of 3: 303 ms per loop
sage: lgp = [p._gp_() for p in l]
sage: timeit("for p in lgp: _ = p.poliscyclo()")
25 loops, best of 3: 34.8 ms per loop

As an aside, isn't it a problem for a method named is_cyclotomic() to return an integer rather than a boolean? Of course, the certificate is interesting to have, and one can still test with is_cyclotomic() as before. It is just that I expect a method is_...() to return a boolean. Maybe I am wrong...

You are right. The best option would be:

def is_cyclotomic(self, certificate=False):
    ...
    return ans if certificate else bool(ans)

A compromise can be:

• add the keyword certificate (if it is True we use self._gp_().polisyclco())
• do your initial suggested change
• add a TODO to the documentation to switch from gp to pari (and open a ticket for it). We will then consider the fact of making it the default.

Vincent

[bgrenet]

comment:13

Replying to @videlec:

A compromise can be:

• add the keyword certificate (if it is True we use self._gp_().polisyclco())
• do your initial suggested change
• add a TODO to the documentation to switch from gp to pari (and open a ticket for it). We will then consider the fact of making it the default.

I like this proposition, I'll do it.

[sagetrac-git]

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

053a30c

Add certificate keyword

[sagetrac-git]

Changed commit from f79eeaa to 053a30c

[videlec]

comment:15

Just typographical comment. The TODO block must be indented. In other words

.. TODO::

    As soon as #XYZ is finished we should do blih and blah.

Have a look to Docstring section of the developer manual.

It is _gp_() and not __gp__() (i.e. one underscore).

I would move the ._gp_().poliscyclo() below the tests for characteristic.

Vincent

[bgrenet]

comment:16

I pushed too fast, and noticed my mistake on __gp__(). Actually, there are other ones. I'm correcting this.

Thanks anyway for the feedback, and especially for the documentation!

[sagetrac-git]

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

c7b5f68

Correcting mistakes + docstring formatting

[sagetrac-git]

Changed commit from 053a30c to c7b5f68

[videlec]

comment:19

One more thing for the doc: for trac tickets you should use :trac:`17730`. That will create automatically a link in the compiled documentation.

More seriously, we definitely do not want to test irreducibility at the beginning. You can possibly add this test in between step 2 and 3. That would prevent some useless factorization. But if you do so, I want to see some benchmarks first.

Vincent

[bgrenet]

comment:21

Replying to @videlec:

One more thing for the doc: for trac tickets you should use :trac:`17730`. That will create automatically a link in the compiled documentation.

OK.

More seriously, we definitely do not want to test irreducibility at the beginning. You can possibly add this test in between step 2 and 3. That would prevent some useless factorization. But if you do so, I want to see some benchmarks first.

I do not understand your remark: irreducibility was tested as the first step before, and moving this test later is precisely the purpose of my ticket.

[videlec]

comment:22

Replying to @bgrenet:

More seriously, we definitely do not want to test irreducibility at the beginning. You can possibly add this test in between step 2 and 3. That would prevent some useless factorization. But if you do so, I want to see some benchmarks first.

I do not understand your remark: irreducibility was tested as the first step before, and moving this test later is precisely the purpose of my ticket.

Right ;-) I will try a bit more optimization on my side.

[sagetrac-git]

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

1b01fad

improved formatting (trac ticket)

[sagetrac-git]

Changed commit from c7b5f68 to 1b01fad

[videlec]

comment:25

Actually, the factor at the end is overkill as we only want a square root. I did have a look at the method .is_square() and it was dirty enough to be slow. So I just make it simpler/clearer. The speedup is very small and concerns mostly positive cases.

I am happy with your commits. Just review mine.

Vincent

---

New commits:

90ac73a

trac #17695: simplify factor step + faster .sqrt()

[videlec]

Changed branch from u/bruno/slightly_fasten_is_cyclotomic__ to u/vdelecroix/17695

[videlec]

Changed commit from 1b01fad to 90ac73a

[bgrenet]

Changed branch from u/vdelecroix/17695 to u/bruno/17695

[bgrenet]

comment:27

Replying to @videlec:

Actually, the factor at the end is overkill as we only want a square root. I did have a look at the method .is_square() and it was dirty enough to be slow. So I just make it simpler/clearer. The speedup is very small and concerns mostly positive cases.

I am happy with your commits. Just review mine.

I am fine with your commit too, but for one typo that I've corrected in the docstring of is_square():

 def is_square(self, root=False):
"""
Returns whether or not polynomial is square. If the optional
- argument ``root`` is set to ``True, then also returns the square root
+ argument ``root`` is set to ``True``, then also returns the square root
(or ``None``, if the polynomial is not square).

I do not know whether it is acceptable for me to set "positive review". I let you do so if you consider everything's OK. Thanks for the patient feedback!

---

New commits:

85d8ffc

Review commit

[bgrenet]

Changed commit from 90ac73a to 85d8ffc

[fchapoton]

comment:28

Looks good to me. I have just made a few typographical changes, so I allow myself to put this into positive review.

---

New commits:

a2a6f02

trac #17695 some doc and code formatting details

[fchapoton]

Changed commit from 85d8ffc to a2a6f02

[fchapoton]

Changed branch from u/bruno/17695 to public/ticket/17695

[videlec]

comment:29

Ho great! Nice catch!

Vincent

[mezzarobba]

Reviewer: Frédéric Chapoton

[vbraun]

Changed branch from public/ticket/17695 to a2a6f02