HiveBrain v1.2.0
Get Started
← Back to all entries
patternMinor

Untouchable Numbers

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
untouchablenumbersstackoverflow

Problem

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:_  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  [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.

Solution

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.

Code Snippets

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]
import Data.List

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

Context

StackExchange Code Review Q#99066, answer score: 4

Revisions (0)

No revisions yet.