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 nth prime,

module Main (main) where
import SievePrimes (primesUpTo)
import System.Environment (getArgs)
printNthPrime :: Int > IO ()
printNthPrime n = print (n, primesUpTo (estim n) !! (n1))
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) ?
