testing prime speed

MainCLI.hs

module Main where

import System.Environment(getArgs)

import qualified Prime

main :: IO ()
main = do
  a <- getArgs
  Prime.main (a !! 0) (read (a !! 1))

prime.c

// same algorithm and data structure as prime_array
// gcc -o prime -std=c99 -O3 -g prime.c
// ./prime 500000  3.03s user 0.00s system 99% cpu 3.033 total
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char** argv)
{
    int iMax = atoi(argv[1]);
    int* arr = malloc(sizeof(int) * (iMax + 1));
    int iN = 0;
    int n2 = 4;
    int iMT = 0;
    int x = 2;
    while (iN <= iMax) {
        if (x >= n2) {
            ++iMT;
            n2 = arr[iMT] * arr[iMT];
        }
        int prime = 1;
        for (int i = 0; i < iMT; ++i) {
            if (x % arr[i] == 0) {
                prime = 0;
                break;
            }
        }
        if (prime) {
            arr[iN++] = x;
        }
        ++x;
    }
    printf("%d\n", arr[iMax]);
    free(arr);
}

Prime.hs

{-# LANGUAGE NoMonomorphismRestriction, BangPatterns, RankNTypes #-}
module Prime where

-- compile: ghc --make -main-is Prime.main Prime.hs -O2
-- $time ./Prime i 500000
-- ./Prime i 500000  4.35s user 0.00s system 99% cpu 4.453 total
-- ./Prime s 500000  7.69s user 0.03s system 99% cpu 7.726 total (!!!)
-- lookup OPT in comments about various optimization points

import qualified Data.Vector.Generic.Mutable as M
import qualified Data.Vector.Unboxed         as U
import qualified Control.Monad.ST as ST
import Control.Monad.Primitive (PrimMonad,PrimState)



divides :: Integral a => a -> a -> Bool
divides n x = (x `mod` n) == 0
{-# INLINE divides #-}

prime_dumb :: Integral a => [a]
prime_dumb = f [2..]
  where
    f (x : xs) = x : f (filter (not . divides x) xs)

-- this is better algorithm compared to previous. It tests divisions only up to sqrt(n)
-- not optimized at all, but surprisingly fast -- less than 2 times slower than the next one and less than 3 times slower than C version
prime_sq :: Integral a => (a -> Bool) -> [a]
prime_sq err = res
  where
    res = f 4 res (takeWhile (not . err) [2..] ++ error "overflow")
    f n2 ~(p0 : pt@(p1 : _)) xs@(x : _) | x >= n2 = f (p1 * p1) pt (filter (not . divides p0) xs)
    f n2 ps (x : xt) = x : f n2 ps xt

errNo :: b -> Bool
errNo = const False

errBound :: (Bounded a, Eq a) => a -> Bool
errBound = \x -> x == maxBound


-- the same algorithm as previous, but use unboxed array
-- OPT: turning on specialization makes it SLOWER (?!) same also for iter
-- {- SPECIALIZE prime_array :: Int -> GHC.ST.ST s (Data.Array.ST.STUArray s Int Int) -}
-- {- SPECIALIZE prime_array :: Int -> IO (Data.Array.IO.IOUArray Int Int) -}
-- prime_array :: (Integral e, Data.Array.IO.MArray a e m) => Int -> m (a Int e)

prime_array :: (Integral a, PrimMonad m, M.MVector v a) 
            => Int -> m (v (PrimState m) a)
{-# INLINE prime_array #-}
prime_array iMax =
  do arr <- M.new (iMax + 1)
     fill arr 0 4 0 2 -- OPT: here was list [2 ..] with processor, switching to explicit iteration gave ~0.5s
     return arr
  where
    fill _ iN _ _ _ | iN > iMax = return ()
    fill arr iN n2 iMT x | x >= n2 = do
      nN0 <- M.unsafeRead arr iMT
      nN1 <- M.unsafeRead arr (iMT + 1)
      fill arr iN (nN1 * nN1) (iMT + 1) (x + 1)
    fill arr iN n2 iMT x = do
      isPrime <- check arr x iMT
              -- OPT: this fix-based one is 0.5s slower
              {- fix (\next i -> case i of
                                  _ | i == iMT -> return True
                                  _ -> do
                                        nT <- unsafeRead arr i
                                        if nT `divides` x then return False else next (i + 1)) 0 -}
      if isPrime
        then do
              M.write arr iN x
              fill arr (iN + 1) n2 iMT (x + 1)
        else fill arr iN n2 iMT (x + 1)
-- -}

check :: (Integral a, PrimMonad m, M.MVector v a) 
      => v (PrimState m) a -> a -> Int -> m Bool
{-# INLINE check #-}
check arr x iMT = iter 0
  where
    {-# INLINE iter #-}
    iter i | i == iMT = pure True
    iter i = do
      nT <- M.unsafeRead arr i
      if nT `divides` x
        then pure False
        else iter (i + 1)


main :: String -> Int -> IO ()
main func k =
  case func of
    "d" -> print $ prime_dumb !! k
    "s" -> print (prime_sq errBound !! k :: Int)
    "a" -> do let v = U.create (prime_array k)
              print (v U.! k :: Int)
    "i" -> do v <- U.unsafeFreeze =<< prime_array k
              print (v U.! k :: Int)

prime.js

function main(iMax)
{
    var arr = [];
    var iN = 0;
    var n2 = 4;
    var iMT = 0;
    var x = 2;
    while (iN <= iMax) {
        if (x >= n2) {
            ++iMT;
            n2 = arr[iMT] * arr[iMT];
        }
        var prime = 1;
        for (var i = 0; i < iMT; ++i) {
            if (x % arr[i] == 0) {
                prime = 0;
                break;
            }
        }
        if (prime) {
            arr[iN++] = x;
        }
        ++x;
    }
    console.log(arr[iMax]);
}

main(process.argv[2]);

at Intel(R) Core(TM)2 Duo CPU 1200MHz:
max@prime-speed$time ./prime 500000
7368791
./prime 500000 2.13s user 0.00s system 99% cpu 2.138 total
max@prime-speed$time ./Prime i 500000
7368791
./Prime i 500000 2.28s user 0.01s system 99% cpu 2.301 total
max@prime-speed$time ./Prime s 500000
7368791
./Prime s 500000 5.30s user 0.06s system 99% cpu 5.365 total
max@prime-speed$time ./Prime a 500000
7368791
./Prime a 500000 29.41s user 0.34s system 99% cpu 29.826 total

ghc 7.4.1

unsafePerformIO (print ("divides", n, x)) seq

Для этого есть Debug.Trace.traceShow

есть Debug.Trace.traceShow

спасибо, буду знать

gcc (Gentoo Hardened 5.4.0-r4 p1.7, pie-0.6.5) 5.4.0
The Glorious Glasgow Haskell Compilation System, version 8.0.2
Intel(R) Core(TM) i5-3320M CPU @ 2.60GHz

test@l29ah-x201 ~/tmp $ time ./prime 500000
7368791

real	0m0.732s
user	0m0.730s
sys	0m0.000s
test@l29ah-x201 ~/tmp $ time ./Prime i 500000
7368791

real	0m3.617s
user	0m3.603s
sys	0m0.000s
test@l29ah-x201 ~/tmp $ time ./Prime s 500000
7368791

real	0m6.074s
user	0m6.030s
sys	0m0.023s
test@l29ah-x201 ~/tmp $ time ./Prime a 500000
7368791

real	0m3.644s
user	0m3.630s
sys	0m0.000s