Multiple updates in ST-Monad

Question

I want to learn using the ST-Monad. Therefore I want to rewrite a bit of code computing for every integer – up to a limit – the list of all its proper divisors. The result should be an array and the entry of index 'n' should be the list of it's proper divisors.

It is done by computing for each Integer 'n' a list 'l' of its multiples and adding for each multiple 'm' out of 'l' at the index 'm' it's divisor 'n' to the list.

Here's the code I want to modify:

properDivisorsOf' :: forall a. (Integral a, Ix a) => a -> Array a [a]
properDivisorsOf' limit =
  let generate :: (Integral a, Ix a) => a -> Array a [a] -> Array a [a]
      generate n acc
        | n > (limit `div` 2) = acc
        | otherwise           =
              let acc' = acc // [(i, n : (acc ! i)) | i <- [2*n, 3*n .. limit]]
              in  generate (n + 1) acc'

  in generate 1 (array (1, limit) [(i, [])| i <- [1..limit]])

And that's the way how I try it:

properDivisorsOf :: forall a. (Integral a, Ix a) => a -> Array a [a]
properDivisorsOf limit =
  let result :: ST s (STArray s a [a])
      result = newArray (1, limit) [] -- In the beginning for every number: no divisors known (empty list) 

      update (index, divisor) = do
        l <- readArray result index -- extracting list of divisors of number 'index'
        let u = divisor : l
        writeArray result index u   -- and adding 'divisor' to the list

      content :: [(a, a)]
      content = do
        n <- [1 .. (limit `div` 2)]
        multiple <- [2*n, 3*n .. limit]
        return (multiple, n)

      doUpdate = map update content -- update result for all multiples (in content)

in runSTArray result

Unfortunatley it is not compiling and the error message doesn't say anything to me. I have two questions:

1. Why doesn't it compile? How can I extract the entry correctly?
2. How would an experienced Haskell-Programm solve this problem under the restriction that he has to work with the ST-Monad (for efficiency purposes)

EDIT: Compiler message

    Couldn't match expected type `[t0]'
            with actual type `STArray i0 a [a]'
    In the second argument of `(:)', namely `l'
    In the expression: divisor : l
    In an equation for `u': u = divisor : l

Failed, modules loaded: none.

Answer 1

Solution

2. How would an experienced Haskell-Programm solve this problem under the restriction that he has to work with the ST-Monad (for efficiency purposes)?

I'm by no means a experienced Haskell programmer, but I would use the following code the imperative code below, but here's the direct transition from your code:

properDivisorsOf :: forall a. (Integral a, Ix a) => a -> Array a [a]
properDivisorsOf limit =
  runSTArray $ do
    result <- newArray (1, limit) []
    mapM_ (update result) content
    return result
  where
      content :: [(a, a)]
      content = do
        n <- [1 .. (limit `div` 2)]
        multiple <- [2*n, 3*n .. limit]
        return (multiple, n)

      update arr (index, divisor) = do
        l <- readArray arr index -- extracting list of divisors of number 'index'
        let u = divisor : l
        writeArray arr index u   -- and adding 'divisor' to the list

Comparing non-ST vs ST:

non-ST variant

Your original function:

main = print $ properDivisorsOf' 100000

$ .\Original.exe +RTS -s
Original.exe: out of memory

We exchange 100000 by 10000:

     221,476,488 bytes allocated in the heap
      21,566,328 bytes copied during GC
     172,813,068 bytes maximum residency (9 sample(s))
       4,434,480 bytes maximum slop
             210 MB total memory in use (5 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0       378 colls,     0 par    0.41s    0.43s     0.0011s    0.0024s
  Gen  1         9 colls,     0 par    0.36s    0.37s     0.0409s    0.1723s

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.28s  (  0.60s elapsed)
  GC      time    0.77s  (  0.80s elapsed)
  EXIT    time    0.00s  (  0.02s elapsed)
  Total   time    1.05s  (  1.42s elapsed)

  %GC     time      73.1%  (56.4% elapsed)

  Alloc rate    787,471,957 bytes per MUT second

  Productivity  26.9% of total user, 19.8% of total elapsed

Even though the program is pretty fast (1s after all), the memory footprint of 210MB is quite noticeable. Also, the memory needed grows squre, for a limit of 20000 you would need around 800 MB, for 20000 around 1.9GB.

ST variant

$ .\ST.exe +RTS -s > $null
     300,728,400 bytes allocated in the heap
      89,696,288 bytes copied during GC
      13,628,272 bytes maximum residency (10 sample(s))
         292,972 bytes maximum slop
              38 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0       565 colls,     0 par    0.14s    0.14s     0.0002s    0.0008s
  Gen  1        10 colls,     0 par    0.09s    0.08s     0.0076s    0.0223s

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.11s  (  0.16s elapsed)
  GC      time    0.23s  (  0.21s elapsed)
  EXIT    time    0.00s  (  0.00s elapsed)
  Total   time    0.34s  (  0.38s elapsed)

  %GC     time      68.2%  (56.6% elapsed)

  Alloc rate    2,749,516,800 bytes per MUT second

  Productivity  31.8% of total user, 28.9% of total elapsed

Only 38 MB, which is ~17% of your original implementation, but calculates 10 times more values (10000 vs 100000).

Imperative variant

If you throw content away, you'll end up with the following program:

import Data.Array (Array(..), Ix)
import Data.Array.ST (newArray, readArray, writeArray, runSTArray)
import Control.Monad (forM_)

properDivisorsOf :: (Integral a, Ix a) => a -> Array a [a]
properDivisorsOf limit =
  runSTArray $ do
    result <- newArray (1, limit) []
    forM_ [1.. (limit `div`2)] $ \n -> do
      forM_ [2*n, 3*n .. limit] $ \index -> do
        l <- readArray result index
        writeArray result index (n:l)
    return result

main = print $ properDivisorsOf 100000

$ .\Imperative.exe +RTS -s > $null
     237,323,392 bytes allocated in the heap
      63,304,856 bytes copied during GC
      13,628,276 bytes maximum residency (7 sample(s))
         138,636 bytes maximum slop
              34 MB total memory in use (0 MB lost due to fragmentation)

                                    Tot time (elapsed)  Avg pause  Max pause
  Gen  0       447 colls,     0 par    0.12s    0.09s     0.0002s    0.0008s
  Gen  1         7 colls,     0 par    0.05s    0.06s     0.0087s    0.0224s

  INIT    time    0.00s  (  0.00s elapsed)
  MUT     time    0.11s  (  0.18s elapsed)
  GC      time    0.17s  (  0.16s elapsed)
  EXIT    time    0.00s  (  0.00s elapsed)
  Total   time    0.30s  (  0.34s elapsed)

  %GC     time      57.9%  (45.9% elapsed)

  Alloc rate    2,169,813,869 bytes per MUT second

  Productivity  42.1% of total user, 36.9% of total elapsed

Why doesn't it compile?

1. Why doesn't it compile? How can I extract the entry correctly?

I a sense, you extract correctly (see above, that's almost the same code you used), but the inferred type for update is wrong, since your usage isn't correct. For example update should work in the ST monad, and therefore you would use it with mapM, not map. Also, there are some other things off, for example you define doUpdate but never use it.