Three postings here (1, 2, 3) have discussed what happens when you form all pairwise sums and differences from a finite set of integers. The number of differences almost always exceeds the number of sums—a fact that lends special interest to the occasional exceptional sets, with More Sums Than Differences (MSTD).

A question left up in the air when I last wrote on this subject was the size of the smallest MSTD set. The smallest known set was {0, 2, 3, 4, 7, 11, 12, 14}, which has 26 sums but only 25 differences. This set has eight elements. Could there be an MSTD set with seven or fewer elements? The question has now been answered in the negative by Peter V. Hegarty of the Chalmers University of Technology and Göteborg University. He proves there is no smaller set, and furthermore that {0, 2, 3, 4, 7, 11, 12, 14} is the only MSTD set of size eight. (Apart from other eight-element sets generated from {0, 2, 3, 4, 7, 11, 12, 14} by affine transformations.)

Read all about it at the arXiv.

I’ve received the following note from Barry Cipra, bit-player’s Bureau Chief in Northfield, Minnesota, (where all hail broke loose yesterday):

Your latest postings [here and here] have motivated me to idle away some time with a variant of the problem(s) you discussed. To wit, I got to wondering what happens if you restrict yourself to taking sums of distinct pairs and positive differences of numbers from a (finite) set S — i.e., x+y and x–y with x > y in S.

It was easy, after a while, to see that the sum set is usually larger than the difference set. In particular, for the set {0, 1, 2, …, m}, the sum set is {1, 2, …, 2m–1}, while the difference set is {1, 2, …, m–1}. It seems, or at least seemed, obvious that the sum set is never smaller than the difference set, but a proof keeps slipping through my (fumbling) fingers, which suggests that either it isn’t true, or the proof is complicated, or the slow decline of my mathematical abilities has steepened and accelerated. Intuitively, each duplicate sum is associated with two duplicate differences: If a+d = b+c with a < b < c < d, then d–b = c–a and d–c = b–a. The problem, of course, is that you can have duplicate duplicates.

It seemed worth looking at how the sum/difference ratio tends toward 2 as one “grows” a set S by randomly adding numbers between 1 and some large upper bound, say m = 10,000. I wrote, ran, debugged, ran, debugged,…., and finally ran successfully a little (Basic) program to do this, and tried it with various values of m, printing out the ratio as each new (random) element got added. The number of elements required to reach any given ratio r (less than 2) seems to be proportional to sqrt(m). In particular, the number of elements to reach r=1.5 seems to be 2*sqrt(m). It would be most gratifying if this turned out to be asymptotically exact.

Barry also points out a connection with the birthday “paradox.” Here on Earth, any gathering of 23 or more people has a better-than-even chance of including at least two people with the same birthday. On Mars, where the year has 687 days, you need to assemble 32 people to have a fifty-fifty chance of such a coincidence. On Saturn, where people can have any of 1,055 birthdays, the minimum size of the group is 39. (I would add that on Pluto, with a year of 14,158 Plutonian days, the minimum likely size for a double-birthday party is 141 people. However, there can’t be that many people on Pluto, since it’s only a dwarf planet.) Anyway, the point of all this is that the number of people needed for a birthday coincidence grows in proportion to the square root of the number of possible birthdays. The calculation is similar when it comes to the number of elements in a set drawn at random from {0, 1, 2, …, m}, except that the probability we’re interested in is not that of finding two identical elements but rather that of having two pairs of elements with the same sum or difference.

Barry’s version of the problem has a satisfying symmetry. For both sums and differences, we’re looking at the upper triangle of a matrix, excluding the diagonal and everything below it. If there are no duplicates, or coincidences, then the sum set and the difference set are necessarily of the same size.

Inspired by Barry’s example, I have idled away a few of my own declining hours, and gone off on a bit of a tangent.

In an arbitrary finite set of non-negative integers, suppose there are n elements, and the largest element is m. If m is much larger than n, coincidences are unlikely among both the sums and the differences; in the limiting case, they have probability zero. If m = n–1, then the set must be the arithmetic progression {0, 1, 2, …, m}, and the number of coincidences is maximized. Barry’s algorithm for adding randomly chosen numbers to a set, one at a time, can be viewed as a way of exploring the transition between these two extremes. We can ask, at each stage of the process: What is the probability that a randomly chosen number, distinct from all the existing elements of the set, can be adjoined to the set without generating a coincidence?

Here is a somewhat different question. Instead of choosing numbers at random to add to the set, suppose we always insert the smallest number that can be adjoined to the set without creating a coincidence. In other words, we build the set from the bottom up. What does the resulting set look like?

First consider just the constraints imposed by the sum set, ignoring coincidences among differences. Starting from the empty set {}, we can obviously go on to {0} and then to {0, 1}. Next we can insert 2, because the three sums 0+1 = 1, 0+2 = 2 and 1+2 = 3 are all distinct. But we cannot add 3 to the set, because that would introduce the coincidence 0+3 = 1+2. The next permissible element is 4; the set {0, 1, 2, 4} produces the sum set {1, 2, 3, 4, 5, 6} without coincidences. Going on from there, both 5 and 6 are excluded, but 7 is okay. Working with pencil and paper, I extended this set as far as {0, 1, 2, 4, 7, 12, 20, 29}. This gave me enough to go consult the Oracle of AT&T. I found the sequence listed as No. A010067 in Neil Sloane’s Encyclopedia of Integer Sequences. The sequence is described as follows: “A B₂ sequence: a(n) = least value such that sequence increases and pairwise sums of distinct elements are all distinct.” This was just what I had been looking for, and the success boosted my confidence that at least I had done the arithmetic right.

Now I tried the same method with the differences. Again the first few sets that yield no coincidences are {], {0} and {0, 1}, but in the subtractive case {0, 1, 2} fails, because 2–1 = 1–0. On the other hand, {0, 1, 3} is without conflicts. The series continues {0, 1, 3, 7, 12, 20, 30}. It was time to visit the Oracle again. This sequence is No. A025582, but the identifying tag line is not quite what I was expecting: “A B₂ sequence: a(n) = least value such that sequence increases and pairwise sums of elements are all distinct.” I had generated the series by looking at differences of distinct elements, but it turns out that sums of not-necessarily-distinct elements yield the same series. Curious.

Let’s consider the element-by-element differences between these two sequences, which I’m going to call S(n) for the sums of distinct pairs and D(n) for the differences of distinct pairs:

   n        D(n)       S(n)      D(n)-S(n)
   0          0          0            0
   1          1          1            0
   2          3          2            1
   3          7          4            3
   4         12          7            5
   5         20         12            8
   6         30         20           10
   7         44         29           15
   8         65         38           27
   9         80         52           28
  10         96         73           23
  11        122         94           28
  12        147        127           20
  13        181        151           30
  14        203        181           22
  15        251        211           40
  16        289        257           32
  17        360        315           45
  18        400        373           27
  19        474        412           62
  20        564        475           89

From the first few entries in this table one might well guess that the sums will always be lagging behind the differences, so that the D(n) – S(n) sequence grows monotonically. But then there’s a reversal at n = 10, and after that it’s not clear exactly what the trend will be. The obvious question is: Can S(n) ever be greater than D(n)? This time the Oracle can’t help us; a Sequence Server search for {0, 0, 1, 3, 5, 8, 10, 15, 27, 28} comes up empty.

I tried producing a few more terms of S(n) and D(n) and their differences. Struggling with a grotesquely inefficient program, I invested an hour of cpu time, which got me as far as n = 60:

   n        D(n)       S(n)      D(n)-S(n)
   0          0          0            0
   1          1          1            0
   2          3          2            1
   3          7          4            3
   4         12          7            5
   5         20         12            8
   6         30         20           10
   7         44         29           15
   8         65         38           27
   9         80         52           28
  10         96         73           23
  11        122         94           28
  12        147        127           20
  13        181        151           30
  14        203        181           22
  15        251        211           40
  16        289        257           32
  17        360        315           45
  18        400        373           27
  19        474        412           62
  20        564        475           89
  21        592        530           62
  22        661        545          116
  23        774        607          167
  24        821        716          105
  25        915        797          118
  26        969        861          108
  27       1015        964           51
  28       1158       1059           99
  29       1311       1160          151
  30       1394       1306           88
  31       1522       1385          137
  32       1571       1434          137
  33       1820       1555          265
  34       1895       1721          174
  35       2028       1833          195
  36       2253       1933          320
  37       2378       2057          321
  38       2509       2260          249
  39       2779       2496          283
  40       2924       2698          226
  41       3154       2873          281
  42       3353       3060          293
  43       3590       3196          394
  44       3796       3331          465
  45       3997       3628          369
  46       4296       3711          585
  47       4432       3867          565
  48       4778       4139          639
  49       4850       4446          404
  50       5122       4639          483
  51       5242       5021          221
  52       5297       5064          233
  53       5750       5322          428
  54       5997       5613          384
  55       6373       6003          370
  56       6800       6273          527
  57       6924       6493          431
  58       7459       6641          818
  59       7546       6979          567
  60       7788       7275          513

Would you care to place a bet on the further behavior of this series? At n = 60, D(n) – S(n) looks like it’s still growing, on the average, but it’s also still fluctuating. Will any of the dips in the curve ever plunge deep enough to bring the balance below zero? Or will the overall upward trend prevail in the long run? Maybe it’s like π(x) – Li(x) in number theory. This function—the difference between the actual number of primes up to x and an estimate proposed by C. F. Gauss—is negative for all values of x anyone has ever calculated explicitly. Nevertheless, there’s a proof that π(x) – Li(x) does eventually become positive; indeed, it changes sign infinitely often. But the first such zero crossing is somewhere out in the remotest exurbs of the number line.

I debated: Is it worth the bother of writing a better program to add a few more lines to the table? If D(n) – S(n) is positive for all values of n up to n = 60, what is the probability of discovering a negative value at n = 61, or even n = 600? I have no idea how to answer that question mathematically, but psychologically the prospects seemed pretty dim. Nevertheless, I wrote the program. Here are the next four entries in the table:

   n        D(n)       S(n)      D(n)-S(n)
  61       8219       7587          632
  62       8502       8017          485
  63       8729       8373          356
  64       8941       9071         -130

There it is: S(64) is greater than D(64). And I soon learned that this excursion into negative territory is not a one-of-a-kind freak occurrence. In the interval between n = 0 and n = 500, D(n) – S(n) is positive for 301 values, negative for 198 values and 0 for 2 values. It seems a reasonable conjecture that positive and negative values occur equally often. Here’s a graph of D(n) – S(n) for n = 0 to n = 250:

Final note: Although I have just today blundered down this interesting path, I am assuredly not the first to pass this way. I gather that I need to do some reading on Sidon sequences. Are there also connections with Golomb rulers?

Kevin O’Bryant, whose work on sets that have more sums than differences was mentioned in this recent post, writes as follows:

Here’s a related problem that Mel Nathanson and myself (with Ruzsa and a few students) have also been thinking about. For a finite set A of integers, let S = {a + b : a,b \in A} be the sumset, and let T = {2a + b : a,b \in A}. Can it happen that S is larger than T? The answer is yes (we have a proof) but we don’t have a single example. Our proof gives a set that would have at least e²⁰⁰⁰ elements, so although our proof is “constructive”, it really isn’t in any practical sense. My professors thought “constructive” referred to a proof that didn’t use the axiom of choice, my generation thinks of it as an algorithm that runs in polynomial time, my students will have to worry about which polynomial. :)

Note that this problem, like the original MSTD problem, has the pleasant property of affine invariance: You can scale or translate the elements of the set by any constant, and the size of the S and T sets remains unchanged. In the case of an arithmetic progression, such as the n-element set {0, 1, 2,…, n–1}, the sumset S has 2n–1 elements, whereas the set T has 3n–2 elements. Thus |T| is greater than |S| by about n. For a random sampling of sets other than arithmetic progressions, it appears that |T| usually exceeds |S| by an even larger margin. I have no idea where to look for that elusive example with |T| < |S|.

Take a set of integers, say {0, 2, 5, 8, 11}, and write down all the numbers that can be represented as sums of two elements drawn from this set. For our example the answer is {0, 2, 4, 5, 7, 8, 10, 11, 13, 16, 19, 22}. Now construct the corresponding set of pairwise differences: {–11, –9, –8, –6, –5, –3, –2, 0, 2, 3, 5, 6, 8, 9, 11}. Note that there are only 12 distinct sums but 15 differences.

Let’s try it with another set, {5, 8, 17, 26, 41}:

sums: {10, 13, 16, 22, 25, 31, 34, 43, 46, 49, 52, 58, 67, 82}
diffs: {–36, –33, –24, –21, –18, –15, –12, –9, –3, 0, 3, 9, 12, 15, 18, 21, 24, 33, 36}

Again the differences outnumber the sums, this time by a margin of 19 to 14.

It’s not hard to see why differences tend to be more numerous: Addition is commutative but subtraction isn’t. Thus the sums 5+8 and 8+5 both yield the single result 13, whereas 5–8 and 8–5 produce two distinct differences, –3 and +3. It is rumored that someone other than John Horton Conway once conjectured that the number of sums formed in this way never exceeds the number of differences. But the conjecture is false. A counterexample is the set {0, 2, 3, 4, 7, 11, 12, 14}, which has 26 distinct pairwise sums but only 25 differences:

sums: {0, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 21, 22, 23, 24, 25, 26, 28}
diffs: {–14, –12, –11, –10, –9, –8, –7, –5, –4, –3, –2, –1, 0, 1, 2, 3, 4, 5, 7, 8, 9, 10, 11, 12, 14}

Melvyn B. Nathanson of Lehman College—the Bronx campus of the City University of New York—has lately called attention to such anomalous sets of integers, which he identifies by the abbreviation MSTD (more sums than differences). He has lots of questions. Why do such sets exist? Where are they found? How many are there? What is their structure?

This past April Nathanson discussed MSTD sets in a talk titled “Problems in Additive Number Theory” at the University of Montreal; the talk is available on the arXiv as math.NT/0604340. In June Nathanson delivered a follow-up talk, “Sets with More Sums than Differences,” at the SIAM Conference on Discrete Mathematics in Victoria, British Columbia; that paper was released last week on the arXiv as math.NT/0608148. Meanwhile Kevin O’Bryant of the College of Staten Island (another CUNY unit) has addressed somewhat different aspects of the MSTD problem in a paper titled “Many Sets Have More Sums than Differences” (math.NT/0608131).

The appeal of a problem like this one is that it seems to get a lot of mileage out of the simplest mathematics: adding, subtracting and counting—operations that most of us know how to do. As the papers of Nathanson and Bryant show, the math is not all so trivial, and yet an amateur like me can still hope to have some fun with these questions. I’ve been toying with MSTD sets for the past week or so.

First, a few preliminaries. A set, as defined for this discussion, is a collection of items without duplicates. For example, {1, 3} is a two-element set. There are four ways to add these elements in pairs—1+1, 1+3, 3+1 and 3+3—but two of those summations yield the same result, and so the “sumset” has just three elements, {1, 4, 6}. The order of the elements in a set has no significance, but for convenience I’ll always list them in ascending sequence. All the sets discussed here are finite.

For a clearer understanding of set sums and differences, it helps to write down an example in matrix format:

A set of n elements has n² pairwise sums and differences, but they are not all distinct. In the case of the sums, the matrix is symmetric, and so everything above the main diagonal is duplicated below it. Thus the maximum possible number of unique sums is given by counting the entries along the diagonal plus those in either the upper or the lower triangle, but not both. This number is equal to n(n+1)/2; for the example given here, where n=4, the maximum size of the sumset is 10. However, for the specific set shown here the maximum is not attained, because of a few “coincidences”: 4 appears as both 0+4 and 2+2, and 6 arises as both 2+4 and 3+3. Thus the sumset has only eight elements: {0, 2, 3, 4, 5, 6, 7, 8}.

The difference matrix is antisymmetric, and so the elements of both the upper and the lower triangles need to be counted. On the other hand, the diagonal of the difference matrix is all zeros. The maximum number of distinct differences is n(n–1)+1. Again, though, the maximum is not reached in this example; coincidences reduce the size of the “diffset” from 13 to 9. Still, the differences outnumber the sums, and so {0, 2, 3, 4} is not an MSTD set.

The smallest possible sumset or diffset has 2n–1 elements. (You might want to work out why.) It’s easy to construct a set that attains this minimum: Just choose elements in an arithmetic progression. For example, the set {0, 2, 4, 6} has the sumset {0, 2, 4, 6, 8, 10, 12} and the diffset {–6, –4, –2, 0, 2, 4, 6}, both of size 7. It’s also straightforward to build a set that generates the largest possible sumset and diffset; the trick is to make each element more than twice as large as the next smaller element, as in the set {0, 1, 3, 7}. This structure eliminates all coincidental duplicates in both the sums and the differences.

In the search for MSTD sets we don’t have to look at all possible sets of integers. It turns out that both the number of sums and the number of differences generated by a set remain unchanged if you add a constant to each member of the set. Likewise, multiplying each element by a constant also leaves the number of sums and differences invariant. In other words, you can transform each element x into ax+b (an affine transformation) without altering the size of the sumset or the diffset. This property is important because it means we can represent any MSTD set in a canonical form. We can shift it along the number line until its smallest element is 0, and we can shrink it down to its smallest possible span of integers by dividing out any factors that are common to all the nonzero elements. For example, the set {5, 8, 17, 26, 41} mentioned above has the canonical form {0, 1, 4, 7, 12}. Both of these sets have 19 differences and 14 sums.

Now for some questions.

What is the smallest MSTD set? The smallest known set is the example {0, 2, 3, 4, 7, 11, 12, 14} that I have already introduced. It has eight elements, and in the canonical representation the largest element is 14. There is one other known eight-element example, {0, 2, 3, 7, 10, 11, 12, 14}. A few seconds of computing is all it takes to show that there is no smaller eight-element MSTD set. But is there an example with fewer than eight elements? I don’t think so, but I can’t prove it. A brute-force search rules out any MSTD set with seven or fewer elements where the largest element is less than 81. Imre J. Ruzsa of the Mathematical Institute of the Hungarian Academy of Sciences claims that any MSTD set must have at least seven elements.

How many MSTD sets are there? In one sense, this is a very easy question. Given the affine invariance of MSTD sets, if just one such set exists, then we can generate infinitely many of them by translation and dilation. But most people would agree that these are all just copies of the same set in disguise. What we really want to know is the number of MSTD sets when they are all reduced to canonical form. Nathanson has shown that in this scheme of reckoning, too, the number of sets is infinite. He gives a formula for generating infinite families of MSTD sets. Starting with the example {0, 2, 3, 4, 7, 11, 12, 14}, the formula yields a sequence of progressively larger sets that Nathanson proves must all have more sums than differences: {0, 2, 3, 4, 7, 11, 15, 16, 18}, then {0, 2, 3, 4, 7, 11, 15, 19, 20, 22}, then {0, 2, 3, 4, 7, 11, 15, 19, 23, 24, 26}, and so on. The question then arises, are all MSTD sets members of such infinite families, or are there also “sporadic” MSTD sets?

How rare are MSTD sets? Having already established that there are infinitely many MSTD sets, it might seem that they can’t be very rare, but that’s not necessarily true. The primes are also infinite, yet they are vanishingly rare. Take the ratio of the number of primes less than N to the number of integers less than N; as N goes to infinity, the ratio goes to zero. MSTD sets could be rare in a similar sense. O’Bryant shows that within a certain infinite series of integer sets, the probability of finding an MSTD set is greater than zero. Does that result hold also for integer sets in general?

By how much can the number of sums exceed the number of differences? Let’s define the discrepancy Δ of a set as the number of differences minus the number of sums. For all “ordinary” sets, Δ ≥ 0. For MSTD sets, Δ is negative. In all the MSTD examples I’ve shown so far, Δ = –1, or in other words the number of sums is just 1 greater than the number of differences. I’ve been able to find lots of sets with Δ = –2; the smallest, with 11 elements, is {0, 1, 2, 4, 5, 9, 12, 13, 14, 16, 17}, which has 35 sums and only 33 differences. I’ve also stumbled upon a few sets with Δ = –3, such as the 16-element set {0, 1, 2, 4, 5, 9, 12, 13, 14, 16, 17, 21, 24, 25, 26, 28}, which has 56 sums and 53 differences. I’m guessing there is no lower bound on Δ. (Update: As I was smoothing out the last details of this report, I discovered a 1973 paper by Sherman K. Stein. He proves that the ratio of the number of differences to the number of sums can be made arbitrarily large or small.)

How are MSTD sets distributed among all integer sets? The trouble with even asking a question like this is that we can’t look at all integer sets. If we try to answer the question statistically, by choosing a representative sample of sets, then we have to wade into the messy business of deciding what sort of sample is representative. For lack of a better idea, I have tried the following approach. Assuming that all sets are in canonical form, I classify them by two parameters: n, the number of elements, and m, the largest element. Then for any given values of n and m I can either examine all sets (if n and m are small enough) or generate a random sample. Note that m cannot be less than n–1. If m = n–1, then there is only one possible set, namely the counting sequence 0, 1, 2,…, m. As m increases for a fixed value of n, so does the number of possible sets and the average gap between the elements. Now we can ask how Δ varies as a function of m and n. In the case of m = n–1, the answer can be given unequivocally: Δ = 0, because the set is an arithmetic progression, and both the number of sums and the number of differences is 2n–1. When m is much greater than n, Δ is almost surely positive. The reason is that the elements of the set are widely dispersed, and a coincidence in which different pairs of elements yield the same sum or the same difference is unlikely. In almost all sets with m >> n, the number of sums takes its maximum value n(n+1)/2 and the number of differences is also at its maximum, n(n–1)+1. If you do the subtraction, you’ll find that Δ increases in proportion to n².

The interesting region would appear to be the middle ground, where m is not too much larger than n. Here is a graph showing the frequency of Δ values for n = 10 and all values of m between 9 and 27.

As m increases, the distribution grows wider and shifts to the right—toward more-positive values of Δ. (Note that the graph is based on a complete enumeration of all sets, not on a statistical sampling.)

MSTD sets are so rare that in the graph above their frequency is indistinguishable from zero. Below we look exclusively at the frequency of MSTD sets as a function of m for three values of n.

It appears that MSTD sets are most common (or maybe one should say least rare) at the smallest value of m where such sets first appear. But this impression is somewhat misleading. As the graph below reveals, the absolute number of MSTD sets increases as a function of m (or, in the case of n = 8, remains constant). Although it’s true that the proportion of MSTD sets declines, that’s only because the total number of integer sets grows exponentially with m.

Thus if you are wandering aimlessly in the universe of integer sets, the number of sets with more sums than differences increases with m, but the probability that a randomly encountered member of the population has the MSTD property falls to zero as m increases.

Finally, where did this problem come from? The locus classicus appears to be an unpublished 1967 edition of a list of unsolved problems compiled by Hallard T. Croft of the University of Cambridge. Other authors, citing the Croft list, attribute the conjecture that differences always outnumber sums to John Horton Conway. Nathanson writes, however: “I asked Conway about this at the Logic Conference in Memory of Stanley Tennenbaum at the CUNY Graduate Center on April 7, 2006. He said that he had actually found a counterexample to the conjecture, and that this is recorded in unpublished notes of Croft.” The mention of Croft refers to the same 1967 list that others cite in attributing the conjecture to Conway. I have not seen this document; if anyone can send me a copy, I would be most grateful. Here are a few more slightly less obscure references:

Marica, John. 1969. On a conjecture of Conway. Canadian Mathematical Bulletin 12:233–234.

Ruzsa, Imre Z. 1984. Sets of sums and differences. In Séminaire de Théorie des Nombres, Paris, 1982–83, pp. 267–273. Boston: Birkhäuser.

Ruzsa, I. Z. 1992. On the number of sums and differences. Acta Mathematica Hungarica 59:439–447.

Stein, Sherman K. 1973. The cardinalities of A+A and A–A. Canadian Mathematical Bulletin 16:343–345.

At the library the other day I was perusing the Collected Rantings and Ravings of Edsger W. Dijkstra [Note 1]. While leafing through the pages in search of something else [Note 2], I stumbled across “An Exercise for Dr. R. M. Burstall” [Note 3], a brief essay in the form of an open letter, written in 1976. The “exercise” involves a little recursive function of the natural numbers, which Dijkstra named fusc. I’m a collector of such trifles, so I immediately perked up. The fusc function, though 30 years old, was new to me. Here’s how Dijkstra defined it [Note 4]:

fusc(1) = 1
fusc(2n) = fusc(n)
fusc(2n+1) = fusc(n) + fusc(n+1)

And here’s the Lisp version [Note 5] I quickly knocked together:

(defun fusc (n)
  (cond ((= n 1) 1)
        ((evenp n) (fusc (/ n 2)))
        ((oddp n) (+ (fusc (/ (- n 1) 2))
                     (fusc (/ (+ n 1) 2))))))

In the letter to Burstall, Dijkstra had this to say about fusc:

(The function fusc is of mild interest on account of the following property: with f1 = fusc(n1) and f2 = fusc(n2) the following two statements hold for n1 ≠ n2: “if there exists an N such that n1 + n2 = 2^N, then f1 and f2 are relatively prime” and “if f1 and f2 are relatively prime, then there exist an n1, an n2, and an N, such that n1 + n2 = 2^N”. In the above recursive definition, this is no longer obvious, at least not to me; hence its name.)

It looked pretty fuscular to me too. Neither the definition nor Dijkstra’s explanation gave me an immediate sense of what the function was meant to calculate, so I ran the program and tabulated a few argument-value pairs:

n       =  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
fusc(n) =  1 1 2 1 3 2 3 1 4  3  5  2  5  3  4  1  5  4  7  3

The first pattern I was able to spot here was in the values of fusc(n) for n = 2, 4, 8, 16…. Whenever n is a power of 2, fusc(n) = 1. It’s easy to see why: When n = 2^k, the function merely invokes itself on successive values of n/2 until eventually n = 1, at which point the recursion terminates. It’s also clear that only for n = 1 and n = 2^k can fusc have the value 1. I noted a similar pattern in the intervals between values of n for which fusc(n) = 2, namely the values n = 3, 6, 12, 24…. Interesting, but not all that illuminating. And the distribution of 3s in the sequence looks less regular, with values of 3 at n = 5, 7, 10, 14, 20….

In another attempt to figure out what was going on, I sketched the tree of recursive function calls the program must make when it is started on a specific value of n:

The general pattern seemed to make sense; in particular, each odd value of n gave rise to two recursive calls (on (n–1)/2 and (n+1)/2), whereas each even n had only one descendant in the call tree, on n/2. And I could see that the eight 1s forming the leafy fringe of this tree are summed up to produce the returned value fusc(21) = 8. But there were still plenty of questions. For example, I could not predict the next term in the fusc series.

If I had persisted, perhaps I would have worked all this out on my own, but persist I did not. What I did was succumb to temptation. I cheated. I looked up the series in Neil J. A. Sloane’s Encyclopedia of Integer Sequences. In 3.547 seconds flat I had a detailed report listing 91 terms of the series and discussing its properties, applications and history. The history, it turns out, extends back long before Dijkstra [Note 6].

Sloane’s EIS also revealed some remarkably conspicuous patterns in the series that I had somehow missed entirely. Consider the subsequence extracted by starting at fusc(2) and continuing with every second value: 1, 1, 2, 1, 3, 2…. It’s the original sequence all over again! This also works for every fourth element starting at fusc(4), for every eighth element starting at fusc(8), and so on. The sequence has a nested, self-similar structure:

1 1 2 1 3 2 3 1 4 3 5 2 5 3 4 1 5 4 7 3 8 5 7 2 7
  1   1   2   1   3   2   3   1   4   3   5   2
      1       1       2       1       3       2
              1               1               2

Another distinctive pattern involved consecutive pairs of elements from the series. If you haven’t already discovered it for yourself, you might want to try again before reading on. (Assuming that you, too, haven’t already looked everything up in the EIS.)

Here’s the pattern: Take each pair of adjacent numbers and form a fraction: 1/1, 1/2, 2/1, 1/3, 3/2, 2/3, 3/1, 1/4, 4/3, 3/5, 5/2, 2/5, 5/3, 3/4, 4/1…. If you keep going in this way, you will eventually generate every rational number exactly once, all of them reduced to lowest terms [Note 7].

The EIS identifies the output of the fusc function as “Stern’s diatomic sequence” and refers to an 1858 paper by M. A. Stern [Note 8]. Seeing that reference rang a bell for me. I recognized both the name of the author and the title of the paper. And several of the other references also looked familiar. And then, as I scanned farther down the list, I had the disconcerting experience of finding one of my own articles mentioned. Not that I’m complaining about this. It’s always a fine thing to have your work cited—to see your name up in lights—but it’s a little worrisome, all the same, when you’re struggling to understand a new idea and you’re referred back to yourself for the answer.

What I knew of Stern was the Stern-Brocot tree [Note 9], which also generates all the rationals, but by a somewhat different construction. Or at least it appears different. The basic idea is to take any two fractions and create a new intermediate fraction, the mediant, by adding the numerators and denominators. For example, 1/2 and 2/3 form a mediant of 3/5. Applying this one rule systematically produces an infinite tree of rational numbers:

All the ratios from the fusc sequence appear here, row by row, although their order within each row is none too obvious. Even less obvious (to me) is the connection between the fusc recursion and the algorithm for building the Stern-Brocot tree. The fusc process seems to be focused on the distinction between even and odd, whereas the Stern-Brocot algorithm is all about a weird kind of averaging. Yet it can’t be coincidence that both operations yield the same set of numbers.

In an effort to better understand what’s happening in the fusc procedure, I created an “instrumented” version of the program. The fusc recursion can terminate only when n is reduced to 1, and the value returned by the function is a count of how many times this has happened. I wanted to see how many times the program triggers the other two clauses—those for cases where n is even and those where n is odd and greater than 1. Here are some results:

n      1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
ones   1  1  2  1  3  2  3  1  4  3  5  2  5  3  4  1  5  4  7  3
evens  0  1  1  2  2  2  3  3  4  3  4  3  5  4  6  4  7  5  7  4
odds   0  0  1  0  2  1  2  0  3  2  4  1  4  2  3  0  4  3  6  2

The row labeled ones records the number of times the (= n 1) clause is executed; the evens row give the number of times the (evenp n) predicate is satisfied, and odds counts the number of times the (oddp n) predicate in the third clause fires [Note 10]. As expected, the numbers in the ones row correspond exactly to the values of fusc(n). Also note that whenever n is an exact power of 2, the entry in the evens row is that power, and the entry in the odds row is 0.

I want to call attention to still another pattern. Observe that odds(n) never exceeds evens(n)—at least within the range of n shown here—and that the two quantities are equal only in a few cases, namely n = 1, 3, 5, 11, 21. Is there something special about these numbers? Here are some more values of n for which odds(n) = evens(n)—in fact, all such values less than 1 million:

1 3 5 11 21 43 85 171 341 683 1365 2731 5461 10923 21845 43691
87381 174763 349525 699051

Before spoiling all our fun by turning to the magic of the EIS, consider what these numbers signify. In a certain sense, they are the oddest of all numbers: When decomposed by the fusc recursion, each of these numbers runs through the odd branch of the circuit more often than any other numbers of similar size. But why do these particular numbers, and no others, have that property? I’m not sure I know the answer to that question, but I have a way of thinking about it that I find amusing. Suppose I ask you to construct a series of natural numbers with the following properties: Every number is odd, and every number is double the preceding number in the series. This task is obviously impossible: No integer that is twice another integer can be odd. But the series of “oddest” numbers listed above represents a kind of best-effort attempt to comply with the rules. All the numbers in the series are certifiably odd, and the ratio of successive terms converges on 2 in the limit of large n.

There’s more to be said. Each of these “oddest” numbers lies between two successive powers of 2. Let us index the oddest numbers as M_k, with M₀=1, M₁=3, and so on. Then for any k, 2^k < M_k < 2^k+1. Furthermore, as k increases, the ratio of M_k to 2^k+1 approaches 2/3, and the ratio of M_k to 2^k approaches 4/3. Why those ratios? Well, where would you expect the “oddest” numbers to be found—at the halfway point between successive powers of 2? Of course not: That’s the line of even division.

The series I have been calling the oddest numbers is sequence A001045 in the EIS, where it is identified as the Jacobsthal sequence [Note 11]. It turns out the Jacobsthal numbers are close kin of the Fibonacci numbers; they are generated by the recursion J_n = J_n–1 + 2J_n–2, with J₀=0 and J₁=1. And they turn up in quite a variety of places. For example J_n gives “the number of ways to tie a necktie using n+2 turns.” Closer to home, a note from Amarnath Murthy points out that the sums of consecutive terms of the Jacobsthal series yield the powers of 2 in order. (I hadn’t noticed.)

In recent weeks I have been chastised (in the friendliest way!) for writing blog entries that leave no room for readers to contribute thoughts of their own. Well, there’s no shortage of loose ends in this story. Here are a few questions for which I would like to know answers.

1. As noted above, the series of n for which fusc(n)=1 begins 1, 2, 4, 8, 16, 32, 64, and the corresponding series for fusc(n)=2 is 3, 6, 12, 24, 48, 96, 192. But the distribution of other fusc values is not so easy to describe. The positions where fusc(n)=3 begin 5, 7, 10, 14, 20, 28, 40, 56, 80; this series is identified in the EIS as a set of numbers whose “binary expansion is 1x100…0 where x = 0 or 1.” Does that property explain why the numbers appear in this context? For fusc(n)=4, the series of n is 9, 15, 18, 30, 36, 60, 72, which are numbers whose “binary expansion is 1xx100…0 where xx = 00 or 11.” Huh?

2. In the tabulation of evens and odds above, I have a million reasons to believe that the number of odds never exceeds the number of evens, but I don’t have a proof. Am I missing something obvious?

3. The action of the fusc function depends critically on the distinction between odd and even, or in other words on residues modulo 2. We can write analogous functions that classify integers modulo 3, 4, 5, and so on. For example, here’s fus3 in Lisp:

(defun fus3 (n)
  (cond ((= n 1) 1)
        ((= n 2) 2)
        ((= (mod n 3) 0) (fus3 (/ n 3)))
        ((= (mod n 3) 1)
         (+ (fus3 (/ (- n 1) 3))
            (fus3 (/ (+ n 2) 3))))
        ((= (mod n 3) 2)
         (+ (fus3 (/ (- n 2) 3))
            (fus3 (/ (+ n 1) 3))))))

And here’s a generic function, fusk [Note 12], that performs the equivalent calculation modulo k. (Interestingly, it’s simpler than the mod-3 version.)

(defun fusk (n k)
  (cond ((< n k) n)
        ((= (mod n k) 0) (fusk (/ n k) k))
        (t (let ((r (mod n k)))
             (+ (fusk (/ (- n r) k) k)
                (fusk (/ (+ n (- k r)) k) k))))))

Here’s some output from a few of these variant functions (fus2 is of course the standard fusc):

n      1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
fus2   1  1  2  1  3  2  3  1  4  3  5  2  5  3  4  1  5  4  7  3
fus3   1  2  1  3  3  2  3  3  1  4  4  3  6  6  3  5  5  2  5  5
fus4   1  2  3  1  3  3  3  2  5  5  5  3  4  4  4  1  4  4  4  3
fus5   1  2  3  4  1  3  3  3  3  2  5  5  5  5  3  7  7  7  7  4

You won’t find any of the series other than fus2 in the EIS. There seems to be a pattern of repetition emerging in the higher-order series, but what is it exactly, and what does it mean?

Notes:

Note 1: No, that’s not really the title of the book. It’s Selected Writings on Computing: A Personal Perspective, Springer-Verlag, 1982. But my title gives a clearer impression of the contents. Dijkstra was a major figure in the early development of computer science, and by the time of his death in 2002 he was a senior statesman of the field. Nevertheless, his preferred role was always that of the small boy who tells the truth about the emperor’s nakedness. All of the material from the book, and much else as well, is now available online at the E. W. Dijkstra Archive at the University of Texas at Austin.

Note 2: I was looking for “How Do We Tell Truths that Might Hurt?,” a 1975 memo that serves as a fine exemplar of Dijkstra at his most vituperative. I was looking it up because the version published a few years later in SIGPLAN Notices is slightly bowdlerized, and I wanted to see what the fig leaves covered. The answer was not very titillating; the unmentionable word in Dijkstra’s text turned out to be “IBM.” Dijkstra’s original version is here; the SIGPLAN Notices redaction (fee or subscription needed) is here.

Note 3: Burstall did not take up the challenge.

Note 4: Please allow me some ranting and raving of my own. Am I the only one who chafes at this inside-out style of function definition? The statement fusc(2n+1) = fusc(n) + fusc(n+1) describes a relation, but it’s not an algorithm for calculating anything. When we are given “2n+1″ as a parameter of a function, we still have to figure out what n must be before we can start computing. This layer of obfuscation (to borrow a term) is an invitation to error and misinterpretation. The present instance is a notably treacherous case, since the sum mentioned in the definition, fusc(n) + fusc(n+1), happens to be close to the right answer, but not quite on the mark.

Note 5: For the benefit of the illisperate, here’s the same algorithm in pseudo-Pascal:

function fusc(n): integer;
begin
    if n=1 then fusc := 1;
    elseif even(n) then fusc := fusc(n/2);
    elseif odd(n) then fusc := fusc((n-1)/2) + fusc((n+1)/2)
end;

The (odd n) Lisp predicate and odd(n) Pascal predicate are not strictly needed, since any natural number that’s not even had better be odd; but I’ve included the predicates to make the logic of the procedure a little clearer.

Note 6: A few pages further along in the Rantings and Ravings we come to “More About the Function ‘fusc’ (A Sequel to EWD570),” where we learn that Dijkstra too eventually looked it up in the Encyclopedia of Integer Sequences, which in those days was a book.

Note 7: Don’t take my word for it. The proof is given by Neil Calkin and Herbert S. Wilf in “Recounting the Rationals,” 2000, American Mathematical Monthly 107:360–363.

Note 8: Stern, M. A. 1858. Ueber eine zahlentheoretische Funktion. Journal für die Reine und angewandte Mathematik, 55:193–220. Moritz Abraham Stern was Gauss’s first doctoral student, and eventually succeeded him in the chair of mathematics at Göttingen. In the opening paragraphs of the paper, Stern gives credit for the basic idea to Gotthold Eisenstein (a fact pointed out to me by Donald E. Knuth).

Note 9: Brocot was Achille Brocot, a French clockmaker, whose interest in this number-theoretical function was highly pragmatic: He used it to help calculate gear ratios.

Note 10: I tried both the evens and the odds series at the Encyclopedia of Integer Sequences. No matches.

Note 11: Question: Who was (is) Jacobsthal? The EIS listing offers no clue, and the three or four references I followed also fail to identify Jacobsthal or cite any works by an author of that name. On the other hand, I did come across the following tidbit in Eric W. Weisstein’s MathWorld:

Amazingly, when interpreted in binary, the Jacobsthal numbers J(n+2) give the nth iteration of applying the rule 28 cellular automaton to initial conditions consisting of a single black cell (E. W. Weisstein, Apr. 12, 2006).

Note 12: I realize I am straying perilously close to the sort of thing that will get me blacklisted by parental-control filters.

Note n: You’d think that footnotes would be a breeze in HTML. What could be more hypertexty? Would it were so.

I am the child of a bookkeeper, and I’ve inherited the habit of double-checking receipts and balancing accounts. My friends make fun of me when I carefully note down the dime that I put in a parking meter, but lately I’ve been fretting over even smaller sums—charges that come to less than a penny. It’s all about sales tax.

For the benefit of those who live outside the U.S. (or in New Hampshire) I should explain how American sales tax works. The price marked on an item in a store excludes the tax, which is added at the cash register when you pay for your purchases. Where I live, in North Carolina, the tax rate for most merchandise is 7 percent. Thus if I buy a magazine with a cover price of $3.95, I’ll actually pay $3.95 + (0.07 × $3.95), or $4.2265. Fractions of a cent are rounded to the nearest penny, making the final price $4.23. The rounding is what I’ve been worrying about.

If rounding is done fairly, one might expect that the tax would be rounded up and rounded down with roughly equal frequency, and everything would come out even in the end. But for a long time I’ve been noticing that the sales tax on my purchases is almost always rounded up; very seldom, it seems, does my shopping basket produce a total for which the tax amount needs to be rounded down. Is this bias real, or do I just imagine that numbers are conspiring against me? As I was gathering up papers for a different kind of tax—the annual income-tax ritual—I decided to find out.

For any tax rate that’s an integer percentage, all possible rounding situations are covered by considering prices modulo $1. In other words, if the tax on $0.xy rounds in a certain way, then $1.xy will get exactly the same treatment, and likewise $2.xy, $3.xy, and so on. Thus for purposes of tax rounding, we can pretend that all prices lie between $0.00 and $0.99. In the sales-tax tables (PDF) issued by the state of North Carolina, the tax on 50 of these amounts is rounded up; in 49 cases it is rounded down; the one remaining case needs no rounding in either direction. This rounding protocol introduces a very slight upward bias, amounting to $0.00005 when averaged over all possible prices. Apart from this tiny asymmetry, the system looks like it ought to be fair if purchase prices are uniformly distributed among the 100 possibilities.

The problem, of course, is that the distribution of prices is far from uniform. I took a fat envelope full of sales receipts from 2005 and recorded the prices of all items subject to the 7 percent North Carolina tax. I was able to identify 471 items. Here is the distribution of their prices, modulo $1:

Red bars are prices for which the tax is rounded up, blue bars represent prices whose tax rounds down; the tax on $0.00 (black bar) is exact. The source of the rounding bias that I’ve suspected is immediately obvious: Almost two-thirds of the items I purchased over the course of the year had a price (modulo $1) of 99 cents, which happens to be an amount for which the tax rounds upward. I’ve always heard that the popularity of 99-cent pricing has a psychological premise: Shaving that penny off makes things seem cheaper and therefore encourages sales. Could there be another motive as well?

Before I start accusing shopkeepers of a nefarious plot to mulct American consumers, there’s a complicating factor to consider. Sales tax is not calculated separately on each item purchased; instead, when you dump your basketful of goods at the checkout counter, the prices are totalled, and the tax is calculated only once, on the sum. (In accounting lingo, the calculation is done “per invoice” not “per item.”) It turns out that my 471 taxable purchases were grouped into 195 invoices. (On a typical trip to the store I bought 2.4 items.) When I reanalyze my receipts on this basis, the 99-cent effect is somewhat softened:

Nevertheless, the invoice totals are still strongly clustered in the region between $0.93 and $0.99, where tax amounts are all rounded upward.

With these data in hand I can now calculate the total “roundage,” which I’m going to define as N(T_exact â€“ T_rounded), where N is the number of transactions at a given invoice amount. With this convention, the roundage is positive when it favors the merchant and negative when it favors me.

An interesting statistical question is whether this pattern offers any evidence of a deliberate policy of setting prices in order to maximize roundage for the merchant. Obviously the prevalence of prices just under a dollar has this effect, but since prices in that range have other possible justifications, they can’t support any firm conclusion. Are the sharp upward spikes at 79 cents and 50 cents a result of careful strategizing, or are they just random accidents? I have not tried to answer that question, and I don’t think I have enough data to do so. It should be noted that if you exclude the region of the graph between $0.93 and $0.99, the net roundage of the remaining transactions is slightly in my favor.

How much money are we talking about here? For calendar year 2005, based on the 195 transactions I was able to document, the skewed rounding of sales tax cost me 13.43 cents. This deficit is probably not the biggest hole in my pocket. On the other hand, if my case is typical, and if we extrapolate my loss to the entire U.S. population, that comes to $40 or $50 million slipping through the cracks.

Where is the insidious excess of roundage going to? That’s the question. If merchants are required to remit to the state the exact amount of tax actually collected, then the bias in roundage merely raises the effective tax rate a tad. No big deal. As far as I can tell, however, that’s not what happens—at least not here in North Carolina. The form on which sales tax is reported asks for total receipts (exclusive of tax collected) and then multiplies this amount by the appropriate percentage. In this way all of a merchant’s sales over an entire month or quarter are lumped together before any rounding is done, and any bias in the roundage of individual transactions will not be taken into account. The form does include a line for “excess collections,” and so a conscientious merchant could keep track of the roundage and report it there. But the instructions that accompany the form do not mention this possibility, and my supposition is that the excess roundage usually stays with the retailer.

We can fight back, though. Using my sample of 471 item prices, I ran a Monte Carlo simulation to estimate the average roundage as a function of the number of items purchased on each shopping trip:

The model assumes that I buy exactly k items, selected at random (with replacement) from the set of 471 prices in the 2005 data set. The sales-tax roundage on this shopping basket is calculated, and then is multiplied by 471/k, to give the total expected roundage for the year. The graph shows averages of 10⁶ repetitions of this process for each value of k from 1 through 25. The lesson is clear: If I bought 9 or 10 items every time I went to the store, I would come out slightly ahead in the roundage game.

However, I have a better solution to propose, one that could make the system more resistant to manipulation by either the seller or the buyer. Any integer tax rate brings the curse of commensurability: a pattern of roundage that repeats every dollar, rewarding strategies such as setting all prices to $x.99. Commerce might be more interesting with an irrational tax rate. For example, if the tax percentage were not 7 but rather the square root of 50 (equal to 7.0710678118654755…), then it would be a little harder to rig prices so that the tax consistently rounds upward. Besides, it would help with a problem discussed in my February 20th post: Pundits could no longer claim that algebra has no role in ordinary life. Every time you buy a cup of coffee, you’d have to solve the quadratic equation x² â€“ 50 = 0.

Further notes and questions.

The Monte Carlo model discussed above chooses k items at random. If you could plan a year’s worth of shopping in advance, you could list all N items you intend to buy and search for the best possible partitioning of these goods into k-item batches. How hard a computational problem is that? The brute-force approach (trying all possible partitionings) is exponential in N, but some such problems turn out to be easier than they look. What if you remove the constraint that each batch must have exactly k items?

Playing the game from the point of view of the vendor, what pricing policy will maximize roundage gains, given a buyer who is determined (at all costs!) to minimize them? Under a fixed, integer tax rate, is there any set of prices that guarantees a win for the merchant, no matter how the shopper assembles purchases into market baskets? Setting all prices to $x.00 enforces a trivial tie. Can the seller do better than that?

Instead of an irrational tax rate, we might try giving the dollar a prime number of cents. The obvious choice is 101 (in which case we might call them not cents but centunos). How would this affect roundage calculations?

The state of Ohio has recently changed its sales-tax regulations, allowing merchants to calculate tax either on a per-invoice or a per-item basis. The per-item option foils all the market-basket strategies for manipulating roundage, and it gives the merchant the potential of earning up to an extra $0.005 per item sold. It will be interesting to see whether pricing patterns change in Ohio.

Among children of a certain age, everyone has a best friend—and exactly one. Ideally, the best-friend relationship is symmetric: If I am your best friend, then you are my best friend, too. But symmetry is not guaranteed, and it can happen that I like you best, but you have someone else you like better than me. Sad, but life is like that sometimes.

We can model best-friendship geometrically by letting distance—or, rather, nearness—stand for intensity of affection. Sprinkle a bunch of points at random on a plane, and then draw an arrow from each point to its nearest neighbor, which we take to be the point’s best friend. When the best friends are mutual, a bidirectional arrow links the two points. In other cases, a chain of arrows points from a to b to c and so on. Here is a best-friend graph constructed in just this way:

There are 100 dots in the diagram, which have spontaneously formed 30 constellations, or connected clusters. Some 40 of the dots represent disconsolate children who feel an attachment for someone who’s attached to someone else. (Note that some of the dots are so close that the arrows are obscured.)

Questions.

1. In general, what proportion of the children in this model are stuck in unrequited best-friendships?
2. How many clusters form, on average?
3. How do these quantities vary as a function of the overall number of children?
4. Adding a new child to the class or the neighborhood could disrupt some existing friendships: If a new point is inserted at a random position, how many existing bonds are likely to be broken or rearranged?
5. What about the geometry and topology of the model? Would it make a difference if the points were plotted on the surface of a torus? Would the results be different in one dimension (with all the points arrayed along a line) or in three dimensions (with the points distributed throughout a volume of space)?

Note that the best-friend problem is not equivalent to the well-known stable-marriage problem (on which there is an extensive literature). In the stable-marriage situation, matchings are always two-by-two: If I like you best but you prefer someone else, then I simply have to find another partner.

For a rather different perspective on the mathematics of friendship (and also enmity) see “Dynamics of Social Balance on Networks,” by T. Antal, P. L. Krapivsky, and S. Redner.