Difference between revisions of "Prime numbers miscellaneous"
(→Oneliners: +1) 
m (→Oneliners) 

Line 199:  Line 199:  
. (, ([2,3],[4..]))  (xs \\ [p*p, p*p+p..]) 
. (, ([2,3],[4..]))  (xs \\ [p*p, p*p+p..]) 

−  +   segmentwise 

_Y (\ps> 2 : [n  (r:q:_, px) < (zip . tails . (2:) . map (^2) <*> inits) ps, 
_Y (\ps> 2 : [n  (r:q:_, px) < (zip . tails . (2:) . map (^2) <*> inits) ps, 

n < [r+1..q1], all ((> 0).rem n) px]) 
n < [r+1..q1], all ((> 0).rem n) px]) 
Revision as of 18:01, 21 June 2019
For the context to this, please see Prime numbers.
Contents
Implicit Heap
The following is an original implicit heap implementation for the sieve of
Eratosthenes, kept here for historical record. Also, it implements more sophisticated, lazier scheduling. The Prime_numbers#Tree merging with Wheel section simplifies it, removing the People a
structure altogether, and improves upon it by using a wheel optimization.
See also the message threads Re: "nocoding" functional data structures via lazyness for more about how merging ordered lists amounts to creating an implicit heap and Re: Code and Perf. Data for Prime Finders for an explanation of the People a
structure that makes it work.
data People a = VIP a (People a)  Crowd [a]
mergeP :: Ord a => People a > People a > People a
mergeP (VIP x xt) ys = VIP x $ mergeP xt ys
mergeP (Crowd xs) (Crowd ys) = Crowd $ merge xs ys
mergeP xs@(Crowd (x:xt)) ys@(VIP y yt) = case compare x y of
LT > VIP x $ mergeP (Crowd xt) ys
EQ > VIP x $ mergeP (Crowd xt) yt
GT > VIP y $ mergeP xs yt
merge :: Ord a => [a] > [a] > [a]
merge xs@(x:xt) ys@(y:yt) = case compare x y of
LT > x : merge xt ys
EQ > x : merge xt yt
GT > y : merge xs yt
diff xs@(x:xt) ys@(y:yt) = case compare x y of
LT > x : diff xt ys
EQ > diff xt yt
GT > diff xs yt
foldTree :: (a > a > a) > [a] > a
foldTree f ~(x:xs) = x `f` foldTree f (pairs xs)
where pairs ~(x: ~(y:ys)) = f x y : pairs ys
primes, nonprimes :: [Integer]
primes = 2:3:diff [5,7..] nonprimes
nonprimes = serve . foldTree mergeP . map multiples $ tail primes
where
multiples p = vip [p*p,p*p+2*p..]
vip (x:xs) = VIP x $ Crowd xs
serve (VIP x xs) = x:serve xs
serve (Crowd xs) = xs
nonprimes
effectively implements a heap, exploiting lazy evaluation.
Prime Wheels
The idea of only testing odd numbers can be extended further. For instance, it is a useful fact that every prime number other than 2 and 3 must be of the form or . Thus, we only need to test these numbers:
primes :: [Integer]
primes = 2:3:prs
where
1:p:candidates = [6*k+r  k < [0..], r < [1,5]]
prs = p : filter isPrime candidates
isPrime n = all (not . divides n)
$ takeWhile (\p > p*p <= n) prs
divides n p = n `mod` p == 0
Here, prs
is the list of primes greater than 3 and isPrime
does not test for divisibility by 2 or 3 because the candidates
by construction don't have these numbers as factors. We also need to exclude 1 from the candidates and mark the next one as prime to start the recursion.
Such a scheme to generate candidate numbers first that avoid a given set of primes as divisors is called a prime wheel. Imagine that you had a wheel of circumference 6 to be rolled along the number line. With spikes positioned 1 and 5 units around the circumference, rolling the wheel will prick holes exactly in those positions on the line whose numbers are not divisible by 2 and 3.
A wheel can be represented by its circumference and the spiked positions.
data Wheel = Wheel Integer [Integer]
We prick out numbers by rolling the wheel.
roll (Wheel n rs) = [n*k+r  k < [0..], r < rs]
The smallest wheel is the unit wheel with one spike, it will prick out every number.
w0 = Wheel 1 [1]
We can create a larger wheel by rolling a smaller wheel of circumference n
along a rim of circumference p*n
while excluding spike positions at multiples of p
.
nextSize (Wheel n rs) p =
Wheel (p*n) [r2  k < [0..(p1)], r < rs,
let r2 = n*k+r, r2 `mod` p /= 0]
Combining both, we can make wheels that prick out numbers that avoid a given list ds
of divisors.
mkWheel ds = foldl nextSize w0 ds
Now, we can generate prime numbers with a wheel that for instance avoids all multiples of 2, 3, 5 and 7.
primes :: [Integer]
primes = small ++ large
where
1:p:candidates = roll $ mkWheel small
small = [2,3,5,7]
large = p : filter isPrime candidates
isPrime n = all (not . divides n)
$ takeWhile (\p > p*p <= n) large
divides n p = n `mod` p == 0
It's a pretty big wheel with a circumference of 210 and allows us to calculate the first 10000 primes in convenient time.
A fixed size wheel is fine, but adapting the wheel size while generating prime numbers quickly becomes impractical, because the circumference grows very fast, as primorial, but the returns quickly diminish, the improvement being just (p1)/p
. See Euler's Sieve, or the functional pearl titled Lazy wheel sieves and spirals of primes for more.
Using IntSet for a traditional sieve
module Sieve where
import qualified Data.IntSet as I
 findNext  finds the next member of an IntSet.
findNext c is  I.member c is = c
 c > I.findMax is = error "Ooops. No next number in set."
 otherwise = findNext (c+1) is
 mark  delete all multiples of n from n*n to the end of the set
mark n is = is I.\\ (I.fromAscList (takeWhile (<=end) (map (n*) [n..])))
where
end = I.findMax is
 primes  gives all primes up to n
primes n = worker 2 (I.fromAscList [2..n])
where
worker x is
 (x*x) > n = is
 otherwise = worker (findNext (x+1) is) (mark x is)
(doesn't look like it runs very efficiently).
Oneliners
Produce an unbounded (for the most part) list of primes.
Using
(\\) = Data.List.Ordered.minus  _not_ Data.List.\\
_Y g = g (_Y g)
fix g = x where x = g x  from Data.Function
f <$> (a,b) = (a, f b)  from Data.Functor
unfold f a  (xs,b) < f a = xs ++ unfold f b
 unfold f = concat . Data.List.unfoldr (Just . f)
Here goes:
[n  n<[2..], product [1..n1] `rem` n == n1]  Wilson's theorem
Data.List.nubBy (((>1).).gcd) [2..]  the shortest one
[n  n<[2..], all ((> 0).rem n) [2..n1]]  trial division
[n  n<[2..], []==[i  i<[2..n1], rem n i==0]]  implementing `all`
[n  n<[2..], []==[i  i<[2..n1], last [n, ni..0]==0]]  and `rem`
 is this still a trial division?
[n  n<[2..], []==[i  i<[2..n1], j<[0, i..n], j==n]]
 or a "forgetful" ... sieve of Eratosthenes?
[n  n<[2..], []==[i  i<[2..n1], j<[i*i, i*i+i..n], j==n]]
[n  let p n = []==[i  i<[2..n1], rem n i==0], n<[2..], p n]
 testing only by primes should be faster, right?
[n  let p n = []==[i  i<[2..n1], p i && rem n i==0], n<[2..], p n]
 (someone actually wrote this, though in C)
fix $ map head . scanl (\\) [2..] . map (\p > [p, p+p..])  executable spec
fix $ map head . scanl ((\\).tail) [2..] . map (\p > [p*p, p*p+p..])
_Y $ (2:) . concat . snd . mapAccumL (\(x:t) ms > (t \\ ms, [x]))
[3..] . map (\p > [p*p, p*p+p..])
_Y $ (2:) . concat . snd
. mapAccumL (\xs ms > (\(h,t) > (t \\ ms, h)) $
span (< head ms) xs) [3..] . map (\p > [p*p, p*p+p..])
[n  n<[2..], all ((> 0).rem n) [2..floor.sqrt.fromIntegral$n]]
2 : [n  n<[3,5..], all ((> 0).rem n) [3,5..floor.sqrt.fromIntegral$n]]
 optimal trial division
fix (\xs> 2 : [n  n<[3..], all ((> 0).rem n) $ takeWhile ((<= n).(^2)) xs])
2 : fix (\xs> 3 : [n  n < [5,7..],  the wheel: 2,
foldr (\p r> p*p>n  (rem n p>0 && r)) True xs])
2:3 : fix (\xs> 5 : [n  n < scanl (+) 7 (cycle [4,2]),  3,
foldr (\p r> p*p>n  (rem n p>0 && r)) True xs])
2:3:5 : fix (\xs> 7 : [n  n < scanl (+) 11 (cycle [2,4,2,4,6,2,6,4]),  5
foldr (\p r> p*p>n  (rem n p>0 && r)) True xs])
foldr (\x r> x : filter ((> 0).(`rem`x)) r) [] [2..]  bottomup
nub . map head . scanl (\xs x> filter ((> 0).(`rem`x)) xs) [2..] $ [2..]
map head . iterate (\(x:xs)> filter ((> 0).(`rem`x)) xs) $ [2..]
unfold (\(x:xs)> ([x], filter ((> 0).(`rem`x)) xs)) [2..]  topdown
fix $ \ps > 2 : [ n  (p, px) < zip ps (inits ps), n < take 1  sideways
[n  n < [p+1..], and [mod n p > 0  p < px]]]
fix $ concatMap (fst.snd)  by spans between primes' squares
. iterate (\(p:t,(h,xs)) > (t,span (< head t^2) [y  y<xs, rem y p>0]))
. (, ([2,3],[4..]))  (xs \\ [p*p, p*p+p..])
 segmentwise
_Y (\ps> 2 : [n  (r:q:_, px) < (zip . tails . (2:) . map (^2) <*> inits) ps,
n < [r+1..q1], all ((> 0).rem n) px])
 n < [r+1..q1] Data.List.\\ [m  p < px, m < [0,p..q1]]])
 n < foldl (\\) [r+1..q1] [[0,p..q1]  p < px]])
 a sieve finds primes by eliminating multiples, isn't it?
 (found this one in the wild, too)
[n  n<[2..], not $ elem n [j*k  j<[2..n1], k<[2..n1]]]
[n  n<[2..], not $ elem n [j*k  j<[2..n`div`2], k<[2..n`div`j]]]
zipWith (flip (!!)) [0..] . scanl1 (\\)  APLstyle
. scanl1 (zipWith (+)) $ repeat [2..]
tail . unfold (\(a:b:t) > (:t) . (\\ b) <$> span (< head b) a)
. scanl1 (zipWith (+) . tail) $ tails [1..]
fix (\sv > dropWhile ((> 0).snd) >>> \((p,_):xs) >  counting by 1s
p : sv (zipWith (fmap . max) (cycle $ [0  _<[2..p]] ++ [1]) xs))
. (`zip` repeat 0) $ [2..]
2 : fix ((3:) . unfold (\(x,u:us)> ([xu==0], (x+2,us)))  counting by 2s
. (,) 5 . ([0,0] ++)
. foldi (\(x:xs) ys > let n=div (head ys  x) 2  1 in
x:take n xs ++ zipWith max ys (drop n xs)) []
. map (\p> (p*p :) . tail . cycle $ 1 : replicate (p1) 0) )
primesTo n = foldl (\a i> a \\ [i*i, i*i+2*i..n]) (2:[3,5..n])
[3,5..(floor . sqrt . fromIntegral) n]
primesTo n = 2 : foldr (\r z> if (head r^2) <= n then head r : z else r) []
(fix $ \rs> [3,5..n] : [t \\ [p*p, p*p+2*p..]  (p:t)<rs])
primesTo = foldr (\n ps> foldl (\a p>  with repeated square root
a \\ [p*p, p*p+p..n]) [2..n] ps)
[] . takeWhile (> 1) . iterate (floor . sqrt . fromIntegral)
concatMap (\(a:b:_)> drop (length a) b)  its dual, with repeated
. tails . scanl (\ps n> foldl (\a p>  squaring
a \\ [p*p, p*p+p..n]) [2..n] ps) [] $ iterate (^2) 2
2 : unfold (\(xs,p:ps)> let (h,t)=span (< p*p) xs in
(h, (filter ((> 0).(`rem`p)) t, ps))) ([3,5..],[3,5..])
2 : fix (\xs> 3 : unfold (\(xs,p:ps)> let (h,t)=span (< p*p) xs in
(h, ((\\ [p*p, p*p+2*p..]) t, ps))) ([5,7..], xs))
2 : _Y (\ps> concatMap snd $ iterate (\((fs:ft, x, p:t),_) >
((ft,p*p+2,t), [x  x < [x, x+2 .. p*p2],
all ((/= 0).rem x) fs])) ((inits ps, 5, ps), [3]) )
2 : _Y (\ps> concatMap snd $ iterate (\((fs:ft, x, p:t),_) >
((ft,p*p+2,t), minus [x, x+2 .. p*p2]
$ foldi union [] [[o, o+2*i .. p*p2]  i < fs,
let o=x+mod(ix)(2*i)])) ((inits ps, 5, ps), [3]) )
let { sieve (x:xs) = x : [n  n < sieve xs, rem n x > 0] }  bottomup
in sieve [2..]
let { sieve (x:xs) = x : sieve [n  n < xs, rem n x > 0] }  topdown 
in sieve [2..]  Turner's sieve
fix $ \xs> let { sieve xs (p:ps) = let (h,t)=span (< p*p) xs in  postponed
h ++ sieve (filter ((> 0).(`rem`p)) t) ps }
in 2 : 3 : sieve [5,7..] (tail xs)
fix $ \xs> let { sieve xs (p:ps) = let (h,t)=span (< p*p) xs in
h ++ sieve (t \\ [p*p, p*p+2*p..]) ps }
in 2 : 3 : sieve [5,7..] (tail xs)
 producing only unique composites
fix $ \xs> 2 : minus [3..] (foldr (\(p:ps) r> fix $ ((p*p) :) .
merge r . map (p*) . merge ps) [] $ tails xs)
2 : minus [3,5..] (foldi (\(x:xs)> (x:).union xs) []  sieving by all odds
$ map (\x>[x*x, x*x+2*x..]) [3,5..])
unfold (\a@(1:p:_) > ([p], a \\ map (*p) a)) [1..]  Euler's sieve
unfold (\(p:xs) > ([p], xs \\ map (*p) (p:xs))) [2..]
unfold (\(p:xs) > ([p], xs \\ [p, p+p..])) [2..]  the basic sieve
unfold (\xs@(p:_) > (\\ [p, p+p..]) <$> splitAt 1 xs) [2..]
fix $ (2:) . unfold (\(p:ps,xs) > (ps ,) .  postponed sieve
(\\ [p*p, p*p+p..]) <$> span (< p*p) xs) . (, [3..])
 Richard Bird's combined sieve
fix $ (2:) . minus [3..] . foldr (\(x:xs)> (x:) . union xs) []
. map (\p> [p*p, p*p+p..])
2 : _Y ( (3:)  unbounded treeshaped folding
. minus [5,7..]
. unionAll  ~= foldi (\(x:xs) ys> x:union xs ys) []
. map (\p> [p*p, p*p+2*p..]) )
[2,3,5,7] ++ _Y ( (11:)  with a wheel
. minus (scanl (+) 13 $ tail wh11)
. unionAll
. map (\p> map (p*) . dropWhile (< p) $
scanl (+) (p  rem (p11) 210) wh11) )
[2,3,5,7] ++ _Y ( (11:)  same, but 1.4x
. minus (scanl (+) 13 $ tail wh11)  slower than
. unionAll  the above
. map (\(w,p)> scanl (\c d> c + p*d) (p*p) w)
. isectBy (compare . snd)
(tails wh11 `zip` scanl (+) 11 wh11) )
wh11 = 2:4:2:4:6:2:6:4:2:4:6:6:2:6:4:2:6:4:6:8:4:2:4:2:  210wheel,
4:8:6:4:6:2:4:6:2:6:6:4:2:4:6:2:6:4:2:4:2:10:2:10:wh11  from 11
 cycle $ zipWith () =<< tail $ [i  i < [11..221], gcd i 210 == 1]
foldi
is an infinitely rightdeepening tree folding function found here. minus
and union
of course are on the main page here (and merge
a simpler, duplicatesignoring variant of union
), such that xs `minus` ys == xs Data.List.\\ ys
and xs `union` ys == nub . sort $ xs ++ ys
for any finite, increasing lists xs
and ys
.
Some definitions use functions from the dataordlist
package, and the 2357 wheel wh11
; they might be leaking space unless minus
is fused with its input (into gaps
/gapsW
from the main page).
A Tale of Sieves
Old stuff, just put together in one place to better see the whole picture at once.
{# LANGUAGE ViewPatterns #}
import Data.List ( (\\), tails, inits )
import Data.List.Ordered ( minus, unionAll )
import Data.Array.Unboxed
primes = sieve [2..]
where
sieve (p:t) = [p] ++ sieve [n  n < t, rem n p > 0]
primes = 2 : sieve primes [3..]
where
sieve (p:ps) (span (< p*p) > (h,t)) =
h ++ sieve ps [n  n < t, rem n p > 0]
primes = 2 : [n  (r:q:_, px) < (zip . tails . (2:) . map (^2) <*> inits) primes,
n < [r+1..q1], all ((> 0).rem n) px]
{ n < [r+1..q1] \\ concat [ [0,p..q]  p < px]]
n < [r+1..q1] `minus` unionAll
[ dropWhile (<= r) [0,p..q1]  p < px]]
(n,True) < assocs ( accumArray (\_ _ > False) True (r+1, q1)
[(m,())  p < px,
let s = (r+p)`div`p*p,
m < [s,s+p..q1]] :: UArray Int Bool )] }
primes = map head $ scanl minus [2..] [[p, p+p..]  p < primes]
 = fix $ map head . scanl minus [2..] . map (\p > [p, p+p..])
primes = 2 : 3 : [5,7..] `minus` unionAll [[p*p, p*p+2*p..]  p < tail primes]
