Posted on

Problem

I had a challenge a while back from a friend to try and produce the fastest possible program to compute `A(3,16)`, where `A` is the Ackermann function.

He wrote it in Java, and it took ~4.4 seconds. The below (very ugly) Haskell program that I wrote took ~1.7 seconds:

``````{-# LANGUAGE BangPatterns #-}

import           Data.Map.Strict (Map)
import qualified Data.Map.Strict as Map

-- Main method:

main :: IO ()
main = print \$ fst \$ ack (I 3 16) Map.empty

-- Definitions:

data I = I {-# UNPACK #-} !Int {-# UNPACK #-} !Int
deriving (Eq, Ord)

ack :: I -> Map I Int -> (Int, Map I Int)
ack   (I 0 n) r  = (n + 1, r)

ack i@(I m 0) r1 = maybe bak fun (Map.lookup i r1)
where
(!val, !r2) = ack (I (m-1) 1) r1
bak       = (val, Map.insert i val r2)
fun v     = (v, r1)

ack i@(I m n) r1 = maybe bak2 fun2 (Map.lookup call2 r2)
where
call1       = (I m (n - 1))
(!val1, !re1) = ack call1 r1
bak1        = (val1, Map.insert call1 val1 re1)
fun1 v      = (v, r1)
(!v1, !r2)    = maybe bak1 fun1 (Map.lookup call1 r1)
call2       = (I (m - 1) v1)
(!val2, !re2) = ack call2 r2
bak2        = (val2, Map.insert call2 val2 re2)
fun2 v      = (v, r2)
``````

This is used on `Int`s specifically because `A(3,16)` is 524285, which within `Int`s range.

How can I make this function more memory efficient, and faster for greater Ackermann calls?

Also, the syntax is disgusting. How can I make this more readable (maybe using monads?) without having horrible performance consequences?

Solution

Using Data.Memocombinators

Have a look at Data.Memocombinators module. It offers combinators for memoizing functions which is essentially what you are doing with the `Map`.

Here is the example from the documentation on how to use it to create a memoizing fibonacci function:

``````import Data.MemoCombinators

fib = Memo.integral fib'
where
fib' 0 = 0
fib' 1 = 1
fib' x = fib (x-1) + fib (x-2)
^^^         ^^^
``````

Things to note:

• `fib` is the memoized version of `fib'` which is just a helper function
• `fib'` calls `fib` for any recursive calls

The other thing which will help is to keep in mind that the type `Memo a` represents a a combinator which is able to memoize a function whose argument type is `a`.

So:

• `integral` is a combinator which can memoize functions taking an `Int`
• `pair integral integral` is a combinator which memoize functions taking an `(Int,Int)`

And thus to memoize the Ackerman function:

``````ack = (pair integral integral) ack'
where ack' (0,n) = n+1
ack' (m,0) = ack (m-1,1)
ack' (m,n) = ack (m-1, ack (m, n-1))
``````

Again, note how `ack` is defined as the memoized version of `ack'` and `ack'` calls `ack` for recursive cases.

Using Array memoization

Faster results can be obtained by using arrays to memoize the functions `ack(1,.)`, `ack(2,.)` and `ack(3,.)`:

``````import Data.Array
import Data.Ix

-- memoize the function f using arrays
arrayMemo bnds f = g
where g i = arr ! i
arr = array bnds [ (i,f arr i) | i <- range bnds ]

a0 n = n+1

a1 = arrayMemo (0,530000) f
where f arr 0 = a0 1
f arr n = a0 (arr ! (n-1))

a2 = arrayMemo (0,270000) f
where f arr 0 = a1 1
f arr n = a1 (arr ! (n-1))

a3 = arrayMemo (0,16) f
where f arr 0 = a2 1
f arr n = a2 (arr ! (n-1))
``````

The only drawback is that explicit bounds have to be determined for each level function. This approach computes `a3 16` in about half a second.