[Haskell-cafe] Re: FASTER primes

[Will Ness]

Daniel Fischer <daniel.is.fischer <at> web.de> writes:


> roll       = scanl (+)
> wheel      = 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:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wheel
> wheel11 = res
>       where
>         snms = scanl (+) 11 (take 47 wheel)
>         nums = tail $ scanl (+) 11 (take 529 wheel)
>         cops = nums `minus` map (*11) snms
>         diffs = zipWith (-) (tail cops) cops
>         res = foldr (:) res diffs
> wheel13 = res
>       where
>         snms = take 480 $ scanl (+) 13 wheel11
>         nums = take (480*13+1) . tail $ scanl (+) 13 wheel11
>         cops = nums `minus` map (*13) snms
>         diffs = zipWith (-) (tail cops) cops
>         res = foldr (:) res diffs
> 


BTW have you seen my take on the "faithful Euler's sieve"? It shows another way 
to look at the wheels, which for me was really the first time I really 
understood what's going on there.

It also makes for easier wheel extention (IMO):


>   euler ks@(p:rs) = p : euler (rs `minus` map (*p) ks)
>   primes = euler [2..]


The essence of Euler's sieve is the wheel: after each step we're left with 
lists of non-multiples of the preceding primes:


    primes  = euler $ listFrom [2] 1
     = 2:euler ( listFrom [3] 1 `minus` map(2*) (listFrom [2] 1)) )
                 listFrom [3,4] 2 `minus` listFrom [4] 2
                       -- listFrom [3] 2 --
     = 2:3:euler (listFrom [5] 2 `minus` map(3*) (listFrom [3] 2))
                  listFrom [5,7,9] 6 `minus` listFrom [9] 6
                       -- listFrom [5,7] 6 --
     = 2:3:5:euler (listFrom [7,11] 6 `minus` listFrom [25,35] 30)
                      [7,11, 13,17, 19,23, 25,29, 31,35] 30
                    -- listFrom [7,11,13,17,19,23,29,31] 30 --
     = .....

where

  listFrom xs by = concat $ iterate (map (+ by)) xs

so

  -- startRoll = ([2],1)

  nextRoll r@(xs@(x:xt),b) 
             = ( (x,r') , r')
     where
           ys@(y:_) = xt ++ [x + b]
           r' = (xs',b')
           b' = x*b
           xs' = takeWhile (< y + b') (listFrom ys b) 
                 `minus` map (x*) xs

  rolls = unfoldr (Just . nextRoll) ([2],1)

  nthWheel n = let (ps,rs) = unzip $ take n rolls
                   (x:xs,b) = last rs
               in ((ps, x), zipWith (-) (xs++[x+b]) (x:xs))

{-
 *Main> mapM print $ take 4 rolls
 (2,([3],2))
 (3,([5,7],6))
 (5,([7,11,13,17,19,23,29,31],30))
 (7,([11,13,17,19,23,29,31,37,41,43,47,53,59,61,67,71,73,79,83,89,97,
      101,103,107,109,113,121,127,131,137,139,143,149,151,157,163,167,169,
      173,179,181,187,191,193,197,199,209,211],210))

 *Main> nthWheel 3
 (([2,3,5],7),[4,2,4,2,4,6,2,6])

 *Main> nthWheel 4
 (([2,3,5,7],11),[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,4,6,2,6,6,4,2,4,6,2,6,4,2,4,2,10,2,10])
-}


Coincidentally, the function specified  by


 eulerPrimes n = let (ps,rs) = unzip $ take n rolls
                     (qs@(q:_),b) = last rs
                 in ps ++ takeWhile (< q^2) qs


can be used to write the specialized nthEulerPrime etc., whose complexity seems 
to be about O(n^1.5).


Maybe this reinvents some spiraling wheels or somethn'. :)