Ticket #3925: trac-3925-primegen.patch

module_list.py

@@ -414 +414 @@
               #depends = [SAGE_ROOT + 'local/include/ratpoints.h'],
               libraries = ["ratpoints", "gmp"]),
     
+    Extension('sage.libs.primegen.primegen',
+              sources = ['sage/libs/primegen/primegen.pyx'],
+              libraries = [ 'primegen' ],
+              language="c",
+              include_dirs=debian_include_dirs + [SAGE_ROOT +'/local/include']
+              ),
+    
     Extension('sage.libs.singular.singular',
               sources = ['sage/libs/singular/singular.pyx'],
               libraries = ['m', 'readline', 'singular', 'givaro', 'gmpxx', 'gmp'],

new file sage/libs/primegen/primegen.pxd

@@ +1 @@
+include "../../ext/cdefs.pxi"
+include "primegen.pxi"
+
+cdef class PrimeGen:
+    cdef primegen pg
+    cdef int done

new file sage/libs/primegen/primegen.pxi

@@ +1 @@
+cdef extern from "primegen/primegen.h":
+    # TODO: Are there any standard types we can use for this?
+    ctypedef unsigned long long uint64
+
+    ctypedef struct primegen:
+        pass
+
+    void primegen_init(primegen* pg)
+    uint64 primegen_next(primegen* pg)
+    uint64 primegen_peek(primegen* pg)
+    uint64 primegen_count(primegen* pg, uint64 to)
+    void primegen_skipto(primegen* pg, uint64 to)

new file sage/libs/primegen/primegen.pyx

@@ +1 @@
+from sage.rings.integer import Integer
+
+# TODO: check what the highest guaranteed valid prime returned by primegen is
+#       (Currently set to (10^15).previous_prime() per the documentation.)
+
+# (10^15).previous_prime()
+max_prime = 999999999999989
+
+cdef class PrimeGen:
+    """
+    A wrapper class around D.J. Bernstein's primegen library
+    that uses the Sieve of Atkin to generate primes efficiently.
+
+    See http://cr.yp.to/primegen.html
+    """
+
+    def __new__(self, start=2):
+        """
+        Construct a new primegen wrapper and initialize it to the given start.
+
+        EXAMPLES::
+
+            sage: from sage.libs.primegen.primegen import PrimeGen
+            sage: pg = PrimeGen()
+            sage: pg.next()
+            2
+
+            sage: pg = PrimeGen(37)
+            sage: pg.next()
+            37
+        """
+        self.reset(start=start)
+
+    def next(self):
+        """
+        Return the next prime. If the end of the valid range for
+        primegen has been reached, it will raise a RuntimeError.
+
+        EXAMPLES::
+
+            sage: from sage.libs.primegen.primegen import PrimeGen
+            sage: pg = PrimeGen(10)
+            sage: pg.next()
+            11
+            sage: pg.next()
+            13
+        """
+        if self.done:
+            raise RuntimeError("End of valid primegen range reached.")
+        p = primegen_next(&self.pg)
+        self.done = (p == max_prime)
+        return Integer(p)
+
+
+    def peek(self):
+        """
+        Return the next prime without advancing its state. If the end of
+        the valid range for primegen has been reached, it will return None.
+
+        EXAMPLES::
+
+            sage: from sage.libs.primegen.primegen import PrimeGen
+            sage: pg = PrimeGen(100)
+            sage: pg.peek()
+            101
+            sage: pg.next()
+            101
+        """
+ 
+        if self.done:
+            return None
+        return Integer(primegen_peek(&self.pg))
+
+    def count(self, uint64 to):
+        """
+        Count primes smaller than to starting at the current location.
+
+        EXAMPLES::
+
+            sage: from sage.libs.primegen.primegen import PrimeGen
+            sage: pg = PrimeGen()
+            sage: pg.count(7)  # 2, 3, 5
+            3
+            sage: pg.count(10) # 7
+            1
+
+            sage: pg = PrimeGen()
+            sage: pg.count(1000000000)
+            50847534
+
+            sage: pg = PrimeGen()
+            sage: pg.count(2000000000) == prime_pi(2000000000)
+            True
+        """
+        if to > max_prime:
+            raise RuntimeError("Counting past end of valid primegen range.")
+        return Integer(primegen_count(&self.pg, to))
+
+    def skipto(self, uint64 to):
+        """
+        Advance primegen state so the next prime returned will be
+        the smallest prime after to. This function has no effect if the
+        current primegen location is already after to.
+
+        EXAMPLES::
+
+            sage: from sage.libs.primegen.primegen import PrimeGen
+            sage: pg = PrimeGen()
+            sage: pg.skipto(13)
+            sage: pg.next()
+            13
+            sage: pg.skipto(10)
+            sage: pg.next()
+            17
+        """
+        if to > max_prime:
+            raise RuntimeError("Skipping past end of valid primegen range.")
+        primegen_skipto(&self.pg, to)
+
+    def reset(self, start=2):
+        """
+        Reinitialize the prime generator to restart at the given value.
+
+        EXAMPLES::
+
+            sage: from sage.libs.primegen.primegen import PrimeGen
+            sage: pg = PrimeGen(1000)
+            sage: pg.reset()
+            sage: pg.next()
+            2
+            sage: pg.reset(100)
+            sage: pg.next()
+            101
+        """
+        primegen_init(&self.pg)
+        self.done = False
+        self.skipto(start)

sage/sets/primes.py

@@ -18 +18 @@
 #*****************************************************************************
 
 from sage.rings.all import Integer, ZZ, infinity
+from sage.libs.primegen.primegen import PrimeGen
 
 from set import Set_generic
 
@@ -40 +41 @@
 
             sage: P = Primes()
         """
-        pass
+        self.pg = PrimeGen()
+
 
     def cardinality(self):
         """
@@ -104 +106 @@
             sage: iter(P).next()
             2
         """
-        p = Integer(2)
+        p = Integer(1)
+        self.pg.reset()
+        while self.pg.peek():
+            p = self.pg.next()
+            yield p
+        p = p.next_prime()
         while True:
             yield p
             p = p.next_prime()

setup.py

@@ -790 +790 @@
                      'sage.libs.ntl',
                      'sage.libs.flint',
                      'sage.libs.pari',
+                     'sage.libs.primegen',
                      'sage.libs.singular',
                      'sage.libs.symmetrica',
                      'sage.libs.cremona',