Solving hard problems

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

Numerical (in)stability of recurrence relations

Posted on 24 August 2026 by John

The previous post gave several examples of three-term recurrence relations for special functions. These relations can be computationally useful, but they have to be applied carefully.

Several years ago I wrote a post on stable and unstable recurrences. In that post I show that the stability of the recurrence relation for Bessel functions produces depends on which kind of Bessel function and which direction the recurrence is applied.

In the forward direction, computing higher order values from lower order values, works well for Bessel functions of the second kind Y_n but not for Bessel functions of the first kind J_n. In the reverse direction, the recurrence is stable for J_n but not for Y_n.

I didn’t explain in that post why this is. In this post I will.

Second order linear difference equations have two independent solutions, just like second order linear differential equations. For both kinds of equations, all solutions are linear combinations of the two solutions. Suppose one solution grows with n and the other decays. You may want to compute the decaying solution, but in doing so you might pick up a small component of the growing solution due to rounding error. This post illustrates this phenomena for differential equations, and this post illustrates it for difference equations.

When you look at a plot of Bessel functions in a text book, you’ll probably see a few plots of J_n(x) andY_n(x) for a few small values of n. The functions seem to behave roughly the same way, like sine and cosine. And that’s true, as functions of x.

But it’s not true for J_n(x) andY_n(x) as functions of n for fixed x. As n increases, J_n(x) decays to zero and Y_n(x) goes off to −∞.

That’s the source of numerical instability. And there will be similar instability problems for other recurrences where the ratios of the two independent solutions goes to zero or infinity as a function of n.

There are techniques for computing the solution that does not diverse, the so-called minimal solution, such as Miller’s algorithm mentioned here.

Three-term recurrences

Posted on 24 August 2026 by John

There many examples of families of functions where each function can be computed as a linear combination of the two previous terms

$f_{n+1}(x) = a(x) f_n(x) + b(x) f_{n-1}(x)$

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:

$\begin{align*} J_{\nu+1}(x) &= \frac{2\nu}{x}\,J_\nu(x) - J_{\nu-1}(x) \\ Y_{\nu+1}(x) &= \frac{2\nu}{x}\,Y_\nu(x) - Y_{\nu-1}(x) \end{align*}$

Modified Bessel functions of the first and second kind:

$\begin{align*} I_{\nu+1}(x) &= I_{\nu-1}(x) - \frac{2\nu}{x}\,I_\nu(x) \\ K_{\nu+1}(x) &= K_{\nu-1}(x) + \frac{2\nu}{x}\,K_\nu(x) \end{align*}$

Chebyshev polynomials of the first and second kind:

$\begin{align*} T_{n+1}(x) &= 2x\,T_n(x) - T_{n-1}(x) \\ U_{n+1}(x) &= 2x\,U_n(x) - U_{n-1}(x) \end{align*}$

Hermite polynomials (physicists’ convention):

$H_{n+1}(x) = 2x\,H_n(x) - 2n\,H_{n-1}(x)$

Legendre polynomials:

$P_{n+1}(x) = \frac{2n+1}{n+1}\,x\,P_n(x) - \frac{n}{n+1}\,P_{n-1}(x)$

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

The von Mises-Fisher distribution

Posted on 24 August 2026 by John

Probability density function must integrate to 1, and so if you know a density function up to a constant, the constant is determined.

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

$f(x \mid \mu, \kappa) = \frac{\exp(\kappa \cos(x - \mu))}{2\pi I_0(\kappa)}$

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

$f(\mathbf{x} \mid \boldsymbol{\mu}, \kappa) = C_{p}(\kappa) \exp \left( {\kappa \boldsymbol{\mu}^\mathsf{T} \mathbf{x} } \right)$

where the normalization constant C_p(κ) is

$C_{p}(\kappa)=\frac {\kappa^{p/2-1}} {(2\pi)^{p/2}I_{p/2-1}(\kappa)}$

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.

What exactly is modified about a modified Bessel function?

Posted on 23 August 2026 by John

Special functions often have arcane names that not very helpful without some context. The previous post goes into some reasons for this. This post will expand on a point at the end of the post about “modified” functions.

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,

$I_n(x) = i^{-n} J_n(ix)$

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

$\begin{align*} \sin(ix) &= i \sinh(x) \\ \cos(ix) &= \phantom{i} \cosh(x) \\ \tan(ix) &= \phantom{i} \tanh(x) \end{align*}$

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

$I_n(x) = i^{-n} J_n(ix)$

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

$I_\nu(x) = \exp(\nu\pi i/2) J_n(ix)$

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

$x^2 y'' + x y' + (x^2 - \nu^2) y = 0$

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

$x^2 y^{\prime \prime} + x y^{\prime} - (x^2 + \nu^2) y = 0$

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.

Why special function terminology is arcane

Posted on 23 August 2026 by John

Special functions are special because they’re useful. They can also be shrouded in arcane terminology. These two facts are related.

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.

The difference orbit inclination makes

Posted on 22 August 2026 by John

Suppose you wanted to find the distance between Earth and Mars over time. To first approximation, both planets orbit the sun in elliptic orbits in the same plane.

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).

Coming soon

Posted on 22 August 2026 by John

There’s a pizza shop near my home with a sign out front that says “Coming Soon.” When I drove by it this morning I thought about how you would model the time until an event happens that is “coming soon.”

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 = ∞.

How would you know whether an ancient culture had zero?

Posted on 21 August 2026 by John

A few weeks ago I wrote about the number system used in labeling spreadsheet columns. Labels run from A through Z, then AA through AZ, etc. This looks a lot like base 26, but it’s not quite the same. It has no analog of zero. If Z were like zero, Y would be followed by AZ. The Excel labeling system is not base 26, but what’s called bijective base 26.

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.

AI-generated ASCII diagrams

Posted on 20 August 2026 by John

I like AI-generated ASCII diagrams. Because nobody would ask AI to generate ASCII diagrams, and so, it’s congruous. I like incongruity [1].

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)     |
   +-------------+                    +-------------+

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

Big little hexagon

Posted on 18 August 2026 by John

A new paper just came out, The Maximum-Area Small Polygon Problem. The paper solves the problem of finding, for each n, the n-gon with diameter 1 and maximum area.

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.

Solving hard problems

Pure math perspective

Applied math perspective

Other modified functions

Related posts