[Haskell-cafe] Re: FASTER primes
Will Ness
will_n48 at yahoo.com
Fri Jan 8 15:05:18 EST 2010
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'. :)
More information about the Haskell-Cafe
mailing list