[Haskell-cafe] Re: FASTER primes

[Will Ness]

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

> 
> 
> Am Dienstag 05 Januar 2010 14:49:58 schrieb Will Ness:
> > > ... There are two attempts to eliminate 45.
> >
> > I would say there are two requests to not have 45 in the output.
> >
> Thers are many possible ways to phrase it.
> 
> >
> > 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))
> >
> >

Re-write:

 euler s = head s:euler (tail s `minus` map(head s*) s)
 primes  = euler [2..]

 primes  = euler $ rollFrom [2] 1
  = 2:euler ( rollFrom [3] 1 `minus` map(2*) (rollFrom [2] 1)) )
              rollFrom [3,4] 2 `minus` rollFrom [4] 2
                    -- rollFrom [3] 2 --
  = 2:3:euler (rollFrom [5] 2 `minus` map(3*) (rollFrom [3] 2))
               rollFrom [5,7,9] 6 `minus` rollFrom [9] 6
                    -- rollFrom [5,7] 6 --
  = 2:3:5:euler (rollFrom [7,11] 6 `minus` rollFrom [25,35] 30)
                   [7,11, 13,17, 19,23, 25,29, 31,35] 30
                 -- rollFrom [7,11,13,17,19,23,29,31] 30 --
  = .....

   where
     rollOnce (x:xs) by = (x, tail xs ++ [x+by])
     rollFrom xs by = concat $ iterate (map (+ by)) (xs)
     multRoll xs@(x:_) by p = takeWhile (< (x+p*by)) $ rollFrom xs by


so, reifying, we get


 data Roll a = Roll [a] a

 rollOnce (Roll (x:xs) by) 
  = (x,Roll (xs ++ [x+by]) by)

 rollFrom (Roll xs by) 
  = concat $ iterate (map (+ by)) (xs)

 multRoll r@(Roll (x:_) by) p 
  = Roll (takeWhile (< (x+p*by)) $ rollFrom r) (by*p)

 primes  = euler $ Roll [2] 1
 euler r@(Roll xs _)
        = x:euler (Roll (mxs `minus` map (x*) xs)  mby)
  where  
   (x,r') = rollOnce r
   (Roll mxs mby) = multRoll r' x


There's much extra primes pre-calculated inside the Roll, of course (upto p^2 in
fact, for p = primes!!n ), so this needs to be somehow tapped into, by writing a
specialized

 nthPrime n = ....

to be called instead of (!! n), and also

 primesUpTo n = ....


This calculated primes !! 1000 == 7927 in 3 seconds for me, interpreted, on my
old slowish laptop. It's supposed to have all the primes upto 7927^2 = 62837329
inside the Roll (if I'm not mistaken - or is it?). That's about 3.72 millionth
prime, according to WolframAlpha. (nah, it cant be that much). But it is surely
not horribly slow.

Is this, in fact, the wheels' "spiral"?


> 
> >
> > not just primes - all the composites not yet removed, too.
> 
> Between p and p^2, there are only primes left, fortunately.

but (map (*p) ks) needs access to the original, non-culled numbers - primes,
composites and all. (?)



> 
> > So it can't even be implemented on shared mutable storage if we
> > want to delay the removal (as we must, in unbounded case) -
> 
> Yes. And it's not nice in the bounded case either.
> 
> 
> 
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe <at> haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>