Pack primes numbers by slices of 30

In any given 30 successive integers greater than 5, there are at most
8 prime numbers. Use this to tightly pack the result of the Eratostene
sieve. This lets us store more prime numbers than before in less space.

basis/math/primes/erato/erato-docs.factor

@@ -3,10 +3,8 @@ IN: math.primes.erato
 
 HELP: sieve
 { $values { "n" "the greatest odd number to consider" } { "arr" "a bit array" } }
-{ $description "Return a bit array containing a primality bit for every odd number between 3 and " { $snippet "n" } " (inclusive). " { $snippet ">index" } " can be used to retrieve the index of an odd number to be tested." } ;
+{ $description "Apply Eratostene sieve up to " { $snippet "n" } ". Primality can then be tested using " { $link sieve } "." } ;
 
-HELP: >index
-{ $values { "n" "an odd number" } { "i" "the corresponding index" } }
-{ $description "Retrieve the index corresponding to the odd number on the stack." } ;
-
-{ sieve >index } related-words
+HELP: marked-prime?
+{ $values { "n" "an integer" } { "arr" "a byte array returned by " { $link sieve } } { "?" "a boolean" } }
+{ $description "Check whether a number between 3 and the limit given to " { $link sieve } " has been marked as a prime number."} ;

basis/math/primes/erato/erato-tests.factor

@@ -1,3 +1,10 @@
-USING: bit-arrays math.primes.erato tools.test ;
+USING: byte-arrays math math.bitwise math.primes.erato sequences tools.test ;
 
-[ ?{ t t t f t t f t t f t f f t } ] [ 29 sieve ] unit-test
+[ B{ 255 251 247 126 } ] [ 100 sieve ] unit-test
+[ 1 100 sieve marked-prime? ] [ bounds-error? ] must-fail-with
+[ 120 100 sieve marked-prime? ] [ bounds-error? ] must-fail-with
+[ f ] [ 119 100 sieve marked-prime? ] unit-test
+[ t ] [ 113 100 sieve marked-prime? ] unit-test
+
+! There are 25997 primes below 300000. 1 must be removed and 3 5 7 added.
+[ 25997 ] [ 299999 sieve [ bit-count ] sigma 2 + ] unit-test

basis/math/primes/erato/erato.factor

@@ -1,25 +1,41 @@
 ! Copyright (C) 2009 Samuel Tardieu.
 ! See http://factorcode.org/license.txt for BSD license.
-USING: bit-arrays kernel math math.functions math.ranges sequences ;
+USING: arrays byte-arrays kernel math math.bitwise math.functions math.order
+math.ranges sequences sequences.private ;
 IN: math.primes.erato
 
-: >index ( n -- i )
-    3 - 2 /i ; inline
+<PRIVATE
 
-: index> ( i -- n )
-    2 * 3 + ; inline
+CONSTANT: masks B{ 0 128 0 0 0 0 0 64 0 0 0 32 0 16 0 0 0 8 0 4 0 0 0 2 0 0 0 0 0 1 }
 
-: mark-multiples ( i arr -- )
-    [ index> [ sq >index ] keep ] dip
-    [ length 1 - swap <range> f swap ] keep
-    [ set-nth ] curry with each ;
+: bit-pos ( n -- byte/f mask/f )
+    30 /mod masks nth-unsafe dup zero? [ 2drop f f ] when ;
 
-: maybe-mark-multiples ( i arr -- )
-    2dup nth [ mark-multiples ] [ 2drop ] if ;
+: marked-unsafe? ( n arr -- ? )
+    [ bit-pos ] dip swap [ [ nth-unsafe ] [ bitand zero? not ] bi* ] [ 2drop f ] if* ;
 
-: init-sieve ( n -- arr )
-    >index 1 + <bit-array> dup set-bits ;
+: unmark ( n arr -- )
+    [ bit-pos swap ] dip
+    over [ [ swap unmask ] change-nth-unsafe ] [ 3drop ] if ;
+
+: upper-bound ( arr -- n ) length 30 * 1 - ;
+
+: unmark-multiples ( i arr -- )
+    2dup marked-unsafe? [
+        [ [ dup sq ] [ upper-bound ] bi* rot <range> ] keep
+        [ unmark ] curry each
+    ] [
+        2drop
+    ] if ;
+
+: init-sieve ( n -- arr ) 29 + 30 /i 255 <array> >byte-array ;
+
+PRIVATE>
 
 : sieve ( n -- arr )
-    [ init-sieve ] [ sqrt >index [0,b] ] bi
-    over [ maybe-mark-multiples ] curry each ; foldable
+    init-sieve [ 2 swap upper-bound sqrt [a,b] ] keep
+    [ [ unmark-multiples ] curry each ] keep ;
+
+: marked-prime? ( n arr -- ? )
+    2dup upper-bound 2 swap between? [ bounds-error ] unless
+    over { 2 3 5 } member? [ 2drop t ] [ marked-unsafe? ] if ;

basis/math/primes/primes-tests.factor

@@ -10,6 +10,9 @@ IN: math.primes.tests
 { { 4999963 4999999 5000011 5000077 5000081 } }
 [ 4999962 5000082 primes-between >array ] unit-test
 
+{ { 8999981 8999993 9000011 9000041 } }
+[ 8999980 9000045 primes-between >array ] unit-test
+
 [ 2 ] [ 1 next-prime ] unit-test
 [ 3 ] [ 2 next-prime ] unit-test
 [ 5 ] [ 3 next-prime ] unit-test

basis/math/primes/primes.factor

@@ -1,31 +1,31 @@
 ! Copyright (C) 2007-2009 Samuel Tardieu.
 ! See http://factorcode.org/license.txt for BSD license.
 USING: combinators kernel math math.bitwise math.functions
-math.order math.primes.erato math.primes.miller-rabin
-math.ranges random sequences sets fry ;
+math.order math.primes.erato math.primes.erato.private
+math.primes.miller-rabin math.ranges literals random sequences sets ;
 IN: math.primes
 
 <PRIVATE
 
-: look-in-bitmap ( n -- ? ) >index 4999999 sieve nth ;
+: look-in-bitmap ( n -- ? ) $[ 8999999 sieve ] marked-unsafe? ; inline
 
-: really-prime? ( n -- ? )
-    dup 5000000 < [ look-in-bitmap ] [ miller-rabin ] if ; foldable
+: (prime?) ( n -- ? )
+    dup 8999999 <= [ look-in-bitmap ] [ miller-rabin ] if ;
 
 PRIVATE>
 
 : prime? ( n -- ? )
     {
-        { [ dup 2 < ] [ drop f ] }
+        { [ dup 7 < ] [ { 2 3 5 } member? ] }
         { [ dup even? ] [ 2 = ] }
-        [ really-prime? ]
+        [ (prime?) ]
     } cond ; foldable
 
 : next-prime ( n -- p )
     dup 2 < [
         drop 2
     ] [
-        next-odd [ dup really-prime? ] [ 2 + ] until
+        next-odd [ dup prime? ] [ 2 + ] until
     ] if ; foldable
 
 : primes-between ( low high -- seq )
@@ -65,5 +65,5 @@ ERROR: too-few-primes n numbits ;
 
 : unique-primes ( n numbits -- seq )
     2dup 2^ estimated-primes > [ too-few-primes ] when
-    2dup '[ _ random-prime ] replicate
+    2dup [ random-prime ] curry replicate
     dup all-unique? [ 2nip ] [ drop unique-primes ] if ;