I tried to learn how the STArray works, but I couldn't. (Doc is poor, or at least the one I found).
Any way, I have the next algorithm, but it uses a lot of !!, which is slow. How can I convert it to use the STArray monad?
-- The Algorithm prints the primes present in [1 .. n]
main :: IO ()
main = print $ primesUpTo 100
type Nat = Int
primesUpTo :: Nat -> [Nat]
primesUpTo n = primesUpToAux n 2 [1]
primesUpToAux :: Nat -> Nat -> [Nat] -> [Nat]
primesUpToAux n current primes =
if current > n
then primes
else primesUpToAux n (current + 1) newAcum
where newAcum = case isPrime current primes of
True -> primes++[current]
False -> primes
isPrime :: Nat -> [Nat] -> Bool
isPrime 1 _ = True
isPrime 2 _ = True
isPrime x neededPrimes = isPrimeAux x neededPrimes 1
isPrimeAux x neededPrimes currentPrimeIndex =
if sqrtOfX < currentPrime
then True
else if mod x currentPrime == 0
then False
else isPrimeAux x neededPrimes (currentPrimeIndex + 1)
where
sqrtOfX = sqrtNat x
currentPrime = neededPrimes !! currentPrimeIndex
sqrtNat :: Nat -> Nat
sqrtNat = floor . sqrt . fromIntegral
Edit
Oops, the !! wasn't the problem; in the next version of the algorithm (below) I've removed the use of !!; also, I fixed 1 being a prime, which is not, as pointed by @pedrorodrigues
main :: IO ()
main = print $ primesUpTo 20000
type Nat = Int
primesUpTo :: Nat -> [Nat]
primesUpTo n = primesUpToAux n 1 []
primesUpToAux :: Nat -> Nat -> [Nat] -> [Nat]
primesUpToAux n current primesAcum =
if current > n
then primesAcum
else primesUpToAux n (current + 1) newPrimesAcum
where newPrimesAcum = case isPrime current primesAcum of
True -> primesAcum++[current]
False -> primesAcum
isPrime :: Nat -> [Nat] -> Bool
isPrime 1 _ = False
isPrime 2 _ = True
isPrime x neededPrimes =
if sqrtOfX < currentPrime
then True
else if mod x currentPrime == 0
then False
else isPrime x restOfPrimes
where
sqrtOfX = sqrtNat x
currentPrime:restOfPrimes = neededPrimes
sqrtNat :: Nat -> Nat
sqrtNat = floor . sqrt . fromIntegral
Now this question is about 2 questions really:
1.- How to transform this algorithm to use arrays instead of lists? (Is for the sake of learning how to handle state and arrays in Haskell) Which somebody already answered in the comments, but pointing to a not very good explained example.
2.- How to eliminate the concatenation of lists every time a new prime is found?
True -> primesAcum++[current]
Here's a more or less direct translation of your code into working with an unboxed array of integers:
import Control.Monad
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Control.Arrow
main :: IO ()
main = print . (length &&& last) . primesUpTo $ 1299709
primesUpTo :: Int -> [Int]
primesUpTo = takeWhile (> 0) . elems . primesUpToUA
primesUpToUA :: Int -> UArray Int Int
primesUpToUA n = runSTUArray $ do
let k = floor( 1.25506 * fromIntegral n / log (fromIntegral n)) -- WP formula
ar <- newArray (0,k+1) 0 -- all zeroes initially, two extra spaces
let
loop c i | c > n = return ar -- upper limit is reached
| otherwise = do -- `i` primes currently in the array
b <- isPrime 0 i c -- is `c` a prime?
if b then do { writeArray ar i c ; loop (c+1) (i+1) }
else loop (c+1) i
isPrime j i c | j >= i = return True -- `i` primes
| otherwise = do -- currently in the array
p <- readArray ar j
if p*p > c then return True
else if rem c p == 0 then return False
else isPrime (j+1) i c
loop 2 0
This is more or less self-explanatory, when you read it slowly, one statement at a time.
Using arrays, there's no problems with list concatenation, as there are no lists. We use the array of primes as we're adding new items to it.
Of course you can re-write your list-based code to behave better; the simplest re-write is
ps :: [Int]
ps = 2 : [i | i <- [3..],
and [rem i p > 0 | p <- takeWhile ((<=i).(^2)) ps]]
primesTo n = takeWhile (<= n) ps
The key is to switch from recursive thinking to corecursive thinking - not how to add at the end, explicitly, but to define how a list is to be produced — and let the lazy semantics take care of the rest.