I took your Haskell program as a base and have refactored it into a version that is about the same speed as your original C++ program. Will follow up with details when I have a little more time.
On Tue, Dec 18, 2012 at 12:42 PM, Branimir Maksimovic <bm...@hotmail.com>wrote: > This time I have tried fasta benchmark since current entries does not > display correct output. > Program is copy of mine > http://benchmarksgame.alioth.debian.org/u64q/program.php?test=fasta&lang=gpp&id=1 > c++ benchmark, but unfortunately executes more than twice time. > > Seems to me that culprit is in function random as I have tested rest of > code > and didn't found speed related problems. > > bmaxa@maxa:~/shootout/fasta$ time ./fastahs 25000000 > /dev/null > > real 0m5.262s > user 0m5.228s > sys 0m0.020s > > bmaxa@maxa:~/shootout/fasta$ time ./fastacpp 25000000 > /dev/null > > real 0m2.075s > user 0m2.056s > sys 0m0.012s > > Since I am planning to contribute program, perhaps someone can > see a problem to speed it up at least around 3.5 secs which is > speed of bench that display incorrect result (in 7.6.1). > > Program follows: > > {-# LANGUAGE BangPatterns #-} > {- The Computer Language Benchmarks Game > > http://shootout.alioth.debian.org/ > > contributed by Branimir Maksimovic > -} > > import System.Environment > import System.IO.Unsafe > > import Data.IORef > import Data.Array.Unboxed > import Data.Array.Storable > import Data.Array.Base > import Data.Word > > import Foreign.Ptr > import Foreign.C.Types > > type A = UArray Int Word8 > type B = StorableArray Int Word8 > type C = (UArray Int Word8,UArray Int Double) > > foreign import ccall unsafe "stdio.h" > puts :: Ptr a -> IO () > foreign import ccall unsafe "string.h" > strlen :: Ptr a -> IO CInt > > main :: IO () > main = do > n <- getArgs >>= readIO.head > > let !a = (listArray (0,(length alu)-1) > $ map (fromIntegral. fromEnum) alu:: A) > make "ONE" "Homo sapiens alu" (n*2) $ Main.repeat a (length alu) > make "TWO" "IUB ambiguity codes" (n*3) $ random iub > make "THREE" "Homo sapiens frequency" (n*5) $ random homosapiens > > make :: String -> String -> Int -> IO Word8 -> IO () > {-# INLINE make #-} > make id desc n f = do > let lst = ">" ++ id ++ " " ++ desc > a <- (newListArray (0,length lst) > $ map (fromIntegral. fromEnum) lst:: IO B) > unsafeWrite a (length lst) 0 > pr a > make' n 0 > where > make' :: Int -> Int -> IO () > make' !n !i = do > let line = (unsafePerformIO $ > newArray (0,60) 0 :: B) > if n > 0 > then do > !c <- f > unsafeWrite line i c > if i+1 >= 60 > then do > pr line > make' (n-1) 0 > else > make' (n-1) (i+1) > else do > unsafeWrite line i 0 > l <- len line > if l /= 0 > then pr line > else return () > > pr :: B -> IO () > pr line = withStorableArray line (\ptr -> puts ptr) > len :: B -> IO CInt > len line = withStorableArray line (\ptr -> strlen ptr) > > repeat :: A -> Int -> IO Word8 > repeat xs !n = do > let v = unsafePerformIO $ newIORef 0 > !i <- readIORef v > if i+1 >= n > then writeIORef v 0 > else writeIORef v (i+1) > return $ xs `unsafeAt` i > > random :: C -> IO Word8 > random (a,b) = do > !rnd <- rand > let > find :: Int -> IO Word8 > find !i = > let > !c = a `unsafeAt` i > !p = b `unsafeAt` i > in if p >= rnd > then return c > else find (i+1) > find 0 > > rand :: IO Double > {-# INLINE rand #-} > rand = do > !seed <- readIORef last > let > newseed = (seed * ia + ic) `rem` im > newran = fromIntegral newseed * rimd > rimd = 1.0 / (fromIntegral im) > im, ia, ic :: Int > im = 139968 > ia = 3877 > ic = 29573 > writeIORef last newseed > return newran > where > last = unsafePerformIO $ newIORef 42 > > alu :: [Char] > alu = > "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG\ > \GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA\ > \CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT\ > \ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA\ > \GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG\ > \AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC\ > \AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA" > > mkCum :: [(Char,Double)] -> [(Word8,Double)] > mkCum lst = map (\(c,p) -> ((fromIntegral.fromEnum) c,p)) $ > scanl1 (\(_,p) (c',p') -> (c', p+p')) lst > > homosapiens, iub :: C > > iub' = mkCum [('a',0.27),('c',0.12),('g',0.12),('t',0.27),('B',0.02) > ,('D',0.02),('H',0.02),('K',0.02),('M',0.02),('N',0.02) > ,('R',0.02),('S',0.02),('V',0.02),('W',0.02),('Y',0.02)] > > homosapiens' = mkCum [('a',0.3029549426680),('c',0.1979883004921) > ,('g',0.1975473066391),('t',0.3015094502008)] > > iub = (listArray (0, (length iub')-1) $ map fst iub', > listArray (0, (length iub')-1) $ map snd iub') > > homosapiens = (listArray (0, (length homosapiens')-1) $ map fst > homosapiens', > listArray (0, (length homosapiens')-1) $ map snd > homosapiens') > > > _______________________________________________ > Haskell-Cafe mailing list > Haskell-Cafe@haskell.org > http://www.haskell.org/mailman/listinfo/haskell-cafe > >
_______________________________________________ Haskell-Cafe mailing list Haskell-Cafe@haskell.org http://www.haskell.org/mailman/listinfo/haskell-cafe