faster is_prime

[videlec]

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.

Depends on #16997

CC: @nathanncohen

Component: number theory

Author: Vincent Delecroix, Jeroen Demeyer

Branch: 1f8abd9

Reviewer: Jeroen Demeyer, Vincent Delecroix

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

[videlec]

New commits:

dd44a7e

trac #16878: faster is_prime

[videlec]

Branch: u/vdelecroix/16878

[videlec]

Commit: dd44a7e

[videlec]

Description changed:

--- 
+++ 
@@ -15,7 +15,7 @@
 10000 loops, best of 3: 4.19 µs per loop

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

 sage: timeit("121.is_prime()", number=10000)

[jdemeyer]

comment:3

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.

[jdemeyer]

Reviewer: Jeroen Demeyer

[jdemeyer]

comment:4

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()

[jdemeyer]

comment:5

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

cpdef _pari_(self)

accomplish the same thing?

[videlec]

comment:6

Thanks for the remark. One question: are we sure that 2⁶⁴ 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

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from dd44a7e to d7dc389

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

37fc8e8	Upgrade to PARI-2.7.1
5db54b6	Trac 15767: reviewer patch
1d103da	Trac 15767: fix doctest in Ser()
62ca821	Trac 15767: explain application of Sturm's theorem
bf435f8	Merge tag '6.4.beta1' into ticket/15767
6e9f686	trac #16878: faster is_prime
d7dc389	trac #16878: review

[videlec]

Dependencies: #15767

[jdemeyer]

comment:9

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

[jdemeyer]

comment:10

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

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from d7dc389 to 175d6b5

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

175d6b5

trac #16878: review 2

[jdemeyer]

comment:13

This is bikeshedding I know, but wouldn't

mpz_ui_pow_ui(PARI_PSEUDOPRIME_LIMIT, 2, 64)

be easier to understand?

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 175d6b5 to d55dc67

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

d55dc67

trac #16878: bikeshedding

[videlec]

comment:15

I don't mind...

---

New commits:

d55dc67

trac #16878: bikeshedding

[videlec]

Changed commit from 175d6b5 to d55dc67

[videlec]

comment:16

ping! (applies cleanly on 6.4.beta2)

[jdemeyer]

Description changed:

--- 
+++ 
@@ -23,3 +23,5 @@
 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.

[jdemeyer]

comment:17

Big reviewer patch coming up...

[jdemeyer]

Description changed:

--- 
+++ 
@@ -24,4 +24,4 @@
 10000 loops, best of 3: 730 ns per loop

-We also modify is_prime_power() to return False for 1. +We also modify is_prime_power() to return False for 1 or non-integral rationals like 1/2.

[jdemeyer]

Description changed:

--- 
+++ 
@@ -24,4 +24,4 @@
 10000 loops, best of 3: 730 ns per loop

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

[jdemeyer]

Changed branch from u/vdelecroix/16878 to u/jdemeyer/ticket/16878

[jdemeyer]

Changed author from Vincent Delecroix to Vincent Delecroix, Jeroen Demeyer

[jdemeyer]

comment:21

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:

7f61321

Reorganize is_prime_power() and next_prime() logic

[jdemeyer]

Changed commit from d55dc67 to 7f61321

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 7f61321 to 5937af5

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

5937af5

Sloane considers 1 to be a prime power

[videlec]

Changed commit from 5937af5 to f0d716f

[videlec]

comment:23

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 pⁿ (following pari convention). I think it is worth it!

Vincent

---

New commits:

f0d716f

trac #16878: certificate flag in is_prime_power

[videlec]

Changed branch from u/jdemeyer/ticket/16878 to u/vdelecroix/16878

[videlec]

comment:24

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

[jdemeyer]

comment:25

Replying to @videlec:

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 pⁿ (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?

[jdemeyer]

comment:26

Replying to @videlec:

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.

[videlec]

comment:27

Replying to @jdemeyer:

Replying to @videlec:

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 pⁿ (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

[jdemeyer]

Changed branch from u/vdelecroix/16878 to u/jdemeyer/ticket/16878

[jdemeyer]

Changed commit from f0d716f to bfd69d3

[jdemeyer]

comment:29

One more reviewer patch...

---

New commits:

bfd69d3

is_prime_power() reviewer fixes

[videlec]

comment:30

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

Vincent

[jdemeyer]

Changed reviewer from Jeroen Demeyer to Jeroen Demeyer, Vincent Delecroix

[vbraun]

Changed branch from u/jdemeyer/ticket/16878 to bfd69d3

[vbraun]

Changed commit from bfd69d3 to none

[vbraun]

comment:33

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

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

[videlec]

comment:34

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

[videlec]

comment:35

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

[videlec]

Changed branch from bfd69d3 to u/vdelecroix/16878