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

[Daniel Fischer]

Am Mittwoch 30 Dezember 2009 01:23:32 schrieb Will Ness:
> Daniel Fischer <daniel.is.fischer <at> web.de> writes:
> > Gee, seems my mail server reads your posts very thoroughly today :)
>
> I hope it's not a bad thing. :)
>

It means, twenty minutes after I replied to the previous, I got your hours old post which 
mentions points I could have included in the aforementioned reply.

Not really bad, but somewhat inconvenient.

> > Am Dienstag 29 Dezember 2009 14:58:10 schrieb Will Ness:
> > > If I realistically needed primes generated in a real life setting, I'd
> > > probably had to use some C for that.
> >
> > If you need the utmost speed, then probably yes. If you can do with a
> > little less, my STUArray bitsieves take about 35% longer than the
> > equivalent C code and are roughly eight times faster than ONeillPrimes.
> > I can usually live well with that.
>
>
> Wow! That's fast! (that's the code from haskellwiki's primes page, right?)
>

No, it's my own code. Nothing elaborate, just sieving numbers 6k±1, twice as fast as the 
haskellwiki code (here) and uses only 1/3 the memory. For the record:

----------------------------------------------------------------------
{-# LANGUAGE BangPatterns #-}
module SievePrimes where

import Data.Array.Base (unsafeRead, unsafeWrite, unsafeAt)
import Data.Array.ST
import Data.Array.Unboxed
import Control.Monad (when)
import Data.Bits

primesUpTo :: Integer -> [Integer]
primesUpTo bound
    = 2:3:[fromIntegral (3*i + (if testBit i 0 then 4 else 5))
            | i <- [0 .. bd], par `unsafeAt` i]
      where
        bd0 = 2*(bound `div` 6) - case bound `mod` 6 of { 0 -> 2; 5 -> 0; _ -> 1 }
        bd = fromInteger bd0
        sr = floor (sqrt (fromIntegral bound))
        sbd = 2*(sr `div` 6) - case sr `mod` 6 of { 0 -> 2; 5 -> 0; _ -> 1 }
        par :: UArray Int Bool
        par = runSTUArray $ do
                ar <- newArray (0,bd) True
                let tick i s1 s2
                        | bd < i    = return ()
                        | otherwise = do
                            p <- unsafeRead ar i
                            when p (unsafeWrite ar i False)
                            tick (i+s1) s2 s1
                    sift i
                        | i > sbd   = return ar
                        | otherwise = do
                            p <- unsafeRead ar i
                            when p (let (!start,!s1,!s2) = if even i then 
(i*(3*i+10)+7,2*i+3,4*i+7) else (i*(3*i+8)+4,4*i+5,2*i+3) in tick start s1 s2)
                            sift (i+1)
                sift 0
----------------------------------------------------------------------

The index magic is admittedly not obvious, but if you take that on faith, the rest is 
pretty clear. To find the n-th prime,

----------------------------------------------------------------------
module Main (main) where

import SievePrimes (primesUpTo)
import System.Environment (getArgs)

printNthPrime :: Int -> IO ()
printNthPrime n = print (n, primesUpTo (estim n) !! (n-1))

main = do
    args <- getArgs
    printNthPrime $ read $ head args

estim :: Int -> Integer
estim k
    | n < 13    = 37
    | n < 70    = 349
    | otherwise = round x
      where
        n = fromIntegral k :: Double
        x = n * (log n + log (log n) - 1 + 1.05 * log (log n) / log n)
----------------------------------------------------------------------

If I don't construct the prime list but count directly on the array, it's ~6% faster.

> > > If OTOH we're talking about a tutorial
> > > code that is to be as efficient as possible without loosing it clarity,
> > > being a reflection of essentials of the problem, then any overly
> > > complicated advanced Haskell wouldn't be my choice either.
> >
> > +1
> > Though perhaps we view mutable array code differently. In my view, it's
> > neither advanced nor complicated.
>
> convoluted, then. Not using "higher level" concepts, like map and foldr, :)
> or head.until isSingleton (pairwise merge).map wrap , that kind of thing.
> :)
>
> BTW I think a really smart compiler should just get a specification, like
> Turner's sieve, and just derive a good code from that by itself.

Go ahead and write one. I would love such a beast.

>
> Another example would be
>
>   qq n m = sum $ take n $ cycle [1..m]
>
> which should really get compiled into just a math formula, IMO. Now _that_
> I would call a good compiler.

Dream on, dream on, with hope in your heart.

> That way I really won't have to learn how to use STUArray`s you see.
>
> I've seen this question asked a lot, what should be a good programming
> language?
>
> IMO, English (plus math where needed, and maybe some sketches by hand). :)
>

Maybe you'd like http://en.wikipedia.org/wiki/Shakespeare_(programming_language) ?

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://www.haskell.org/pipermail/haskell-cafe/attachments/20091229/1750823b/attachment.html