factor/basis/math/primes/erato/erato.factor

! Copyright (C) 2009 Samuel Tardieu.
! See http://factorcode.org/license.txt for BSD license.
USING: arrays byte-arrays kernel math math.bitwise math.functions math.order
math.ranges sequences sequences.private ;
IN: math.primes.erato

<PRIVATE

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 }

: bit-pos ( n -- byte/f mask/f )
    30 /mod masks nth-unsafe dup zero? [ 2drop f f ] when ;

: marked-unsafe? ( n arr -- ? )
    [ bit-pos ] dip swap [ [ nth-unsafe ] [ bitand zero? not ] bi* ] [ 2drop f ] if* ;

: 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 [ 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 ;