[Haskell-cafe] Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)

[apfelmus]

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 cleverly constructed heap doesn't force
the multiples of prime p until it starts considering candidate primes >=
p*p.

Also, large parts of the filtered lists have to be kept in memory since
subsequent primes need to filter them further. Using a fresh throw-away
wheel for every prime eliminates this memory eater.


> until it blows up on garbage collection

Maybe the major reason for this is that many parts of shape of the
ternary tree are forced too early. In other words, irrefutable pattern
matches in  triple , group  and  root  should do the job. Also, using
explicitly infinite lists is probably better.

    -- an explicitly infinite list
  data Stream a = a :> Stream a

  triples ~(x :> ~(y :> ~(z :> zs))) = mergeA x y z : triples zs
  group   ~(x :> ~(y :> ys)) = x :> y :> group (triples ys)

  root (A x xt) yss                 = x :> (root xt yss)
  root (B xs)  ~(ys :> ~(zs :> zss) = root (mergeB xs ys zs) zss

I'm not sure whether it's really necessary to establish a perfectly
balanced tree with  root . Maybe it's better to have the smaller prime
multiples near the top of the tree since they yield the minimum more
often. In other words, I'd simply use

  venturi ~(xs :> ~(ys :> yss)) = mergeA xs ys (venturi $ triples yss)

> Thanks to apfelmus for various helpful remarks that lead me to think of
> this approach.
λ(^_^)λ

Regards,
apfelmus