[Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
apfelmus quantentunnel.de> writes: > > Dave Bayer wrote: > > What I'm calling a "venturi" > > > >venturi :: Ord a => [[a]] -> [a] > > > > merges an infinite list of infinite lists into one list, under the > > assumption that each list, and the heads of the lists, are in > > increasing order. > > > > I wrote this as an experiment in coding data structures in Haskell's > > lazy evaluation model, rather than as explicit data. The majority of the > > work done by this code is done by "merge"; the multiples of each prime > > percolate up through a tournament consisting of a balanced tree of > > suspended "merge" function calls. In order to create an infinite lazy > > balanced tree of lists, the datatype > > > >data List a = A a (List a) | B [a] > > > > is used as scaffolding. One thinks of the root of the infinite tree as > > starting at the leftmost child, and crawling up the left spine as > > necessary. > > After some pondering, the List a data structure for merging is really > ingenious! :) Here's a try to explain how it works: > > The task is to merge a number of sorted and infinite lists when it's > known that their heads are in increasing order. In particular, we want > to write > > primes = (2:) $ diff [3..] $ venturi $ map multiple primes > > Thus, we have a (maybe infinite) list > > xss = [xs1, xs2, xs3, ...] > > of infinite lists with the following properties > > all sorted xss > sorted (map head xss) > > where sorted is a function that returns True if the argument is a > sorted list. A first try to implement the merging function is > > venturi xss = foldr1 merge xss > = xs1 `merge` (xs2 `merge` (xs3 `merge` ... > > where merge is the standard function to merge to sorted lists. > > However, there are two problems. The first problem is that this doesn't > work for infinite lists since merge is strict in both arguments. But > the property head xs1 < head xs2 < head xs3 < ... we missed to exploit > yet can now be used in the following way > > venturi xss = foldr1 merge' xss > > merge' (x:xt) ys = x : merge xt ys > > In other words, merge' is biased towards the left element > > merge' (x:_|_) _|_ = x : _|_ > > which is correct since we know that (head xs < head ys). > > The second problem is that we want the calls to merge to be arranged > as a balanced binary tree since that gives an efficient heap. It's not > so difficult to find a good shape for the infinite tree, the real > problem is to adapt merge' to this situation since it's not associative: > > .. > > The problem is that the second form doesn't realize that y is also > smaller than the third argument. In other words, the second form has to > treat more than one element as "privileged", namely x1,x2,... and y. > This can be done with the aforementioned list data structure > > data People a = VIP a (People a) | Crowd [a] > > The people (VIPs and crowd) are assumed to be _sorted_. Now, we can > start to implement > > merge' :: Ord a => People a -> People a -> People a Hi, ... replying to a two-years-old post here, :) :) and after consulting the full "VIP" version in haskellwiki/Prime_Numers#Implicit_Heap ... It is indeed the major problem with the merged multiples removing code (similar one to Richard Bird's code from Melissa O'Neill's JFP article) - the linear nature of foldr, requiring an (:: a->b->b) merge function. To make it freely composable to rearrange the list into arbitrary form tree it must indeed be type uniform (:: a->a->a) first, and associative second. The structure of the folded tree should be chosen to better suit the primes multiples production. I guestimate the total cost as Sum (1/p)*d, where p is a generating prime at the leaf, and d the leaf's depth, i.e. the amount of merge nodes its produced multiple must pass on its way to the top. The structure used in your VIP code, 1+(2+(4+(8+...))), can actually be improved upon with another, (2+4)+( (4+8)+( (8+16)+...)), for which the estimated cost is about 10%-12% lower. This can be expressed concisely as the following: primes :: () -> [Integer] primes () = 2:primes' where primes'= [3,5] ++ drop 2 [3,5..] `minus` comps mults = map (\p-> fromList [p*p,p*p+2*p..]) $ primes' (comps,_) = tfold mergeSP (pairwise mergeSP mults) fromList (x:xs) = ([x],xs) tfold f (a: ~(b: ~(c:xs))) = (a `f` (b `f` c)) `f` tfold f (pairwise f xs) pairwise f (x:y:ys) = f x y : pairwise f ys mergeSP (a,b) ~(c,d) = let (bc,b') = spMerge b c in (a ++ bc, merge b' d) where spMerge u@(x:xs) w@(y:ys) = case compare x y of LT -> (x:c,d) where (c,d) = spMerge xs w EQ -> (x:c,d) where (c,d) = spMerge xs ys GT -> (y:c,d) where (c,d) = spMerge u ys spMerge u [] = ([], u) spMerge [] w = ([], w) with ''merge'' and ''minus'' defined in the usual way. Its run times ar
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Hello, On Wednesday 25 July 2007 01:42, Thorkil Naur wrote: > Hello Melissa, > > On Tuesday 24 July 2007 19:09, Melissa O'Neill wrote: > > ... > > (See ONeillPrimes.hs in http://www.cs.hmc.edu/~oneill/code/haskell- > > primes.zip for the complete code. I've also added Thorkil Naur's > > code from March, which is also a good performer, > > Do you have detailed measurements that you wish to share? I would be most > interested, I assure you. > > > although its another > > case where most readers would find a detailed explanation of the code > > instructive.) > > I'll do my very best to provide such an explanation within, say, the next > couple of weeks. > ... And now that time has come, so brace yourselves. For your convenience, my "code from March" is thorkilnaur.dk/~tn/T64_20070303_1819.tar.gz See also a preliminary description in http://www.haskell.org/pipermail/haskell-cafe/2007-March/023095.html. The new explanation is here: thorkilnaur.dk/~tn/Haskell/EratoS/EratoS2.txt Best regards Thorkil ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Hello Melissa, On Tuesday 24 July 2007 19:09, Melissa O'Neill wrote: > apfelmus wrote: > > After some pondering, the List a data structure for merging is > > really > > ingenious! :) Here's a try to explain how it works: > > Thanks apfelmus! A detailed explanation of this code is really > helpful for anyone trying to understand what is going on. The VIP/ > Crowd analogy is also very nice. > > > I think that this approach has the potential to outperform O'Neill's > > (old?) code because it doesn't incorporates another prime number to > > the > > sieve mechanism until it really has to. I mean the following: in the > > > > 1st, 2nd, 3rd, 4th, ... step, > > > > O'Neill's code adds the multiples > > > > 2*2, 3*3, 5*5, 7*7, ... > > > > to the priority queue and uses them to sieve for potential prime > > numbers > > from then on. > > Yeah, that's (only) what the code in my paper does -- it's good > enough for explicative purposes, but all recent versions have used a > slightly augmented priority queue. It's a priority queue coupled > with a "feeder list" that provides entries to add to the queue (in > increasing order). They are only admitted to the heap data structure > only once when the root of the heap "gets there". > > The two most important bits are: > > type HybridQ k v = (PriorityQ k v, [(k,v)]) > > -- postRemoveHQ is called when the min element of the heap > has changed > postRemoveHQ :: Ord k => HybridQ k v -> HybridQ k v > postRemoveHQ mq@(pq, []) = mq > postRemoveHQ mq@(pq, (qk,qv) : qs) > | qk < minKeyPQ pq = (insertPQ qk qv pq, qs) > | otherwise= mq > > > (See ONeillPrimes.hs in http://www.cs.hmc.edu/~oneill/code/haskell- > primes.zip for the complete code. I've also added Thorkil Naur's > code from March, which is also a good performer, Do you have detailed measurements that you wish to share? I would be most interested, I assure you. > although its another > case where most readers would find a detailed explanation of the code > instructive.) I'll do my very best to provide such an explanation within, say, the next couple of weeks. > > > the approach here adds 5*5=25 to the heap only after considering > > the 9th prime 23. > > Yep, that's what mine does too. > > Best Regards, > > Melissa. > > ___ > Haskell-Cafe mailing list > Haskell-Cafe@haskell.org > http://www.haskell.org/mailman/listinfo/haskell-cafe > Thanks a lot and the best regards Thorkil ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
apfelmus wrote: After some pondering, the List a data structure for merging is really ingenious! :) Here's a try to explain how it works: Thanks apfelmus! A detailed explanation of this code is really helpful for anyone trying to understand what is going on. The VIP/ Crowd analogy is also very nice. I think that this approach has the potential to outperform O'Neill's (old?) code because it doesn't incorporates another prime number to the sieve mechanism until it really has to. I mean the following: in the 1st, 2nd, 3rd, 4th, ... step, O'Neill's code adds the multiples 2*2, 3*3, 5*5, 7*7, ... to the priority queue and uses them to sieve for potential prime numbers from then on. Yeah, that's (only) what the code in my paper does -- it's good enough for explicative purposes, but all recent versions have used a slightly augmented priority queue. It's a priority queue coupled with a "feeder list" that provides entries to add to the queue (in increasing order). They are only admitted to the heap data structure only once when the root of the heap "gets there". The two most important bits are: type HybridQ k v = (PriorityQ k v, [(k,v)]) -- postRemoveHQ is called when the min element of the heap has changed postRemoveHQ :: Ord k => HybridQ k v -> HybridQ k v postRemoveHQ mq@(pq, []) = mq postRemoveHQ mq@(pq, (qk,qv) : qs) | qk < minKeyPQ pq = (insertPQ qk qv pq, qs) | otherwise= mq (See ONeillPrimes.hs in http://www.cs.hmc.edu/~oneill/code/haskell- primes.zip for the complete code. I've also added Thorkil Naur's code from March, which is also a good performer, although its another case where most readers would find a detailed explanation of the code instructive.) the approach here adds 5*5=25 to the heap only after considering the 9th prime 23. Yep, that's what mine does too. Best Regards, Melissa. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Dave Bayer wrote: > What I'm calling a "venturi" > >venturi :: Ord a => [[a]] -> [a] > > merges an infinite list of infinite lists into one list, under the > assumption that each list, and the heads of the lists, are in > increasing order. > > I wrote this as an experiment in coding data structures in Haskell's > lazy evaluation model, rather than as explicit data. The majority of the > work done by this code is done by "merge"; the multiples of each prime > percolate up through a tournament consisting of a balanced tree of > suspended "merge" function calls. In order to create an infinite lazy > balanced tree of lists, the datatype > >data List a = A a (List a) | B [a] > > is used as scaffolding. One thinks of the root of the infinite tree as > starting at the leftmost child, and crawling up the left spine as > necessary. After some pondering, the List a data structure for merging is really ingenious! :) Here's a try to explain how it works: The task is to merge a number of sorted and infinite lists when it's known that their heads are in increasing order. In particular, we want to write primes = (2:) $ diff [3..] $ venturi $ map multiple primes Thus, we have a (maybe infinite) list xss = [xs1, xs2, xs3, ...] of infinite lists with the following properties all sorted xss sorted (map head xss) where sorted is a function that returns True if the argument is a sorted list. A first try to implement the merging function is venturi xss = foldr1 merge xss = xs1 `merge` (xs2 `merge` (xs3 `merge` ... where merge is the standard function to merge to sorted lists. However, there are two problems. The first problem is that this doesn't work for infinite lists since merge is strict in both arguments. But the property head xs1 < head xs2 < head xs3 < ... we missed to exploit yet can now be used in the following way venturi xss = foldr1 merge' xss merge' (x:xt) ys = x : merge xt ys In other words, merge' is biased towards the left element merge' (x:_|_) _|_ = x : _|_ which is correct since we know that (head xs < head ys). The second problem is that we want the calls to merge to be arranged as a balanced binary tree since that gives an efficient heap. It's not so difficult to find a good shape for the infinite tree, the real problem is to adapt merge' to this situation since it's not associative: merge' xs (merge' (y:_|_) _|_) = merge' xs (y:_|_) = x1:x2:..:y:_|_ =/= merge' (merge' xs (y:_|_)) _|_ = merge' (x1:x2:...:y:_|_) _|_ = x1:_|_ The problem is that the second form doesn't realize that y is also smaller than the third argument. In other words, the second form has to treat more than one element as "privileged", namely x1,x2,... and y. This can be done with the aforementioned list data structure data People a = VIP a (People a) | Crowd [a] The people (VIPs and crowd) are assumed to be _sorted_. Now, we can start to implement merge' :: Ord a => People a -> People a -> People a The first case is merge' (VIP x xt) ys = VIP x (merge' xt ys) In other words, the invariant is that every VIP on the left of a merge' is guaranteed to be smaller than anyone on the right and thus will be served first. The next case merge' (Crowd xs) (Crowd ys) = Crowd (merge xs ys) is clear since it doesn't involve the invariant. What about the last case merge' (Crowd xs) (VIP y yt) = ?? Here, someone from the crowd xs may be smaller than y. But should we return a crowd or a VIP? The crucial answer is to always return a VIP merge' xs@(Crowd (x:xt)) ys@(VIP y yt) | x <= y= VIP x (merge' (Crowd xt) ys) | otherwise = VIP y (merge' xs yt) because doing otherwise would turn a VIP into a member of some crowd. But turning x into a VIP is no problem since that doesn't violated the invariant. Now merge' is associative and everything works as we want. I think that this approach has the potential to outperform O'Neill's (old?) code because it doesn't incorporates another prime number to the sieve mechanism until it really has to. I mean the following: in the 1st, 2nd, 3rd, 4th, ... step, O'Neill's code adds the multiples 2*3, 3*3, 5*5, 7*7, ... to the priority queue and uses them to sieve for potential prime numbers from then on. But the approach here adds 5*5=25 to the heap only after considering the 9th prime 23. > >> trim p = let f m x = mod x m /= 0 in filter (f p) > > lurks in the prime sieve code, but it is only used with primes of size > up to the square root of the largest output prime. I tried more > thoughtful alternatives, and they all slowed down the sieve. Sometimes > dumb is beautiful. I still think that filtering the wheel when generating multiples is wrong. In fact, most of the algorithmic work would be done here if there wasn't the lucky coincidence that "it is only used with primes of size up to the square root of the largest output prime". You are saved to some extend by the fact that the cle
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Neil Mitchell wrote: -O3 is slower than -O2 and -O in ghc? If you want "fast code" then specify -O2, not -O3. Oops. That ought to have been -O2. But from what I can tell, -O3 is harmless (at least in this case). Both invocations generate identical executables, at least for these examples on my machine. Melissa. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Hi But for the current version of my code, there is still a bit of a performance gap between our two methods. Here are the stats I get (ghc -O3, 2.4GHz x86): Are you aware that -O3 is slower than -O2 and -O in ghc? If you want "fast code" then specify -O2, not -O3. Thanks Neil ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Dave Bayer wrote: Here is another prime sieve. It's great to know people are still having fun with this stuff... I've added your implementation to the zipfile available at http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip (FWIW, I added specializations for Int and Integer and also primesToNth and primesToLimit). It is about half the length of the fastest contributed code (ONeillPrimes) I think you must be comparing it against an older version of my code (the zip file above now contains the most up to date version). That older version actually contained two slightly different copies of the algorithm. The more recent version doesn't. FWIW, here are the statistics I see for lines of code (ignoring comments, blank lines, and compiler directives): ONeillPrimes: 91 lines, 750 words, 4111 chars, 75628 bytes in .o file BayerPrimes: 72 lines, 604 words, 2649 chars, 74420 bytes in .o file So, I'd say the difference is at best 25% in source size and 2% in final binary size. But in reality, a big chunk of my code is a general purpose heap/ priority-queue implementation. If priority queue operations were provided as a standard part of Haskell in the same way that lists and BSTs are, the statistics would be: ONeillPrimes: 47 lines, 331 words, 2039 chars and nearly as fast Your results are great! It actually beats early versions of my method, before I made a couple of tweaks to improve performance. But for the current version of my code, there is still a bit of a performance gap between our two methods. Here are the stats I get (ghc -O3, 2.4GHz x86): 1*10^6 | 2*10^6 | 3*10^6 | 4*10^6 | 5*10^6 | 6*10^6 ---+++++- ONeillPrimes1.36 | 3.08 | 4.98 | 6.98 | 9.05 | 11.21 ONeillPrimes* 1.35 | 3.07 | 4.94 | 6.93 | 8.99 | 11.14 BayerPrimes 2.18 | 4.49 | 8.99 | 11.18 | 16.60 | 25.77 The "*" version is one that uses ``calcPrimes()'' rather than ``primes'' to provide its list of primes, and thereby avoids needless remembering of the whole list of primes and all the memory that entails. until it blows up on garbage collection: I think that is the biggest issue with many of the other prime generators. At a glance (just looking at RSS with Unix's top command), your space usage seems like its about double mine. And when I use ``calcPrimes()'' rather than ``primes'' I barely need any space at all (O(sqrt(N)) at the point where I've calculated N primes. The difference there is striking -- a couple of MB vs hundreds. Anyway, fun stuff... Melissa. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
It appears that at least on gmane, my posts to this thread ended up as singletons, breaking the thread. Here are the posts: http://article.gmane.org/gmane.comp.lang.haskell.cafe/26426 http://article.gmane.org/gmane.comp.lang.haskell.cafe/26466 ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
As an exercise, trying to understand the beautiful paper
Stream Fusion
From Lists to Streams to Nothing at All
Duncan Coutts, Roman Leshchinskiy and Don Stewart
http://www.cse.unsw.edu.au/~dons/papers/CLS07.html
http://www.cse.unsw.edu.au/~dons/streams.html
I recoded my prime sieve using a pared down version of their Stream
datatype; this is the simplest version I could write that achieves a
significant speedup.
My reaction to their paper was, if streams are better internally than
lists, why not code directly in streams? Lists enjoy a serious
notational advantage in Haskell, but one could imagine a language
where the list notation was reserved for stream semantics.
My sieve was spending half its time in "merge", so I made only the
changes necessary to convert "merge" to use streams. My streams are
infinite, and "merge" can be written to not use Skip, so Step goes away.
Even though "nextx" and "nexty" only have one case now, using case
statements is significantly faster than using let or where clauses.
I'm imagining that I read about this somewhere, but if I did, it
didn't sink in until I was tuning this code. I don't know if this is
related to fusion optimization, or a general effect.
The timings are
[Integer] -O2 1*10^6 2*10^6 3*10^6 4*10^6 5*10^6
-
ONeillPrimes | 3.338 | 7.320 | 11.911 | 18.225 | 21.785
StreamPrimes | 3.867 | 8.405 | 13.656 | 21.542 | 37.640
BayerPrimes | 3.960 | 8.940 | 18.528 | 33.221 | 38.568
Here is the code:
{-# OPTIONS_GHC -fglasgow-exts #-}
module StreamPrimes (primes) where
-- stream code
data Stream a = forall s. Stream (s -> (a,s)) s
data AStream a = A a (AStream a) | B (Stream a)
stream :: [a] -> Stream a
stream xs = Stream next xs
where
next [] = undefined
next (x:xt) = (x,xt)
astream :: [a] -> AStream a
astream [] = undefined
astream (x:xt) = A x $ B $ stream xt
merge :: Ord a => Stream a -> Stream a -> Stream a
merge (Stream nextx vs) (Stream nexty ws) =
Stream next (vt,ws,Left v)
where
(v,vt) = nextx vs
next (xs,ys,Left x) =
case nexty ys of
(y,yt) ->
if x < y
then (x,(xs,yt,Right y))
else (y,(xs,yt,Left x))
next (xs,ys,Right y) =
case nextx xs of
(x,xt) ->
if x < y
then (x,(xt,ys,Right y))
else (y,(xt,ys,Left x))
mergeA :: Ord a => AStream a -> AStream a -> AStream a
mergeA (A x xt) ys = A x (mergeA xt ys)
mergeA (B xs) ys = mergeB xs ys
mergeB :: Ord a => Stream a -> AStream a -> AStream a
mergeB s@(Stream next xs) ys@(A y yt) =
case next xs of
(x,xt) ->
if x < y
then A x (mergeB (Stream next xt) ys)
else A y (mergeB s yt)
mergeB xs (B ys) = B $ merge xs ys
-- Code for venturi :: Ord a => [[a]] -> [a]
root :: Ord a => AStream a -> [AStream a] -> [a]
root (A x xt) yss= x : (root xt yss)
root (B xs) (ys:yst) = root (mergeB xs ys) yst
root _ _ = undefined
pair :: Ord a => [AStream a] -> [AStream a]
pair (x:y:xt) = mergeA x y : (pair xt)
pair _ = undefined
group :: Ord a => [AStream a] -> [AStream a]
group (x:xt) = x : (group $ pair xt)
group _ = undefined
venturi :: Ord a => [[a]] -> [a]
venturi (x:xt) = root (astream x) $ group $ map astream xt
venturi _ = undefined
-- Code for primes :: Integral a => [a]
diff :: Ord a => [a] -> [a] -> [a]
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)
diff _ _ = undefined
trim :: Integral a => a -> [a] -> [a]
trim p = let f m x = mod x m /= 0 in filter (f p)
seed :: Integral a => [a]
seed = [2,3,5,7,11,13,17]
wheel :: Integral a => [a]
wheel = drop 1 [ m*j + k | j <- [0..], k <- ws ]
where m = foldr1 (*) seed
ws = foldr trim [1..m] seed
multiples :: Integral a => [a] -> [[a]]
multiples ws = map fst $ tail $ iterate g ([], ws)
where g (_,ps@(p:pt)) = ([ m*p | m <- ps ], trim p pt)
g _ = undefined
primes :: Integral a => [a]
primes = seed ++ (diff wheel $ venturi $ multiples wheel)
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Here is another prime sieve. It is about half the length of the fastest contributed code (ONeillPrimes) and nearly as fast until it blows up on garbage collection: % cat ONeillPrimes.hs | grep -v "^--" | wc 18511056306 % cat BayerPrimes.hs | grep -v "^--" | wc 85 5662418 [Integer] -O 1*10^6 2*10^6 3*10^6 4*10^6 5*10^6 - ONeillPrimes | 3.555 | 7.798 | 12.622 | 18.927 | 23.529 BayerPrimes | 3.999 | 8.895 | 18.003 | 22.977 | 38.053 I wrote this as an experiment in coding data structures in Haskell's lazy evaluation model, rather than as explicit data. The majority of the work done by this code is done by "merge"; the multiples of each prime percolate up through a tournament consisting of a balanced tree of suspended "merge" function calls. In order to create an infinite lazy balanced tree of lists, the datatype data List a = A a (List a) | B [a] is used as scaffolding. One thinks of the root of the infinite tree as starting at the leftmost child, and crawling up the left spine as necessary. What I'm calling a "venturi" venturi :: Ord a => [[a]] -> [a] merges an infinite list of infinite lists into one list, under the assumption that each list, and the heads of the lists, are in increasing order. This could be a generally useful function. If one can think of a better way to write "venturi", swapping in your code would in particular yield a faster prime sieve. I found that a tertiary merge tree was faster than a binary merge tree, because this leads to fewer suspensions. One can speed this code up a bit by interleaving strict and lazy calls, but I prefer to leave the code short and readable. It is bizarre that trim p = let f m x = mod x m /= 0 in filter (f p) lurks in the prime sieve code, but it is only used with primes of size up to the square root of the largest output prime. I tried more thoughtful alternatives, and they all slowed down the sieve. Sometimes dumb is beautiful. Thanks to apfelmus for various helpful remarks that lead me to think of this approach. Here is the code: module BayerPrimes (venturi,primes) where -- Code for venturi :: Ord a => [[a]] -> [a] merge :: Ord a => [a] -> [a] -> [a] -> [a] merge xs@(x:xt) ys@(y:yt) zs@(z:zt) | x <= y = if x <= z then x : (merge xt ys zs) else z : (merge xs ys zt) | otherwise = if y <= z then y : (merge xs yt zs) else z : (merge xs ys zt) merge _ _ _ = undefined data List a = A a (List a) | B [a] mergeA :: Ord a => List a -> List a -> List a -> List a mergeA (A x xt) ys zs = A x (mergeA xt ys zs) mergeA (B xs) ys zs = mergeB xs ys zs mergeB :: Ord a => [a] -> List a -> List a -> List a mergeB xs@(x:xt) ys@(A y yt) zs = case compare x y of LT -> A x (mergeB xt ys zs) EQ -> A x (mergeB xt yt zs) GT -> A y (mergeB xs yt zs) mergeB xs (B ys) zs = mergeC xs ys zs mergeB _ _ _ = undefined mergeC :: Ord a => [a] -> [a] -> List a -> List a mergeC xs@(x:xt) ys@(y:yt) zs@(A z zt) | x < y = if x < z then A x (mergeC xt ys zs) else A z (mergeC xs ys zt) | otherwise = if y < z then A y (mergeC xs yt zs) else A z (mergeC xs ys zt) mergeC xs ys (B zs) = B $ merge xs ys zs mergeC _ _ _ = undefined root :: Ord a => List a -> [List a] -> [a] root (A x xt) yss = x : (root xt yss) root (B xs) (ys:zs:yst) = root (mergeB xs ys zs) yst root _ _ = undefined wrap :: [a] -> List a wrap [] = B [] wrap (x:xt) = A x $ B xt triple :: Ord a => [List a] -> [List a] triple (x:y:z:xs) = mergeA x y z : (triple xs) triple _ = undefined group :: Ord a => [List a] -> [List a] group (x:y:xt) = x : y : (group $ triple xt) group _ = undefined venturi :: Ord a => [[a]] -> [a] venturi (x:xt) = root (wrap x) $ group $ map wrap xt venturi _ = undefined -- Code for primes :: Integral a => [a] diff :: Ord a => [a] -> [a] -> [a] 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) diff _ _ = undefined trim :: Integral a => a -> [a] -> [a] trim p = let f m x = mod x m /= 0 in filter (f p) seed :: Integral a => [a] seed = [2,3,5,7,11,13,17] wheel :: Integral a => [a] wheel = drop 1 [ m*j + k | j <- [0..], k <- ws ] where m = foldr1 (*) seed ws = foldr trim [1..m] seed multiples :: Integral a => [a] -> [[a]] multiples ws = map fst $ tail $ iterate g ([], ws) where g (_,ps@(p:pt)) = ([ m*p | m <- ps ], trim p pt) g _ = undefined primes :: Integral a => [a] primes = seed ++ (diff wheel $ venturi $ multiples wheel) ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
G'day all. This one is pretty elegant. A Pritchard sieve is actually an Eratosthenes sieve with the loops reversed. Unfortunately, it's a bit slower. Maybe someone else can speed it up a bit. mergeRemove :: [Integer] -> [Integer] -> [Integer] mergeRemove [] ys = [] mergeRemove xs [] = xs mergeRemove xs'@(x:xs) ys'@(y:ys) = case compare x y of LT -> x : mergeRemove xs ys' EQ -> mergeRemove xs ys GT -> mergeRemove xs' ys pritchardSieve :: Integer -> [Integer] pritchardSieve n | n <= 16 = takeWhile (<=n) [2,3,5,7,11,13] | otherwise = removeComposites [2..n] (sieve [2..n`div`2]) where removeComposites ps [] = ps removeComposites ps (cs@(c:_):css) = removeComposites' ps where removeComposites' [] = [] removeComposites' (p:ps) | p < c = p : removeComposites' ps | otherwise = removeComposites (mergeRemove ps cs) css pjs = pritchardSieve sn sn = isqrt n sieve [] = [] sieve (f:fs) = composites pjs : sieve fs where composites [] = [] composites (p:ps) | pf > n || f `mod` p == 0 = [pf] | otherwise = pf : composites ps where pf = p*f Cheers, Andrew Bromage ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Here's another program you can add. It's fairly short and efficient. -- Lennart import System (getArgs) infixr :> data StreamInt = !Int :> StreamInt (!>) :: StreamInt -> Int -> Int (x :> _) !> 0 = x (_ :> xs) !> n = xs !> (n-1) -- By replacing lprimes on the next line by '5 :> gen 7 4 2' this algorithm -- runs in very little space, but is somewhat slower. primes = 2 :> 3 :> lprimes where isPrime (p:>ps) n = n `rem` p /= 0 && (p*p > n || isPrime ps n) lprimes = 5 :> gen 7 4 2 gen n a b = if isPrime lprimes n then n :> gen (n+a) b a else gen (n+a) b a printNthPrime n = print (n, primes !> (n-1)) main = do args <- getArgs printNthPrime $ read $ head args On Feb 25, 2007, at 12:51 , Melissa O'Neill wrote: For those enjoying the fun with prime finding, I've updated the source at http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip I've tweaked my code a little to improve its space behavior when finding primes up to some limit, added an up-to-limit version of the Naive Primes algorithm, and added Oleg's prime finding code too. I also got a chance to look at space usage more generally. I won't reproduce a table here, but the conclusions were more-or-less what you'd expect. The "unlimited list" algorithms used O(n) space to find n primes (except for Runciman's algorithm, which appeared to be much worse), and the "primes up to a limit" algorithms used O (sqrt(n)) space to find the nth prime. Both of these are better than the classic C algorithm, which uses O (n log n) space to find the nth prime. For example, heap profiling shows that my own O(sqrt(n)) algorithm uses only 91200 bytes to find the 10^7th prime, whereas the classic C algorithm needs at least 11214043 bytes for its array -- a factor of more than 100 different, and one that gets worse for larger n. Lennart Augustsson wrote: Another weird thing is that much of the Haskell code seems to work with Integer whereas the C code uses int. Originally, I was comparing Haskell with Haskell, and for that purpose I wanted to have a level playing field, so going with Integer everywhere made sense. That doesn't seem fair. Actually, to the extent that any of the comparisons are "fair", I think this one is too. After all, typical Haskell code uses Integer and typical C code uses int. I could use arrays in my Haskell code and never use laziness, but when I program in Haskell, I'm not trying to exactly recreate C programs, but rather write their Haskell equivalents. For example, to me, producing a lazy list was essential for a true Haskell feel. For some people, the "Haskell feel" also includes treating the language as a declarative specification language where brevity is everything -- but for me, other things (like fundamental algorithmic efficiency and faithfulness to the core ideas that make the Sieve of Eratosthenes an *efficient* algorithm) are universal and ought to be common to both C and Haskell versions. But to allow a better comparison with C, I've added a run for an Int version of my algorithm. With that change, my code is closer to the speed of the C code. More interestingly, for larger n, I seem to be narrowing the gap. At 10^6, my code runs nearly 30 times slower than the classic C version, but at 10^8, I'm only about 20 times slower. This is especially interesting to me there was some (reasonable looking) speculation from apfelmus several days ago, that suggested that my use of a priority queue incurred an extra log(n) overhead, from which you would expect a worse asymptotic complexity, not equivalent or better. Melissa. Enc. (best viewed with a fixed-width font) -- Time (in seconds) for Number of Primes Algorithm 10^310^4 10^5 10^6 10^7 10^8 -- C-Sieve 0.00 0.00 0.01 0.29 5.1288.24 O'Neill (#3) 0.01 0.04 0.55 8.34122.62 1779.18 O'Neill (#2) 0.01 0.06 0.9513.85194.96 2699.61 O'Neill (#1) 0.01 0.07 1.0715.95230.11 - Bromage 0.02 0.39 6.50 142.85 - - "sieve" (#3) 0.01 0.25 7.28 213.19 - - Naive (#2)0.02 0.5914.70 386.40 - - Naive (#1)0.32 0.6616.04 419.22 - - Runciman 0.02 0.7429.25- - - Reinke0.04 1.2141.00- - - Zilibowitz0.02 2.50 368.33- - - Gale (#1) 0.12 17.99-- - - "sieve" (#1) 0.16 32.59-- - - "sieve" (#2) 0.01 32.76--
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
For those enjoying the fun with prime finding, I've updated the source at http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip I've tweaked my code a little to improve its space behavior when finding primes up to some limit, added an up-to-limit version of the Naive Primes algorithm, and added Oleg's prime finding code too. I also got a chance to look at space usage more generally. I won't reproduce a table here, but the conclusions were more-or-less what you'd expect. The "unlimited list" algorithms used O(n) space to find n primes (except for Runciman's algorithm, which appeared to be much worse), and the "primes up to a limit" algorithms used O(sqrt (n)) space to find the nth prime. Both of these are better than the classic C algorithm, which uses O(n log n) space to find the nth prime. For example, heap profiling shows that my own O(sqrt(n)) algorithm uses only 91200 bytes to find the 10^7th prime, whereas the classic C algorithm needs at least 11214043 bytes for its array -- a factor of more than 100 different, and one that gets worse for larger n. Lennart Augustsson wrote: Another weird thing is that much of the Haskell code seems to work with Integer whereas the C code uses int. Originally, I was comparing Haskell with Haskell, and for that purpose I wanted to have a level playing field, so going with Integer everywhere made sense. That doesn't seem fair. Actually, to the extent that any of the comparisons are "fair", I think this one is too. After all, typical Haskell code uses Integer and typical C code uses int. I could use arrays in my Haskell code and never use laziness, but when I program in Haskell, I'm not trying to exactly recreate C programs, but rather write their Haskell equivalents. For example, to me, producing a lazy list was essential for a true Haskell feel. For some people, the "Haskell feel" also includes treating the language as a declarative specification language where brevity is everything -- but for me, other things (like fundamental algorithmic efficiency and faithfulness to the core ideas that make the Sieve of Eratosthenes an *efficient* algorithm) are universal and ought to be common to both C and Haskell versions. But to allow a better comparison with C, I've added a run for an Int version of my algorithm. With that change, my code is closer to the speed of the C code. More interestingly, for larger n, I seem to be narrowing the gap. At 10^6, my code runs nearly 30 times slower than the classic C version, but at 10^8, I'm only about 20 times slower. This is especially interesting to me there was some (reasonable looking) speculation from apfelmus several days ago, that suggested that my use of a priority queue incurred an extra log(n) overhead, from which you would expect a worse asymptotic complexity, not equivalent or better. Melissa. Enc. (best viewed with a fixed-width font) -- Time (in seconds) for Number of Primes Algorithm 10^310^4 10^5 10^6 10^7 10^8 -- C-Sieve 0.00 0.00 0.01 0.29 5.1288.24 O'Neill (#3) 0.01 0.04 0.55 8.34122.62 1779.18 O'Neill (#2) 0.01 0.06 0.9513.85194.96 2699.61 O'Neill (#1) 0.01 0.07 1.0715.95230.11 - Bromage 0.02 0.39 6.50 142.85 - - "sieve" (#3) 0.01 0.25 7.28 213.19 - - Naive (#2)0.02 0.5914.70 386.40 - - Naive (#1)0.32 0.6616.04 419.22 - - Runciman 0.02 0.7429.25- - - Reinke0.04 1.2141.00- - - Zilibowitz0.02 2.50 368.33- - - Gale (#1) 0.12 17.99-- - - "sieve" (#1) 0.16 32.59-- - - "sieve" (#2) 0.01 32.76-- - - Oleg 0.18 68.40-- - - Gale (#2) 1.36268.65-- - - -- - The dashes in the table mean "I gave up waiting" (i.e., > 500 seconds) - "sieve" (#1) is the classic example we're all familiar with - "sieve" (#2) is the classic example, but sieving a list without multiples of 2,3,5, or 7 -- notice how it makes no real difference - "sieve" (#3) is the classic example, but generating a lazy-but- finite list (see below) - O'Neill (#1) is basically the algorithm of mine discussed in http:// www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf, with a few minor tweaks - O'Neill (#2) is a variant of that algorithm that
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
OK. Another weird thing is that much of the Haskell code seems to work with Integer whereas the C code uses int. That doesn't seem fair. -- Lennart On Feb 25, 2007, at 02:40 , Melissa O'Neill wrote: Someone asked if I'd include a "classic C" version of the Sieve in my comparisons. Having done so, Lennart wrote (slightly rephrased): How did you compare the C version with the Haskell versions? The Haskell programs produce the Nth prime, whereas the C code produces the last prime less than M. True. But since I have to know what M is to find the Nth prime, it's easy enough to ask the C code to produce the right prime. To make the C code to what the Haskell code does you need to set some upper bound that is related to the prime number distribution. I see no trace of this in your code. The Haskell versions that go up to a limit do this, so I could easily have written code to do it -- it's not hard, but has no real bearing on the time complexity of the code, so I didn't bother. You could argue that it's cheating to tell it so blatantly when to stop, but I hate the C code I'd found enough that I didn't really want to touch it any more than I had to. A much more legitimate complaint about the comparison with the C code is actually on space usage. It uses much more space than some of the algorithms it's competing with. More about that in an upcoming message. Melissa. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
oneill: > Someone asked if I'd include a "classic C" version of the Sieve in my > comparisons. Having done so, Lennart wrote (slightly rephrased): > >How did you compare the C version with the Haskell versions? The > >Haskell programs produce the Nth prime, whereas the C code produces > >the last prime less than M. I've taken the liberty of adding the benchmark programs to the nobench suite, http://www.cse.unsw.edu.au/~dons/code/nobench/spectral/primes2007/ so they'll be run across a range of haskell compilers. -- Don ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Someone asked if I'd include a "classic C" version of the Sieve in my comparisons. Having done so, Lennart wrote (slightly rephrased): How did you compare the C version with the Haskell versions? The Haskell programs produce the Nth prime, whereas the C code produces the last prime less than M. True. But since I have to know what M is to find the Nth prime, it's easy enough to ask the C code to produce the right prime. To make the C code to what the Haskell code does you need to set some upper bound that is related to the prime number distribution. I see no trace of this in your code. The Haskell versions that go up to a limit do this, so I could easily have written code to do it -- it's not hard, but has no real bearing on the time complexity of the code, so I didn't bother. You could argue that it's cheating to tell it so blatantly when to stop, but I hate the C code I'd found enough that I didn't really want to touch it any more than I had to. A much more legitimate complaint about the comparison with the C code is actually on space usage. It uses much more space than some of the algorithms it's competing with. More about that in an upcoming message. Melissa. ___ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
How did you compare the C version with the Haskell versions?
The Haskell programs produce the Nth prime, whereas the C code
produces the last prime less than N.
To make the C code to what the Haskell code does you need to set some
upper bound that is related to the prime number distribution. I see
no trace of this in your code.
-- Lennart
On Feb 23, 2007, at 05:27 , Melissa O'Neill wrote:
Bulat Ziganshin asked:
but how it looks compared with classic C implementation of sieve
algorithm?
It's still worse. I Googled for a "typical" implementation and
added it to the collection. The best Haskell implementation is
still about two orders of magnitude slower, but remember that the
Haskell versions we'd looked at so far are able to incrementally
produce a list of primes of arbitrary length.
Andrew Bromage wrote:
Just to fill out the implementations:
http://andrew.bromage.org/darcs/numbertheory/
Math/Prime.hs has an implementation of the Atkin-Bernstein sieve.
Cool, thanks. When I ran your code trying to find the 10,000th
prime, I got
AtkinSieveTest: Ix{Integer}.index: Index (36213) out of range
((0,36212))
but that went away when I made your array one bigger.
Here's the updated table...
--
Time (in seconds) for Number of Primes
Algorithm 10^310^4 10^510^6 10^7 10^8
--
C-Sieve 0.00 0.00 0.01 0.29 5.1288.24
O'Neill (#2) 0.01 0.09 1.4522.41393.28-
O'Neill (#1) 0.01 0.14 2.9347.08- -
Bromage 0.02 0.39 6.50 142.85- -
"sieve" (#3) 0.01 0.25 7.28 213.19- -
Naive 0.32 0.6616.04 419.22- -
Runciman 0.02 0.7429.25 - - -
Reinke0.04 1.2141.00 - - -
Gale (#1) 0.12 17.99- - - -
"sieve" (#1) 0.16 32.59- - - -
"sieve" (#2) 0.01 32.76- - - -
Gale (#2) 1.36268.65- - - -
--
Notes:
- Bromage is Andrew Bromage's implementation of the Atkin-Bernstein
sieve. Like O'Neill (#2) and "sieve" (#3), asks for some upper
limit on the number of primes it generates. Unlike O'Neill (#2)
and "sieve" (#3), it uses arrays, and the upper limit causes a
large initial array allocation. Also, unlike the other Haskell
algorithms, it does not produce a lazy list; no output is produced
until sieving is complete
- C-Sieve is a "typical" simple implementation of the sieve in C
found with Google; it skips multiples of 2 and uses a bit array.
Also, obviously, it doesn't produce incremental output.
I've also updated the zip file of implementations at
http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip
Enjoy,
Melissa.
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
Re: [Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
G'day all.
Quoting Melissa O'Neill <[EMAIL PROTECTED]>:
> Cool, thanks. When I ran your code trying to find the 10,000th
> prime, I got
>AtkinSieveTest: Ix{Integer}.index: Index (36213) out of range
> ((0,36212))
> but that went away when I made your array one bigger.
Fixed, thanks.
Cheers,
Andrew Bromage
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
[Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
Bulat Ziganshin asked:
but how it looks compared with classic C implementation of sieve
algorithm?
It's still worse. I Googled for a "typical" implementation and added
it to the collection. The best Haskell implementation is still about
two orders of magnitude slower, but remember that the Haskell
versions we'd looked at so far are able to incrementally produce a
list of primes of arbitrary length.
Andrew Bromage wrote:
Just to fill out the implementations:
http://andrew.bromage.org/darcs/numbertheory/
Math/Prime.hs has an implementation of the Atkin-Bernstein sieve.
Cool, thanks. When I ran your code trying to find the 10,000th
prime, I got
AtkinSieveTest: Ix{Integer}.index: Index (36213) out of range
((0,36212))
but that went away when I made your array one bigger.
Here's the updated table...
--
Time (in seconds) for Number of Primes
Algorithm 10^310^4 10^510^6 10^7 10^8
--
C-Sieve 0.00 0.00 0.01 0.29 5.1288.24
O'Neill (#2) 0.01 0.09 1.4522.41393.28-
O'Neill (#1) 0.01 0.14 2.9347.08- -
Bromage 0.02 0.39 6.50 142.85- -
"sieve" (#3) 0.01 0.25 7.28 213.19- -
Naive 0.32 0.6616.04 419.22- -
Runciman 0.02 0.7429.25 - - -
Reinke0.04 1.2141.00 - - -
Gale (#1) 0.12 17.99- - - -
"sieve" (#1) 0.16 32.59- - - -
"sieve" (#2) 0.01 32.76- - - -
Gale (#2) 1.36268.65- - - -
--
Notes:
- Bromage is Andrew Bromage's implementation of the Atkin-Bernstein
sieve. Like O'Neill (#2) and "sieve" (#3), asks for some upper limit
on the number of primes it generates. Unlike O'Neill (#2) and
"sieve" (#3), it uses arrays, and the upper limit causes a large
initial array allocation. Also, unlike the other Haskell algorithms,
it does not produce a lazy list; no output is produced until sieving
is complete
- C-Sieve is a "typical" simple implementation of the sieve in C
found with Google; it skips multiples of 2 and uses a bit array.
Also, obviously, it doesn't produce incremental output.
I've also updated the zip file of implementations at
http://www.cs.hmc.edu/~oneill/code/haskell-primes.zip
Enjoy,
Melissa.
___
Haskell-Cafe mailing list
Haskell-Cafe@haskell.org
http://www.haskell.org/mailman/listinfo/haskell-cafe
