Skip to content

Commit

Permalink
Browse files Browse the repository at this point in the history
  • Loading branch information
mkoeppe committed Sep 19, 2021
2 parents 2273487 + 41df930 commit 90b6f3d
Show file tree
Hide file tree
Showing 3 changed files with 64 additions and 33 deletions.
2 changes: 2 additions & 0 deletions src/sage/libs/pari/convert_sage.pxd
Original file line number Diff line number Diff line change
Expand Up @@ -11,3 +11,5 @@ cpdef Gen new_gen_from_rational(Rational self)

cpdef pari_is_prime(Integer p)
cpdef pari_is_prime_power(Integer q, bint get_data)
cpdef unsigned long pari_maxprime()
cpdef list pari_prime_range(long c_start, long c_stop, bint py_ints=*)
49 changes: 48 additions & 1 deletion src/sage/libs/pari/convert_sage.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,8 @@ from cypari2.gen cimport objtogen
from cypari2.stack cimport new_gen
from .convert_gmp cimport INT_to_mpz, new_gen_from_mpz_t, new_gen_from_mpq_t, INTFRAC_to_mpq

from sage.libs.gmp.mpz cimport mpz_fits_slong_p, mpz_sgn, mpz_get_ui, mpz_set, mpz_set_si
from sage.ext.stdsage cimport PY_NEW
from sage.libs.gmp.mpz cimport mpz_fits_slong_p, mpz_sgn, mpz_get_ui, mpz_set, mpz_set_si, mpz_set_ui
from sage.libs.gmp.mpq cimport mpq_denref, mpq_numref
from sage.rings.integer cimport smallInteger
from sage.rings.all import RealField, ComplexField, QuadraticField
Expand Down Expand Up @@ -527,3 +528,49 @@ cpdef pari_is_prime_power(Integer q, bint get_data):
return (smallInteger(p), smallInteger(n)) if get_data else True
else:
return (q, smallInteger(0)) if get_data else False


cpdef unsigned long pari_maxprime():
"""
Return to which limit PARI has computed the primes.
EXAMPLES::
sage: from sage.libs.pari.convert_sage import pari_maxprime
sage: a = pari_maxprime()
sage: res = prime_range(2, 2*a)
sage: b = pari_maxprime()
sage: b >= 2*a
True
"""
return maxprime()


cpdef list pari_prime_range(long c_start, long c_stop, bint py_ints=False):
"""
Return a list of all primes between ``start`` and ``stop - 1``, inclusive.
.. SEEALSO::
:func:`sage.rings.fast_arith.prime_range`
TESTS::
sage: from sage.libs.pari.convert_sage import pari_prime_range
sage: pari_prime_range(2, 19)
[2, 3, 5, 7, 11, 13, 17]
"""
cdef long p = 0
cdef byteptr pari_prime_ptr = diffptr
res = []
while p < c_start:
NEXT_PRIME_VIADIFF(p, pari_prime_ptr)
while p < c_stop:
if py_ints:
res.append(p)
else:
z = <Integer>PY_NEW(Integer)
mpz_set_ui(z.value, p)
res.append(z)
NEXT_PRIME_VIADIFF(p, pari_prime_ptr)
return res
46 changes: 14 additions & 32 deletions src/sage/rings/fast_arith.pyx
Original file line number Diff line number Diff line change
Expand Up @@ -36,13 +36,7 @@ Basic arithmetic with C integers
# The int definitions

from libc.math cimport sqrt
from sage.libs.gmp.mpz cimport mpz_set_ui

from sage.ext.stdsage cimport PY_NEW

from cypari2.paridecl cimport *
from cypari2.gen cimport Gen as pari_gen
from sage.libs.pari.all import pari
from sage.rings.integer cimport Integer

cpdef prime_range(start, stop=None, algorithm=None, bint py_ints=False):
Expand Down Expand Up @@ -152,8 +146,6 @@ cpdef prime_range(start, stop=None, algorithm=None, bint py_ints=False):
- Robert Bradshaw (speedup using Pari prime table, py_ints option)
"""
cdef Integer z
cdef long c_start, c_stop, p
cdef byteptr pari_prime_ptr
# input to pari.init_primes cannot be greater than 436273290 (hardcoded bound)
DEF init_primes_max = 436273290
DEF small_prime_max = 436273009 # a prime < init_primes_max (preferably the largest)
Expand Down Expand Up @@ -187,42 +179,32 @@ cpdef prime_range(start, stop=None, algorithm=None, bint py_ints=False):
algorithm = "pari_isprime"

if algorithm == "pari_primes":
from sage.libs.pari.convert_sage import pari_maxprime, pari_prime_range
from sage.libs.pari import pari

if max(start, stop or 0) > small_prime_max:
raise ValueError('algorithm "pari_primes" is limited to primes larger than'
+ ' {}'.format(small_prime_max - 1))

if stop is None:
# In this case, "start" is really stop
c_start = 1
c_stop = start
stop = start
start = 1
else:
c_start = start
c_stop = stop
if c_start < 1:
c_start = 1
if c_stop <= c_start:
start = start
stop = stop
if start < 1:
start = 1
if stop <= start:
return []

if maxprime() < c_stop:
if pari_maxprime() < stop:
# Adding prime_gap_bound should be sufficient to guarantee an
# additional prime, given that c_stop <= small_prime_max.
pari.init_primes(min(c_stop + prime_gap_bound, init_primes_max))
assert maxprime() >= c_stop

pari_prime_ptr = diffptr
p = 0
res = []
while p < c_start:
NEXT_PRIME_VIADIFF(p, pari_prime_ptr)
while p < c_stop:
if py_ints:
res.append(p)
else:
z = <Integer>PY_NEW(Integer)
mpz_set_ui(z.value, p)
res.append(z)
NEXT_PRIME_VIADIFF(p, pari_prime_ptr)
pari.init_primes(min(stop + prime_gap_bound, init_primes_max))
assert pari_maxprime() >= stop

res = pari_prime_range(start, stop, py_ints)

elif (algorithm == "pari_isprime") or (algorithm == "pari_primes"):
from sage.arith.all import primes
Expand Down

0 comments on commit 90b6f3d

Please sign in to comment.