Opened 6 years ago

Closed 6 years ago

#16878 closed enhancement (fixed)

faster is_prime

Reported by: vdelecroix Owned by:
Priority: major Milestone: sage-6.6
Component: number theory Keywords:
Cc: nochen Merged in:
Authors: Vincent Delecroix, Jeroen Demeyer Reviewers: Jeroen Demeyer, Vincent Delecroix
Report Upstream: N/A Work issues:
Branch: 1f8abd9 (Commits) Commit: 1f8abd999dcf4935f92b2d83aad01684397e62d1
Dependencies: #16997 Stopgaps:

Description (last modified by jdemeyer)

Right now to test if a Sage integer is prime it is faster to call prime_range rather than .is_prime()...

sage: timeit("bool(prime_range(121,122))", number=10000)
10000 loops, best of 3: 1.09 µs per loop
sage:  timeit("bool(prime_range(1009,1010))", number=10000)
10000 loops, best of 3: 1.2 µs per loop

versus

sage: timeit("121.is_prime()", number=10000)
10000 loops, best of 3: 5.3 µs per loop
sage: timeit("1009.is_prime()", number=10000)
10000 loops, best of 3: 4.19 µs per loop

The patch does some tiny modifications in integer.pyx and we get

sage: timeit("121.is_prime()", number=10000)
10000 loops, best of 3: 812 ns per loop
sage: timeit("1009.is_prime()", number=10000)
10000 loops, best of 3: 730 ns per loop

We also modify is_prime_power() to return False for 1 and raise an error for non-integral rationals like 1/2.

See also this sage-devel discussion.

Change History (70)

comment:1 Changed 6 years ago by vdelecroix

  • Branch set to u/vdelecroix/16878
  • Commit set to dd44a7e8d05a2b3978c679c363a2729bf0cf5e93
  • Status changed from new to needs_review

New commits:

dd44a7etrac #16878: faster is_prime

comment:2 Changed 6 years ago by vdelecroix

  • Description modified (diff)

comment:3 follow-up: Changed 6 years ago by jdemeyer

  • Reviewers set to Jeroen Demeyer
  • Status changed from needs_review to needs_work

I would prefer not to include stuff from pari/pari.h directly, but instead use the declarations from src/sage/libs/pari/decl.pxi. And please rebase this over #15767.

comment:4 Changed 6 years ago by jdemeyer

Can you also change the value of PARI_PSEUDOPRIME_LIMIT (see PARI docs) to 2^64 and use mpz_fits_ui() before calling mpz_get_ui()

comment:5 Changed 6 years ago by jdemeyer

Instead of changing _pari_() to _pari_c() everywhere, wouldn't

cpdef _pari_(self)

accomplish the same thing?

comment:6 in reply to: ↑ 3 Changed 6 years ago by vdelecroix

Thanks for the remark. One question: are we sure that 264 fits into an unsigned long on any platform? Otherwise we can not safely call uisprime and uisprimepower. If not I will use a double check

if mpz_cmp(PARI_PSEUDOPRIME_LIMIT, v) > 0 and mpz_cmp_ui(v, ULONG_MAX) < 0:
     # use uisprime or uisprimepower

Vincent

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

comment:7 Changed 6 years ago by git

  • Commit changed from dd44a7e8d05a2b3978c679c363a2729bf0cf5e93 to d7dc38918903a97af742a90050ce87c58f6f222a

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

37fc8e8Upgrade to PARI-2.7.1
5db54b6Trac 15767: reviewer patch
1d103daTrac 15767: fix doctest in Ser()
62ca821Trac 15767: explain application of Sturm's theorem
bf435f8Merge tag '6.4.beta1' into ticket/15767
6e9f686trac #16878: faster is_prime
d7dc389trac #16878: review

comment:8 Changed 6 years ago by vdelecroix

  • Dependencies set to #15767
  • Status changed from needs_work to needs_review

comment:9 Changed 6 years ago by jdemeyer

  • Status changed from needs_review to needs_work

10000000000000000 != 2^64 :-)

I don't think we should lower PARI_PSEUDOPRIME_LIMIT because it's also used in places where it's not bounded by ULONG_MAX, see next_prime() for example.

Personally, I would assume that ULONG_MAX < 2^64 and just use mpz_fits_uint_p().

comment:10 Changed 6 years ago by jdemeyer

Also, you wrote "set set" instead of "set".

comment:11 Changed 6 years ago by git

  • Commit changed from d7dc38918903a97af742a90050ce87c58f6f222a to 175d6b5fd339989df12ec9c1f26d3b0d7e6dd82c

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

175d6b5trac #16878: review 2

comment:12 Changed 6 years ago by vdelecroix

  • Status changed from needs_work to needs_review

comment:13 Changed 6 years ago by jdemeyer

This is bikeshedding I know, but wouldn't

mpz_ui_pow_ui(PARI_PSEUDOPRIME_LIMIT, 2, 64)

be easier to understand?

comment:14 Changed 6 years ago by git

  • Commit changed from 175d6b5fd339989df12ec9c1f26d3b0d7e6dd82c to d55dc6737936077eeee903687207f7af54588439

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

d55dc67trac #16878: bikeshedding

comment:15 Changed 6 years ago by vdelecroix

  • Commit changed from 175d6b5fd339989df12ec9c1f26d3b0d7e6dd82c to d55dc6737936077eeee903687207f7af54588439

I don't mind...


New commits:

d55dc67trac #16878: bikeshedding

comment:16 Changed 6 years ago by vdelecroix

ping! (applies cleanly on 6.4.beta2)

comment:17 Changed 6 years ago by jdemeyer

  • Description modified (diff)

Big reviewer patch coming up...

comment:18 Changed 6 years ago by jdemeyer

  • Description modified (diff)

comment:19 Changed 6 years ago by jdemeyer

  • Description modified (diff)

comment:20 Changed 6 years ago by jdemeyer

  • Branch changed from u/vdelecroix/16878 to u/jdemeyer/ticket/16878
  • Created changed from 08/25/14 10:31:53 to 08/25/14 10:31:53
  • Modified changed from 08/29/14 10:34:32 to 08/29/14 10:34:32

comment:21 Changed 6 years ago by jdemeyer

  • Authors changed from Vincent Delecroix to Vincent Delecroix, Jeroen Demeyer
  • Commit changed from d55dc6737936077eeee903687207f7af54588439 to 7f6132188b7ddafd97630aaf721875f9d60a7106

As you can see, I made quite a lot of additional changes, mostly trying to make the code cleaner. I am currently running doctests.

Can you review this extra commit?


New commits:

7f61321Reorganize is_prime_power() and next_prime() logic

comment:22 Changed 6 years ago by git

  • Commit changed from 7f6132188b7ddafd97630aaf721875f9d60a7106 to 5937af54afa0562007d710c495ecd0f1160f7498

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

5937af5Sloane considers 1 to be a prime power

comment:23 follow-up: Changed 6 years ago by vdelecroix

  • Branch changed from u/jdemeyer/ticket/16878 to u/vdelecroix/16878
  • Commit changed from 5937af54afa0562007d710c495ecd0f1160f7498 to f0d716fd34029f65d0be0aa45d5cb1783f8da2f5

Great! I especially like that now isprimepower is accessible.

In my last commit, I added the keyword certificate to is_prime_power in order to get back the pair (n,p) such that the integer is pn (following pari convention). I think it is worth it!

Vincent


New commits:

f0d716ftrac #16878: certificate flag in is_prime_power

comment:24 follow-up: Changed 6 years ago by vdelecroix

Note the treatment of 1 is different inis_prime_powers and prime_powers... either we fix it in that ticket or it can wait #16880.

Vincent

comment:25 in reply to: ↑ 23 ; follow-up: Changed 6 years ago by jdemeyer

Replying to vdelecroix:

Great! I especially like that now isprimepower is accessible.

In my last commit, I added the keyword certificate to is_prime_power in order to get back the pair (n,p) such that the integer is pn (following pari convention). I think it is worth it!

I like the idea, but I'm not sure that I like the interface with the certificate keyword.

How about a separate method write_as_prime_power() (raising ArithmeticError if it's not a prime power) or something?

comment:26 in reply to: ↑ 24 Changed 6 years ago by jdemeyer

Replying to vdelecroix:

Note the treatment of 1 is different inis_prime_powers and prime_powers... either we fix it in that ticket or it can wait #16880.

Let's ignore that for now and try to merge this ticket first.

comment:27 in reply to: ↑ 25 Changed 6 years ago by vdelecroix

Replying to jdemeyer:

Replying to vdelecroix:

Great! I especially like that now isprimepower is accessible.

In my last commit, I added the keyword certificate to is_prime_power in order to get back the pair (n,p) such that the integer is pn (following pari convention). I think it is worth it!

I like the idea, but I'm not sure that I like the interface with the certificate keyword.

How about a separate method write_as_prime_power() (raising ArithmeticError if it's not a prime power) or something?

Note that there is already factor that does the job. Moreover, having two methods (is_prime_power and write_as_prime_power) with the exact same code looks bad as well. From my implementation, it is clear that the certificate keyword only modifies the output.

My design choice also comes from the fact that this certifcate keyword is already in use (see sage.graphs or sage.combinat.designs). I am not particularly happy with it but I find it better than having two functions.

On the other hand, treating a Python exception takes time, and this is_prime_power is the kind of function that I might use in time critical code.

Vincent

comment:28 Changed 6 years ago by jdemeyer

  • Branch changed from u/vdelecroix/16878 to u/jdemeyer/ticket/16878
  • Modified changed from 08/29/14 14:20:49 to 08/29/14 14:20:49

comment:29 Changed 6 years ago by jdemeyer

  • Commit changed from f0d716fd34029f65d0be0aa45d5cb1783f8da2f5 to bfd69d3b188a8ec36ab373b4e958b54c74b64d2b

One more reviewer patch...


New commits:

bfd69d3is_prime_power() reviewer fixes

comment:30 Changed 6 years ago by vdelecroix

I am happy with the current version (at commit bfd69d3). Do you agree to set to positive review?

Vincent

comment:31 Changed 6 years ago by jdemeyer

  • Reviewers changed from Jeroen Demeyer to Jeroen Demeyer, Vincent Delecroix
  • Status changed from needs_review to positive_review

comment:32 Changed 6 years ago by vbraun

  • Branch changed from u/jdemeyer/ticket/16878 to bfd69d3b188a8ec36ab373b4e958b54c74b64d2b
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:33 Changed 6 years ago by vbraun

  • Commit bfd69d3b188a8ec36ab373b4e958b54c74b64d2b deleted
  • Resolution fixed deleted
  • Status changed from closed to new

This takes about a second before this ticket and really long after this ticket:

sage: with proof.WithProof('arithmetic', False):  
....:    sage.rings.arith.is_prime_power((10^1000 + 453)^2)

In fact, it get so slow that the ARM buildbot times out in

sage -t --long src/sage/rings/finite_rings/constructor.py
    Timed out
**********************************************************************
...
sage: k = FiniteField(10^1000 + 453, proof=False) ## line 230 ##
sage: k = FiniteField((10^1000 + 453)^2, 'a', proof=False)      # long time -- about 5 seconds ## line 231 ##

**********************************************************************

comment:34 follow-up: Changed 6 years ago by vdelecroix

very strange indeed

sage: sage.rings.arith.is_prime_power((10^999 + 453)^2)   # very fast
False
sage: sage.rings.arith.is_prime_power((10^1001 + 453)^2)  # very fast
False
sage: sage.rings.arith.is_prime_power((10^1000 + 453)^2)
... <still waiting> ...

comment:35 in reply to: ↑ 34 Changed 6 years ago by vdelecroix

I see... seems to be because we do not care about the flag proof=True or proof=False anymore in is_prime_power. We should still rely on is_pseudoprime_small_power when proof=False.

And by the way, in is_pseudoprime_small_power the choice of the keyword to return (n,p) is get_data. I guess we should change the certificate in is_prime_power to get_data as well.

Vincent

comment:36 Changed 6 years ago by vdelecroix

  • Branch changed from bfd69d3b188a8ec36ab373b4e958b54c74b64d2b to u/vdelecroix/16878
  • Commit set to e6f88c5f9adef669ede74e3c161e19fe37203276
  • Status changed from new to needs_review

New commits:

6e9f686trac #16878: faster is_prime
d7dc389trac #16878: review
175d6b5trac #16878: review 2
d55dc67trac #16878: bikeshedding
7f61321Reorganize is_prime_power() and next_prime() logic
5937af5Sloane considers 1 to be a prime power
f0d716ftrac #16878: certificate flag in is_prime_power
bfd69d3is_prime_power() reviewer fixes
e6f88c5trac #16878: test is_pseudoprime_small_power + clean

comment:37 follow-up: Changed 6 years ago by jdemeyer

In is_pseudoprime_small_power, the fixed exponent bound of 1024 is clearly not acceptable.

The PARI implementation of isprimepower() is as follows: first do trial division by small primes up to 100. If we do not find a divisor, then find out if the number is a perfect power. This is now a lot cheaper since the exponent can be at most log(n)/log(p), where n is the number to check and p the smallest prime that we didn't check in the trial division.

But overall, I would prefer PARI to add a new function ispseudoprimepower() such that Sage can interface it.

comment:38 in reply to: ↑ 37 Changed 6 years ago by vdelecroix

Replying to jdemeyer:

In is_pseudoprime_small_power, the fixed exponent bound of 1024 is clearly not acceptable.

The PARI implementation of isprimepower() is as follows: first do trial division by small primes up to 100. If we do not find a divisor, then find out if the number is a perfect power. This is now a lot cheaper since the exponent can be at most log(n)/log(p), where n is the number to check and p the smallest prime that we didn't check in the trial division.

But overall, I would prefer PARI to add a new function ispseudoprimepower() such that Sage can interface it.

Me too. I saw that you already submit the question http://pari.math.u-bordeaux.fr/archives/pari-dev-1409/msg00011.html...

On the other hand, we need a reasonable fix in order to get this ticket in.

Vincent

comment:39 follow-up: Changed 6 years ago by vbraun

Jeroen, did you end up writing a pari patch? Maybe we can use the current is_pseudoprime_small_power until there is a better implementation?

comment:40 in reply to: ↑ 39 Changed 6 years ago by jdemeyer

Replying to vbraun:

Jeroen, did you end up writing a pari patch?

Not yet, but it would be trivial to do that. I will add it to #16997.

comment:41 Changed 6 years ago by jdemeyer

On the other hand, perhaps better add the patch here...

comment:42 Changed 6 years ago by jdemeyer

...or a new ticket.

comment:43 Changed 6 years ago by jdemeyer

Turns out adding the PARI patch to #16997 was the easiest (since there would have been merge conflicts otherwise). It's now done.

comment:44 Changed 6 years ago by jdemeyer

  • Status changed from needs_review to needs_work

I have now a version of #16997 which is ready for review. For this ticket, I would wait until that ticket is settled and then have a look again at this.

comment:45 Changed 6 years ago by vdelecroix

  • Dependencies changed from #15767 to #15767, #16997

comment:46 Changed 6 years ago by vdelecroix

Hi,

Good news: Karim recently committed a ispseudoprimepower in pari (commit 0f20474).

Vincent

comment:47 Changed 6 years ago by jdemeyer

  • Dependencies changed from #15767, #16997 to #16997, #17852

comment:48 Changed 6 years ago by jdemeyer

  • Branch changed from u/vdelecroix/16878 to u/jdemeyer/16878

comment:49 Changed 6 years ago by jdemeyer

  • Commit changed from e6f88c5f9adef669ede74e3c161e19fe37203276 to 9751ee8ccf4887970ed5cea1e6f52447b633e2c3
  • Dependencies changed from #16997, #17852 to #16997
  • Status changed from needs_work to needs_review

New commits:

ade4f75Upgrade PARI to git master
9751ee8Trac #16878: faster is_prime

comment:50 Changed 6 years ago by jdemeyer

  • Milestone changed from sage-6.4 to sage-6.6
  • Status changed from needs_review to needs_work

comment:51 Changed 6 years ago by vdelecroix

(I guess you did it because of conflicts with #17852)

comment:52 follow-up: Changed 6 years ago by jdemeyer

I'm doing various fixes, like remove is_pseudoprime_small_power.

I think it would be better to rebase #17852 on this ticket, since this ticket removes lots of things that #17852 fixes.

comment:53 in reply to: ↑ 52 Changed 6 years ago by vdelecroix

Replying to jdemeyer:

I'm doing various fixes, like remove is_pseudoprime_small_power.

I think it would be better to rebase #17852 on this ticket, since this ticket removes lots of things that #17852 fixes.

Fine for me.

comment:54 Changed 6 years ago by jdemeyer

I would also like to switch around the get_data result: when given p^n as input, I would find it more natural to return (p, n) instead of (n, p). What do you think?

comment:55 Changed 6 years ago by vdelecroix

gp is doing the other way around: it returns k as the main argument since it is compatible with the boolean answer. But doing (p,n) will be consistent with Integer.perfect_power so I would go for it.

comment:56 Changed 6 years ago by git

  • Commit changed from 9751ee8ccf4887970ed5cea1e6f52447b633e2c3 to 96f1f5037d18eb3b2e4c8e0e2a9c7a5cdac0d643

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

96f1f50Trac #16878: further fixes to is_prime() and friends

comment:57 Changed 6 years ago by jdemeyer

As far as I'm concerned, this is the final version. If you agree, then you can set this to positive_review.

I haven't yet run all doctests, doing that right now...

comment:58 Changed 6 years ago by jdemeyer

  • Status changed from needs_work to needs_review

comment:59 Changed 6 years ago by git

  • Commit changed from 96f1f5037d18eb3b2e4c8e0e2a9c7a5cdac0d643 to 1f8abd999dcf4935f92b2d83aad01684397e62d1

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

1f8abd9Doctest fix

comment:60 Changed 6 years ago by jdemeyer

Now all doctests pass.

comment:61 follow-up: Changed 6 years ago by vdelecroix

  • Status changed from needs_review to needs_info
  1. Now arith.is_prime_power and arith.is_pseudoprime_power do not handle negative powers anymore (i.e. rational numbers p^-n). This is different from previous versions.
  1. In Integer.is_prime_power why did you use a long for n and p instead of an unsigned long. Is pari uisprimepower not able to handle up to ULONG_MAX?

comment:62 in reply to: ↑ 61 ; follow-up: Changed 6 years ago by jdemeyer

Replying to vdelecroix:

  1. Now arith.is_prime_power and arith.is_pseudoprime_power do not handle negative powers anymore (i.e. rational numbers p^-n). This is different from previous versions.

Yes, this was discussed on sage-devel and the former behaviour was considered to be a bug => no deprecation needed

  1. In Integer.is_prime_power why did you use a long for n and p instead of an unsigned long. Is pari uisprimepower not able to handle up to ULONG_MAX?

I assume that pari can handle it, but smallInteger() for the result takes a long.

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

Replying to jdemeyer:

Replying to vdelecroix:

  1. Now arith.is_prime_power and arith.is_pseudoprime_power do not handle negative powers anymore (i.e. rational numbers p^-n). This is different from previous versions.

Yes, this was discussed on sage-devel and the former behaviour was considered to be a bug => no deprecation needed

True, I forgot about here and it is this thread.

  1. In Integer.is_prime_power why did you use a long for n and p instead of an unsigned long. Is pari uisprimepower not able to handle up to ULONG_MAX?

I assume that pari can handle it, but smallInteger() for the result takes a long.

I see... I do not like the fact that we avoid pari uisprimepower for the range between LONG_MAX and ULONG_MAX. Especially because otherwise we import the proof module which takes some non negligible amount of time. We can either get rid of smallInteger (that would be bad if either p or n would be below 256) or add a smallunsignedInteger

cdef inline Integer smallunsignedInteger(unsigned long value):
    """
    This is the fastest way to create a (likely) small positive Integer.
    """
    cdef Integer z
    if value <= small_pool_max:
        return <Integer>small_pool[value - small_pool_min]
    else:
        z = PY_NEW(Integer)
        mpz_set_ui(z.value, value)
        return z

What do you think?

comment:64 in reply to: ↑ 63 Changed 6 years ago by jdemeyer

To be honest, I think introducing smallunsignedInteger is just overkill. I don't think it makes a big difference if the limit for a fast is_prime_power is 2^63 or 2^64.

comment:65 Changed 6 years ago by jdemeyer

  • Description modified (diff)

comment:66 Changed 6 years ago by jdemeyer

  • Status changed from needs_info to needs_review

comment:67 follow-up: Changed 6 years ago by vdelecroix

  • Status changed from needs_review to needs_work

The buildbot is complaining about a digit

sage -t --long --warn-long 35.8 src/sage/libs/pari/gen.pyx
**********************************************************************
File "src/sage/libs/pari/gen.pyx", line 6237, in sage.libs.pari.gen.gen.ellheight
Failed example:
    e.ellheight([1,0], precision=128).sage()
Expected:
    0.47671165934373953737948605888465305945902294217            
Got:
    0.47671165934373953737948605888465305945902294218
**********************************************************************

Might be related to the doctest failure in pari update in #16997.

comment:68 in reply to: ↑ 67 Changed 6 years ago by jdemeyer

  • Status changed from needs_work to needs_review

Replying to vdelecroix:

The buildbot is complaining about a digit

Fixed by #16997.

Last edited 6 years ago by jdemeyer (previous) (diff)

comment:69 Changed 6 years ago by vdelecroix

  • Status changed from needs_review to positive_review

Merges cleanly with #16997 and all tests pass! Cool.

Vincent

comment:70 Changed 6 years ago by vbraun

  • Branch changed from u/jdemeyer/16878 to 1f8abd999dcf4935f92b2d83aad01684397e62d1
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.