basis/math/primes/erato/erato.factor

   1 ! Copyright (C) 2009 Samuel Tardieu.
   2 ! See http://factorcode.org/license.txt for BSD license.
   3 USING: arrays byte-arrays kernel math math.bitwise math.functions math.order
   4 math.ranges sequences sequences.private ;
   5 IN: math.primes.erato
   6
   7 <PRIVATE
   8
   9 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 }
  10
  11 : bit-pos ( n -- byte/f mask/f )
  12     30 /mod masks nth-unsafe [ drop f f ] when-zero ;
  13
  14 : marked-unsafe? ( n arr -- ? )
  15     [ bit-pos ] dip swap [ [ nth-unsafe ] [ bitand zero? not ] bi* ] [ 2drop f ] if* ;
  16
  17 : unmark ( n arr -- )
  18     [ bit-pos swap ] dip
  19     over [ [ swap unmask ] change-nth-unsafe ] [ 3drop ] if ;
  20
  21 : upper-bound ( arr -- n ) length 30 * 1 - ;
  22
  23 : unmark-multiples ( i arr -- )
  24     2dup marked-unsafe? [
  25         [ [ dup sq ] [ upper-bound ] bi* rot <range> ] keep
  26         [ unmark ] curry each
  27     ] [
  28         2drop
  29     ] if ;
  30
  31 : init-sieve ( n -- arr ) 30 /i 1 + 255 <array> >byte-array ;
  32
  33 PRIVATE>
  34
  35 : sieve ( n -- arr )
  36     init-sieve [ 2 swap upper-bound sqrt [a,b] ] keep
  37     [ [ unmark-multiples ] curry each ] keep ;
  38
  39 : marked-prime? ( n arr -- ? )
  40     2dup upper-bound 2 swap between? [ bounds-error ] unless
  41     over { 2 3 5 } member? [ 2drop t ] [ marked-unsafe? ] if ;