How to improve this algorithm by 1)use arrays, 2)avoid list concatenation (lazy lists?)?

Question

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]

Answer 1

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.