The oddest numbers

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.