Solving the Pairs of Percentages Puzzle

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.

Exploring the space⌗

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:

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.

Profiling⌗

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:

with 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.

with 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.

with 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):

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.

Parallelization⌗

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

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.

Expression⌗

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.

Mapping to a simpler domain⌗

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.

The 9-pair⌗

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).

Numbers and 9-pairs⌗

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.

Generalization: elements and element pairs⌗

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.

Combinatorics on that domain⌗

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

Unmemoized timing 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)

Memoization 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]

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$.

Memoization 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.

Conclusions⌗

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.

Update (11/11/2016)⌗

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} + (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.

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