Opened 6 years ago
Closed 6 years ago
#16878 closed enhancement (fixed)
faster is_prime
Reported by:  vdelecroix  Owned by:  

Priority:  major  Milestone:  sage6.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 )
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 nonintegral rationals like 1/2
.
See also this sagedevel discussion.
Change History (70)
comment:1 Changed 6 years ago by
 Branch set to u/vdelecroix/16878
 Commit set to dd44a7e8d05a2b3978c679c363a2729bf0cf5e93
 Status changed from new to needs_review
comment:2 Changed 6 years ago by
 Description modified (diff)
comment:3 followup: ↓ 6 Changed 6 years ago by
 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
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
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
Thanks for the remark. One question: are we sure that 2^{64} 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
comment:7 Changed 6 years ago by
 Commit changed from dd44a7e8d05a2b3978c679c363a2729bf0cf5e93 to d7dc38918903a97af742a90050ce87c58f6f222a
Branch pushed to git repo; I updated commit sha1. This was a forced push. New commits:
37fc8e8  Upgrade to PARI2.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

comment:8 Changed 6 years ago by
 Dependencies set to #15767
 Status changed from needs_work to needs_review
comment:9 Changed 6 years ago by
 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
Also, you wrote "set set" instead of "set".
comment:11 Changed 6 years ago by
 Commit changed from d7dc38918903a97af742a90050ce87c58f6f222a to 175d6b5fd339989df12ec9c1f26d3b0d7e6dd82c
Branch pushed to git repo; I updated commit sha1. New commits:
175d6b5  trac #16878: review 2

comment:12 Changed 6 years ago by
 Status changed from needs_work to needs_review
comment:13 Changed 6 years ago by
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
 Commit changed from 175d6b5fd339989df12ec9c1f26d3b0d7e6dd82c to d55dc6737936077eeee903687207f7af54588439
Branch pushed to git repo; I updated commit sha1. New commits:
d55dc67  trac #16878: bikeshedding

comment:15 Changed 6 years ago by
 Commit changed from 175d6b5fd339989df12ec9c1f26d3b0d7e6dd82c to d55dc6737936077eeee903687207f7af54588439
comment:16 Changed 6 years ago by
ping! (applies cleanly on 6.4.beta2)
comment:18 Changed 6 years ago by
 Description modified (diff)
comment:19 Changed 6 years ago by
 Description modified (diff)
comment:20 Changed 6 years ago by
 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
 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:
7f61321  Reorganize is_prime_power() and next_prime() logic

comment:22 Changed 6 years ago by
 Commit changed from 7f6132188b7ddafd97630aaf721875f9d60a7106 to 5937af54afa0562007d710c495ecd0f1160f7498
Branch pushed to git repo; I updated commit sha1. New commits:
5937af5  Sloane considers 1 to be a prime power

comment:23 followup: ↓ 25 Changed 6 years ago by
 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 p^{n} (following pari convention). I think it is worth it!
Vincent
New commits:
f0d716f  trac #16878: certificate flag in is_prime_power

comment:24 followup: ↓ 26 Changed 6 years ago by
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 ; followup: ↓ 27 Changed 6 years ago by
Replying to vdelecroix:
Great! I especially like that now
isprimepower
is accessible.In my last commit, I added the keyword
certificate
tois_prime_power
in order to get back the pair (n,p) such that the integer is p^{n} (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
Replying to vdelecroix:
Note the treatment of 1 is different in
is_prime_powers
andprime_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
Replying to jdemeyer:
Replying to vdelecroix:
Great! I especially like that now
isprimepower
is accessible.In my last commit, I added the keyword
certificate
tois_prime_power
in order to get back the pair (n,p) such that the integer is p^{n} (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()
(raisingArithmeticError
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
 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
 Commit changed from f0d716fd34029f65d0be0aa45d5cb1783f8da2f5 to bfd69d3b188a8ec36ab373b4e958b54c74b64d2b
comment:30 Changed 6 years ago by
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
 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
 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
 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 followup: ↓ 35 Changed 6 years ago by
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
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
 Branch changed from bfd69d3b188a8ec36ab373b4e958b54c74b64d2b to u/vdelecroix/16878
 Commit set to e6f88c5f9adef669ede74e3c161e19fe37203276
 Status changed from new to needs_review
New commits:
6e9f686  trac #16878: faster is_prime

d7dc389  trac #16878: review

175d6b5  trac #16878: review 2

d55dc67  trac #16878: bikeshedding

7f61321  Reorganize is_prime_power() and next_prime() logic

5937af5  Sloane considers 1 to be a prime power

f0d716f  trac #16878: certificate flag in is_prime_power

bfd69d3  is_prime_power() reviewer fixes

e6f88c5  trac #16878: test is_pseudoprime_small_power + clean

comment:37 followup: ↓ 38 Changed 6 years ago by
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
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 mostlog(n)/log(p)
, wheren
is the number to check andp
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.ubordeaux.fr/archives/paridev1409/msg00011.html...
On the other hand, we need a reasonable fix in order to get this ticket in.
Vincent
comment:39 followup: ↓ 40 Changed 6 years ago by
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
comment:41 Changed 6 years ago by
On the other hand, perhaps better add the patch here...
comment:42 Changed 6 years ago by
...or a new ticket.
comment:43 Changed 6 years ago by
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
 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
 Dependencies changed from #15767 to #15767, #16997
comment:46 Changed 6 years ago by
comment:47 Changed 6 years ago by
 Dependencies changed from #15767, #16997 to #16997, #17852
comment:48 Changed 6 years ago by
 Branch changed from u/vdelecroix/16878 to u/jdemeyer/16878
comment:49 Changed 6 years ago by
 Commit changed from e6f88c5f9adef669ede74e3c161e19fe37203276 to 9751ee8ccf4887970ed5cea1e6f52447b633e2c3
 Dependencies changed from #16997, #17852 to #16997
 Status changed from needs_work to needs_review
comment:50 Changed 6 years ago by
 Milestone changed from sage6.4 to sage6.6
 Status changed from needs_review to needs_work
comment:51 Changed 6 years ago by
(I guess you did it because of conflicts with #17852)
comment:52 followup: ↓ 53 Changed 6 years ago by
comment:53 in reply to: ↑ 52 Changed 6 years ago by
comment:54 Changed 6 years ago by
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
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
 Commit changed from 9751ee8ccf4887970ed5cea1e6f52447b633e2c3 to 96f1f5037d18eb3b2e4c8e0e2a9c7a5cdac0d643
Branch pushed to git repo; I updated commit sha1. New commits:
96f1f50  Trac #16878: further fixes to is_prime() and friends

comment:57 Changed 6 years ago by
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
 Status changed from needs_work to needs_review
comment:59 Changed 6 years ago by
 Commit changed from 96f1f5037d18eb3b2e4c8e0e2a9c7a5cdac0d643 to 1f8abd999dcf4935f92b2d83aad01684397e62d1
Branch pushed to git repo; I updated commit sha1. New commits:
1f8abd9  Doctest fix

comment:60 Changed 6 years ago by
Now all doctests pass.
comment:61 followup: ↓ 62 Changed 6 years ago by
 Status changed from needs_review to needs_info
 Now
arith.is_prime_power
andarith.is_pseudoprime_power
do not handle negative powers anymore (i.e. rational numbersp^n
). This is different from previous versions.
 In
Integer.is_prime_power
why did you use along
forn
andp
instead of anunsigned long
. Is pariuisprimepower
not able to handle up toULONG_MAX
?
comment:62 in reply to: ↑ 61 ; followup: ↓ 63 Changed 6 years ago by
Replying to vdelecroix:
 Now
arith.is_prime_power
andarith.is_pseudoprime_power
do not handle negative powers anymore (i.e. rational numbersp^n
). This is different from previous versions.
Yes, this was discussed on sagedevel
and the former behaviour was considered to be a bug => no deprecation needed
 In
Integer.is_prime_power
why did you use along
forn
andp
instead of anunsigned long
. Is pariuisprimepower
not able to handle up toULONG_MAX
?
I assume that pari
can handle it, but smallInteger()
for the result takes a long
.
comment:63 in reply to: ↑ 62 ; followup: ↓ 64 Changed 6 years ago by
Replying to jdemeyer:
Replying to vdelecroix:
 Now
arith.is_prime_power
andarith.is_pseudoprime_power
do not handle negative powers anymore (i.e. rational numbersp^n
). This is different from previous versions.Yes, this was discussed on
sagedevel
and the former behaviour was considered to be a bug => no deprecation needed
True, I forgot about here and it is this thread.
 In
Integer.is_prime_power
why did you use along
forn
andp
instead of anunsigned long
. Is pariuisprimepower
not able to handle up toULONG_MAX
?I assume that
pari
can handle it, butsmallInteger()
for the result takes along
.
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
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
 Description modified (diff)
comment:66 Changed 6 years ago by
 Status changed from needs_info to needs_review
comment:67 followup: ↓ 68 Changed 6 years ago by
 Status changed from needs_review to needs_work
The buildbot is complaining about a digit
sage t long warnlong 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
 Status changed from needs_work to needs_review
comment:69 Changed 6 years ago by
 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
 Branch changed from u/jdemeyer/16878 to 1f8abd999dcf4935f92b2d83aad01684397e62d1
 Resolution set to fixed
 Status changed from positive_review to closed
New commits:
trac #16878: faster is_prime