[Haskell-cafe] Re: FASTER primes

[Will Ness]

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

> 
> 
> Am Montag 04 Januar 2010 19:16:32 schrieb Will Ness:
> > Daniel Fischer <daniel.is.fischer <at> web.de> writes:
> > > Am Montag 04 Januar 2010 13:25:47 schrieb Will Ness:
> > > > Euler's sieve is
> > > >
> > > >  sieve (p:ps) xs = h ++ sieve ps (t `minus` map (p*) [p,p+2..])
> > > >                       where (h,t) = span (< p*p) xs
> > >
> > > Not quite. That's a variant of the postponed filters, it crosses off e.g.
> > > 45 twice, once as 3*15 and once as 5*9 (okay, it has already been removed
> > > by the first, so let's say 45 appears in two lists of numbers to be
> > > removed if present).
> >
> > there won't be any such. whatever is removed first, is just skipped second
> > (and third, etc). 
> 
> ((45:(offer 47 when demanded)) `minus` (45:(next will be 51 when demanded)))
>  `minus` (45:(next will be 55 when demanded))
> 
> So there are two attempts to tell the generator to not output 45. To the
> second, it answers "I already knew that", but the request is made
> nevertheless.

yes, of course. 

> ... There are two attempts to eliminate 45.

I would say there are two requests to not have 45 in the output.

> > I don't see any problem here. As Melissa (and yourself, I think) have
> > shown, double hits on multiples are few and far between.
> 
> It's not a problem, it just means it's not Euler's sieve, because that
> attempts to eliminate each composite exactly once.

yes I see now. My bad. Wasn't reading that wikipedia page with enough attention 
to detail. It uses the modified (culled, so far) numbers to find out the next 
multiples to be eliminated, not just the prime itself like my code does.

You solution is indeed, exponential:

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

  
  primes 
   = 2:euler (as@[3,5..])
       3:euler (bs@(tail as `minus` map (3*) as))
         5:euler (cs@(tail bs `minus` map (5*) bs))


There are two separate look-back pointers to /as/ in /bs/, and there are two 
such in /cs/, to /bs/. The book-keeping machinery explodes.



> > Also, it uses no filters, i.e. no individual number divisibility testing.
> > The "filters" refers first of all to testing an individual number to decide
> > whether to keep it in or not.
> 
> Umm, the postponed filters says keep everything until p^2, then eliminate
> (filter out, remove) multiples of p in the remainder, after that, pick next
> prime.
> That's precisely what the above does. It doesn't do the filtering out by
> divisibility testing but by minus (hence more efficiently). I would say
> that's a variant of the postponed filters.
> 

Filter is usually (as in Haskell's 'filter') is about testing individual 
elements by a predicate function. There is of course a filtering effect in two 
lists elts' comparison that 'minus' performs, so that's debatable. Even the PQ 
code performs "filtering" in this wider sense.


> > > Euler's sieve is never attempting to remove a number more than once,
> > > that's
> >
> > How's that possible?
> 
> http://en.wikipedia.org/wiki/Sieve_of_Eratosthenes#Euler.27s_Sieve
> 
> C) The number after the previous prime is also a prime. *Multiply each 
> number /that's left/
> in the list starting from this prime by this prime and discard the products*.


yes. Wasn't paying attention to that, more to the "intent" of it.

There's of course enourmous vagueness in how exactly it is to be performed, in 
the unbounded case, which you uncovered here.

 
> > It can't have foresight, right?
> >
> 
> But it has :) By only trying to eliminates products of the prime p 
> currently under consideration *with numbers (>= p) /which have not/
> /yet been eliminated/ from the list*, it is known in advance that all these
> products are still in the list.


missed that.


> When p is under consideration, the list contains (apart from the primes < p) 
precisely the numbers whose smallest prime factor is >= p.
> 
> > > the outstanding feature of it. Unfortunately, that also makes it hard to
> > > implement efficiently. The problem is that you can't forget the primes
> > > between p and p^2, you must remember them for the removal of multiples of
> > > p later on.


not just primes - all the composites not yet removed, too. So it can't even be 
implemented on shared mutable storage if we want to delay the removal (as we 
must, in unbounded case) - the composites will get removed so their multiples 
must actually all be produced first! 


> > The more pressing problem with that code is its linear structure of course
> > which gets addressed by the tree-folding merge etc.
> >
> 
> Which unfortunately introduces its own space problem :(



> Try a different minus:
> 
> xxs <at> (x:xs) `minus` yys <at> (y:ys)
>    = case compare x y of
>        LT -> x : xs `minus` yys
>        EQ -> xs `minus` ys
>        GT -> error ("trying to remove " ++ show y ++ " a second time")
> 
> Your code is not. It is, however, much faster.

I understand now. Thanks!


 
> > > Its performance is really horrible though.


exponential, empyrically as well.



> It just occured to me that the accumulation list is just
>   (takeWhile (< head cs) (p:ps)), so we could improve it as
> 
> euler (p:ps) cs = h ++ euler ps (t `minus` comps)
>       where
>         (h,t) = span (< p*p) cs
>         comps = map (*p) (takeWhile (< head cs) (p:ps) ++ cs)
> 
> Still Bad.

yes, _both_ /t/ and /comps/ have their back-pointers into cs. 


> > No point in improving
> > anything until the linear structure isn't turned into a tree.
> >
> You're a tree hugger? :D
> No, the point isn't linear vs. tree, it's that Euler's sieve is a) great for
> theoretical reasoning (look at a proof of Euler's product formula for the
> (Riemann) Zeta-function to appreciate it) b) reasonably good for sieving
> with a list of numbers written on paper, but c) not good to implement on a
> computer. For the computer, keeping track of which numbers are left is so
> much more complicated than just ticking off numbers a couple of times that
> it isn't even worth attempting it.


     
> 
> 
> 
> 
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe <at> haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>