I’m beginning to learn F# and the best way for me to go about that is to start applying it. There are some interesting computational problems over at Project Euler . This post is the continuation of a series of posts each involving a separate Project Euler problem. This post involves solutions to Project Euler problem #92.

A number chain is created by continuously adding the square of the digits in a number to form a new number until it has been seen before.

For example,

44 → 32 → 13 → 10 → 1 → 1
85 → 89 → 145 → 42 → 20 → 4 → 16 → 37 → 58 → 89

Therefore any chain that arrives at 1 or 89 will become stuck in an endless loop. What is most amazing is that EVERY starting number will eventually arrive at 1 or 89.

How many starting numbers below ten million will arrive at 89?

First define some constants we’ll use in the various versions of the programs:

let NumDigits = 7
let limit numDigits = int(10.0**(float)numDigits) - 1 // limit 7 = 9999999 |

My first attempt utilizes this memoize function from a previous post.

let memoize (func:'a->'b) =
let cache = Dictionary<_, _>()
fun n ->
let (found, result) = cache.TryGetValue(n)
if found then
result
else
let result = func n
cache.[n] <- result
result |

My first attempt converted each number into a string, sliced that string into characters, and converted each individual character to a single-digit number and summed the squares of the digits to generate the next number in the sequence. Using that generator, it then determined whether each starting number terminated at 89 using a memoized function.

let digitize : string -> seq<int> = fun str -> seq {
for i = 0 to str.Length-1 do
yield int(str.[i]) - int('0')}
let generateNext1 =
let toString obj = obj.ToString()
let square x = x * x
toString >> digitize >> Seq.map square >> Seq.sum |> memoize
let rec terminatesAt1 : (int -> int) -> int -> int =
let f generateNext start =
match generateNext start with
| 1 -> 1
| 89 -> 89
| x -> terminatesAt1 generateNext x
memoize f
let numTerminatedAt1 terminatesAt terminal numDigits=
{1..limit numDigits}
|> Seq.map (terminatesAt >> (=) terminal)
|> Seq.countBy id
|> Seq.find (fun (bool, count) -> bool)
|> snd
let answer1 = lazy (numTerminatedAt1 (terminatesAt1 generateNext1) 89 NumDigits) |

Even being memoized, this performed pretty badly, taking 140 seconds to produce the answer. Within the `generateNext1`

function, every number gets converted into a string which produces a lot of garbage. There is a way to calculate the digits of a number without converting it into a string and slicing the sub-digits out of it:

let generateNext2 =
let square x = x * x
let rec digits x = seq {
let rem = x % 10
yield rem
let dividend = x / 10
if dividend > 0 then yield! digits dividend
}
digits >> Seq.map square >> Seq.sum |> memoize
let answer2 = lazy (numTerminatedAt1 (terminatesAt1 generateNext2) 89 NumDigits) |

Which completes in 65 seconds. However, 65 seconds is still very slow. We can further optimize that function by not generating a sequence for each and every number:

let generateNext3 x =
let mutable sum = 0
let mutable dividend = x
while dividend > 0 do
let rem = dividend % 10
sum <- sum + rem * rem
dividend <- dividend / 10
sum
let answer3 = lazy (numTerminatedAt1 (terminatesAt1 generateNext3) 89 NumDigits) |

`generateNext3`

performs much better coming in at 44 seconds.

Next we can focus on enumerating all 10 million starting numbers. Instead of memoizing, what about storing the results in an array of 10 million elements, and prepopulating element 1 and 89?

let numTerminatedAt2 terminal numDigits=
let array : bool option array = Array.create ((limit numDigits)+1) None
match terminal with
| 1 -> array.[1] <- Some true; array.[89] <- Some false
| 89 -> array.[1] <- Some false; array.[89] <- Some true
| _ -> failwith "1 and 89 are the only valid terminals."
let rec terminatesAt start =
match array.[start] with
| Some bool -> bool
| None ->
let result = generateNext3 start |> terminatesAt
array.[start] <- Some result
result
for i = 1 to limit numDigits do terminatesAt i |> ignore
let mapBoolOptionToSum x =
match x with
| None -> 0
| Some false -> 0
| Some true -> 1
Array.sumBy mapBoolOptionToSum array
let answer4 = lazy (numTerminatedAt2 89 NumDigits) |

That got us to 5.5 seconds. What if we reduced the garbage generated by replacing the Option array with an array of ints?

let numTerminatedAt3 terminal numDigits=
let array : int array = Array.create ((limit numDigits)+1) 0
array.[1] <- 1
array.[89] <- 89
let rec terminatesAt start =
match array.[start] with
| 0 ->
let result = generateNext3 start |> terminatesAt
array.[start] <- result
result
| x -> x
for i = 1 to limit numDigits do terminatesAt i |> ignore
Array.sumBy (fun x -> if x = terminal then 1 else 0) array
let answer5 = lazy (numTerminatedAt3 89 NumDigits) |

This clocks in at just under 2 seconds! One potential optimization upon this that could be made is to recognize that the biggest result that can be obtained from summing the squares of the digits for all numbers less than 10 million is for 9999999. 9999999 results in 567 (7*9^{2}). So rather than storing all 10 million terminating numbers, we only need to store the first 567.

let generateTerminalArray numDigits =
let limit = numDigits*9*9
let array : int array = Array.create (limit+1) 0
array.[1] <- 1
array.[89] <- 89
let rec terminatesAt start =
match array.[start] with
| 0 ->
let result = generateNext3 start |> terminatesAt
array.[start] <- result
result
| x -> x
for i = 1 to limit do terminatesAt i |> ignore
array
let numTerminatedAt4 terminal numDigits =
let terminals = generateTerminalArray numDigits
let mutable sum = 0
for i = 0 to limit numDigits do
if terminals.[generateNext3 i] = terminal then sum <- sum + 1
sum
let answer6 = lazy (numTerminatedAt4 89 NumDigits) |

This version took 1.5 seconds to compute. Another optimization is to recognize, that for a particular number, all permutations of that number generate the same chain. For example 12 and 21:

12 → 5 → 25 → 29 → 85 → 89
21 → 5 → 25 → 29 → 85 → 89

So if the chain for one particular permutation of a combination of numbers terminates in 89 then so will the chains for all the other permutations for that same combination of numbers. Thus, instead of brute-forcing the problem by calculating the result for each and every *permutation* of numbers, we can generate each *combination* and for every combination which terminates in 89, add to the sum the number of permutations for each combination.

let combinations (numDigits : int) : seq<int array> =
let rec enumerate (combo : int array) (start : int) (digit : int) = seq {
for i = start to 9 do
combo.[i] <- combo.[i] + 1
if digit = numDigits - 1
then yield combo
else yield! enumerate combo i (digit + 1)
combo.[i] <- combo.[i] - 1
}
enumerate (Array.create 10 0) 0 0 |

Then we need a way to generate a sequence of these arrays of digit counts – with each combination produced only once in the sequence.

To simplify the storage and calculation of the digits, each number is stored as an array of counts of each digit. I.e. The number 1223777 is stored as 0121000300 -> 0 zeros, 1 one, 2 twos, 1 three, etc. Furthermore, to reduce garbage and computation time, the combination sequence reuses the array it returns – modifying it slightly each time it is returned. `generateNext4`

is just a modified version of `generateNext3`

changed to operate on these arrays of digit counts.

let generateNext4 (combo : int array) =
let mutable sum = 0
for i = 0 to 9 do
let x = combo.[i]
sum <- sum + x*i*i
sum |

The number of permutations for each combination can be found via the multinomial coefficient:

(n_{1}+n_{2}+n_{3}+…)! |

n_{1}!*n_{2}!*n_{3}!*… |

where n_{1} is the number of 1′s, n_{2} is the number of 2′s, n_{3} is the number of 3′s, etc.

let factorials =
let factorial n = Seq.fold (*) 1 [1 .. n]
// store factorials in an array as a way of memoizing them
Array.init 12 factorial
let multinomialCoefficient (x : int array) : int =
let numerator = Array.sum x |> Array.get factorials
let denominator =
let mutable result = 1
for i = 0 to 9 do
result <- result * (Array.get factorials x.[i])
result
numerator / denominator |

Finally putting it all together:

let numTerminatedAt5 terminal numDigits =
let terminals = generateTerminalArray numDigits
combinations numDigits
|> Seq.filter (generateNext4 >> Array.get terminals >> (=) terminal)
|> Seq.sumBy multinomialCoefficient
let answer7 = lazy (numTerminatedAt5 89 NumDigits) |

This version clocks in at less than 50 milliseconds! But there’s an even better algorithm, mentioned by ceebo on the Project Euler website. Instead of enumerating through all of the combinations, simply enumerate through the array of 567 terminals, and for each one which terminates in 89, count the number of numbers which generate that number. I.e., for each number from 1 to 567 which results in 89, how many numbers, after running through `generateNext`

result in that number?

let `count(digits, n)`

represent the count of numbers with `digit`

digits that result in `n`

when run through `generateNext`

.

```
count(digits, n) = count(digits-1, n-0^2) + count(digits-1, n-1^2) + count(digits-1, n-2^2) + ...
```

With the following base cases:
```
count(0, 0) = 1
count(0, n) = 0 if n > 0
count(digits, n) = 0 if n < 0
```

let singledigits = [|0; 1; 2; 3; 4; 5; 6; 7; 8; 9|]
let rec count digits n =
match digits with
| 0 -> if n = 0 then 1 else 0
| _ ->
let summationTerm i =
let nextN = n - i*i
if nextN < 0 then 0 else count (digits-1) nextN
singledigits |> Array.sumBy summationTerm |

However, due to the recursive nature of this function, it’s much more efficient to memoize it. To avoid create lots of tuple garbage by storing the (digit, n) pairs in a dictionary, I instead memoized this function in a 2-dimensional array with this memoizing function:

let memoizeArray2D (x: int) (y : int) (f : 'a -> 'b -> 'c) : 'a -> 'b -> 'c =
let cache = Array2D.create x y None
fun a b ->
match cache.[a,b] with
| Some value -> value
| None ->
let value = f a b
cache.[a,b] <- Some value
value |

And the `count`

function is modified to use it:

let singledigits = [|0; 1; 2; 3; 4; 5; 6; 7; 8; 9|]
let rec count =
let compute digits n =
match digits with
| 0 -> if n = 0 then 1 else 0
| _ ->
let summationTerm i =
let nextN = n - i*i
if nextN < 0 then 0 else count (digits-1) nextN
singledigits |> Array.sumBy summationTerm
memoizeArray2D (NumDigits+1) (NumDigits*9*9+1) compute |

And finally putting it all together:

let numTerminatedAt6 terminal numDigits =
generateTerminalArray numDigits
|> Array.mapi (fun i term -> if term = terminal then count numDigits i else 0)
|> Array.sum
let answer8 = lazy (numTerminatedAt6 89 NumDigits) |

This clocks in under 10 ms…a 5x improvement over the previous solution. Furthermore, at higher digits, this final function scales better than the previous solution.