5
\$\begingroup\$

I wrote this program in response to a CodeGolf challenge that required generating this sequence of "untouchable" numbers. This sequence can be found on OEIS as A005114.

My initial implementation was prohibitively slow when asked for anything other than the first few terms so, I wrote this implementation with the goal of trading off memory usage for speed by maintaining an internal state representing already calculated values.

import Control.Applicative
import Control.Monad
import Control.Monad.State.Lazy
import Data.List
import qualified Data.Map as M
import System.Environment

main = do 
  n:_ <- getArgs
  print $ evalState (take (read n) <$> untouchables) M.empty

--A state monad containing a map from integers to the sum of their proper divisors
type CalcualtedSumsST = State (M.Map Int Int)  

--Generates an infinite list of untouchable numbers
untouchables :: CalcualtedSumsST [Int]
untouchables = filterM untouchable [1..]

--Tests weather a number is untouchable
untouchable :: Int -> CalcualtedSumsST Bool
untouchable 1 = return False
untouchable n = all (/=n) <$> mapM properDivisorSum [1..(n-1)*(n-1)]

--calculated the sum of the proper divisors of a number and adds it to the state
--or returns the pre-calculated if one is present
properDivisorSum :: Int -> CalcualtedSumsST Int
properDivisorSum n = do
  preCalulated <- get 
  if M.member n preCalulated then
    return $ preCalulated M.! n
  else do
    let dSum = sum $ properDivisors n
    put $ M.insert n dSum preCalulated
    return dSum

--generate a list of all the proper divisors of a number
properDivisors :: Int -> [Int]
properDivisors n = 1:(nub $ do
      --I think this is better than an equivalent list comprehension because 
      --it lets me utilize divMod instead of using div and mod separately
      x <- [2..floor.sqrt.fromIntegral $ n]
      let (q,r) = divMod n x
      if r==0 then [x,q] else [])

I'm primarily looking for review of my use of the State monad as this is an area of Haskell I have very little experience in. This code is still slower than I would like it, struggling to compute the first 50 terms of the sequence.

\$\endgroup\$

1 Answer 1

4
\$\begingroup\$

The State monad

You are using the State monad to memoize the sum of the proper divisors calculation, and there is a pure way (i.e. not using a monad) of doing that.

The easiest way is to use the Data.MemoCombinators module to build a memoizing version of your function:

import qualified Data.MemoCombinators as Memo

divSum n = ...compute the sum of proper divisors of n...

memoDivSum = Memo.integral divSum

isTouchable n = any (\k -> memoDivSum k == n) [1..n*n]

The function memoDivSum works just like divSum except that it caches previously calculated values in a trie.

Update: In practice I've found that a better way to implement isTouchable is to iterate from n*n downwards.

A better algorithm

Note that your algorithm to find the untouchable numbers <= n is actually O(n^3). To determine if k is untouchable is O(k^2) and you're doing that for all the numbers k from 1 through n.

Here is a slightly better approach:

import Data.List

touchableList n = dedup $ sort $ filter (<= n) $ map divSum [1..n*n]
  where dedup [] = []
        dedup (x:xs) = x : dedup ( dropWhile (==x) xs )

It finds the touchable numbers, but they appear in sorted order, so finding the untouchable ones shouldn't be hard and is left as an exercise.

Note you don't need to memoize divSum now because you are only calling it once for each value of k, k from 1 through n*n.

Even faster

The fastest method I can think of gets away from pure functional programming and uses mutable arrays and sieves.

\$\endgroup\$
2
  • \$\begingroup\$ Thanks for the review! I didn't know there was a term for what I was doing let alone that it was already implemented in Haskell. I'll wait a day or so before accepting an answer to give more people time to post. \$\endgroup\$ Aug 6, 2015 at 1:05
  • \$\begingroup\$ dedup is more simply written as dedup = (map head) . group \$\endgroup\$
    – Caridorc
    Aug 30, 2015 at 19:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.