previous_prime, previous_prime_power, next_prime_power

[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

[ncohen]

+1 for "primes"

[git]

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

​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
​175d6b5	trac #16878: review 2
​d55dc67	trac #16878: bikeshedding
​f73bd87	trac #16869: fast prime_powers
​8259a7d	trac #16869: fix prime_powers(0,1)
​6c6cbe2	trac #16880: implement prime_power_range

[git]

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

​e7b9458

trac #16880: prime_power_range

[vdelecroix]

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

Vincent

[git]

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

[vdelecroix]

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

[git]

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

​e3be1ce

trac #16880: prime_powers and prime_power_range

[vdelecroix]

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

Vincent

[jdemeyer]

Needs rebase.

[git]

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

​b365b9e	trac #17852: fix binomial(n,0) for negative n
​faed0c0	trac #17852: remove lcm on rationals (category takes care)
​cbbb06f	trac #17852: make a reference to the symbolic binomial
​1655367	trac #17852: replace powermod_ui by hash in the doc
​d401269	trac #17852: fix output type of arith.binomial
​8fc6b3d	trac #17852: type -> parent
​61b94c2	trac #17852: fix curious failing doctest in multi_polynomial_ideal
​6697e68	trac #17852: be nicer on non-integral domains
​2260680	trac #17852: merge beta-6.6
​25176c5	trac #16880: cleanup for prime and prime powers

[git]

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

​b3cfd3c

trac 16880: fix doctest related to random_primes

[git]

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

​d85a98c	Trac 16880: previous_prime/{next,previous}_prime_power
​db1e869	Trac 16880: make rings.arith uses the new methods

[vdelecroix]

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

[jdemeyer]

1. I don't know if lesser or equal than is proper English. How about less than or equal to?

2. PARI instead of Pari.

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

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

[git]

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

​1e09561

Trac 16880: optimization + Pari -> PARI

[vdelecroix]

Replying to jdemeyer:

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

[jdemeyer]

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.

[git]

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

​35c0211

Trac 16880: optimization

[leif]

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

[vdelecroix]

Replying to leif:

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

Ha right.