steve
/
factor
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

math.primes.erato: faster compressed sieve by two improvements: 

1) fixnum math in unmark-multiples
2) 3 upper sqrt 2 <range> in sieve
db4
John Benediktsson 2015-06-16 08:26:48 -07:00
parent 7e23c12c0f
commit 64db5c5bb4
2 changed files with 27 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
basis/math/primes/erato/erato-tests.factor
View File

@ -1,5 +1,5 @@
USING: kernel byte-arrays sequences tools.test ; USING: fry kernel math math.bitwise math.primes.erato
USING: math math.bitwise math.ranges math.primes.erato ; math.ranges sequences tools.test ;
[ B{ 255 251 247 126 } ] [ 100 sieve ] unit-test [ B{ 255 251 247 126 } ] [ 100 sieve ] unit-test
[ 1 100 sieve marked-prime? ] [ bounds-error? ] must-fail-with [ 1 100 sieve marked-prime? ] [ bounds-error? ] must-fail-with
@ -14,3 +14,11 @@ USING: math math.bitwise math.ranges math.primes.erato ;
! end-point for numbers with all possibilities mod 30. If something ! end-point for numbers with all possibilities mod 30. If something
! were to go wrong, we'd get a bounds-error. ! were to go wrong, we'd get a bounds-error.
[ ] [ 2 100 [a,b] [ dup sieve marked-prime? drop ] each ] unit-test [ ] [ 2 100 [a,b] [ dup sieve marked-prime? drop ] each ] unit-test
{ t } [
{ 2 3 5 7 11 13 } 100 sieve '[ _ marked-prime? ] all?
] unit-test
{ t } [
{ 4 6 8 9 10 12 } 100 sieve '[ _ marked-prime? not ] all?
] unit-test

25
basis/math/primes/erato/erato.factor
View File

@ -1,7 +1,7 @@
! Copyright (C) 2009 Samuel Tardieu. ! Copyright (C) 2009 Samuel Tardieu.
! See http://factorcode.org/license.txt for BSD license. ! See http://factorcode.org/license.txt for BSD license.
USING: kernel kernel.private locals math math.bitwise USING: kernel kernel.private locals math math.bitwise
math.functions math.order math.ranges sequences math.functions math.order math.private math.ranges sequences
sequences.private ; sequences.private ;
IN: math.primes.erato IN: math.primes.erato
@ -15,19 +15,25 @@ CONSTANT: masks
30 /mod masks nth-unsafe 30 /mod masks nth-unsafe
{ maybe{ fixnum } } declare ; inline { maybe{ fixnum } } declare ; inline
: marked-unsafe? ( n sieve -- ? ) :: marked-unsafe? ( n sieve -- ? )
[ bit-pos ] dip swap n bit-pos [
[ [ nth-unsafe ] [ mask zero? not ] bi* ] [ 2drop f ] if* ; inline [ sieve nth-unsafe ] [ mask zero? not ] bi*
] [ drop f ] if* ; inline
: unmark ( n sieve -- ) :: unmark ( n sieve -- )
[ bit-pos swap ] dip pick n bit-pos [
[ [ swap unmask ] change-nth-unsafe ] [ 3drop ] if ; inline swap sieve [ swap unmask ] change-nth-unsafe
] [ drop ] if* ; inline
: upper-bound ( sieve -- n ) length 30 * 1 - ; inline : upper-bound ( sieve -- n ) length 30 * 1 - ; inline
:: unmark-multiples ( i upper sieve -- ) :: unmark-multiples ( i upper sieve -- )
i sieve marked-unsafe? [ i sieve marked-unsafe? [
i sq upper i <range> [ sieve unmark ] each i 2 fixnum*fast :> step
i i fixnum*fast
[ dup upper <= ] [
[ sieve unmark ] [ step fixnum+fast ] bi
] while drop
] when ; inline ] when ; inline
: init-sieve ( n -- sieve ) : init-sieve ( n -- sieve )
@ -38,10 +44,11 @@ PRIVATE>
:: sieve ( n -- sieve ) :: sieve ( n -- sieve )
n integer>fixnum-strict init-sieve :> sieve n integer>fixnum-strict init-sieve :> sieve
sieve upper-bound >fixnum :> upper sieve upper-bound >fixnum :> upper
2 upper sqrt [a,b] 3 upper sqrt 2 <range>
[ upper sieve unmark-multiples ] each [ upper sieve unmark-multiples ] each
sieve ; sieve ;
: marked-prime? ( n sieve -- ? ) : marked-prime? ( n sieve -- ? )
[ integer>fixnum-strict ] dip
2dup upper-bound 2 swap between? [ bounds-error ] unless 2dup upper-bound 2 swap between? [ bounds-error ] unless
over { 2 3 5 } member? [ 2drop t ] [ marked-unsafe? ] if ; over { 2 3 5 } member? [ 2drop t ] [ marked-unsafe? ] if ;