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

29 Dec 2009

      Am Mittwoch 30 Dezember 2009 01:23:32 schrieb Will Ness:
...
Daniel Fischer  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) ?