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 fry kernel kernel.private locals math
   5 math.functions math.private sequences sequences.private ;
   6
   7 IN: math.primes.erato.fast
   8
   9 <PRIVATE
  10
  11 CONSTANT: wheel-2-3-5-7 B{
  12     2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8 4 2 4 2 4 8 6 4 6 2
  13     4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10
  14 }
  15
  16 :: each-prime ( upto sieve quot -- )
  17     11 upto >fixnum '[ dup _ <= ] [
  18         wheel-2-3-5-7 [
  19             over dup 2/ sieve nth-unsafe [ drop ] quot if
  20             fixnum+fast
  21         ] each
  22     ] while drop ; inline
  23
  24 :: mark-multiples ( i upto sieve -- )
  25     i 2 fixnum*fast :> step
  26     i i fixnum*fast upto >fixnum '[ dup _ <= ] [
  27         t over 2/ sieve set-nth-unsafe
  28         step fixnum+fast
  29     ] while drop ; inline
  30
  31 : sieve-bits ( n -- bits )
  32     210 /i 1 + 210 * 2/ 6 + ; inline
  33
  34 PRIVATE>
  35
  36 :: make-sieve ( n -- sieve )
  37     n sieve-bits <bit-array> :> sieve
  38     t 0 sieve set-nth
  39     t 4 sieve set-nth
  40     n sqrt sieve [ n sieve mark-multiples ] each-prime
  41     sieve ; inline
  42
  43 :: sieve ( n -- primes )
  44     V{ 2 3 5 7 } clone :> primes
  45     n dup make-sieve [
  46         dup n <= [ primes push ] [ drop ] if
  47     ] each-prime primes ;
  48
  49 :: sieve-prime? ( i sieve -- prime? )
  50     i even? [ f ] [ i 2/ sieve nth ] if ;