where a and b are functions of x but not on n. This is called a three-term recurrence formula.

It’s amazing how often you can run into three-term recurrence formulas. There are theorems that give conditions for such recurrences to hold, but I haven’t reached the bottom of that rabbit hole [1].

For this post I just want to give examples.

Bessel functions of the first and second kind:

Modified Bessel functions of the first and second kind:

Chebyshev polynomials of the first and second kind:

Hermite polynomials (physicists’ convention):

Legendre polynomials:

[1] See Bochner’s theorem for orthogonal polynomials, the Nikiforov–Uvarov method, and Infeld-Hull factorization.

When you’re looking at a probability density f(x) for the first time, it helps to ignore the normalizing constant. Concentrate on the part of the function involving x and know that the normalizing constant is whatever it has to be. For example, about half of the ink that it takes to write down a beta or chi-squared density is devoted to the normalization constant; the rest of the expression is easier to understand.

This post will do the opposite of the advice above and focus on normalization constants because this ties into the previous post on modified Bessel functions.

The von Mises probability distribution on a circle has two parameters, μ and κ, and its density function is

The normalizing constant is 2π I₀(κ). The factor of 2π is unsurprising for anything defined on a circle. The more interesting part is I₀, the modified Bessel function of order 0.

The von Mises-Fisher distribution is the generalization of the von Mises distribution to a sphere in p dimensions. The density function is

where the normalization constant C_p(κ) is

where I_{p/2 − 1} is the modified Bessel function of order p/2 − 1. The values of x and μ are in bold face because they are now vectors, points on the unit sphere.

When p = 2, we have the “sphere” in two dimensions, i.e. the circle, and the von Mises-Fisher distribution reduces to the von Mises distribution. But where did the cosine go? The inner product of x and μ is the cosine of the angle between the two vectors.

When p = 3, obviously an important special case, the von Mises-Fisher distribution is known as the Fisher distribution. In that case the normalizing constant C₃(κ) can be written without using modified Bessel functions because when ν = ½ + n for an integer n, I_ν(x) is an elementary function.

Things are given their names for reasons. Discovering those reasons may help you understand their motivation and use.

Pure math perspective

For each integer n, the modified Bessel function I_n is essentially the Bessel function J_n evaluated along the imaginary axis. Specifically,

From a certain shallow perspective, that’s the end of the story: modified Bessel functions are modified in the sense that the argument is multiplied by i. And there’s a fiddly constant term up front for no apparent reason.

But of course that’s not the end of the story or else this wouldn’t be worth an entire post.

The equation above is analogous to the relationships between circular and hyperbolic functions

These relationships are interesting because the circular and hyperbolic functions are independently meaningful. If you view these equations merely as definitions you lose their significance. Circular and hyperbolic functions were widely used before Euler discovered the connection between them.

Similarly, there’s a reason the modified Bessel functions were given a name their own. If you were led to Bessel functions and modified Bessel functions separately by different applications, you would regard the equation

as a discovery rather than just a definition. The following section explains why someone would be interested in modified Bessel functions.

Before we move on, I’d like to explain the reason for the term i⁻ⁿ term. In general

for all real ν. The reason for the exp(νπi/2) term is that it makes I_ν(x) real for all real x.

Applied math perspective

Bessel functions often arise from solving problems with radial symmetry. Solving the wave equation in cylindrical coordinates using separation of variables leads to Bessel’s differential equation

and its solutions J_n and Y_n, Bessel functions of the first and second kind.

Solving the heat equation in cylindrical coordinates with separation of variables leads to the modified Bessel equation

and its solutions I_n and K_n, the modified Bessel functions of the first and second kind.

This is the reason behind the complex analysis perspective above: the change of variables sending x to ix changes the sign of the x² term in Bessel’s equation.

Bessel functions describe radially symmetric oscillations, such as the vibrations of a drum head. Modified Bessel functions describe radially symmetric exponential growth or decay [1], such as in the heat in a cylinder.

Other modified functions

Struve functions are closely related to Bessel functions. The (modified) Struve functions also satisfy Bessel’s (modified) differential equation, but with a non-zero right hand side. The modified Struve functions are proportional to the unmodified Struve functions evaluated along the imaginary axis, with a proportionality constant that makes the modified Struve functions real for real arguments.

There’s a similar relationship between the Mathieu functions and modified Mathieu functions. The general pattern is that “modified” in the context of special functions means “evaluated at ix and multiplied by a constant to make the function real for real arguments.”

[1] The functions I_n grow exponentially and the functions K_n decay exponentially. For this reason, A&S didn’t tabulate I_n and K_n per se. Instead it tabulated e^−xI_n and e^xK_n because these functions varied less over their range.

The more widely useful a function is, the more likely it is that the function will be discovered independently multiple times. Independent discoveries lead to varying definitions and notations. For example, there are two widely used definitions of Hermite polynomials, one used in probability and another used in physics, that only differ by a scaling factor. This also explains why there are so many variations on the definitions of the Fourier transform and spherical coordinates.

Special functions were discovered and applied before they were studied systematically. As with most mathematics, practice preceded theory. In hindsight, some names and conventions were less than ideal, at least from the perspective of someone seeking to organize a theory.

Functions can have arcane names for several reasons, one being that their usefulness became apparent long ago. If you’re instinct is that things with strange names are no longer important, you’re instinct might be backward. The strange name may be an indication that something is so important that its usefulness became apparent long ago.

Sometimes special functions have bland, uninformative names because the names stuck before anybody could think of something better. Bob looks into an interesting family of functions [1], then later he finds another interesting family of functions. These become known as “Bob’s functions of the first kind” and “Bob’s functions of the second kind.” These names are quite understandable at the time, though in the future people will want to know what distinguishes the functions, other than the fact that Bob discovered them, and what the groupings have in common other than the order in which Bob found them.

I started this post intending to discuss modified Bessel functions and explain what exactly is modified about them, but my preface became its own post. “Modified” is an example of the bland terminology mentioned above. There are Bessel functions and modified Bessel functions. Without more context, the “modified” term isn’t very informative. But it does provide a clue that there’s some kind of close relationship between the modified and unmodified functions. That’ll be the topic of my next post.

[1] Math education doesn’t place much emphasis on history and motivation. You may have to do some digging to find out why Bob was interested in his functions. What else was Bob known for? Maybe they’re related.

If you wanted to be more accurate, you’d need to take into account the fact that the orbit of Mars is tilted about 1.85° relative to the Earth’s orbit. How much difference does that make?

To simplify things, let’s assume the Earth orbits the sun in a circle of radius 1 and Mars orbits the sun in a circle of radius 1.5. The distance between Earth and Mars over time would be basically sinusoidal.

How much does inclination contribute to this distance? In other words, what is the difference between the distance accounting for the inclination of Mars’ orbit and the distance if we assume the two orbits are in the same plane?

This plot gives the answer.

The effect is not large, about three orders of magnitude smaller than the main effect, but it’s interesting how erratic it is.

The plots were made with the following code.

from numpy import *

R = 1.5
T = R**1.5 # Kepler's third law

def f(t, theta):
    return sqrt(
        (cos(t) - R*cos(t/T)*cos(theta))**2 +
        (sin(t) - R*sin(t/T))**2 +
        (R*sin(theta)*cos(t/T))**2
    )

The first plot graphs f(t, θ) and the second graphs f(t, θ) − f(t, 0).

Suppose I look at the sign one day and guess how many days until the pizza shop will open. When I drive by a week later and guess again, should my guess be smaller? You might argue that the shop will open some day, fixed in time but unknown to me, and so every day I’m one day closer to the eventual opening.

You might model the pizza shop opening like radioactive decay and say that the estimated number of days until it opens is always the same until the day it actually opens.

Now I think this shop has been “coming soon” for over a year. So instead of decreasing, every day I increase my estimate of the time until the shop opens. Something has gone wrong that the owners didn’t expect when they put up the sign.

Maybe the reasonable thing would be for estimated days until opening to decrease over time, but only up to a point. After some point, the longer a business has been “coming soon” the less like that it is coming soon, or coming at all.

This brings up an interesting point about modeling. There are two probability distributions at work: the probability that the shop will eventually open, and the time until opening assuming it eventually opens.

When the sign first goes up saying the business is coming soon, there’s some change that it is in fact not coming. Maybe you’re optimistic and think this probability is small, but it would seem unreasonable to think the probability is zero. That means the expected number of days until opening is always infinite. If there’s a probability ε that the shop never opens, the expected time to opening is

ε × ∞ + (1 − ε) × something = ∞.

If you found fragments of writing from an ancient culture and inferred that five symbols were used as digits, how could you distinguish base 5 from bijective base 5? Suppose you believe these five symbols were digits

★ ︎ ︎ ☗ ︎

but you don’t know in what order. You just see sequences like ︎︎︎ and ★★︎ and believe they’re numbers.

If you noticed that numbers often contain ︎, but ︎ never appears at the beginning of a number, you might infer that ︎ is a zero. But this would take a fairly large sample. If you found only 20 numbers, for example, you could hardly conclude ︎ never appears at the beginning of a number just because it doesn’t come at the beginning of any number you’ve seen.

Now suppose you’ve found writing with more number symbols. Say you’ve found 17 numeric symbols. You might infer that the writing used a base 20 system, because it would be hard to imagine a human culture using base 17. Now imagine you find more fragments and confirmed that indeed there are 20 numeric symbols. Approached as a purely statistical problem, you’d need a very large sample to infer what the digits correspond to and whether they use a base 20 or bijective base 20 system (or some other system).

You’re best hope is to find numbers in some context where you know what number is being represented. If you knew somehow that some symbol corresponds to 20, then you’d know they didn’t use base 20 because base b doesn’t have a single symbol for b.

If you had a huge collection of numbers but no context, which is highly unlikely, you could use Benford’s law to infer the meaning of the number symbols: the most common leading digit is probably 1, the next most common is probably 2, etc. This is interesting to think about, but it seems much more realistic that a number system would be decoded by finding context, such as a list of consecutive numbers or numbers with known meaning.

Aside from the incongruity of using a gazillion-parameter neural network to make 1970’s style ASCII art, ASCII diagrams have some uses. They’re absolutely tiny compared to image files. But more importantly they can be inserted into plain text files, such as source code or markdown. A diagram embedded directly into a source file cannot become separated from the code.

ASCII diagrams are tedious to create, though there are tools to mitigate the tedium. But if an AI can generate the diagram, the tedium goes away.

I was curious how well Claude could create ASCII diagrams, so I tried a few examples. I hope these render well in whatever format you’re reading this post. They look fine for me previewing the post in a browser. I expect they might not turn out so well in an RSS reader.

I asked it to reproduce the graphs from my recently post on the graph imbalance theorem and the first diagram turned out nicely.

                 +-------+                                +-------+
                 |   A   |--------------------------------|   B   |
                 +-------+                                +-------+
                     |                                        |
                     |                                        |
   ------------------|------------------             ---------|---------
   |        |        |        |        |             |        |        |
   |        |        |        |        |             |        |        |
+-----+  +-----+  +-----+  +-----+  +-----+       +-----+  +-----+  +-----+
| a0  |  | a1  |  | a2  |  | a3  |  | a4  |       | b0  |  | b1  |  | b2  |
+-----+  +-----+  +-----+  +-----+  +-----+       +-----+  +-----+  +-----+

The second network is more complicated and so the corresponding ASCII diagram is hard to read.

   +----------------------------------------------------------------+
   |                                                                |
   |+-----------------------------------------------+               |
   ||                                               |               |
   ||+-------------------------------+              |               |
  +-------+       +-------+       +-------+       +-------+       +-------+
  |  R1   |-------|  R2   |-------|  R3   |-------|  R4   |-------|  R5   |
  +-------+       +-------+       +-------+       +-------+       +-------+
      |             | | |          |   |           |   |           |  |  |
     ++             | | |          |   |           |   |           |  |  |
     | +---------------------------+   |           |   |           |  |  |
     | |            ++| |              |           |   |           |  |  |
     | |             || +--------------+           |   |           |  |  |
     | |             || |  +---------------------------------------+  |  |
     | |             |+-|--|-------------+         |   |              |  |
     | |             |  |  |             |  +------+   |              |  |
     | |             |  |  |             |  |  +----------------------+  |
     | |             |  +--|-------------|--|--|-------------+           |
     | |             |  |  |             |  |  |       +-----|--+        |
     | |             |  |  |             |  |  |             |  |  +-----+
     | |             |  |  |             |  |  |             |  |  |
  +-------+         +-------+           +-------+           +-------+
  |  G1   |         |  B1   |           |  B2   |           |  B3   |
  +-------+         +-------+           +-------+           +-------+

For a third example, here is a fairly complicated diagram that nevertheless lends itself to a readable ASCII diagram. It’s a Feistel network diagram for DES encryption.

   +-------------+                    +-------------+
   |   L(i-1)    |                    |   R(i-1)    |--------
   +-------------+                    +-------------+       |
          |                                  |              |
          |                                  |              |
          |                      +-----------------------+  |
          |                      |   E (expand 32->48)   |  |
          |                      +-----------------------+  |
          |                                  |              |
          |                      +-----------------------+  |
          |                      |     XOR with K(i)     |  |
          |                      +-----------------------+  |
          |                                  |              |
          |                      +-----------------------+  |
          |                      |    S-boxes S1..S8     |  |
          |                      +-----------------------+  |
          |                                  |              |
          |                      +-----------------------+  |
          |                      |    P (permutation)    |  |
          |                      +-----------------------+  |
          |                                  |              |
          |                                  |              |
          |            +-------+             |              |
          +------------|  XOR  |-------------+              |
                       +-------+                            |
                           |                                |
          +----------------|--------------------------------+
          |                +------------------+
          |                                   |
   +-------------+                    +-------------+
   |    L(i)     |                    |    R(i)     |
   +-------------+                    +-------------+

Related posts

• The worst tool for the job
• ASCII armor
• Memorizing the ASCII table

[1] See Christian Wolff’s discussion of dogs playing poker in The Accountant (2016).

For odd n, the solution is what you might expect: a regular n-gon. I would expect this to be the solution for even n as well, but it’s not.

In 1974 [1] Ron Graham found a solution for n = 6, a hexagon with unit diameter and area larger than a regular hexagon with unit diameter. Polygons with diameter ≤ 1 are called “small”, and he found the “largest” (i.e. maximum area) small hexagon.

The vertices of Graham’s hexagon are given below.

  A = (0.0000000000,  0.0000000000)
  C = (0.4023506913, -0.5000000000)
  F = (0.9390533483, -0.3437714489)
  B = (1.0000000000,  0.0000000000)
  E = (0.9390533483,  0.3437714489)
  D = (0.4023506913,  0.5000000000)

You can verify that the distance between any pair of vertices is no more than 1 and that the area of Graham’s hexagon is 0.674981.

The area of a regular hexagon of diameter 1 is (3/8)√3 = 0.649519, and the area of Graham’s hexagon is about 3.9% larger.

[1] R. L. Graham. The Largest Small Hexagon. Journal of Combinatorial Theory (A) 18, 165–170 (1975). The paper was submitted February 22, 1974 and published in 1975.

What does the conjecture theorem say? Start with a graph G with no edge between two nodes of the same degree. Then for every edge, calculate the absolute value of the difference of the degree of each end. The imbalance theorem says there exists another graph H whose vertices have degrees corresponding to the differences of degrees in G.

For example, let G be the graph below.

The edges from the top red vertex A to each of the blue vertices around it all have degree difference 5 because A has degree 6 and the vertices a₀ to a₄ have degree 1. The edge between the two red vertices, A and B, has degree difference 2. The remaining vertices have degree difference 3.

So the multiset of degree differences is

{5, 5, 5, 5, 5, 2, 3, 3, 3}

The imbalance theorem says there exists a graph H whose nodes have these degrees. Here is an example of such an H.

Note that in H, the 5 red nodes have degree 5, the single green node has degree 2, and the three blue nodes have degree 3.

More graph posts

• Traveling salesman tour of Texas
• Effective graph resistence
• Visualizing graph spectra

The center of the star is not quite the center of the orbit. The planet moves along an ellipse with the star at one focus of that ellipse.

Let a be the semi-major axis of planet’s orbit, the maximum distance from the center of the ellipse to a point on the ellipse. Then the distance of a focus to the center of the ellipse is ae where e is the eccentricity of the ellipse. This defines eccentricity. The center of earth’s orbit is between three and four solar radii away from the center of the sun [1].

The planet is farthest from the star when it is along the major axis of the ellipse on the opposite side as the star. The distance is then a + ae, the distance to the center plus the distance from the center to the star. On the opposite side of its orbit, the planet is closest to the star. There the distance is a − ae. In summary the maximum distance to the star is

a(1 + e)

and the minimum distance is

a(1 − e).

If you had to guess the average distance between the planet and its star, a would be a good guess since it’s the average of the maximum and minimum distance. And that’s a good approximation, provided e is small. The mean distance over time is

a(1 + ½e²).

See derivation. The average distance is greater than a because the planet moves faster when nearest the star and slower when further from the star.

The relative error in approximating the mean distance by a is then ½e². When e is small, ½e² is very small. For the earth’s orbit, e = 0.01671, and so the approximation is off by around 0.014%.

The eccentricity of Pluto’s orbit is 0.2488, and so in that case the approximation is off by about 3.1%. The eccentricity of a Molniya orbit, used by some Russian satellites, is 0.74 [2]. For such satellites the error in approximating the mean distance to earth as the semimajor axis is around 27%.

Related posts

• Min and max orbital speed
• Roman moon, Greek moon
• From Kepler to Bessel

[1] For earth’s orbit, e = 0.01671, a = 1.496×10¹¹ m, and the sun’s radius is r = 6.957×10⁸ m. And so ea / r = 3.59.

[2] An object in such a highly elliptical orbit will spend a long time at the far side of its orbit, i.e. over Russia. Sort of a poor man’s geostationary orbit.

Starting with a Hadamard matrix H₀ and a Hadamard matrix G, you can construct a sequence of Hadamard matrices by

H_n+1 = G ⊗ H_n

for  positive integers n. This is known as the generalized Sylvester method.

Let p_n be the proportion of 1s in H_n and let q be the proportion of 1s in G. Then you can show that the recurrence holds

p_n+1 = q p_n + (1 − q)(1 − p_n).

You can solve the recurrence to show that

lim_{n → ∞} p_n = ½

and so as the iterations proceed, the ratio of number of 1s to the number of −1s approaches 1.

This doesn’t say anything Hadamard matrices in general, but it does apply to all Hadamard matrices created by repeatedly applying the generalized Sylvester method.

If you set G and H equal to the matrix

then p₀ = q = ¾. Then for n = 1, 2, 3, …, 8 the values of p_n are

0.625
0.5625
0.53125
0.515625
0.5078125
0.50390625
0.501953125
0.5009765625.

What is the probability that a pixel could be repaired if corrupted? That depends on your probability model. We will assume that the probability of each bit being flipped is p and that errors are independent.

(Are errors independent, i.e. if a bit flips, is the next bit more or less likely to flip? That would depend on context.)

It’s straight-forward to calculate the probability that 7 or fewer or fewer bits out of 32 flip; this is the cumulative distribution of a binomial random variable. The following Python code will return the probability of k or fewer successes out of n trials, each with probability of success p:

    from scipy.stats import binom
    print(binom.cdf(k, n, p))

For example, if there is a 10% chance that each bit will flip, there’s a 98.8% chance that 7 or fewer bits out of 32 will flip.

However this only gives a lower bound on the probability of correcting an error. If eight bits flip in transit, we cannot tell with certainty which codeword was sent, but there will be a couple possibilities that stand out. We’ll have to guess, but we’ve narrowed down the possibilities. With even more flipped bits, there’s always a chance of recovering the original data. Still, the lower bound captures most of the probability of recovery.

Now suppose you’re given a desired error recovery rate and have to determine what value of p it can sustain. For example, someone might say they want a 98.8% chance of recovering a pixel correctly, and you could come back and say p must be less than or equal to 0.1. This would be a conservative answer because as discussed above, p = 0.1 gives a pixel recovery probability of something more than 98.8, though it’s messy to calculate how much more.

You could solve for p by trial and error, or you could use some more sophisticated math to compute p directly. Given a probability F, you can solve for p such that the probability of up to k successes out of n trials using the inverse of the regularized incomplete beta function.

from scipy.special import betaincinv
p = 1 - betaincinv(n - k, k + 1, F)

Calculating F given n, k, and p could be a homework exercise in an introductory probability course. Solving for p given F, n, and k either requires some numerical programming or special functions and so would be a more challenging problem.

Related posts

• Golay code used in Voyager
• Vehicle Identification Number checksum
• Credit card checksum

A Hadamard matrix is an orthogonal matrix with all entries equal to ±1. Jacques Hadamard conjectured that there exist Hadamard matrices of order 4n for all positive integers n. It’s necessary that the order be divisible by 4, and Hadamard conjectured that this is sufficient [1].

How could you compactly represent a Hadamard matrix? Since the entries are all either 1 or − 1 each entry could be represented by a single bit, and n² bits could store an n × n Hadamard matrix. But we can do better.

Methodical matrices

If the matrix can be produced by an algorithm, you only need to store the name of the algorithm and the argument to the algorithm. So, for a 1024 × 1024 matrix applied by iterating Sylvester’s algorithm could be stored by saying “Apply Sylvester’s algorithm 10 times” rather than storing a megabyte of data.

Paley’s method can create a Hadamard matrix corresponding to every prime power. So you could determine a Paley type matrix by storing the prime and the exponent.

Next in complexity would be hybrid algorithms, such as start with the Paley method applied to 37⁶ and then apply Sylvester’s method 3 times.

There are more methods of creating Hadamard matrices than Sylvester’s method and Paley’s method, though they’re harder to describe and parameterize.

Sporadic matrices

If a Hadamard matrix cannot be constructed using an algorithm, you can still store the matrix in fewer than n² bits. Since the rows are orthogonal, the last row of the matrix is determined by all the previous rows, up to sign. So you could store a Hadamard matrix using n(n − 1) + 1 bits.

Some Hadamard matrices are symmetric or skew. A symmetric matrix is determined by its diagonal and the elements above the diagonal. So a symmetric Hadamard matrix could be represented by n(n + 1)/2 bits.

A skew Hadamard matrix isn’t quite skew-symmetric. A matrix M is skew symmetric if

M^T = −M.

This implies the diagonal elements are 0, and Hadamard matrices cannot contain 0s. A Hadamard matrix H is called skew if

H + H^T = 2I.

This implies the diagonal elements are all 1s and the elements below the diagonal have the opposite sign of the elements above the diagonal. Since the elements on the diagonal are determined, a skew Hadamard matrix can be sorted using n(n − 1)/2 bits.

Incidentally, there is a conjecture that there exist skew Hadamard matrices of order 4n for all positive n.

[1] There are Hadamard matrices of order 1 and 2, but larger orders must be divisible by 4.

This evening I gave an example, describing how NASA used a Hadamard matrix of order 32 to transmit photos from the Mariner 9 spacecraft in 1971. This post will give another application: sphere packing.

Conway and Sloane [1] give a correspondence between binary codes and sphere packings that they call Construction A. Given an (n, M, d) binary code C, center a sphere on a point x if and only if x is a congruent (mod 2) to codeword in C.

Here (n, M, d) means an error correcting code that encodes M bits of data as strings of n bits, with a minimum Hamming distance between code words of d, i.e. all codewords differ in at least d bits.

The previous post described how to create a (32, 6, 16) code by stacking a Hadamard matrix H of order 32 on top of −H and turning −1’s into 0’s. The analogous construction for a (8, 4, 4) Hadamard code gives E₈, the densest packing in ℝ⁸.

We start with the Hadamard matrix

and obtain the matrix

whose centers form the sphere packing.

This doesn’t look like the E8 sphere packing as it is usually presented, but it’s isomorphic.

[1] J. H. Conway and N. J. A. Sloane. Sphere Packings, Lattices and Groups. Springer. 1999.

The images were encoded for transmission using a code based on Hadamard matrices, specifically a (32, 6, 16) Hadamard code. This means that each 6-bit pixel value was encoded as a 32-bit code word, with all code words differing in at least 16 positions.

The previous post explained a way to construct Hadamard matrices of order 2ⁿ. Use this process to create a 32 × 32 Hadamard matrix H and create a 64 × 32 matrix M by stacking H on top of −H. Then form a matrix M′ by changing all the −1 entries to 0. The rows of M′ are the code words.

For a 6-bit photo pixel value, one of the bits determines whether to read a code word from the top half or bottom half of M′. The other five bits determine which row to choose.

So a pixel is transmitted as a 32-bit codeword c, one of the 64 rows of M′. Ideally c would be received, but possibly some corrupted versions c′ is received with some of bits flipped.

Replace all the 0’s in c′ with −1 to create c″. Now multiply M by c″, thinking of the latter as a column vector. This yields a column vector of length 64. The largest component of this vector corresponds to the row of M′ that was most likely sent.

To see this, suppose there was no corruption: c was transmitted and c was received. Then the product Mc″ has a 32 in the entry corresponding to c and zeros everywhere else. If no more than 7 bits in c were corrupted, the row with the largest entry corresponds to the row that was transmitted.

In practice the product Mc″ can be computed using an algorithm analogous to the FFT using fewer operations than it would take to multiply a general 64 × 32 matrix by a 32 × 1 matrix.

is a Hadamard matrix of order 2. True to Stigler’s law of eponymy, James Joseph Sylvester investigated Hadamard matrices before Jacques Hadamard. Sylvester saw how to bootstrap the example above into more examples. If H is a Hadamard matrix, then the partitioned matrix

Sylvester’s construction can be generalized as follows. If H_m is a Hadamard matrix of order m and H_n is a Hadamard matrix of order n, the the Kronecker product H_m⊗H_n is a Hadamard matrix of order mn. That is, you can form a new Hadamard matrix by taking the matrix H_m and replacing ±1 with the matrix ±H_n.

Let S be the set of all possible Hadamard matrix orders. By the construction above, this set is closed under multiplication. Since 2 is in S, every power of 2 is in S. Hadamard proved that all n ≥ 4 in S are multiples of 4. That is, the condition 4 | n is necessary. He conjectured that it was also sufficient, though that has not been proven.

So the big question is what is the set S. Is there some multiple of 4 not in S? Until that question is answered, what is the smallest multiple of 4 not known to be in S? Hadamard matrices are useful in applications, so constructing Hadamard matrices of various orders is useful even while Hadamard’s conjecture remains open. For example, see the next post for how NASA used Hadamard matrices to transmit photographic images back from Mars.

Paley’s method

Raymond Paley came up with a way of constructing Hadamard matrices of size q + 1 if q is a prime power congruent to 3 mod 4, and of size 2(q + 1) if q is a prime power congruent to 1 mod 4. Let’s see what we can squeeze out of this.

If p is a prime congruent to 1 mod 4, every power of p is also congruent to 1 mod 4, and so there exist Hadamard matrices of order 2(p^k + 1) for every k.

If p is a prime with p = 3 mod 4, then even powers of p are congruent to 1 mod 4 and odd powers of p are congruent to 3 mod 4. So there are Hadamard matrices of order 2(p^2k + 1) and of order p^2k+1 + 1.

Let’s run a script to see what we can learn from this.

from sympy import primerange

s = set()

for p in primerange(20):
    if p % 4 == 1:
        s.update([2*(p**k + 1) for k in range(1, 10)])
    if p % 4 == 3:
        s.update([2*(p**(2*k) + 1) for k in range(1, 6)])
        s.update([p**(2*k + 1) + 1 for k in range(1, 6)])
print(sorted(s)[:20])

This prints

12, 20, 28, 36, 52, 100, 164, 244, 252, 340]

We can add 16 to the list because it’s a power of 2, and we can add 24 because it’s 2 × 12, etc. But there doesn’t seem to be any way to get 44. There is a way to create a Hadamard matrix of order 44, but it doesn’t follow from anything we’ve seen so far.

New records

I have a book published in 1996 that says Hadamard’s conjecture had been verified for n up to 428. Until yesterday, the smallest multiple of 4 for which nobody had found a corresponding Hadamard matrix was 668. Then Levent Alpöge announced that he and his and collaborators found an example of size 668 and filled in all remaining gaps below 2000.

So now the set S is known to contain {1, 2, 4, 8, 12, 16, …, 2000}. It also contains all orders that can be obtained by Paley’s method and other methods. And it contains all products of its elements. But it is not yet known to contain 2004.

Related posts

• Hadamard’s determinant inequality
• Hadamard’s Dirichlet energy example

But then after you use the command line a while, you might begin to notice a pattern. You often have to run a pair of commands with the same last argument, such as when you make a directory then navigate into it.

mkdir /some/long/path/name
cd /some/long/path/name

Maybe someone sees you typing this and offers a helpful suggestion: you can type !$ to repeat the last part of the previous command.

mkdir /some/long/path/name
cd !$

Once you learn this, you’ll see uses for it regularly. Sorta like when you buy a green Toyota Corolla and suddenly become aware how many green Toyota Corollas there are out there.

So then you start to use !$, but you have some doubts: is the syntax !$ or $!? As long as it’s just two cryptic symbols, you may have a hard time remembering.

Although the command line is cryptic, it is also consistent. On the command line, and in Unix [1] culture more generally, ! refers to a command. For example, you can run a command from inside the Vim text editor with the :! command. And $ often refers to the end of something, such as the end of a string in a regex. So you type !, telling the shell that you’re about to refer to part of a string, then type $ to specify that last part.

The characters !$ are an example of an event designator, a sort of mini language for referring to parts of commands. The syntax may seem mysterious, but every bit of it is analogous to syntax elsewhere in Unix culture. For example,

!:s/foo/bar/

replaces the first instance of foo with bar, and

!:gs/foo/bar/

will replace all instances. This is very similar to commands in sed and in vi For more on event designators, see this recent article by Filip Roséen.

There’s a saying that when the student is ready, the teacher will appear [2]. A less woo-woo version of this saying would be that until the student is ready, it doesn’t matter whether a teacher is there. Once you’ve used the command line to see a need for various shortcuts, there’s a much better chance they’ll stick.

Related posts

• Just-in-case versus Just-in-time
• Brace expansion
• Shell variable ~-

[1] Linux has so thoroughly eclipsed Unix that many people say “Linux” when they’re referring to things inherited from Unix that are now part of Linux and other operating systems, as well as various software packages and programming languages.

Here’s a plot of references to Unix versus Linux over time, via Google’s ngram viewer.

[2] When I first heard this, someone told me it was an Asian proverb. It sure sounds like one. But apparently it comes from a 19th century pastiche of eastern and western ideas.

1. Bella
2. Luna
3. Max
4. Charlie
5. Coco
6. Lola
7. Rocky
8. Milo
9. Teddy
10. Lucy

I wondered if the name frequencies might fit a power-law distribution. They do not, but they follow a log-normal distribution remarkably well.

Related posts

• Metabolism and power laws
• Why isn’t everything normally distributed?
• Power laws and the generalized CLT

Manual cryptography

If you told someone about RSA encryption, for example, you’d lose them right after you said “First you find a couple 1000-digit primes.” But there’s no need for using 1000-digit primes if no attacker has a computer. You could use 100 digit primes. Could you use 10 digit primes? If you chose primes just big enough to make the method unbreakable by hand, could someone implement it by hand?

Kirchoff’s principle says the strength of an encryption method should depend only on keeping the key secret, not the method. If you could keep the method secret, RSA would be unbreakable because nobody thought of anything like it before computers. But to make our thought experiment more interesting, let’s suppose that an enemy has also traveled back in time. If you tell your side about RSA, he can tell his side about it as well. So we’re back to Kirchoff’s principle.

An encryption method combinining substitution and permutation would have been practical to carry our manually. The ADFGVX cipher from 1918 was a start in this direction. That idea could been extended further, with a larger substitution set and longer permutations, and with more than one round of substitution and permutation, approaching what would be come the approach used in modern symmetric encryption. Such a method might have been manually implementable without being manually breakable.

Mechanized cryptography

World War II was a time of transition from manual cryptography to computerized cryptography. Encryption machines were attacked by cryptanalysis machines, though these machines were general-purpose computers. If you could implement a symmetric encryption method like AES in a mechanical device, no mechanical device could break it.

You could use something like DES, simpler than AES but still unbreakable at the time. DES is considered broken because now you could throw enough compute power at it to break it by brute force, but that would not be possible with only mechanical devices.

My hunch is that the best approach would be stream ciphers. Maybe it would be practical to implement one of these by hand or with the aid of simple calculating machines. Something like PCG, which is not cryptographically secure today [1], would have been then, though I don’t know how practical it would have been to carry out PCG, say, in the 1940s.

More pre-computer cryptography

• Vigenère cipher
• Playfair cipher
• Straddling checkerboard

[1] In 2020, Charles Bouillaguet, Florette Martinez, and Julia Sauvage were able to break PCG using 20,000 CPU-hours. See their paper Practical seed-recovery for the PCG Pseudo-Random Number Generator. IACR Transactions on Symmetric Cryptology. ISSN 2519-173X, Vol. 2020, No. 3, pp. 175–196.

Many of the exercises below ask a modern reader to find and/or to correct errors in the literature of bygone days. The point is not to gloat over how smart we are in the 21st century; the point is rather to understand that even the pioneers of a subject can stumble. One good way to learn that a set of ideas is not really as simple as it might seem to today’s computer scientists and mathematicians is to observe that some of the world’s leading thinkers had to struggle with the concepts when they were new.

You might also do the opposite, especially in a mathematics class: look at what differential equation a set of functions satisfy in order to understand those functions.

Bessel functions came out of solving differential equations from astronomy. But then they turned out to be useful, not just in other areas of science, but in pure mathematics as well. So there are Bessel function users who came to the functions first and haven’t seen the differential equation they came from.

You can learn a lot about Bessel functions, and other special functions, by looking at their defining differential equation even if you’re not directly interested in the differential equation or the physical problem that motivated it.

Bessel functions satisfy

x² y″ + x y′ + (x² − n²) y = 0.

You can tell a lot about Bessel functions just by inspecting this equation without solving it. If we divide by x² and write the equation in the form

y″ + (p(x)/x) y′ + (q(x)/x²) y = 0

then p(x) = 1 and q(x) = x² − n². The indicial equation

r(r − 1) + p(0) r + q(0) = 0

reduces to

r² = n²

and so r = ±n. That alone tells us there are two solutions, one analytic at zero and one singular at zero. These are J_n and Y_n respectively. It also tells us the behavior of these functions as x goes to zero and as x goes to infinity. To find out more, look up “method of Frobinius.”

According to Claude,

When X-rays are diffracted by a helical structure, the resulting diffraction pattern breaks into a series of horizontal “layer lines.” Cochran, Crick, and Vand showed mathematically that the diffracted amplitude on the n-th layer line is proportional to a Bessel function of the first kind, order n:

J_n(2π r R)

where r is the radius of the helix and R is the distance out from the center (the meridian) in the diffraction pattern.

The citation for this paragraph is a paper from 1952 [2] that amazingly is behind a paywall.

Update: Evan pointed out something fishy in the expression above: “Isn’t there a dimensional problem with J(2π r R)? Maybe 2πr/R instead?”

This was Claude’s response:

• r is the radius of the helix, a real-space length (say, in Å).
• R is the radial coordinate in reciprocal space — it has units of inverse length (Å⁻¹), because that’s how X-ray diffraction/Fourier space works (spatial frequency, not distance).

So Evan and Claude were both right. There is an implicit division happening, but crystallographers do write r R.

Related posts

• Mr. Bessel’s eponymous functions
• From Kepler to Bessel
• Bessel zero spacing

[1] Watson and Crick didn’t “discover DNA” as is commonly said. DNA was discovered in 1878. Watson and Crick discovered the structure of DNA in 1953.

[2] Cochran, W., Crick, F. H. C., & Vand, V. (1952). “The Structure of Synthetic Polypeptides. I. The Transform of Atoms on a Helix.” Acta Crystallographica, 5(5), 581–586.

If you want to compute the sine or cosine of an angle x you could start by reducing x mod 2π since that would not change the result. However, accurately reducing a number mod 2π is not trivial; that’s why range reduction is an area of algorithm development.

Range reduction mod π/2

Even better would be to reduce x mod π/2. Reducing to a smaller range means that power series method, and other methods such as rational approximation, will be more efficient.

So suppose you can find an integer k such that

x − k π/2 = y

where 0 ≤ y ≤ π/2. Then sin(x) is ±sin(y) or ±cos(y), depending on k mod 4 equals 0, 1, 2, or 3.

from math import *

def reduced_sin(x, k):
     match k % 4:
        case 0: return sin(x)
        case 1: return cos(x)
        case 2: return -sin(x)
        case 3: return -cos(x)

Naive range reduction

Now let’s set x = 500. Then k = 318 because that’s the multiple of π/2 we need to subtract to bring x into range, and the sine of x should be the negative of the sine of the reduced value y because 318 = 2 mod 4.

The following code computes sin(x) with naive range reduction

def naive_sin(x):
    k = floor(x / (pi/2))
    y = x % (pi/2)
    return reduced_sin(y, k)

and when x = 500 the error is on the order of 1.7 × 10⁻¹⁴.

Better range reduction

The value of k above is fine, but we’d like to calculate y more accurately. The following code is much better.

def Cody_Waite_sin(x):
    C1 = 1686629713 / 2**30
    C2 = 4701928774853425 / 2**86

    k = floor(x / (pi/2))
    y = (x - k*C1) - k*C2
    return reduced_sin(y, k)

This will compute sin(500) to full machine precision. What kind of magic is this?

The trick is that the exact value of C1 + C2 equals π/2 to more precision than is possible in a single float [1]. You can confirm, with bc or some other extended precision software, that the difference between C1 + C2 and π/2 is roughly 2⁻⁸⁸, while the limit of float precision is 2⁻⁵².

If we compute

y = x - k*(C1 + C2)

then we’re doing the same calculation as naive_sin and will get the same error. But if we compute

y = (x - k*C1) - k*C2

we will get a more accurate result, provided x isn’t too large.

You can use the following code to play around and see how large x can be before errors start to creep in. For small enough x, like 500, the Cody and Waite sine returns full precision. For larger x it’s better than naive sine but does not return full precision. And for large enough x it completely breaks down.

def compare(x):
    y0 = naive_sin(x) 
    y1 = Cody_Waite_sin(x)
    y2 = sin(x)
    print("Naive error:     ", y2 - y0)
    print("Cody Waite error:", y2 - y1)

Now this may seem circular since we’re using math.sin as our gold standard. However, this function is calling the sine function on your CPU, which is using sophisticated range reduction to compute its result accurately down to the last bit, assuming you run the code on a computer that’s less than 40 years old.

The Cody and Waite algorithm is inadequate for large x, but it’s a good place to begin studying range reduction. It shows there are clever ways of squeezing out more precision than seems possible.

[1] The numerator n₁ of C1 is ⌊2³⁰ π/2⌋. The numerator n₂ of C2 is the solution to

2⁸⁶⁻³⁰ n₁ + n₂ = ⌊2⁸⁶ π/2⌋.

When I type an apostrophe ' on my laptop, it becomes ’ on my phone. And when I type 's on my phone, it becomes 痴 on my laptop.

Apparently the phone turns the apostrophe (U+0027) into a right single quote (U+2019), then bungles bytes in the UTF-8 encoding of U+2019 as three Windows-1252 characters. The bytes E28099_hex are interpreted as â (E2_hex), € (80_hex), and ™ (99_hex).

When I type 's on my phone, it is encoded as two Windows-1252 characters 92_hex and 73_hex. Then by the time the text appears on my laptop, the bytes 9273_hex are interpreted as a Shift-JIS encoding of the CJK character 痴 (U+75F4).

Here’s Python code to reproduce the problem.

def mojibake(s: str, n: int, bad_encoding: str = 'cp1252') -> str:
    for _ in range(n):
        s = s.encode('utf-8').decode(bad_encoding, errors='replace')
    return s

print(mojibake("’s", 1))
print(mojibake("’s", 1, 'shift_jis'))

The code won’t corrupt ASCII text. The problem started with an ASCII character being replaced by a similar non-ASCII character. Text containing a non-ASCII character gets more corrupted with each round.

Power series are an excellent way to calculate functions near the center of the series, such as computing sine for small angles. But the further you get from the center, the less useful power series are.

Let’s suppose you want to calculate cos(200) using the power series for cosine. The nth term of that series is

(−1)ⁿ x²ⁿ / (2n)!

This is an alternating series, and so the error in truncating the series after n terms is bounded by the size of the n+1 term, if you’ve gone far enough out in the series that the terms are monotonically decreasing in absolute value.

To calculate cos(200) to machine precision, i.e. with an error of less than 2⁻⁵², we’d need to sum the series up to n where

| 200²ⁿ⁺² / (2n + 2)! | < 2⁻⁵²

Actually, that will ensure that the absolute error is small enough, but not that the relative error is small enough; if the value of cos(200) is small, we’d need more terms. Let’s ignore that and assume we’re only concerned with absolute error.

Turns out we’d need 287 terms. That’s a lot of terms. But you might say “That’s fine. I’m not in a hurry, and it’s just more work for the computer, not for me.” OK, so let’s try.

from math import *

s = 0
for n in range(288):
    s += (-1)**n * 200**(2*n) / factorial(2*n)
print(s)

This prints -3.6840358571084123e+67. You may suspect the answer is incorrect since values of cosine are on the order of 1, not on the order of 10⁶⁷. Something went spectacularly bad. On closer inspection, it’s remarkable the code didn’t crash.

If you changed 200 to 200.0 above, the code would crash. Calculating 200.0**(2*n) overflows when n = 67. But when we calculate 200**(2*n), the result is an integer. And we’re dividing by factorial(2*n), which is also an integer. Both of these integers become too large to fit in a float, but their ratio has a maximum value of around 10⁸⁰, smaller than the maximum float, which is on the order of 10³⁰⁸.

When we don’t overflow, we have a different problem: catastrophic cancellation. You can’t calculate a number between −1 and 1 as an alternating sum of numbers as large as 10⁸⁰. You’d need more than 80 + 16 = 96 decimal places of precision to compute the sum accurately, and floating point only gives you between 15 and 16 decimal places of precision.

So how would you calculate cos(200)? The first step would be to use some sort of range reduction on 200. You could reduce 200 mod 2π to get a smaller number to work with.

>>> from math import cos, pi
>>> x = 200 % (2*pi)
>>> x
5.221255477432827

Using a power series to compute the cosine of 5.221255477432827 is feasible, but not optimal. There’s also another problem: the naive range reduction above loses some precision.

>>> cos(x)
0.48718767500701254
>>> cos(x) - cos(200)
6.661338147750939e-15

The error is small, but it’s still an order of magnitude larger than machine precision. You can’t simply reduce n mod 2π with ordinary float division because the integer part of n / 2π pushes some digits of precision off the right end. I intend to write about how range reduction works in future posts.

Update: See this post for a simple range reduction algorithm that is fine for values of x such as 200, but not adequate for much larger values.

In this post I’ll use n = 200! as my example rather than 1000! because this value of N is larger than the largest representable floating point number but small enough to be more convenient to work with.

Suppose someone calculates 200! for you:

78865786736479050355236321393218506229513597768717326329474253324435\
94499634033429203042840119846239041772121389196388302576427902426371\
05061926624952829931113462857270763317237396988943922445621451664240\
25403329186413122742829485327752424240757390324032125740557956866022\
60319041703240623517008587961789222227896237038973747200000000000000\
00000000000000000000000000000000000

You could now calculate log(n) using

n = 7.886578673647905 × 10³⁷⁴

and so

log(n) = log(7.886578673647905 × 10³⁷⁴)
= log(7.886578673647905) + 374 log(10) = 863.2319871924055.

The key thing that makes this possible is that the least significant digits of n only affect the least significant digits of log(n). In the calculation above I kept the first 16 digits of n. Python couldn’t make use of any more digits, and had no need of any more digits, in order to produce the logarithm to machine precision.

Cosine doesn’t work that way. The cosine of n depends on the remainder when n is divided by 2π, and that remainder depends on every single digit of n. I’ll illustrate that below.

Using bc -l and setting the scale to 400, I can calculated n then calculate

cos(n + 10ⁱ)

for i running from 0 to 374, tweaking each digit one at a time. (Except when a digit is a 9 and the addition results in a carry.)

    n = 1
    for (i = 1; i <= 200; i++) n *= i
    scale = 400
    for (i = 1; i <= 374; i++) {
        x = c(n+10^i)
        scale = 16
        print x/1, "\n"
        scale = 400
    }

Here’s what a plot of the results look like.

The value of cos(n) is about −0.985, but the values above are all over the map. We can look at the range by projecting all the points over to the left edge then rotating a quarter turn:

The remarkable thing about this image is that there are a few gaps, i.e. a few values the cosine does not take on.

Here’s a more sophisticated way to look at it. The sequence 10ⁱ mod 2π is dense in [0, 2π], and so by going far enough out in the sequence, we can find a value that shifts the phase of n by any desired amount within any given tolerance.

Every digit in n matters, and changing any digit can change the value of cosine to be essentially any value. You cannot calculate the cosine of an enormous number without using some kind of extended precision arithmetic. There are clever range reduction algorithms that minimize the amount of extended arithmetic necessary, but extended arithmetic cannot be completely eliminated.

>>> from math import log, factorial
>>> log(factorial(1000))
5912.128178488163

If you don’t find this unexpected, note that if you replace math.log with numpy.log the code will fail [1]. Functions like natural logarithm operate on real numbers. Real numbers are represented as floating point numbers in programming languages, and 1000! factorial is too large to represent as a standard floating point number. (More on that here.)

In this post I’d like to look at how you might calculate log(1000!) with less capable software, and even without software.

One approach would be to sum the logarithms of the numbers 1 through 1000. This will give essentially the same result as above, with a little difference in the last couple decimal places due to rounding error.

If you have a way to calculate 1000! but not a way to cast it to a floating point number, you could do this manually.

>>> s = str(factorial(1000))
>>> s[:16]
'4023872600770937'
>>> len(s)
2568

This tells us 1000! = 4.023872600770937 × 10²⁵⁶⁷. Therefore

log(1000!) = log(4.023872600770937) + 2567 log(10)

which only requires working with numbers of modest size.

Calculating by hand

Now suppose it’s 1964. You don’t have a computer, or even a calculator, but you do have a copy of the recently published Handbook of Mathematical Functions by Abramowitz and Stegun (A&S). You turn to Table 6.6 “Factorials for large arguments.” This has values of factorial for 100, 200, 300, …, 1000, so you can simply look up your answer to 20 decimal places.

That was too easy; I didn’t expect that to be there when I started writing this post. If you wanted to compute log(950!), for example, you’d have to work harder. You could find A&S equation 6.1.41 (Stirling’s series) which says

So how would you use this formula to calculate log(1000!)? Since n! = Γ(n + 1), you set z = 1001.

You’d need to decide how many terms you need to use. Assuming the error is on the order of the first term you leave out, you’d reason that you could probably stop with the 1/12z term because the next term is between 10⁻¹¹ and 10⁻¹².

You find Table 4.2 has natural logarithms, but not for 1001. You can look up log(1.001), however, and at the bottom of the same page is log(10) to 16 decimal places, and you can find log(10) to 24 decimal places in Table 1.1. So you calculate

log(1001) = log(1.001 × 10³) = log(1.001) + 3 log(10).

You can find log(2) and log(π) in Table 1.1, and average them to find ½ log(2π).

Here’s Python code to simulate the hand calculations.

log2     = 0.6931_47180_55994_53094_172321 # Table 1.1
log10    = 2.3025_85092_99404_56840_179915 # Table 1.1
logpi    = 1.1447_29885_84940_01741_43427  # Table 1.1
log1_001 = 0.00099_95003_330835            # Table 4.2

z = 1001
logz = log1_001 + 3*log10
s = (z - 0.5)*logz - z + (log2 + logpi)/2 + 1/(12*z)

print(s)

This result differs from the one at the top of the post only in the last decimal place.

Related posts

Doing calculations with tables is not as simple as “just look it up.” It takes a bit of skill.

• Using a table of logarithms
• Using a trig table
• How much precision can you squeeze out of a table?

[1] The code will also fail if you replace math.log with math.cos. Both logarithm and cosine return moderate sized real numbers when given enormous inputs like 1000!, so representing the output as a float is not the problem. But logarithms of huge numbers can be computed with ordinary precision functions, as above. But computing the cosine of a huge number requires extended precision.

Update: The next post expands on why computing the cosine of a large number is more difficult than computing the log.

The code computes logarithms for integers larger than the largest representable float. For example, the largest float is on the order of 10³⁰⁸, and yet the following code works.

>>> import math
>>> math.log10(10**400)
400.0

The log, log2, and log10 functions have some code inside that handles large integers specially. It doesn’t simply convert the integers to floats before taking the logarithm. If it did, it would overflow. If you replace math with numpy above, the code will fail. NumPy’s implementation of logarithms is more what I would expect.

While playing around with this I also noticed that you can define floats larger than the largest float without warnings.

>>> math.log(1e308)
709.1962086421661
>>> math.log(1e309)
inf

This isn’t a feature of math.log but of how Python handles scientific notation. The expression 1e308 is the floating point representation of 10³⁰⁸. It is a float, not an int.

>>> type(1e308)
<class 'float'>

The expression 1e309 is also a float. But since it’s larger than is possible for a float, Python interprets it as inf. The code

math.log(1e309)

returns inf based on the reasoning that log(∞) = ∞.

That explains the following behavior:

>>> 1e309 == 1e310
True

The expressions 1e309 and 1e310 are equal because both are alternate ways of writing inf.

The number of rooted trees with n nodes is the same as the number of ways to configure n − 1 non-overlapping circles. Not only are the counts the same, there is a natural correspondence between the trees and the circles. It’s not obvious that there should be such a correspondence, with the right notation the correspondence is sort of a pun.

The standard way to represent unlabeled trees is as a multiset of their children. We use a multiset, not a set, because some elements will be repeated. We represent a leaf as a pair of parentheses: ().

There is only one rooted tree with one node: ().

There is only one rooted tree with one two nodes: (()). Here the outer parentheses represent the root node and the inner parentheses represent its child.

There are two rooted trees with three nodes, and we can represent them as ((())) and ((),()). The first is the straight line tree: a node that has a single child node that has a single child node. The second is a node that branches to two nodes. (Here’s where we need multisets.)

The four rooted trees with four nodes can be represented as (((()))), ((((),())), ((),(())), and ((),(),(),()).

Here are the nine rooted trees with five nodes:

((((()))))
((((),())))
(((),(())))
(((),(),()))
((()),(()))
((),((())))
((),((),()))
((),(),(()))
((),(),(),())

The correspondence with non-overlapping circles removes the outer parentheses then joins the rest to form circles, with nested parentheses corresponding to concentric circles. A more geometric way to see the correspondence is to start at the bottom of the tree, replace leaves with circles, then work your way up circling connected components.