From 41d804ddbd00b7041b97f98894a6b385a9ea6d3c Mon Sep 17 00:00:00 2001 From: Samuel Tardieu Date: Wed, 24 Jun 2009 13:04:20 +0200 Subject: [PATCH] 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 | 10 ++--- basis/math/primes/erato/erato-tests.factor | 11 +++++- basis/math/primes/erato/erato.factor | 46 +++++++++++++++------- basis/math/primes/primes-tests.factor | 3 ++ basis/math/primes/primes.factor | 18 ++++----- 5 files changed, 56 insertions(+), 32 deletions(-) diff --git a/basis/math/primes/erato/erato-docs.factor b/basis/math/primes/erato/erato-docs.factor index b12ea45052..1e32818fe3 100644 --- a/basis/math/primes/erato/erato-docs.factor +++ b/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."} ; diff --git a/basis/math/primes/erato/erato-tests.factor b/basis/math/primes/erato/erato-tests.factor index 917824c9c1..e78e5210f9 100644 --- a/basis/math/primes/erato/erato-tests.factor +++ b/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 \ No newline at end of file diff --git a/basis/math/primes/erato/erato.factor b/basis/math/primes/erato/erato.factor index 70a9c10ff5..673f9c97cd 100644 --- a/basis/math/primes/erato/erato.factor +++ b/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 + ( 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 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 + 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 ] keep + [ unmark ] curry each + ] [ + 2drop + ] if ; + +: init-sieve ( n -- arr ) 29 + 30 /i 255 >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 diff --git a/basis/math/primes/primes-tests.factor b/basis/math/primes/primes-tests.factor index 6580f0780e..a950395bf4 100644 --- a/basis/math/primes/primes-tests.factor +++ b/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 diff --git a/basis/math/primes/primes.factor b/basis/math/primes/primes.factor index e3985fc600..ea8c60508d 100644 --- a/basis/math/primes/primes.factor +++ b/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 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 ;