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."} ;
-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
\ No newline at end of file
! 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 ;
\ No newline at end of file
{ { 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
! 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 )
: 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 ;