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) ?