the jsomers.net blog.

[This post turned into a kind of narrative about solving a particular problem, which I'm happy with, but if you're interested in my general remarks about the online encyclopedia of integer sequences (the original intention of this post) jump here.]

Over the weekend I had a chance to work on a few Project Euler problems. I think of them as little "mental workouts," in that they provide a somewhat frictionless (controlled) environment in which to practice math and hacking. And like physical workouts, PE problems are satisfying in direct proportion to their difficulty.

Problem 106, reproduced here, has been one of my all-time favorites (of the 111 I've solved):

Let S(A) represent the sum of elements in set A of size n. We shall call it a special sum set if for any two non-empty disjoint subsets, B and C, the following properties are true:

1. S(B) ≠ S(C); that is, sums of subsets cannot be equal.
2. If B contains more elements than C then S(B) > S(C).

For this problem we shall assume that a given set contains n strictly increasing elements and it already satisfies the second rule.

Surprisingly, out of the 25 possible subset pairs that can be obtained from a set for which n = 4, only 1 of these pairs need to be tested for equality (first rule). Similarly, when n = 7, only 70 out of the 966 subset pairs need to be tested.

For n = 12, how many of the 261625 subset pairs that can be obtained need to be tested for equality?

(If you're interested in the rest of this post, I encourage you to read the description again, at least until you have some idea of what it's about; it's a bit of a mouthful the first time through. And if you're interested in solving the problem, by all means do it -- but be sure to come back when you've finished.)

The first thing that struck me about this question is that it doesn't feel well-defined. What does it mean for us to "need" to test a pair? Isn't that a judgment call?

My first step, then, was to print out the 25 subset pairs where n = 4; the idea would be to go through them and cross out the ones that "obviously" couldn't be equal, and see if there was a pattern to it.

I immediately discovered one rule: pairs of singletons can't possibly have equal "sums." Since the sets are disjoint by assumption, any member of one can't be part of the other -- which, when your sets only have one member each, implies a strict inequality. Simple.

The next rule, too, follows directly from a fact mentioned in the problem description, namely, that property #2 is satisfied. So we can throw out any pair containing sets of different sizes, because these can't be equal (by assumption).

Now it turns out that when you account for those two rules, you're only left with three pairs -- out of which only one "needs" to be tested, according to the problem description. They are:

[1, 2] [3, 4]
[2, 3] [1, 4]
[1, 3] [2, 4]

Can you see the rule?

Well, in the first and third cases, set A is strictly dominated by set B: if you move through both sets simultaneously, proceeding from left to right, the number you choose from set B is in every case greater than its partner from set A. So set B unequivocally has the larger sum.

With that in mind, it looks like we have all of the rules for telling us exactly how many pairs we need to check. But before we jump all the way up to n = 12, it makes sense to code up a little script on n = 7 to make sure our rules eliminate all but 70 of the 966 possible subset pairs (as we're told in the problem description). Here's a little code that does this (using a fast powerset() function pasted to the Python mailing list):

s = [1, 2, 3, 4, 5, 6, 7]
subsets = [list(A) for A in list(powerset(s))]

ct = 0
# Enumerate all possible pairs, then filter:
for i, s in enumerate(subsets[1:]):
    for j, t in enumerate(subsets[i:]):
        if len((s + t)) == len(set(s + t)): # Disjoint
            if len(s) == len(t) and len(s) > 1: # Rules 1 and 2
                diffs = [a[1] - a[0] for a in zip(s, t)]
                if len(filter(lambda x: x <= 0, diffs)) > 0: # Rule 3
                    ct += 1
print ct

It prints 70, which verifies that our rules are sound.

Unfortunately, though, we can't just add five elements to s and be done with it. When n = 12, the power set of s contains 2¹² elements, which means that if we "enumerated all possible pairs," as in the script above, we'd run through our loop (2¹²)² times, which would take something like 9 hours to complete on a computer like mine.

So the problem comes down to our being able to find (a) the number of pairs P_x of disjoint, equally-sized subsets for sizes x = 2, ..., 6 (or 12/2), and (b) subtract from each of these P_x the number of pairs that are "strictly dominated," in the sense above.

It took me a long time to figure out how to even approach (a). What I finally did was to make a table, for n = 4 and x = 2, of all the possible subset pairs; I'd then cross out the ones that weren't disjoint and try to find a pattern. So I started to draw something like this:

         (1, 2) (1, 3) ... (6, 7)
        +------+------+---+------
(1, 2)  |
(1, 3)  |
 ...    |
(6, 7)  |

which was going to be filled with X's in every spot where two sets contained the same element. But I stopped way early, because in the process of filling out my first row, I realized that the pair (1, 2) had essentially "blotted out" exactly two possibilities (1 and 2) from my list, leaving 3, 4, 5, 6, and 7. Which meant that if I was going to ensure "disjointness," all I had to do was figure out the number of ways to select 2 elements out of the 5 remaining -- which is just "5 choose 2", or, 5!/2!(5-2)!.

To figure out my total P₂, then, I just had to multiply by the number of rows, which turns out to be "7 choose 2", or as a Python function, nCr(7, 2).

I checked my work using the lists I had printed out earlier, and it turned out that I was always off by a factor of two: I was double-counting. In terms of the table, I had forgotten that you only had to look above (or below) the diagonal.

With that done, all that was left was to subtract out the number of strictly dominated pairs.

But this was a lot more difficult than it sounds.

To restate the problem concretely, let's say N = 7 and x = 2. What we're after is the number of ways of selecting two pairs of elements from our 7-element list such that each of the guys from pair A has a partner in pair B that is strictly greater than it. Visually, where our pair A is represented by stars below the numbers and pair B by lines above, what we want is something like this:

      _   _
1 2 3 4 5 6 7
  *     *

as opposed to this:

    _   _
1 2 3 4 5 6 7
      *   *

because in this second case, 6 (from pair A) isn't dominated by anyone from pair B.

My first stab at a systematic rule here was to try fixing two stars, like so:

1 2 3 4 5 6 7
  *   *

and to see how many ways there were to arrange two lines so that they strictly dominated the stars. Working with a pencil, I started drawing, erasing, and tallying.

This was a brutal mess. There were simply too many possibilities, and I wasn't at all clear about what would happen if I increased n or x; I couldn't find a system, and was even having trouble moving the lines systematically without losing track. Arranging the stars differently just muddled things up even more.

I went to dinner and tried to think a bit while eating. One idea I had was to write, above each number, the different star-numbers that would be dominated if I chose that number (i.e., if I drew a line above it) -- sort of like how people solve Sudoku puzzles by listing all the possibilities and propagating constraints.

I don't think that was a good idea, but it may have got me thinking, because when I got back to work to try it out I came up with something much better. It was the clearest "insight" of the day.

What I realized was that n doesn't matter, at least in one important sense. Say x = 2 and n = 1,000,000. For any arrangement of stars and lines, there are all kinds of ways to shift things around; probably millions of them. But for the most part, none of these moves will change the essential dominance relationship -- the ordering of lines and stars relative to one another. There are really only a fixed number of arrangements that matter, and increasing n just adds empty space.

In fact, we can systematically investigate every arrangement-that-matters by letting n = 4, which is small enough to work out by hand. Then, and here's the trick, anytime we increase n, all we need to think about is the number of ways of choosing 4 numbers from that set.

Let's say, just to be explicit about it, that there are 2 ways of arranging two stars and two lines around the list 1 2 3 4 such that the lines dominate the stars. (There are.) Then, if we let n = 7, or 12, or 2,000, all we need to do is calculate the number of ways of embedding our 4-member list in there -- given by nCr(n, 4) -- and multiply that number by 2. Done.

We do have a bit more work to do, though, because all we know precisely at this point is what happens when x = 2, and we need to consider the cases where x = 3, 4, 5, and 6, because those are all the possibilities when n = 12. (Singletons aren't allowed, and anything bigger than 6 runs out of space.)

I was able to calculate these values explicitly for x = 2, 3, 4, and roughly for x = 4 (I wasn't sure about the last one, because it required lots of drawing, erasing, and tallying). So I had a series that looked like 2:2, 3:5, 4:14, 5:~43. I could have tried to figure it out for x = 6, but I saw that the series was growing quickly and I wasn't too keen on making hundreds of tallies. But I also wasn't sure that I could figure out a formula; I assumed it would be some tricky recursive combinatorial expression, and frankly, I had a much simpler idea:

* * *

I went to the online encyclopedia of integer sequences, otherwise called "OEIS" or "Sloane's."

What it is is a repository of integer sequences (obviously) to which anyone can contribute. For each entry, there is

• A list of the first twenty or so terms, along with links to more. There's even a link to an audio file where you can hear the terms spoken out loud.
• Comments -- usually one or many "interpretations" of the list, assuming it has some number-theoretic, algebraic, or geometric significance.
• Heavy-duty (i.e. peer-reviewed) references of papers that use the sequence in some way; occasionally this includes the paper that "discovered" it.
• Links to online explanations or discussions.
• Formulas and computer code (usually Maple and Mathematica) for producing it.
• Cross references to related sequences.

The way I think about it is less as an encyclopedia of sequences per se than it is a kind of projection of mathematics onto sequence-space. The point being that an integer sequence, because it often has so many mathematical interpretations, acts in some sense like the intersection of those interpretations. There are of course many ways in which two or more mathematical concepts are linked, but a shared integer sequence is among the most useful, precisely because it can be browsed and searched like an encyclopedia.

It also helps that this encyclopedia is online -- and therefore hyperlinked -- and that its contributors and editors are mostly professional mathematicians or computer scientists. It makes for an incredible resource, one more modest than Wikipedia or IMDb but in some ways more interesting.

The chief difference, I think, is that one can more readily imagine someone using the OEIS to do original research. It's as much a rich set of raw data as it is a compendium of what's already known.

Of course it's possible to be misled by a particular entry, which is why one still has to do the math before taking a connection too seriously -- for instance, I bet that nearly any random sequence of five or six strictly increasing numbers less than 100 will show up on Sloane's; so it might not be prudent to, say, walk into an elevator, see which numbers are highlighted, look them up online, and assume you've found some connection between your building and the Stolarsky array read by antidiagonals. But it's a rich resource for those who know what to do with it.

* * *

So that's where I went to find the next member of my series (x = 6). Since I wasn't sure of my value for 5 (43), I just omitted it and searched for "2, 5, 14", which brought up the "Catalan numbers" as the first entry.

As you may gather from the length of that page, the Catalan sequence is quite important and admits to a whole slew of interpretations. One of these in particular jumps out to me as straightforwardly relevant to the problem at hand:

Ways of joining 2n points on a circle to form n nonintersecting chords.

All you have to do is think of "chords crossing" as "a violation of strict dominance," and the correspondence should be clear. In my mind it's a beautiful visualization of the idea.

In any case, it turns out that the 5th member of our series is going to be 132. Which means we're ready to solve the problem, because we now have a formula for determining the number of strictly dominated pairs given an array length n and subset pairs each of size x. It's simply the number of ways of embedding 2x elements in n, i.e., nCr(n, 2x), multiplied by what turns out to be the xth Catalan number for that embedded array of 2x points.

Setting n = 12 and ranging over the possible values for x, 2 through 6, gives our answer. Here's the code:

def fact(n):
    "factorial"
    return reduce(lambda x, y: x * y, range(1, n + 1))

def nCr(n, r):
    "n choose r"
    if n == r:
        return 1
    else:
        return fact(n) / (fact(r) * fact(n - r))

def x_with_x(n, x):
    "number of *disjoint* x-subset pairs one can form with an array of length n"
    if n - x >= x:
        return (nCr(n, x) * nCr(n - x, x)) / 2
    else:
        return 0

catalans = {2:2, 3:5, 4:14, 5:42, 6:132}

def dominated(n, x):
    "number of 'strictly dominated' x_with_x guys -- these need not be checked"
    return nCr(n, 2 * x) * catalans[x]

def need_to_check(n):
    "number of disjoint subset pairs we need to check for array length n"
    tot = 0
    for x in range(2, n // 2 + 1):
        tot += x_with_x(n, x) - dominated(n, x)
    return tot

print need_to_check(12)

It takes .27 seconds to run on my machine.