Home · Blog · Games · Resume

2016-10-29

- Exploring the space
- Profiling
- Parallelization
- Expression
- Generalization: elements and element pairs
- Conclusions
- Update (11/11/2016)

Here’s a math puzzle inspired by the election. I noticed a while back that the
two candidates were at 81.5% and 18.5%, respectively. **The numbers 185 and 815
are anagrams, and add to 1000 ($10^3$). How many such pairs of positive integers
add to $10^3$?** I’ve been doing a bit of Haskell recently, so that’s my
language of choice.

The first step is to define a function which takes an integer pair and decides if the pair is valid (that is, if the two integers are anagrams or not.)

```
import Data.List
valid :: Int -> Int -> Bool
valid a b = list a == list b
where list = sort . filter (/=0) . digits
digits n = map (\x -> read [x] :: Int) (show n)
```

Let’s use this `valid`

function to check a range of numbers `(a,b)`

where
$a+b=1000$, and `valid a b == True`

.

```
import Data.List
λ> [(a,b) | a <- [1..500], let b = 1000 - a, valid a b]
[(95,905),(185,815),(275,725),(365,635),(455,545),(500,500)]
```

This is pretty cool. Let’s fit this into a function $f(n)$, where

- $f(n)$ is the power $10^n$ we’re matching against,
`pairs n`

returns the list of potential pairs,`validPairs n = filter(\(a,b) -> valid a b) $ pairs n`

, and`f n = length $ validPairs n`

, which is what we’re ultimately trying to find an expression for.

```
pairs :: Int -> [(Int, Int)]
pairs pow = [(a, (10^pow)-a) | a <- candidates]
where candidates pow = [0,5..(div (10^pow) 2)]
validPairs = filter(\(a,b) -> valid a b) . pairs
f = length . validPairs
```

If we run `validPairs 4`

, we see a cool pattern.

```
λ> validPairs 4
[(995,99005), (1895,98105), (1985,98015), ... (49055,50945), (49505,50495), (50000,50000)]
λ> f 4
6
λ> f 5
141
λ> f 6
...
```

Before we go further, a few notes on this: for performance reasons, we only
compute pairs `(a,b)`

, not `(b,a)`

. Because the last pair computed is always
something like `(50,50)`

, `(500,500)`

, etc., the total number of pairs is
actually $2f(n) - 1$. We’ll keep this in mind for later.

Here’s what we have so far:

`n` |
1 | 2 | 3 | 4 | 5 | 6 |
---|---|---|---|---|---|---|

`f(n)` |
1 | 1 | 6 | 6 | 141 | 141 |

`2*f(n)-1` |
1 | 1 | 11 | 11 | 281 | 281 |

Computing $f(6)$ took a really long time, so let’s eventually wrap this in a
`.hs`

file, compile, and execute it with timing.

But before we go further, let’s check the Online Encyclopedia of Integer
Sequences to see if any of the sequences
`(1,1,6,6,141,141)`

, `(1,6,141)`

, or `(1,11,281)`

appear at all.
It turns out `(1,6,141)`

are the first few digits of
A241015, which is the

Number of pairs of endofunctions $f$, $g$ on $[n]$ satisfying $g(g(g(f(i)))) = f(i)$ for all $i$ in $[n]$.

…whatever that means.

If our pattern matches theirs, we might expect the next number in the sequence to be is 6184. We’ll find out in a bit, but first let’s bring runtime down a bit.

OK, it’s time to wrap this in a file, compile, and execute it with timing.

```
--puzzle.hs <num> (L | _)
import System.Environment
import Data.List
import Data.Char
valid :: Int -> Int -> Bool
valid a b = list a == list b
where list = sort . filter (/=0) . digits
digits n = map (\x -> read [x] :: Int) (show n)
pairs :: Int -> [(Int, Int)]
pairs pow = [(a, (10^pow)-a) | a <- candidates pow]
where candidates pow = [0,5..(div (10^pow) 2)]
validPairs = filter(\(a,b) -> valid a b) . pairs
f = length . validPairs
main :: IO ()
main = do
[num, action] <- getArgs
let n = digitToInt (head num)
let shouldList = (head action == 'L')
if shouldList
then print $ validPairs n
else print $ f n
```

Just for kicks, `./puzzle`

takes two arguments; the first is the number $n$, and
the second is either `L`

, to print the full list of matches, or `_`

, to print
the value $f(n)$.

We have a few ways to run this. One is to open up a `ghci`

shell, import the
file, and run `validPairs <num>`

or `f <num>`

right there in the shell:

`ghci`

```
$ ghci
GHCi, version 7.10.3: http://www.haskell.org/ghc/ :? for help
λ> :l puzzle
[1 of 1] Compiling Main ( puzzle.hs, interpreted )
Ok, modules loaded: Main.
λ> validPairs 3
[(95,905),(185,815),(275,725),(365,635),(455,545),(500,500)]
λ> f 3
6
```

This is great for testing, but we don’t get too many stats. Let’s use a
real-time thing like `runhaskell`

.

`runhaskell`

```
$ runhaskell puzzle.hs 3 L
[(95,905),(185,815),(275,725),(365,635),(455,545),(500,500)]
$ runhaskell puzzle.hs 3 _
6
```

To add stats, we use the `+RTS`

runtime system flag with the `-s`

summary flag.

```
$ runhaskell puzzle.hs 3 _ +RTS -s
6
101,696 bytes allocated in the heap
3,464 bytes copied during GC
68,912 bytes maximum residency (1 sample(s))
13,008 bytes maximum slop
1 MB total memory in use (0 MB lost due to fragmentation)
Tot time (elapsed) Avg pause Max pause
Gen 0 0 colls, 0 par 0.000s 0.000s 0.0000s 0.0000s
Gen 1 1 colls, 0 par 0.000s 0.000s 0.0001s 0.0001s
INIT time 0.000s ( 0.000s elapsed)
MUT time 0.001s ( 0.189s elapsed)
GC time 0.000s ( 0.000s elapsed)
EXIT time 0.000s ( 0.000s elapsed)
Total time 0.006s ( 0.189s elapsed)
%GC time 2.5% (0.1% elapsed)
Alloc rate 178,562,837 bytes per MUT second
Productivity 96.1% of total user, 2.9% of total elapsed
```

Very cool. As a final step, we’ll properly compile and optimize it with `ghc`

.

`ghc`

```
$ ghc -Odph puzzle.hs -rtsopts
[1 of 1] Compiling Main ( puzzle.hs, puzzle.o )
Linking puzzle ...
$ ./puzzle 3 _
6
```

We can run this with `+RTS -s`

too. `$ ./puzzle 3 _ +RTS -s`

generates the following table for `f(1)`

through `f(7)`

:

n | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
---|---|---|---|---|---|---|---|

`f(n)` |
1 | 1 | 6 | 6 | 141 | 141 | 5591 |

time | 0.001s | 0.001s | 0.003s | 0.017s | 0.190s | 2.301s | 26.711s |

Let’s see if this matches our sequence from earlier. We now have the next value, putting our sequence at (1,6,141,5591). This is not in the OEIS.. So we’ll have to find an expression ourselves, and check with the scripts.

Unfortunately, our timing statisics are not optimistic. At this rate, $f(8)$ will take around five minutes to calculate. It would be a good idea to parallelize.

We can use the `Control.Parallel`

library, (see the
docs
for more) which gives us the handy `par`

and `pseq`

functions, which force the
compiler to split evaluations into separate
threads and then wait to recombine. Much of the code has to be
modified slightly. I run this on a four-core processor, so I rewrote this to
split up the list of candidates into four roughly equal parts, evaluate the
validity of each element in each of them, compute the length of that valid
sublist, and recombine afterwards.

We’ll call this new file `puzzle-parallel.hs`

. You’ll notice we drop support for
printing the full list – constructing the sum was easier to parallelize.

```
-- ghc -O2 --make puzzle-parallel.hs -threaded -rtsopts
-- ./puzzle-parallel <n> +RTS -N4 -s
--
import Control.Parallel
import Data.List.Split
... -- all the same imports as before
valid :: Int -> Int -> Bool -- as before
pairs :: Int -> Int -> [(Int, Int)]
pairs pow core = [(x, (10^pow)-x) | x <- candidates]
where candidates = [a+5,a+10..b] -- not as before!
a = (core - 1) * (div (10^pow) 8)
b = (core - 0) * (div (10^pow) 8)
-- validPairs can no longer be point-free, since it takes
-- two arguments: pow and core. There's probably a nice currying fix, though.
validPairs pow core = filter(\(a,b) -> valid a b) $ pairs pow core
f n = s1 `par` s2 `par` s3 `par` s4 `pseq` (s1 + s2 + s3 + s4)
where
s1 = length $ validPairs n 1
s2 = length $ validPairs n 2
s3 = length $ validPairs n 3
s4 = length $ validPairs n 4
main :: IO ()
main = do
[num] <- getArgs
-- we drop support for `./puzzle <num> L`, since computing these
-- large lists, sorting, and printing them is sort of out of
-- scope for the moment.
let n = digitToInt (head num)
print $ f n
```

The file `puzzle-parallel.hs`

can now be compiled to be multithreaded with `ghc`

:

```
$ ghc -O2 --make puzzle-parallel.hs -threaded -rtsopts
```

and executed on four cores with the `+RTS`

flag `-N4`

(or `-N2`

, or however
many cores you have. The above code is written for four cores, but it’s not
hard to modify for any other number of cores.)

```
$ ./puzzle-parallel <n> #unparallelized
$ ./puzzle-parallel <n> +RTS -N4 #force 4 cores
$ ./puzzle-parallel <n> +RTS -N4 -s #with stats
```

n | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
---|---|---|---|---|---|---|---|---|

`f(n)` |
1 | 1 | 6 | 6 | 141 | 141 | 5591 | 5591 |

time w/o parallelization | 0.001s | 0.001s | 0.003s | 0.017s | 0.190s | 2.301s | 26.711s | 501.863s |

time w/ parallelization | 0.001s | 0.001s | 0.003s | 0.021s | 0.133s | 1.567s | 18.489s | 208.285s |

We see around a 2x speedup from this – we expect a 4x from threading, but there is some loss to overhead. This speedup helps more the higher $n$ grows. This is better performance across the board. Now calculating $f(9)$ isn’t quite as impossible.

Rather than delve deeper into optimization, let’s focus on trying to find a closed-form expression for this sequence. We do this first by visual inspection of the types of pairs generated.

```
$ ./puzzle 1 L
[(5,5)]
$ ./puzzle 2 L
[(50,50)]
$ ./puzzle 3 L
[(95,905),(185,815),(275,725),(365,635),(455,545),(500,500)]
$ ./puzzle 4 L
[(950,9050),(1850,8150),(2750,7250),(3650,6350),(4550,5450),(5000,5000)]
$ ./puzzle 5 L
[(995,99005),(1895,98105),(1985,98015),(2795,97205),(2975,97025),(3695,96305),
(3965,96035),(4595,95405),(4955,95045),(5495,94505),(5945,94055),(6395,93605),
(6935,93065),(7295,92705),(7925,92075),(8195,91805),(8915,91085),(9095,90905),
...
(45815,54185),(45905,54095),(46355,53645),(46535,53465),(47255,52745),
(47525,52475),(48155,51845),(48515,51485),(49055,50945),(49505,50495),
(50000,50000)]
```

Let’s introduce some nomenclature here to turn this from a programming puzzle
into a math puzzle. For some $n$, we define $f(n)$ to be the number of possible
pairs, and $S(n)$ to be the set of pairs itself. Thus, `f(n) = length S(n)`

.

Some observations:

- For all even values $n$, we can see
`S(n) = [10*i for i in S(n-1)]`

; that is, every element of`S(even n)`

is an element of`S(odd n)`

multiplied by ten. - Thus, $f(n) = f(n-1)$ for even n.

We can see from the examples above that for sets generated by, say, $n=3$ and $n=4$, only the first two digits really fluctuate. We suspect the problem reduces to a combinatorics problem focused on strings of that length, which we call $m$. So, let’s map the domain of possible $n$ values to a simpler domain:

```
m(n)
^
6 | o o
5 |
4 | o o
3 |
2 | o o
1 +-+-+-+-+-+-+-+--> n
1 2 3 4 5 6 7 8
```

We’ll switch to Python here for readability, and also because our algorithm will be imperative in a moment:

```
def M(n):
return ((n-1)//2)*2
```

This alongside future memoizations will help us avoid having to calculate $f(n)$ when $f(n-1)$ is known.

We first define the concept of a 9-pair. A 9-pair is a set of two positive
numbers which add to 9. 9-pairs are: `(0,9)`

, `(1,8)`

, `(2,7)`

, `(3,6)`

, and
`(4,5)`

.

Here is an example pair from $S(n=5)$: `(48155,51845)`

. We can reduce this to
its tuple of interest: `(4815, 5184)`

. It seems `(1,8)`

and `(4,5)`

are
the 9-pairs which compose the alphabet from which this tuple was generated.

Since `1`

and `8`

are a 9-pair, if `1`

appears in the first number of the
tuple, then we know for sure that:

`1`

appears in the second item of the tuple,`8`

appears in the first item of the tuple,`8`

appears in the second item of the tuple.

Let’s generalize this away from digits and number and into elements, sets, and permutations.

We will call `(a,b)`

a tuple, `a`

and `b`

strings in that tuple, and the
characters in the strings `a`

and `b`

characters.

For a given string `a`

or `b`

, any character in `a`

or `b`

uniquely defines
its counterpart within an element-pair (`i`

,`j`

): if `i in a`

, then `j in b`

. Since the
strings `a`

and `b`

are by definition anagrams, this means `j in a`

as well.
Thus `a`

is a unique ordering of pairs of characters, selected from a known
alphabet of character pairs.

We can reduce this problem to finding the number of unique orderings of a list composed of selecting $x$ element-pairs from a list of $p$ element-pairs. We call this function $G(x,p)$. Because our real problem lies on the domain of digits, there will only eve be five element-pairs (the five 9-pairs), so $p=5$, always.

This is simple (if inefficient) in Python:

```
def G(x,p):
result = 0
# since p=5 always, this domain is always ['1','2','3','4','5']
domain = (str(i) for i in range(p))
itr = itertools.combinations_with_replacement(domain,x)
for i0 in itr:
# turns ['1','2'] into ['1','1~','2','2~'] for proper element-pairs
alphabet = list(i0) + list( map((lambda x: str(x)+"~"), i0) )
#the number of unique permutations in that alphabet of length 2x
result += len(set(list(itertools.permutations(alphabet,2*x))))
return result
```

To get our desired $f(n)$, we can just write:

```
def f(n):
# the same as M(n) above
m = ((n-1)//2)*2
# the answer is recursive -- G(x,5) alone will only return pairs of length x.
# we must also consider pairs of length y<x, right-padded by zeroes.
def inner_f(n):
return (G(n,5) + inner_f(n-1)) if (n > 0) else 0
# divided by two because (185, 815), (815, 185) are the same pair.
# plus one because (500, 500) won't be found by G5(x)
return inner_f(m)/2 + 1
```

`t0`

We can run this with timing:

```
for n in range(1,12):
start = time.time()
localF = f(n)
end = time.time()
print "f({0})\t{1}\t{2}".format(n, localF, end-start)
```

n | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
---|---|---|---|---|---|---|---|---|---|---|

`f(n)` |
1 | 1 | 6 | 6 | 141 | 141 | 5591 | 5591 | 281566 | 281566 |

`t0` |
5.9e-6s | 1.3e-5s | 4.2e-5s | 2.2e-5s | 1.5e-4s | 1.4e-4s | 4.3e-3s | 4.0e-3s | 0.94s | 0.94s |

`t1`

We can improve performance a lot through memoization. Functions like $G(n,5)$ get called repeatedly for the same value of $n$, for example.

```
G5memo = {}
def G5(x):
if x not in G5memo:
G5memo[x] = G(x,5)
return G5memo[x]
```

n | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
---|---|---|---|---|---|---|---|---|---|---|

`t1` |
9.1e-6s | 3.1e-6s | 5.5e-5s | 3.1e-6s | 1.6e-4s | 2.9e-6s | 5.5e-3s | 4.1e-6s | 0.91s | 1.6e-5s |

The most obvious boost is that $f(n)$ where $n$ is even now takes almost no time to compute, since $f(n) = f(n-1)$ for even $n$.

`t2`

We can memoize further. The `alphabet`

of characters fed into
`len(list(set(itr(alphabet))))`

is usually something like `['a','A','a','A']`

or `['b','B','b','B']`

. However, it seems obvious to us that the number of
unique pairs from the first of these two alphabets will be the same as that of
the second. It would be useful to create a non-unique *hash* of some kind,
representing the pair distribution within an alphabet, and check whether
`len(list(set(itr(alphabet))))`

for an equivalent alphabet has already been
calculated.

```
# takes an alphabet like ['1', '1', '2', '3', '1~', '1~', '2~', '3~'] and returns
# \ / | |
# \ / | |
# (2) (1) (1) --> returns [1,1,2] as a list-hash
def makeListHash(alphabet):
d = {}
for i in alphabet[0:len(alphabet)/2]:
try:
d[i] += 1 #if i in d
except:
d[i] = 1
return list(str(x) for x in sorted(d.values()))
# just a wrapper for makeListHash which returns a string ("1,2,2") instead.
def makeStringHash(alphabet):
return ",".join(makeListHash(alphabet))
```

This lets us rewrite $G(x,p)$ as:

```
def G(x,p):
result = 0
alphabet = (str(i) for i in range(p))
itr = itertools.combinations_with_replacement(alphabet,x)
lMemo = {} #inline since memoization is only useful across a single value of x
for i0 in itr:
alphabet = list(i0) + list( map((lambda x: str(x)+"~"), i0) )
key = makeStringHash(alphabet)
if not key in lMemo:
lMemo[key] = len(set(list(itertools.permutations(alphabet,2*x))))
result += lMemo[key]
return result
```

This memoization brings down computation time even further, and lets us calculate $f(11)=f(12)$ for the first time ever.

$n$ | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 |
---|---|---|---|---|---|---|---|---|---|---|---|

$f(n)$ | 1 | 1 | 6 | 6 | 141 | 141 | 5591 | 5591 | 281566 | 281566 | 16397596 |

`t2` |
6.9e-6 | 3.1e-6 | 6.8e-5 | 1.3e-5 | 1.4e-4 | 3.1e-6 | 6.9e-4 | 5.0e-6 | 6.8e-2 | 7.2e-6 | 10.5 |

That’s all! I wasn’t able to further optimize the problem, and I’m not convinced a closed-form solution necessarily exists. That said, we found the number of anagram pairs which sum to $100,000,000,000$, and that’s pretty cool. I also learned a lot about parallelization in Haskell and using the built-in optimization flags in `ghc`

, and a bit of memoization in Python.

It’s been around a week since last I touched this puzzle. My buddy Josh bounced the problem to his boss Piotr, who came up with this phenomenal closed-form solution for $G(x,p=5)$. It’s brilliant and involves reducing the combinatorics problem, searching for and finding the elusive closed-form solution in the literature, and writing a small haskell script to implement it. In short:

The combinatorics section above describes a

sequence $a_n^N := \Sigma \left( \dfrac{n!}{p_1! p_2! … p_N!}\right)^2$, where the sum runs over sets of non-negative integers $p_1,..p_N$ summing to $n$.

This sequence can be found by a recurrence relation, and was performed in
*Sums of squares of binomial coefficients, with applications to Picard-Fuchs
equations* by H. A. Verrill in
2008.

*Table 1* in the above paper describes recurrence relations for $a_n^N$, and
we care about the case where $N=5$: the fourth equation is the object of our
desire. It reads:
\(0 = n^4 a_n^5 - (35n^4 - 70n^2 + 63n^2 - 28n + 5)a^5_{n-1}\)
\(\phantom{0 = } + (n-1)^2(259(n-1)^2 + 26)a^5_{n-2} - (3 \cdot 5)^2(n-1)^2 a^5_{n-3}\)

Piotr implemented this as so: (I’ve changed it the tiniest bit to better perform testing on it.)

```
import System.Environment
import Data.List
import Data.Char
next :: Integer -> Integer -> Integer -> Integer -> Integer
next n prev1 prev2 prev3 =
((35*n*n*n*n - 70*n*n*n + 63*n*n - 28*n + 5) * prev1 -
((n-1)*(n-1)*(259*(n-1)*(n-1) + 26)) * prev2 +
15*15*(n-1)*(n-1)*(n-2)*(n-2) * prev3)
`div` (n*n*n*n)
vals :: [Integer]
vals = 1 : 5 : 45 : (zipWith4 next [3..] (drop 2 vals) (drop 1 vals) vals)
choose :: Integer -> Integer -> Integer
choose n 0 = 1
choose 0 k = 0
choose n k = choose (n-1) (k-1) * n `div` k
finalvals = zipWith (\v n -> v * (choose (n*2) n)) vals [0..]
main :: IO ()
main = do
[num] <- getArgs
let n = read num :: Int
print $ (finalvals !! n)
```

and called with `runhaskell math.hs <num> +RTS -s`

(assuming you save it as `math.hs`

.)

As expected, this has a phenomenal runtime.

$n$ | 1 | 10 | 1e2 | 1e3 | 1e4 | 2e4 | 3e4 | 4e4 | 5e4 |
---|---|---|---|---|---|---|---|---|---|

runtime | 0.183s | 0.183s | 0.180s | 0.198s | 0.542s | 1.232s | 2.032s | 3.203s | 4.643s |

I think that’s about as far as this problem goes. Thanks to Josh and Piotr for doing the actual difficult mathematics.

Email · Twitter · GitHub · LinkedIn · Instagram · RSS