Solving hard problems

We help companies solve hard problems in mathematics, statistics, and computing. Let’s explore how we might work together.

Surprising curves with sine and sn

Posted on 23 June 2018 by John

In the previous post I said that the Jacobi functions are like trig functions. That’s true if you look along the real axis. If you look at the rest of the complex plane you’ll see how they can be very different.

The sine function is periodic along the real axis, but it grows exponentially along the imaginary axis. I mentioned parenthetically last time that the Jacobi functions are not just periodic but doubly periodic. That means that not only are they periodic as you move along the real axis, they’re periodic when you move along any line in the complex plane [1].We’ll illustrate this with some plots.

It will make our code more compact if we first define a function to split a complex number into its real and imaginary parts.

    pair[z_] := {Re[z], Im[z]}

Here’s what it looks like when you plot the real and imaginary parts of the sine function along the line (10 + 0.5i)t.

    
    ParametricPlot[
        pair[ Sin[(0.5 I + 10)t] ], 
        {t, 0, 10}, 
        PlotRange -> All
    ]

plot of sine with complex argument

By contrast, here’s a plot of the sn function along the line 1.25 + it.

    ParametricPlot[
        pair[ JacobiSN[1.25 + I t, 0.5] ], 
        {t, 0, 10}
    ]

It’s a closed loop because sn is periodic in every direction. (By the way, the curve looks like an egg. More posts about egg-shaped curves starting here.)

When you plot sn along various lines you’ll always get closed curves, but you won’t always get such round curves. Here’s an example that doesn’t look like a closed loop because the curve turns around sharply at each end.

 
    ParametricPlot[
        pair[ JacobiCN[ (I + 0.5) t, 0.5] ], 
        {t, 0, 10}
    ]

To show that it really is periodic, we’ll add a vertical time dimension to visualize how the curve is traced out over time.

    triple[z_, t_] = {Re[z], Im[z], t}
    ParametricPlot3D[
        triple[JacobiSN[(I + 0.5) t, 0.5], 0.1 t], 
        {t, 0, 20}
    ]

Here’s another curve plotting sn along a different line.

    ParametricPlot[
        pair[ JacobiSN[(0.9 + I) t, 0.5] ], 
        {t, 0, 55},
        PlotRange -> All, 
        PlotPoints -> 100
    ]

As before, we can add a time dimension to imagine how the curve is traced out.

    ParametricPlot3D[
        triple[JacobiSN[(0.9 + I) t, 0.5], 0.1 t], 
        {t, 0, 150},
        PlotRange -> All,
        PlotPoints -> 200
]

Here’s a variation on the plot above that stretches the vertical axis so you can better see what’s going on.

   
    ParametricPlot3D[
        triple[JacobiSN[(0.9 + I) t, 0.5], 0.2 t], 
        {t, 0, 150},
        PlotRange -> All,
        PlotPoints -> 200
    ]

[1] To be more precise, elliptic functions are periodic in two linearly independent directions, and thus in any direction that’s an integer linear combination of those two directions. So they’re exactly periodic in a countable number of directions almost periodic in every other direction.

System of Jacobi elliptic functions

Posted on 22 June 2018 by John

Jacobi’s elliptic functions are sorta like trig functions. His functions sn and cn have names that reminiscent of sine and cosine for good reason. These functions come up in applications such as the nonlinear pendulum (i.e. when θ is too
large to assume θ is a good enough approximation to sin θ) and in conformal mapping.

I ran across an article [1] yesterday that shows how Jacobi’s three elliptic functions—sn, cn, and dn—could be defined by one dynamical system

$\begin{eqnarray*} x' &=& yz \\ y' &=& -zx \\ z' &=& -k^2 xy \end{eqnarray*}$

with initial conditions x(0) = 0, y(0) = 1, and z(0) = 1.

The parameter k is the modulus. (In Mathematica’s notation below, k² is the modulus.) As k decreases to 0, sn converges to sine, cn to cosine, and dn to 1. As k increases to 1, sn converges tanh, and cn and dn converge to sech. So you could think of k as a knob you turn to go from being more like circular functions (ordinary trig functions) to more like hyperbolic functions.

Since we have a dynamical system, let’s plot the solution, varying the modulus each time.The Jacobi functions are periodic (in fact they’re doubly periodic) and so the plots will be closed loops.

f[t_, m_] = {JacobiSN[t, m], JacobiCN[t, m], JacobiDN[t, m]}
ParametricPlot3D[f[t, 0.707], {t, 0, 10}, AspectRatio -> 1]

phase portrait k = 1/2

ParametricPlot3D[f[t, 0.99], {t, 0, 20}, AspectRatio -> 1]

phase portrait k = 1/2

ParametricPlot3D[f[t, 0.01], {t, 0, 20}, AspectRatio -> 1]

phase portrait k = 1/2

Note that this last plot is nearly flat because the modulus is small and so z is nearly constant. The small modulus also makes the phase portrait nearly circular because x is approximately sine and y is approximately cosine.

[1] Kenneth R. Meyer. Jacobi Elliptic Functions from a Dynamical Systems Point of View. The American Mathematical Monthly, Vol. 108, No. 8 (Oct., 2001), pp. 729-737

Commutative multiplication of triples

Posted on 20 June 2018 by John

The complex numbers make a field out of pairs of real numbers. The quaternions almost make a field out of four-tuples of numbers, though multiplication is not commutatative. Technically, quaternions form a division algebra.

Frobenius’s theorem says only real vector spaces that can be made into division algebras are the real numbers, complex numbers, and quaternions.

So can you multiply triples of real numbers? Sure. Cross product is one way, but there’s no “division” analog to cross product. For example, the cross product of any vector with itself is zero, so cross product has lots of zero divisors, non-zero elements that multiply to zero.

Frobenius tells us we can’t form triples of real numbers into a division algebra, but can we come closer to a division algebra than we get with cross products? It turns out we can! We can define a multiplication that is commutative, has an identity element, and has no zero divisors. But what’s missing? The multiplication doesn’t have the distributive property you’d have in a division algebra.

The multiplication operation is complicated to describe. For details see “A Commutative Multiplication of Number Triplets” by Frank R. Pfaff, The American Mathematical Monthly, Vol. 107, No. 2 (Feb., 2000), pp. 156-162.

Related posts:

Three things about dominoes

Posted on 19 June 2018 by John

Here are three things about dominoes, two easy and one more advanced.

Counting

First, how many pieces are there in a set of dominoes? A domino corresponds to an unordered pair of numbers from 0 to n. The most popular form has n = 6, but there are variations with other values of n. You can show that the number of dominoes is

${ n+1\choose 2} + n + 1$

This is because there are n+1 possible numbers (since blanks are a possibility) and each one is either a double or not. The number of ways to choose two distinct numbers is the binomial coefficient and the number of doubles is n+1.

Another way to look at this is that we are selecting two things from a set of n+1 things with replacement and so the number of possibilities is

$\left({ n+1\choose 2}\right) = {n+2 \choose 2}$

where the symbol on the left is Stanley’s symbol for selection with replacement.

In any case, there are 28 dominoes when n = 6, 55 when n = 9, and 91 when n = 12.

Magic squares

There are a couple ways to make a magic square of sorts from a set of dominoes. To read more about this, see this post.

Domino magic square

Tiling

How many ways can you cover an m by n chess board with dominoes? The answer turns out to be

$\sqrt{\prod_{k=1}^m \prod_{\ell=1}^n \left( 2\cos\left(\frac{\pi k}{m+1} \right) + 2i \cos\left(\frac{\pi \ell}{n+1} \right) \right)}$

See this post for details.

Magical learning

Posted on 16 June 2018 by John

I asked two questions on twitter yesterday. The previous post summarized the results for a question about books that I asked from my personal Twitter account.

This post will summarize the results of a question I asked from @AnalysisFact.

If a genie offered to give you a thorough understanding of one theorem, what theorem would you choose?
— Analysis Fact (@AnalysisFact) June 15, 2018

Here are the results by category. Some are not strictly theorems but rather topics or conjectures.

Computer science

Curry-Howard correspondence
P=NP or not
Collatz Conjecture

Logic

Cohen’s forcing theorem
Godel’s Incompleteness theorem
Continuum hypothesis
Zorn’s lemma

Algebra and number theory

The ABC conjecture (theorem?)
Prime number theorem.
Riemann hypothesis
Fundamental theorem of finite abelian groups
Classification of finite simple groups
Fermat’s last theorem, the unpublished Fermat version

Topology and differential geometry

Thurston’s geometrization conjecture
Gauss Bonnet theorem
de Rham’s theorem
Grothendieck Riemann Roch theorem

Analysis

Banach-Tarski theorem.
Stokes theorem
Carleson-Hunt theorem
The epsilon/delta definition of continuity
Universal approximation theorem

Differential equations

Navier–Stokes existence and smoothness postulate
The relativistic version of the Shrodinger equation
Atiyah-Singer index theorem

Mathematical physics

E = mc²
Noether’s theorem
Liouville’s theorem

Miscellaneous

Existence of general equilibrium prices
Farkas’ lemma
The graph minor theorem
Central limit theorem

Mischievous genie

A couple people picked up the fact that in folk stories, being granted a wish doesn’t usually turn out well.

M. Shah: uh oh. Is this one of those malicious genies that twists language used to make the wish so that you are granted some horrific wish?

Jumpy the Hat: You now understand every single thing about irrational numbers but it’s all wasted because you’re cursed to become NJ Wildberger and you don’t think they exist

M Shah: or you want to thoroughly understand some theorem about Weierstrass’s monster. But little did anyone know that Weierstrass actually did have a pet monster. And it ends up biting your head off because it doesn’t like other things that are continuous.

Books you’d like to have read

Posted on 15 June 2018 by John

I asked on Twitter today for books that people would like to have read, but don’t want to put in the time and effort to read.

What’s a book you would like to have read but don’t want to read?
— John D. Cook (@JohnDCook) June 15, 2018

Here are the responses I got, organized by category.

Literature:

Math, Science, and Software:

History and economics

Religion and philosophy

Misc:

Crucial Conversations

Low-rank matrix perturbations

Posted on 14 June 2018 by John

Here are a couple of linear algebra identities that can be very useful, but aren’t that widely known, somewhere between common knowledge and arcane. Neither result assumes any matrix has low rank, but their most common application, at least in my experience, is in the context of something of low rank added to something of full rank.

Sylvester’s determinant theorem

The first is Sylvester’s determinant theorem. If A is a n by k matrix and B is a k by n matrix, then

$\det(I_n + AB) = \det(I_k + BA)$

where I is an identity matrix whose dimension is given by its subscript. The theorem is true for any k, but it’s especially useful when k is small because the matrices on the right side are of size k. If k = 1, the right side is the determinant of a 1 by 1 matrix, i.e. just a number!

Woodbury matrix inversion lemma

The second is known as the matrix inversion lemma or Woodbury’s matrix identity. It says

$\left(A+UCV \right)^{-1} = A^{-1} - A^{-1}U \left(C^{-1}+VA^{-1}U \right)^{-1} VA^{-1}$

which is a lot to take in at once. We’ll unpack it a little at a time.

First of all, the matrices have whatever properties are necessary for the equation to make sense: A and C are square, invertible matrices, though possibly of different sizes.

Suppose A is n by n. Then U necessarily has n rows. Let k be the number of columns in U. Then C is k by k and V is k by n.

To simplify things a little, let C be the identity.

$\left(A+UV \right)^{-1} = A^{-1} - A^{-1}U \left(I+VA^{-1}U \right)^{-1} VA^{-1}$

Think of k as small, maybe even 1. Then UV is a low-rank perturbation of A. Suppose you have computed the inverse of A, or know something about the inverse of A, and want to know how that inverse changes when you change A by adding a rank k matrix to it.

To simplify things further, assume A is the identity matrix. Now the matrix inversion lemma reduces to something similar to Sylvester’s theorem above.

$\left(I+UV \right)^{-1} = I - U \left(I+VU \right)^{-1} V$

To make things clearer, we’ll add subscripts on the identity matrices as above.

$\left(I_n+UV \right)^{-1} = I_n - U \left(I_k+VU \right)^{-1} V$

If k is small, then the matrix between U and V on the right side is small. If k = 1, it’s just a number, and its inverse is just it’s reciprocal.

Almost prime generators and almost integers

Posted on 13 June 2018 by John

Here are two apparently unrelated things you may have seen before. The first is an observation going back to Euler that the polynomial

produces a long sequence of primes. Namely, the values are prime for n = 1, 2, 3, …, 40.

The second is that the number

$e^{\pi \sqrt{163}}$

is extraordinarily close to an integer. This number is known as Ramanujan’s constant. It differs from its nearest integer by 3 parts in 10³⁰. Ramanujan’s constant equals

262537412640768743.99999999999925…

There is a connection between these two facts: The polynomial

returns primes for n = 1, 2, 3, …, k-1 primes if 4k – 1 is a Heegner number, and

$e^{\pi \sqrt{d}}$

is almost an integer if d is a (large) Heegner number.

Source: The Book of Numbers by Conway and Guy.

Heegner numbers

So what’s a Heegner number and how many are there? An integer d is a Heegner number if the ring generated by appending √-d to the integers has unique factorization. There are nine such numbers:

1, 2, 3, 7, 11, 19, 43, 67, 163.

There’s deeper stuff going on here than I understand—modular forms, the j-function, etc.—so this post won’t explain everything. There’s something unsatisfying about saying something is “almost” an integer without quantifying. There’s a way to be more precise, but we won’t go there. Instead, we’ll just play with the results.

Mathematica computation

First we look at the claim that n² – n + k produces primes for n = 1 through k – 1 if 4k – 1 is a Heegner number. The Heegner numbers of the form 4k + 1 are 2, 3, 5, 11, and 17. The following code shows that the claim is true for these values of k.

k = {2, 3, 5, 11, 17}
claim[x_] := AllTrue[
  Table[n^2 - n + x, {n, x - 1}], 
  PrimeQ
]
AllTrue[k, claim]

This returns True, so the claim is true.

As for exp(π √d) being close to an integer, this apparently only true for the last three Heegner numbers.

h = {1, 2, 3, 7, 11, 19, 43, 67, 163}
For[i = 1, i < 10, i++, 
  Print[
    AccountingForm[
      N[
        Exp[ Pi Sqrt[ h[[i]] ] ], 
        31
      ]
    ]
  ]
]

(The function AccountingForm suppresses scientific notation, making it easier to see where the decimal portion of the number starts.)

Here are the results:

                23.1406926327792
                85.0196952232072
               230.7645883191458
              4071.9320952252610
             33506.1430655924387
            885479.7776801543194
         884736743.9997774660349
      147197952743.9999986624548
262537412640768743.9999999999993

I manually edited the output to align the decimal points and truncate the decimal places beyond that needed to show that the last number is not an integer.

US flag if California splits into three states

Posted on 13 June 2018 by John

There’s a proposal for California to split into three states. If that happens, what would happen to the US flag?

The US flag has had 13 stripes from the beginning, representing the first 13 states. The number of stars has increased over time as the number of states has increased. Currently there are 50 stars, arranged in alternating rows of 6 and 5 stars.

If California breaks into three states, how would we arrange 52 stars? One way would be alternating rows of 7 and 6. Here’s a quick mock up of a possible 52-star flag I just made:

52 star flag

The difference isn’t that noticeable.

What if California [1] were to secede from the US? We’ve had 49 states before, in the brief period between the beginning of Alaskan statehood and the beginning of Hawaiian statehood. During that time [1], the flag had 7 rows of 7 stars, but the rows were staggered as in the current flag, not lined up in columns as in the 48-star flag before it.

***

[1] I don’t know how many people seriously propose any state leaving the US. The last time a state tried the result was the bloodiest war in US history. There’s a group in Vermont that wants their state to withdraw from the US.

Texans talk about seceding from the union, and you’ll occasionally see SECEDE bumper stickers, but it’s a joke. Maybe a few people take it seriously, but certainly not many.

[2] There is a tradition dating back to 1818 that the flag only changes on July 4. So the period of the 49-star flag didn’t exactly coincide with the period of 49 states.

Partition numbers and Ramanujan’s approximation

Posted on 12 June 2018 by John

The partition function p(n) counts the number of ways n unlabeled things can be partitioned into non-empty sets. (Contrast with Bell numbers that count partitions of labeled things.)

There’s no simple expression for p(n), but Ramanujan discovered a fairly simple asymptotic approximation:

$p(n) \sim \frac{1}{4n\sqrt{3}} \exp(\pi \sqrt{2n/3})$

How accurate is this approximation? Here’s a little Matheamtica code to see.

p[n_] := PartitionsP[n]
approx[n_] := Exp[ Pi Sqrt[2 n/3]] / (4 n Sqrt[3])
relativeError[n_] := Abs[p[n] - approx[n]]/p[n]
ListLinePlot[Table[relativeError[n], {n, 100}]]

So for values of n around 100, the approximation is off by about 5%.

Since it’s an asymptotic approximation, the relative error decreases (albeit slowly, apparently) as n increases. The relative error for n = 1,000 is about 1.4% and the relative error for n = 1,000,000 is about 0.044%.

Update: After John Baez commented on the oscillation in the relative error I decided to go back and look at it more carefully. Do the oscillations end or do they just become too small to see?

To answer this, let’s plot the difference in consecutive terms.

relerr[a_, b_] := Abs[a - b]/b
t = Table[relerr[p[n], approx[n]], {n, 300}];
ListLinePlot[Table[ t[[n + 1]] - t[[n]], {n, 60}]]

first differences of relative error

The plot crosses back and forth across the zero line, indicating that the relative error alternately increases and decreases, but only up to a point. Past n = 25 the plot stays below the zero line; the sign changes in the first differences stop.

But now we see that the first differences themselves alternate! We can investigate the alternation in first differences by plotting second differences [1].

ListLinePlot[
    Table[ t[[n + 2]] - 2 t[[n + 1]] + t[[n]], 
    {n, 25, 120}]
]

first differences of relative error

So it appears that the second differences keep crossing the zero line for a lot longer, so far out that it’s hard to see. In fact the second differences become positive and stay positive after n = 120. But the second differences keep alternating, so you could look at third differences …

Related posts: Special numbers

[1] The code does a small algebraic simplification that might make some people scratch their heads. All it does is simplify

(t_n+2 – t_n+1) – (t_n+1 – t_n).

Perl as a better grep

Posted on 12 June 2018 by John

I like Perl’s pattern matching features more than Perl as a programming language. I’d like to take advantage of the former without having to go any deeper than necessary into the latter.

The book Minimal Perl is useful in this regard. It has chapters on Perl as a better grep, a better awk, a better sed, and a better find. While Perl is not easy to learn, it might be easier to learn a minimal subset of Perl than to learn each of the separate utilities it could potentially replace. I wrote about this a few years ago and have been thinking about it again recently.

Here I want to zoom in on Perl as a better grep. What’s the minimum Perl you need to know in order to use Perl to search files the way grep would?

By using Perl as your grep, you get to use Perl’s more extensive pattern matching features. Also, you get to use one regex syntax rather than wondering about the specifics of numerous regex dialects supported across various programs.

Let RE stand for a generic regular expression. To search a file foo.txt for lines containing the pattern RE, you could type

    perl -wln -e "/RE/ and print;" foo.txt

The Perl one-liner above requires more typing than using grep would, but you could wrap this code in a shell script if you’d like.

If you’d like to print lines that don’t match a regex, change the and to or:

    perl -wln -e "/RE/ or print;" foo.txt

By learning just a little Perl you can customize your search results. For example, if you’d like to just print the part of the line that matched the regex, not the entire line, you could modify the code above to

    perl -wln -e "/RE/ and print $&;" foo.txt

because $& is a special variable that holds the result of the latest match.

***

For daily tips on regular expressions, follow @RegexTip on Twitter.

Mathematics of Deep Note

Posted on 12 June 2018 by John

I just finished listening to the latest episode of Twenty Thousand Hertz, the story behind “Deep Note,” the THX logo sound.

There are a couple mathematical details of the sound that I’d like to explore here: random number generation, and especially Pythagorean tuning.

Random number generation

First is that part of the construction of the sound depended on a random number generator. The voices start in a random configuration and slowly reach the target D major chord at the end.

Apparently the random number generator was not seeded in a reproducible way. This was only mentioned toward the end of the show, and a teaser implies that they’ll go more into this in the next episode.

Pythagorean tuning

The other thing to mention is that the final chord is based on Pythagorean tuning, not the more familiar equal temperament.

The lowest note in the final chord is D1. (Here’s an explanation of musical pitch notation.) The other notes are D2, A2, D3, A3, D4, A4, D5, A5, D6, and F#6.

Octaves

Octave frequencies are a ratio of 2:1, so if D1 is tuned to 36 Hz, then D2 is 72 Hz, D3 is 144 Hz, D4 is 288 Hz, D5 is 576 Hz, and D6 is 1152 Hz.

Fifths

In Pythagorean tuning, fifths are in a ratio of 3:2. In equal temperament, a fifth is a ratio of 2^7/12 or 1.4983 [1], a little less than 3/2. So Pythagorean fifths are slightly bigger than equal temperament fifths. (I explain all this here.)

If D2 is 72 Hz, then A2 is 108 Hz. It follows that A3 would be 216 Hz, A4 would be 432 Hz (flatter than the famous A 440), and A5 would be 864 Hz.

Major thirds

The F#6 on top is the most interesting note. Pythagorean tuning is based on fifths being a ratio of 3:2, so how do you get the major third interval for the highest note? By going up by fifths 4 times from D4, i.e. D4 -> A4 -> E5 -> B5 -> F#6.

The frequency of F#6 would be 81/16 of the frequency of D4, or 1458 Hz.

The F#6 on top has a frequency 81/64 that of the D# below it. A Pythagorean major third is a ratio of 81/64 = 1.2656, whereas an equal temperament major third is f 2^4/12 or 1.2599 [2]. Pythagorean tuning makes more of a difference to thirds than it does to fifths.

A Pythagorean major third is sharper than a major third in equal temperament. Some describe Pythagorean major chords as brighter or sweeter than equal temperament chords. That the effect the composer was going for and why he chose Pythagorean tuning.

Detuning

Then after specifying the exact pitches for each note, the composer actually changed the pitches of the highest voices a little to make the chord sound fuller. This makes the three voices on each of the highest notes sound like three voices, not just one voice. Also, the chord shimmers a little bit because the random effects from the beginning of Deep Note never completely stop, they are just diminished over time.

[1] The exponent is 7/12 because a half step is 1/12 of an octave, and a fifth is 7 half steps.

[2] The exponent is 4/12 because a major third is 4 half steps.

Musical score above via THX Ltd on Twitter.

Stirling numbers, including negative arguments

Posted on 7 June 2018 by John

Stirling numbers are something like binomial coefficients. They come in two varieties, imaginatively called the first kind and second kind. Unfortunately it is the second kind that are simpler to describe and that come up more often in applications, so we’ll start there.

Stirling numbers of the second kind

The Stirling number of the second kind

$S_2(n,k) = \left\{ {n \atop k} \right\}$

counts the number of ways to partition a set of n elements into k non-empty subsets. The notation on the left is easier to use inline, and the subscript reminds us that we’re talking about Stirling numbers of the second kind. The notation on the right suggests that we’re dealing with something analogous to binomial coefficients, and the curly braces suggest this thing might have something to do with counting sets.

Since the nth Bell number counts the total number of ways to partition a set of n elements into any number of non-empty subsets, we have

$B_n = \sum_{k=1}^n \left\{ {n \atop k}\right\}$

Another connection to Bell is that S₂(n, k) is the sum of the coefficients in the partial exponential Bell polynomial B_{n, k}.

Stirling numbers of the first kind

The Stirling numbers of the first kind

$S_1(n,k) = \left[ {n \atop k} \right]$

count how many ways to partition a set into cycles rather than subsets. A cycle is a sort of ordered subset. The order of elements matters, but only a circular way. A cycle of size k is a way to place k items evenly around a circle, where two cycles are considered the same if you can rotate one into the other. So, for example, [1, 2, 3] and [2, 3, 1] represent the same cycle, but [1, 2, 3] and [1, 3, 2] represent different cycles.

Since a set with at least three elements can be arranged into multiple cycles, Stirling numbers of the first kind are greater than or equal to Stirling numbers of the second kind, given the same arguments.

Addition identities

We started out by saying Stirling numbers were like binomial coefficients, and here we show that they satisfy addition identities similar to binomial coefficients.

For binomial coefficients we have

${n \choose k} = {n - 1 \choose k} + {n-1 \choose k-1}.$

To see this, imagine we start with the set of the numbers 1 through n. How many ways can we select a subset of k items? We have selections that exclude the number 1 and selections that include the number 1. These correspond to the two terms on the right side of the identity.

The analogous identities for Stirling numbers are

$\left\{{n \atop k} \right\}= k \left\{ {n-1 \atop k} \right\} + \left\{ {n-1 \atop k-1} \right\}$

$\left[ {n \atop k} \right]= (n-1) \left[ {n-1 \atop k} \right] + \left[ {n-1 \atop k-1} \right]$
The combinatorial proofs of these identities are similar to the argument above for binomial coefficients. If you want to partition the numbers 1 through n into k subsets (or cycles), either 1 is in a subset (cycle) by itself or not.

More general arguments

Everything above has implicitly assumed n and k were positive, or at least non-negative, numbers. Let’s look first at how we might handle zero arguments, then negative arguments.

It works out best if we define S₁(0, k) and S₂(0, k) to be 1 if k is 0 and zero otherwise. Also we define S₁(n, 0) and S₂(n, 0) to be 1 if n is 0 and zero otherwise. (See Concrete Mathematics.)

When either n or k is zero the combinatoric interpretation is strained. If we let either be negative, there is no combinatorial interpretation, though it can still be useful to consider negative arguments, much like one does with binomial coefficients.

We can take the addition theorems above, which are theorems for non-negative arguments, and treat them as definitions for negative arguments. When we do, we get the following beautiful dual relationship between Stirling numbers of the first and second kinds:

$\left\{{n \atop k} \right\}= \left[ {-k \atop -n} \right]$

Tetrahedral numbers

Posted on 7 June 2018 by John

Start with the sequence of positive integers:

1, 2, 3, 4, …

Now take partial sums, the nth term of the new series being the sum of the first n terms of the previous series. This gives us the triangular numbers, so called because they count the number of coins at each level of a triangular arrangement:

1, 3, 6, 10, …

If we repeat this again we get the tetrahedral numbers, the number of balls on each level of a tetrahedral stack of balls.

1, 4, 10, 20, …

We can repeat this process and general define T_{n, d}, the nth tetrahedral number of dimension d, recursively. We define T_{n, 1} = n and for d > 1,

$T(n, d) = \sum_{i=1}^n T(i, d-1)$

This is just a formalization of the discussion above.

It turns out there’s a simple expression for tetrahedral number of all dimensions:

$T(n, d) = \binom{n+d-1}{d}$

Here’s a quick proof by induction. The theorem clearly holds when n = 1 or d = 1. Now assume it hold for all n < m and d < m.

$\begin{eqnarray*} T(n, m) &=& \sum_{i=1}^n T(n, m-1) \\ &=& T(n, m-1) + \sum_{i=1}^{n-1} T(n, m-1) \\ &=& T(n, m-1) + T(n-1, m) \\ &=& \binom{n+m-2}{m-1} + \binom{n+m-2}{m} \\ &=& \binom{n+m-1}{m} \end{eqnarray*}$

The last line follows from the binomial addition identity

$\binom{a}{b} = \binom{a-1}{b} + \binom{a-1}{b-1}$

It also turns out that T_{n, d} is the number of ways to select d things from a set of n with replacement.

Related posts:

Bell numbers

Posted on 5 June 2018 by John

The nth Bell number is the number of ways to partition a set of n labeled items. It’s also equal to the following sum.

$B_n = \frac{1}{e}\sum_{k=0}^\infty \frac{k^n}{k!}$

You may have to look at that sum twice to see it correctly. It looks a lot like the sum for eⁿ except the roles of k and n are reversed in the numerator.

The sum reminded me of the equation

$\frac{d}{de}e^x = x e^{x-1}$

that I blogged about a few weeks ago. It’s correct, but you may have to stare at it a little while to see why.

Incidentally, the nth Bell number is the nth moment of a Poisson random variable with mean 1.

There’s also a connection with Bell polynomials. The nth Bell number is the sum of the coefficients of the nth complete Bell polynomial. The nth Bell polynomial is sum over k of the partial exponential Bell polynomials B_n,k. The partial (exponential) Bell polynomials are defined here.

Computing Bell numbers

You can compute Bell numbers in Python with sympy.bell and in Mathematical with BellB. You can compute Bell numbers by hand, or write your own function in a language that doesn’t provide one, with the recursion relation B₀ = 1 and

$B_n = \sum_{k=0}^{n-1}{ n-1 \choose k}B_k$

for n > 0. Bell numbers grow very quickly, faster than exponential, so you’ll need extended precision integers if you want to compute very many Bell numbers.

For theoretical purposes, it is sometimes helpful to work with the exponential generating function of the Bell numbers. The nth Bell number is n! times the coefficient of xⁿ in exp(exp(x) – 1).

$\sum_{n=0}^\infty \frac{B_n}{n!} x^n = \exp(\exp(x) -1)$

An impractical but amusing way to compute Bell numbers would be via simulation, finding the expected value of nth powers of random draws from a Poisson distribution with mean 1.