Ticket #3925: trac-3925-primegen.patch

File trac-3925-primegen.patch, 7.1 KB (added by wjp, 12 years ago)
  • module_list.py

    # HG changeset patch
    # User Willem Jan Palenstijn <wpalenst@math.leidenuniv.nl>
    # Date 1247550421 -7200
    # Node ID dd83478425d73d487e9ea3e28163343828eb073e
    # Parent  ca1f31d6f6bf7b1f8da0cb8973c838ab19921c29
    Add wrapper around D.J. Bernstein's primegen library.
    Based on code by Alex Mendes da Costa.
    
    diff --git a/module_list.py b/module_list.py
    a b  
    414414              #depends = [SAGE_ROOT + 'local/include/ratpoints.h'],
    415415              libraries = ["ratpoints", "gmp"]),
    416416   
     417    Extension('sage.libs.primegen.primegen',
     418              sources = ['sage/libs/primegen/primegen.pyx'],
     419              libraries = [ 'primegen' ],
     420              language="c",
     421              include_dirs=debian_include_dirs + [SAGE_ROOT +'/local/include']
     422              ),
     423   
    417424    Extension('sage.libs.singular.singular',
    418425              sources = ['sage/libs/singular/singular.pyx'],
    419426              libraries = ['m', 'readline', 'singular', 'givaro', 'gmpxx', 'gmp'],
  • new file sage/libs/primegen/primegen.pxd

    diff --git a/sage/libs/primegen/__init__.py b/sage/libs/primegen/__init__.py
    new file mode 100644
    diff --git a/sage/libs/primegen/primegen.pxd b/sage/libs/primegen/primegen.pxd
    new file mode 100644
    - +  
     1include "../../ext/cdefs.pxi"
     2include "primegen.pxi"
     3
     4cdef class PrimeGen:
     5    cdef primegen pg
     6    cdef int done
  • new file sage/libs/primegen/primegen.pxi

    diff --git a/sage/libs/primegen/primegen.pxi b/sage/libs/primegen/primegen.pxi
    new file mode 100644
    - +  
     1cdef extern from "primegen/primegen.h":
     2    # TODO: Are there any standard types we can use for this?
     3    ctypedef unsigned long long uint64
     4
     5    ctypedef struct primegen:
     6        pass
     7
     8    void primegen_init(primegen* pg)
     9    uint64 primegen_next(primegen* pg)
     10    uint64 primegen_peek(primegen* pg)
     11    uint64 primegen_count(primegen* pg, uint64 to)
     12    void primegen_skipto(primegen* pg, uint64 to)
  • new file sage/libs/primegen/primegen.pyx

    diff --git a/sage/libs/primegen/primegen.pyx b/sage/libs/primegen/primegen.pyx
    new file mode 100644
    - +  
     1from sage.rings.integer import Integer
     2
     3# TODO: check what the highest guaranteed valid prime returned by primegen is
     4#       (Currently set to (10^15).previous_prime() per the documentation.)
     5
     6# (10^15).previous_prime()
     7max_prime = 999999999999989
     8
     9cdef class PrimeGen:
     10    """
     11    A wrapper class around D.J. Bernstein's primegen library
     12    that uses the Sieve of Atkin to generate primes efficiently.
     13
     14    See http://cr.yp.to/primegen.html
     15    """
     16
     17    def __new__(self, start=2):
     18        """
     19        Construct a new primegen wrapper and initialize it to the given start.
     20
     21        EXAMPLES::
     22
     23            sage: from sage.libs.primegen.primegen import PrimeGen
     24            sage: pg = PrimeGen()
     25            sage: pg.next()
     26            2
     27
     28            sage: pg = PrimeGen(37)
     29            sage: pg.next()
     30            37
     31        """
     32        self.reset(start=start)
     33
     34    def next(self):
     35        """
     36        Return the next prime. If the end of the valid range for
     37        primegen has been reached, it will raise a RuntimeError.
     38
     39        EXAMPLES::
     40
     41            sage: from sage.libs.primegen.primegen import PrimeGen
     42            sage: pg = PrimeGen(10)
     43            sage: pg.next()
     44            11
     45            sage: pg.next()
     46            13
     47        """
     48        if self.done:
     49            raise RuntimeError("End of valid primegen range reached.")
     50        p = primegen_next(&self.pg)
     51        self.done = (p == max_prime)
     52        return Integer(p)
     53
     54
     55    def peek(self):
     56        """
     57        Return the next prime without advancing its state. If the end of
     58        the valid range for primegen has been reached, it will return None.
     59
     60        EXAMPLES::
     61
     62            sage: from sage.libs.primegen.primegen import PrimeGen
     63            sage: pg = PrimeGen(100)
     64            sage: pg.peek()
     65            101
     66            sage: pg.next()
     67            101
     68        """
     69 
     70        if self.done:
     71            return None
     72        return Integer(primegen_peek(&self.pg))
     73
     74    def count(self, uint64 to):
     75        """
     76        Count primes smaller than to starting at the current location.
     77
     78        EXAMPLES::
     79
     80            sage: from sage.libs.primegen.primegen import PrimeGen
     81            sage: pg = PrimeGen()
     82            sage: pg.count(7)  # 2, 3, 5
     83            3
     84            sage: pg.count(10) # 7
     85            1
     86
     87            sage: pg = PrimeGen()
     88            sage: pg.count(1000000000)
     89            50847534
     90
     91            sage: pg = PrimeGen()
     92            sage: pg.count(2000000000) == prime_pi(2000000000)
     93            True
     94        """
     95        if to > max_prime:
     96            raise RuntimeError("Counting past end of valid primegen range.")
     97        return Integer(primegen_count(&self.pg, to))
     98
     99    def skipto(self, uint64 to):
     100        """
     101        Advance primegen state so the next prime returned will be
     102        the smallest prime after to. This function has no effect if the
     103        current primegen location is already after to.
     104
     105        EXAMPLES::
     106
     107            sage: from sage.libs.primegen.primegen import PrimeGen
     108            sage: pg = PrimeGen()
     109            sage: pg.skipto(13)
     110            sage: pg.next()
     111            13
     112            sage: pg.skipto(10)
     113            sage: pg.next()
     114            17
     115        """
     116        if to > max_prime:
     117            raise RuntimeError("Skipping past end of valid primegen range.")
     118        primegen_skipto(&self.pg, to)
     119
     120    def reset(self, start=2):
     121        """
     122        Reinitialize the prime generator to restart at the given value.
     123
     124        EXAMPLES::
     125
     126            sage: from sage.libs.primegen.primegen import PrimeGen
     127            sage: pg = PrimeGen(1000)
     128            sage: pg.reset()
     129            sage: pg.next()
     130            2
     131            sage: pg.reset(100)
     132            sage: pg.next()
     133            101
     134        """
     135        primegen_init(&self.pg)
     136        self.done = False
     137        self.skipto(start)
  • sage/sets/primes.py

    diff --git a/sage/sets/primes.py b/sage/sets/primes.py
    a b  
    1818#*****************************************************************************
    1919
    2020from sage.rings.all import Integer, ZZ, infinity
     21from sage.libs.primegen.primegen import PrimeGen
    2122
    2223from set import Set_generic
    2324
     
    4041
    4142            sage: P = Primes()
    4243        """
    43         pass
     44        self.pg = PrimeGen()
     45
    4446
    4547    def cardinality(self):
    4648        """
     
    104106            sage: iter(P).next()
    105107            2
    106108        """
    107         p = Integer(2)
     109        p = Integer(1)
     110        self.pg.reset()
     111        while self.pg.peek():
     112            p = self.pg.next()
     113            yield p
     114        p = p.next_prime()
    108115        while True:
    109116            yield p
    110117            p = p.next_prime()
  • setup.py

    diff --git a/setup.py b/setup.py
    a b  
    790790                     'sage.libs.ntl',
    791791                     'sage.libs.flint',
    792792                     'sage.libs.pari',
     793                     'sage.libs.primegen',
    793794                     'sage.libs.singular',
    794795                     'sage.libs.symmetrica',
    795796                     'sage.libs.cremona',