Opened 6 years ago

Closed 6 years ago

Last modified 6 years ago

#16880 closed enhancement (fixed)

previous_prime, previous_prime_power, next_prime_power

Reported by: vdelecroix Owned by:
Priority: major Milestone: sage-6.7
Component: number theory Keywords:
Cc: ncohen, jdemeyer Merged in:
Authors: Vincent Delecroix Reviewers: Jeroen Demeyer
Report Upstream: N/A Work issues:
Branch: 35c0211 (Commits) Commit:
Dependencies: Stopgaps:

Description (last modified by vdelecroix)

In this ticket we add three methods to integers: previous_prime (calling pari precprime), next_prime_power and previous_prime_power. We also make changes in sage.rings.arith to call these (faster) methods.

Timings on the current version (sage-6.7.beta2)

sage: timeit("for i in range(200): next_prime_power(i)", number=5000, repeat=10)
5000 loops, best of 10: 311 µs per loop
sage: a = 33622489
sage: timeit("next_prime_power(a)", number=5000, repeat=10)
5000 loops, best of 10: 53.6 µs per loop

Timings with the branch applied

sage: timeit("for i in range(200): next_prime_power(i)", number=5000, repeat=10)
5000 loops, best of 10: 184 µs per loop
sage: a = 33622489
sage: timeit("next_prime_power(a)", number=5000, repeat=10)
5000 loops, best of 10: 16.9 µs per loop

Follow up: #18298 and #18299

Change History (30)

comment:1 Changed 6 years ago by ncohen

+1 for "primes"

comment:2 Changed 6 years ago by vdelecroix

  • Branch set to u/vdelecroix/16880
  • Description modified (diff)
  • Status changed from new to needs_review

comment:3 Changed 6 years ago by git

  • Commit set to 6c6cbe25b41dece8513ad3527c68415517328710

Branch pushed to git repo; I updated commit sha1. Last 10 new commits:

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
175d6b5trac #16878: review 2
d55dc67trac #16878: bikeshedding
f73bd87trac #16869: fast prime_powers
8259a7dtrac #16869: fix prime_powers(0,1)
6c6cbe2trac #16880: implement prime_power_range

comment:4 Changed 6 years ago by git

  • Commit changed from 6c6cbe25b41dece8513ad3527c68415517328710 to e7b94588ee7ecb82c00afa48f905622f5c374cdd

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

e7b9458trac #16880: prime_power_range

comment:5 Changed 6 years ago by vdelecroix

  • Dependencies #16869, #16878 deleted
  • Description modified (diff)
  • Summary changed from even faster prime_powers to prime_power_range

Rebased on 6.4.beta2, no dependency anymore...

Vincent

comment:6 Changed 6 years ago by git

  • Commit changed from e7b94588ee7ecb82c00afa48f905622f5c374cdd to a3195f7bb4e16cc4eb3f3ae4ef58a4b5d0341ad2

Branch pushed to git repo; I updated commit sha1. This was a forced push. 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
9e04c7atrac #16880: merge #16878 (faster is_prime)
a3195f7trac #16880: prime_power_range

comment:7 Changed 6 years ago by vdelecroix

  • Dependencies set to #16878
  • Description modified (diff)

And with timings

using prime_powers (EDIT: only for commit at a3195f7 and not at e3be1ce)

sage: timeit("prime_powers(1000)", number=4000)
4000 loops, best of 3: 86.7 µs per loop
sage: timeit("prime_powers(5000)", number=2000)
2000 loops, best of 3: 314 µs per loop

using prime_power_range instead

sage: timeit("prime_power_range(1000)", number=4000)
4000 loops, best of 3: 17.6 µs per loop
sage:  timeit("prime_power_range(5000)", number=2000)
2000 loops, best of 3: 72.3 µs per loop

using prime_power_range with py_ints=True

sage: timeit("prime_power_range(1000,py_ints=True)", number=4000)
4000 loops, best of 3: 14.6 µs per loop
sage:  timeit("prime_power_range(5000,py_ints=True)", number=2000)
2000 loops, best of 3: 33.9 µs per loop
Last edited 6 years ago by vdelecroix (previous) (diff)

comment:8 Changed 6 years ago by vdelecroix

  • Description modified (diff)
  • Summary changed from prime_power_range to prime_powers and prime_power_range

comment:9 Changed 6 years ago by git

  • Commit changed from a3195f7bb4e16cc4eb3f3ae4ef58a4b5d0341ad2 to e3be1ce5b13d007baf6982e1b94a2299c68744c1

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

e3be1cetrac #16880: prime_powers and prime_power_range

comment:10 Changed 6 years ago by vdelecroix

It looks more natural to me to follow what primes and prime_range do... this is what I did in e3be1ce.

Vincent

comment:11 Changed 6 years ago by jdemeyer

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

Needs rebase.

comment:12 Changed 6 years ago by jdemeyer

  • Dependencies changed from #16878 to #16878, #17852

comment:13 Changed 6 years ago by git

  • Commit changed from e3be1ce5b13d007baf6982e1b94a2299c68744c1 to 25176c520e538f0a98334c7ffdfe69c14207b32d

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

b365b9etrac #17852: fix binomial(n,0) for negative n
faed0c0trac #17852: remove lcm on rationals (category takes care)
cbbb06ftrac #17852: make a reference to the symbolic binomial
1655367trac #17852: replace powermod_ui by hash in the doc
d401269trac #17852: fix output type of arith.binomial
8fc6b3dtrac #17852: type -> parent
61b94c2trac #17852: fix curious failing doctest in multi_polynomial_ideal
6697e68trac #17852: be nicer on non-integral domains
2260680trac #17852: merge beta-6.6
25176c5trac #16880: cleanup for prime and prime powers

comment:14 Changed 6 years ago by vdelecroix

  • Description modified (diff)
  • Status changed from needs_work to needs_review
  • Summary changed from prime_powers and prime_power_range to further cleanup in arith.py/integer.pyx

comment:15 Changed 6 years ago by git

  • Commit changed from 25176c520e538f0a98334c7ffdfe69c14207b32d to b3cfd3cce04b352ee86b7b84c3880f0c5d5829e9

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

b3cfd3ctrac 16880: fix doctest related to random_primes

comment:16 Changed 6 years ago by vdelecroix

  • Dependencies #16878, #17852 deleted
  • Description modified (diff)
  • Milestone changed from sage-6.6 to sage-6.7
  • Summary changed from further cleanup in arith.py/integer.pyx to previous_prime, previous_prime_power, next_prime_power

comment:17 Changed 6 years ago by git

  • Commit changed from b3cfd3cce04b352ee86b7b84c3880f0c5d5829e9 to db1e869640dbb20d39e565970aaa8c5d58fa8953

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

d85a98cTrac 16880: previous_prime/{next,previous}_prime_power
db1e869Trac 16880: make rings.arith uses the new methods

comment:18 Changed 6 years ago by vdelecroix

  • Description modified (diff)

I split the other parts to different tickets: #18298 and #18299.

comment:19 Changed 6 years ago by vdelecroix

  • Description modified (diff)

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

  1. I don't know if lesser or equal than is proper English. How about less than or equal to?
  1. PARI instead of Pari.
  1. In next_prime_power and previous_prime_power, a small optimization is possible: you can replace
    mpz_set(n.value, self.value)
    mpz_add_ui(n.value, n.value, 1)
    

by

mpz_add_ui(n.value, self.value, 1)

and in next_prime_power, replace

if mpz_cmp_ui(n.value, 2) < 0:
    mpz_set_ui(n.value, 2)
    return n

by

if mpz_cmp_ui(self.value, 2) < 0:
    return smallInteger(2)
  1. Another algorithmic comment for next_prime_power (an analogous comment holds for previous_prime_power): instead of computing the next power of 2 and doing
    while mpz_cmp(m.value, n.value) == 1:
    

I would instead check if mpz_sizeinbase(n.value, 2) increases during the loop. In that case, we know that we skipped a power of 2.

comment:21 Changed 6 years ago by jdemeyer

  • Status changed from needs_review to needs_work

comment:22 Changed 6 years ago by git

  • Commit changed from db1e869640dbb20d39e565970aaa8c5d58fa8953 to 1e0956157ef8f49bc19eb145dc60589f8e59587a

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

1e09561Trac 16880: optimization + Pari -> PARI

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

  • Status changed from needs_work to needs_review

Replying to jdemeyer:

  1. Another algorithmic comment for next_prime_power (an analogous comment holds for previous_prime_power): instead of computing the next power of 2 and doing
    while mpz_cmp(m.value, n.value) == 1:
    

I would instead check if mpz_sizeinbase(n.value, 2) increases during the loop. In that case, we know that we skipped a power of 2.

I used mpz_tstbit instead. Should be the fastest.

Potential micro optimization that would need to use mpn_* functions

  • the call to mpz_sizeinbase can be replaced by mpn_count_leading_zeros and some arithmetic
  • we can avoid a useless bit shift in mpz_tstbit by only considering the & with the most significant limb. Moreover, the allocation of the limbs will not change inside the loop. So we can already use the good pointer to the concerned limb.

I think this is too much for a naive algorithm...

comment:24 Changed 6 years ago by jdemeyer

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

skept -> skipped

This should be moved after the comparison with 2 (and self should be compared with 2):

cdef Integer n = PY_NEW(Integer)

Another simplification that I missed: you can skip the block

    if mpz_even_p(n.value):
    if n.is_prime_power(proof=proof):
    return n
    mpz_add_ui(n.value, n.value, 1)

if you instead just add 2 if self was odd:

mpz_add_ui(n.value, self.value, 1 if mpz_even_p(self.value) else 2)

Then the bit_index should be based on the value of self, such that the check works correctly if self was 2^n - 1.

comment:25 Changed 6 years ago by git

  • Commit changed from 1e0956157ef8f49bc19eb145dc60589f8e59587a to 35c0211766537e6769cd060ed528ad9d77e60068

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

35c0211Trac 16880: optimization

comment:26 Changed 6 years ago by vdelecroix

  • Status changed from needs_work to needs_review

comment:27 Changed 6 years ago by jdemeyer

  • Status changed from needs_review to positive_review

comment:28 Changed 6 years ago by vbraun

  • Branch changed from u/vdelecroix/16880 to 35c0211766537e6769cd060ed528ad9d77e60068
  • Resolution set to fixed
  • Status changed from positive_review to closed

comment:29 follow-up: Changed 6 years ago by leif

  • Commit 35c0211766537e6769cd060ed528ad9d77e60068 deleted

Some copy/paste accidents (could probably be corrected on the follow-ups):

+    def previous_prime(self, proof=None):
+        r"""
+        Returns the previous prime before self.
+
+        This method calls the PARI ``precprime`` function.
+
+        INPUT:
+
+        - ``proof`` - if ``True`` ensure that the returned value is the next
+          prime power and if set to ``False`` uses probabilistic methods
+          (i.e. the result is not guaranteed). By default it uses global
+          configuration variables to determine which alternative to use (see
+          :mod:`proof.arithmetic` or :mod:`sage.structure.proof`).
+    def previous_prime_power(self, proof=None):
+        r"""
+        Return the previous prime power before self.
+
+        INPUT:
+
+        - ``proof`` - if ``True`` ensure that the returned value is the next
+          prime power and if set to ``False`` uses probabilistic methods
+          (i.e. the result is not guaranteed). By default it uses global
+          configuration variables to determine which alternative to use (see
+          :mod:`proof.arithmetic` or :mod:`sage.structure.proof`).

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

Replying to leif:

Some copy/paste accidents (could probably be corrected on the follow-ups):

Ha right.

Note: See TracTickets for help on using tickets.