[email protected]
[Top] [All Lists]

Re: [Haskell-cafe] Re: FASTER primes (was: Re: Code and Perf. Data for P

Subject: Re: [Haskell-cafe] Re: FASTER primes was: Re: Code and Perf. Data for Prime Finders (was: Genuine Eratosthenes sieve)
From: Daniel Fischer
Date: Wed, 30 Dec 2009 04:03:10 +0100

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

_______________________________________________
Haskell-Cafe mailing list
[email protected]
http://www.haskell.org/mailman/listinfo/haskell-cafe
<Prev in Thread] Current Thread [Next in Thread>