Opened 7 years ago

Closed 7 years ago

#17695 closed enhancement (fixed)

Slightly fasten is_cyclotomic()

Reported by: bruno Owned by:
Priority: minor Milestone: sage-6.5
Component: number theory Keywords: cyclotomic polynomials
Cc: chapoton Merged in:
Authors: Bruno Grenet Reviewers: Frédéric Chapoton
Report Upstream: N/A Work issues:
Branch: a2a6f02 (Commits, GitHub, GitLab) Commit: a2a6f026794fb12e1d557d8ca6738c45d5ed937f
Dependencies: Stopgaps:

Status badges

Description (last modified by bruno)

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.

Change History (31)

comment:1 Changed 7 years ago by bruno

  • Branch set to u/bruno/slightly_fasten_is_cyclotomic__

comment:2 Changed 7 years ago by bruno

  • Commit set to f79eeaaa630b1bd966a25a2f7b0f7a5f1dda3378
  • Description modified (diff)
  • Keywords cyclotomic polynomials added
  • Status changed from new to needs_review

comment:3 follow-up: Changed 7 years ago by vdelecroix

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

Last edited 7 years ago by vdelecroix (previous) (diff)

comment:4 Changed 7 years ago by chapoton

  • Cc chapoton added

comment:5 in reply to: ↑ 3 ; follow-up: Changed 7 years ago by bruno

Replying to vdelecroix:

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?

comment:6 in reply to: ↑ 5 ; follow-up: Changed 7 years ago by vdelecroix

Replying to bruno:

Replying to vdelecroix: 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

comment:7 follow-up: Changed 7 years ago by chapoton

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 ?

comment:8 in reply to: ↑ 7 Changed 7 years ago by vdelecroix

Replying to chapoton:

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

Last edited 7 years ago by vdelecroix (previous) (diff)

comment:9 in reply to: ↑ 6 ; follow-up: Changed 7 years ago by bruno

Replying to vdelecroix:

Replying to bruno:

Replying to vdelecroix: 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!

comment:10 in reply to: ↑ 9 ; follow-up: Changed 7 years ago by vdelecroix

Replying to bruno:

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

comment:11 in reply to: ↑ 10 ; follow-up: Changed 7 years ago by bruno

Replying to vdelecroix:

Replying to bruno:

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

comment:12 in reply to: ↑ 11 ; follow-up: Changed 7 years ago by vdelecroix

Replying to bruno:

Replying to vdelecroix:

Replying to bruno:

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

comment:13 in reply to: ↑ 12 Changed 7 years ago by bruno

Replying to vdelecroix:

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.

comment:14 Changed 7 years ago by git

  • Commit changed from f79eeaaa630b1bd966a25a2f7b0f7a5f1dda3378 to 053a30cc574fb122ec8527c53849419ccf871e68

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

053a30cAdd certificate keyword

comment:15 Changed 7 years ago by vdelecroix

  • Status changed from needs_review to needs_work

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

comment:16 Changed 7 years ago by bruno

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!

comment:17 Changed 7 years ago by git

  • Commit changed from 053a30cc574fb122ec8527c53849419ccf871e68 to c7b5f6861c5bbabadd9a9e68cdc8c6d7a5673c50

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

c7b5f68Correcting mistakes + docstring formatting

comment:18 Changed 7 years ago by bruno

  • Status changed from needs_work to needs_review

comment:19 follow-up: Changed 7 years ago by vdelecroix

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

Last edited 7 years ago by vdelecroix (previous) (diff)

comment:20 Changed 7 years ago by vdelecroix

  • Status changed from needs_review to needs_work

comment:21 in reply to: ↑ 19 ; follow-up: Changed 7 years ago by bruno

Replying to vdelecroix:

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.

comment:22 in reply to: ↑ 21 Changed 7 years ago by vdelecroix

Replying to bruno:

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.

comment:23 Changed 7 years ago by git

  • Commit changed from c7b5f6861c5bbabadd9a9e68cdc8c6d7a5673c50 to 1b01fad201f72a254845e763cf8b072d28cde8b8

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

1b01fadimproved formatting (trac ticket)

comment:24 Changed 7 years ago by bruno

  • Status changed from needs_work to needs_review

comment:25 follow-up: Changed 7 years ago by vdelecroix

  • Branch changed from u/bruno/slightly_fasten_is_cyclotomic__ to u/vdelecroix/17695
  • Commit changed from 1b01fad201f72a254845e763cf8b072d28cde8b8 to 90ac73a342710459415a32442d418751e3560be6

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:

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

comment:26 Changed 7 years ago by bruno

  • Branch changed from u/vdelecroix/17695 to u/bruno/17695

comment:27 in reply to: ↑ 25 Changed 7 years ago by bruno

  • Commit changed from 90ac73a342710459415a32442d418751e3560be6 to 85d8ffcdd69d911dcc042e4dd3ab8374f90e4d56

Replying to vdelecroix:

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:

85d8ffcReview commit

comment:28 Changed 7 years ago by chapoton

  • Branch changed from u/bruno/17695 to public/ticket/17695
  • Commit changed from 85d8ffcdd69d911dcc042e4dd3ab8374f90e4d56 to a2a6f026794fb12e1d557d8ca6738c45d5ed937f
  • Status changed from needs_review to positive_review

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


New commits:

a2a6f02trac #17695 some doc and code formatting details

comment:29 Changed 7 years ago by vdelecroix

Ho great! Nice catch!

Vincent

comment:30 Changed 7 years ago by mmezzarobba

  • Reviewers set to Frédéric Chapoton

comment:31 Changed 7 years ago by vbraun

  • Branch changed from public/ticket/17695 to a2a6f026794fb12e1d557d8ca6738c45d5ed937f
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.