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

   1 ! Copyright (C) 2015 John Benediktsson
   2 ! See http://factorcode.org/license.txt for BSD license
   3
   4 USING: bit-arrays kernel literals math math.functions
   5 math.private ranges math.statistics sequences
   6 sequences.private ;
   7
   8 IN: math.primes.erato.fast
   9
  10 <PRIVATE
  11
  12 CONSTANT: wheel-2-3-5-7 $[
  13     11 dup 210 + [a..b] [
  14         { 2 3 5 7 } [ divisor? ] with none?
  15     ] B{ } filter-as differences
  16 ]
  17
  18 :: each-prime ( upto sieve quot -- )
  19     11 upto integer>fixnum-strict '[ dup _ <= ] [
  20         wheel-2-3-5-7 [
  21             over dup 2/ sieve nth-unsafe [ drop ] quot if
  22             fixnum+fast
  23         ] each
  24     ] while drop ; inline
  25
  26 :: mark-multiples ( i upto sieve -- )
  27     i 2 fixnum*fast :> step
  28     i i fixnum*fast upto integer>fixnum-strict '[ dup _ <= ] [
  29         t over 2/ sieve set-nth-unsafe
  30         step fixnum+fast
  31     ] while drop ; inline
  32
  33 : sieve-bits ( n -- bits )
  34     210 /i 1 + 210 * 2/ 6 + ; inline
  35
  36 PRIVATE>
  37
  38 :: make-sieve ( n -- sieve )
  39     n sieve-bits <bit-array> :> sieve
  40     t 0 sieve set-nth
  41     t 4 sieve set-nth
  42     n sqrt >integer sieve
  43     [ n sieve mark-multiples ] each-prime
  44     sieve ; inline
  45
  46 :: sieve ( n -- primes )
  47     V{ 2 3 5 7 } clone :> primes
  48     n dup make-sieve [
  49         dup n <= [ primes push ] [ drop ] if
  50     ] each-prime primes ;
  51
  52 :: marked-prime? ( i sieve -- prime? )
  53     i dup even? [ 2 = ] [ 2/ sieve nth not ] if ;