Most special functions are holonomic. To quantify that statement, I went through the special functions covered in Abramowitz and Stegun. The large majority are holonomic, thought some common functions like the gamma function are not holonomic.

This report goes through the functions in A&S. For those that are holonomic, it gives the differential equation that the function solves. The large majority of these equations are second order, but not all. And the coefficients are nearly always first or second order polynomials, rarely higher order.

The sequence c(n) is also the number of constraints on an n-step Runge-Kutta method; that’s how I became interested in it.

Now the t(n) sequence has been cataloged as OEIS A000081 and OEIS gives the asymptotic estimate of t(n) for large n as

where C = 0.4399… and α = 2.9557….

The cumulative sum of t(n), what I’ve called c(n), is also cataloged in OEIS, sequence number A087803. However, OEIS does not give an asymptotic estimate for this sequence. I’ll give one here.

(Update: After looking closer at the page for A087803 I see that there is an asymptotic formula, the same one derived here.)

The basis for my derivation is to assume the cumulative sum of the asymptotic estimates gives an asymptotic estimate of the cumulative sum. This is justified by the fact that the sequence is increasing rapidly and only the last few terms contribute much relatively to the sum.

The technique illustrated here would be applicable to the cumulative sum of other series whose asymptotic form is known.

Here’s code to visualize the rate of convergence.

import numpy as np
import matplotlib.pyplot as plt

# from https://oeis.org/A000081/b000081.txt
A000081 = [
    0,
    1,
    1,
    2,
    4,
    ...
    51384328351659326880337136395054298255277970,
]  
A087803 = np.cumsum(A000081)

def approx(n):
    C = 0.43992401257102530
    a = 2.95576528565199497
    return C*a**(n+1)*n**(-3/2)/(a - 1)

n = np.arange(len(A087803))
ratio = A087803/approx(n)

plt.plot(n[1:], ratio[1:])
plt.plot(n, 0*n + 1, '--')
plt.xlabel("$n$")
plt.ylabel("exact/approx")
plt.show()

Here’s the plot:

This class of equations seems like a mathematically natural topic, but why is it so important in applications? I did a PhD in differential equations without ever learning why. The theory of second order linear equations with polynomial coefficients is too complicated for undergraduate courses [0] and too well-established for graduate courses [1].

The explanation that I was missing can be found in the first chapter of [2]. The PDEs that are common in physics are separable in various coordinate systems, meaning that in these coordinate systems the PDEs reduce to ODEs. These ODEs either have polynomial coefficients, or there is a change of variables which makes the ODEs have polynomial coefficients.

See this writeup that looks at the Helmholtz and Laplace equations in 11 coordinate systems.

[0] You may see the simplest parts of the theory in a section on solving ODEs with power series. But textbooks don’t go very far for good reasons.

[1] Unfortunately, a lot of really useful topics are left out of the graduate curriculum because they’re too well understood to provide thesis topics. Or the problems that are still open have been open for so long that they’re likely too hard to be cracked by a graduate student.

[2] Gerhard Kristensson. Second Order Differential Equations: Special Functions and their Classification. Springer, 2010.

Counting the number of unlabeled rooted trees [1] with n nodes is a pure mathematics problem. Designing numerical methods for solving differential equations is an applied mathematics problem. And yet the two are closely linked.

Let t(n) be the number of distinct unlabeled rooted trees with n nodes. The diagram below shows that the first few terms of this sequence are 1, 1, 2, and 4.

In an earlier post I showed that designing a 4-stage explicit Runge-Kutta method required solving a system of 8 equations in 10 unknowns, leaving two degrees of freedom in the solutions.

The number of constraints c(s) needed to design an s-stage explicit RK method is equal to the number of rooted trees with up to s nodes:

c(s) = t(1) + t(2) + t(3) + … + t(s)

This is because there is a one-to-one correspondence between constraints on the nth derivative of an RK formula and rooted trees, and an s stage method has to satisfy the constraints of all stages up to s. In the example of the 4th order RK method, we have

c(4) = t(1) + t(2)  + t(3) + t(4) = 1 + 1 + 2 + 4 = 8.

The first few values [2] of t(n) are

1, 1, 2, 4, 9, 20, 48, 115, 286, 719, 1842, 4766, 12486, 32973, …

and so you can see that t(n) grows quickly. In fact, it grows exponentially [3].

However, the number of parameters in an s stage RK method is s(s + 1)/2. The number of equations grows exponentially and the number of variables grows only quadratically, so at some point you have more equations than variables. That’s already the case for s = 5 because you have 17 constraints on 15 variables. The system has a solution because symmetry considerations render some of the equations redundant.

A 10th order RK method requires 17 stages. (See the previous post for why the number of stages exceeds the order when the order is greater than 4.) Designing such a method would require solving over a million equations in 153 variables, and yet it can be done. [4]

Related posts

• Solving the RK4 equations
• Counting permutations with roots
• DNA alignment and Kings

[1] This is a slightly contradictory term. Unlabeled means the we don’t distinguish the nodes. But we do distinguish one node, namely the root.

[2] See OEIS A000081.

[2] Richard Otter proved in 1948 that the number of unlabeled rooted trees with n nodes is asymptotically C αⁿ / n^−3/2 where C = 0.4399… and α = 2.9557…. The cumulative sum is at least this large since Otter’s estimate gives the size of the last term in the sum.

[3] E. Hairer. A Runge-Kutta Method of Order 10. J. Inst. Maths Applics (1978) 21, 47-59

Before going any further, let’s back up and say what we mean by stages and by order.

Stages

The number of stages in an RK method to solve the equation

is the number of evaluations of the function f on the right-hand side. For example, the textbook RK4 method estimates the solution at each step by

where

which requires four stages, i.e. four evaluations of f.

Order

A differential equation solver is said to have order p if the local error, the error after one step of size h, is O(h^{p + 1}). Then after solving an ODE over a period of time T with N = T/h steps, the global error is O(h^p). So, for example, if p = 4, you would expect that cutting your step size h in half would cut your error at T by a factor of 16.

More stages than the order

John C. Butcher proved that an explicit RK method of order p requires s stages where s > p if p > 4.

An important example is the Dormand-Prince method. It is a version of RK that has order 5 and 7 stages. The clever thing about this method is that you can make a 4th order solver out of a subset of its function evaluations.

That means that after you’ve evaluated one step of the 5th order method, you can also evaluate a 4th order method essentially for free. And by comparing them, you can get a sense of the error. If the solutions given by the two methods are substantially different, you have probably taken too big a step and need to back up. If the two solutions essentially agree, you’re probably good to take the next step.

For an explict RK method to have order 5, 6, or 7 you need at least 6, 7, or 9 stages respectively.

These calculations, which are not reproduced in Kutta’s paper (they are however in Huen (1900)), are very tedious.

The calculations are a set of eight constraints that the parameters of a fourth-order Runge-Kutta method must satisfy. I wondered how well Mathematica might have done at assisting Mr. Huen in his “very tedious” calculations if it had been available in 1900.

I go into Runge-Kutta methods in this post. Here I’d like to concentrate on a step in the design of the methods, namely solving the set of equations alluded in the quote above.

The first thing to note is that there are 10 variables and only 8 equations, and so the solution is not fully determined. What we think of as the fourth order Runge-Kutta method is in fact a fourth order Runge-Kutta method.

One could argue that we should have b₂ = b₃ and c₂ = c₃. With these additional equations, the system of equations has a unique solution, and Mathematic finds it easily.

eqs = {
    b1 + b2 + b3 + b4 == 1,
    b2*c2 + b3*c3 + b4*c4 == 1/2,
    b2*c2^2 + b3*c3^2 + b4*c4^2 == 1/3,
    b3*a32*c2 + b4*(a42*c2 + a43*c3) == 1/6,
    b2*c2^3 + b3*c3^3 + b4*c4^3 == 1/4,
    b3*c3*a32*c2 + b4*c4*(a42*c2 + a43*c3) == 1/8,
    b3*a32*c2^2 + b4*(a42*c2^2 + a43*c3^2) == 1/12,
    b4*a43*a32*c2 == 1/24,
    b2 == b3,
    c2 == c3
};

vars = {b1, b2, b3, b4, c2, c3, c4, a32, a42, a43};

solution = Solve[eqs, vars]

This returns the parameters used for the version of Runge-Kutta presented in every textbook.

If you keep the requirement b₂ = b₃ but substitute the requirement 2c₂ = c₃ for c‘s Mathematica will return the coefficients for the so-called Runge-Kutta 3/8 rule. This method has some slight advantages by some criteria.

In 1951 Gill [2] discovered a fourth order Runge-Kutta rule optimized for running in extremely constrained computer hardware. It’s a strange method, with irrational parameters, but one that was a very clever response to the limitations of its time.

Update: See this post for a discussion of the parameters and constriants for higher-ordered RK methods.

Related posts

• Dormand and Prince
• RK4 as a fold
• Stiff differential equations

[1] Hairer, Nørsett, and Wanner. Solving Ordinary Differential Equations I: Nonstiff Problems. Springer-Verlag 1987.

[2] A. Gill. A process for the step-by-step integration of differential equations in an automatic digital computing machine. Proc. Cambridge Philos. Soc., vol 27, pp 95–108.

⌊log₂(n!)⌋ ≥ b

given b. That is, given a number of bits b, find the smallest value of n such that n! ≥ 2^b.

What the code got right

Looking back on the code in that post, there are a few changes I’d like to make. But first of all, I’d like to point out something the post does right: instead of trying to solve

Γ(y) = x

it solves

log Γ(y) = log x.

That’s why the argument to inverse_log_gamma is logarg. That makes the code useful for values of x that would far exceed the maximum floating point value, such as in the calculations for the previous post.

What I’d change

Rounding

The function inverse_factorial from the old post solves finds the closest integer solution. It would be better for it to return the solution without rounding and then let the user round result if they want to. In my calculations in the previous post, I wanted to take the floor, not round.

Newton’s method

The code in the previous post uses the bisection method. This method is very safe, and fast enough for my purposes, but it could be made faster. Newton’s method is faster, but it can be ill-behaved if you don’t start close enough to the solution.

It’s safe to use Newton’s method to invert log Γ for two reasons. First, you can get a good starting point based on Stirling’s approximation. Second, and more importantly, log Γ is convex. Newton’s method will converge from any starting point when applied to a convex function. A little caution is necessary because log Γ is not convex everywhere, but it is convex on the positive real axis.

Another difficulty with Newton’s method is that you need to supply the derivative of the function whose root you’re trying to find. But this isn’t an issue here because the derivative of log Γ is the digamma function, which is implemented in SciPy.

Tolerance

Finally, the previous code used the default tolerance for deciding when to stop refining the solution. The revised method lets the user specify tolerance. It provides a default value, but that default is visible in the function call, not hidden down in SciPy.

Revised code

Here’s the revised code.

from scipy.special import gammaln, digamma
from scipy.optimize import newton

def inverse_log_gamma(logarg, tol=1e-12):
    assert(logarg > 0)    
    x0 = logarg / log(logarg + 1) + 1 if logarg > 1 else 2.0
    def f(z): return gammaln(z) - logarg
    return newton(f, x0, fprime=digamma, tol=tol)

def inverse_factorial(logarg):
    g = inverse_log_gamma(logarg)
    return g - 1 

⌊log₂(52!)⌋ = 225.

Here ⌊x⌋ is x rounded down to the nearest integer.

If we want to store bigger keys, we’re going to need a bigger deck of cards.

Bitcoin

A Bitcoin key has 256 bits, which would require a deck of 58 cards. There is a card game called Zwicker that uses a deck of 58 cards, the usual 52 cards plus six jokers. So you could store a Bitcoin key in the permutation of a Zwicker deck.

You could also use a deck of 52 cards, plus 2 jokers, if you also consider orientation. 30 cards are rotationally symmetric, 22 are not, and neither are jokers. So, including two asymmetric jokers, you could add 24 additional bits. Permutations of a 54 card deck can encode 237 bits, and with 24 orientation bits, this is a total of 261 bits.

RSA

RSA key sizes vary, but 2048 and 3072 are common. A 2048-bit key would require a deck of 301 cards. Casinos often use a shoe of 312 cards, combining six decks of 52 cards, to deal Baccarat or Blackjack. However, casinos combine identical decks. If you were to combine six unique decks, you could store a 2048-bit key.

Storing a 3072-bit key would require a deck of 422 cards. You could make a deck of 432 cards by combining 8 distinguishable packs of 54 cards (52 + 2 jokers).

ML-KEM

ML-KEM is a proposed quantum-resistant replacement for RSA. As with RSA, key sizes for ML-KEM vary, the smallest being ML-KEM-512 with a key size of 1632 bytes, which equals 13056 bits. This would require a deck of 1442 cards. You could combine 28 distinct packs of 52 cards, but that’s unwieldy.

This illustrates one of the difficult trade-offs with post-quantum cryptography: key sizes are much bigger. If you wanted to create a deck of 1442 cards, you’d probably want to make your “cards” something other than standard playing cards. You’d want to use permutations of something else.

Verification

The following Python code verifies the calculations above.

from math import log2, factorial, floor

def capacity(cards):
    return floor(log2(factorial(cards)))

def verify(bits, cards):
    return capacity(cards) >= bits and capacity(cards-1) < bits

print(verify(237, 54))
print(verify(256, 58))
print(verify(2048, 301))
print(verify(3072, 422))
print(verify(1632*8, 1442))

For more on how I came up with the deck sizes, see the next post on computing the inverse factorial.

You could store a 225-bit key as a permutation of 52 cards because

log₂(52!) = 225.581.

But then how would you number permutations so you could go from a number to a particular permutation and later decode the permutation to a number? Is this even practical? For a small number n, you could encode a number k < n by enumerating the first k permutations of a set of n items, and you could decode by enumerating permutations until you find the one you have. But this is completely impractical for large n, such as n = 52.

The process of mapping permutation to an integer is called ranking, and the mapping from an integer to a permutation is called unranking. How efficiently can rankings and unrankings be calculated?

Let n be the number of symbols being permuted. Then there are simple algorithms for ranking and unranking with respect to lexicographical order that have complexity O(n²) and more sophisticated algorithms that have complexity O(n log n). There are also O(n) algorithms that do not preserve lexicographical order.

The Permutations class in SymPy has methods unrank_lex and rank to unrank and rank permutations according to lexicographical order.

The notation the Permutations class uses requires a little explanation. For example, suppose we unrank 2026.

>>> from sympy.combinatorics import Permutation
>>> Permutation.unrank_lex(52, 2026)
Permutation(45, 47, 51, 48, 46, 50)

The output is not a full list of 52 numbers in permuted order; it is only a cycle. The notation refers to the permutation that sends 45 to 47, 47 to 51, …, 50 to 45 and leaves everything else fixed.

If we rank the permutation given above, we get 2026 back.

>>> Permutation.rank(Permutation(45, 47, 51, 48, 46, 50))
2026

Note that we didn’t say how many elements (45, 47, 51, 48, 46, 50) is a permutation of. Because of lexicographical order, the rank would be the same whether we viewed this as a permutation of 52 objects or of more objects.

Now let’s do something larger. Let’s generate a 220-bit number and encode it as a permutation.

>>> n = random.getrandbits(225)
>>> a = Permutation.unrank_lex(52, n)
>>> n
40234719030664563684489051530416964877785781669439875437823431388841
>>> a
Permutation(0, 25, 32, 15, 8, 28)(1, 48, 34, 14, 10, 51, 38, 31, 21, 5, 42, 47, 29, 26, 46, 30, 50, 49, 37, 22, 18, 23)(2, 45, 17, 20, 36, 40, 11, 4, 7, 41, 33, 3, 43, 44, 19, 16, 35, 39, 12, 6, 9)
>>> Permutation.rank(a) == n
True

Now just for fun, let’s display the permutation above applied to a standard (French) deck of 52 cards. As explained here, symbols associated with these cards have a range of Unicode values. By printing these values, we can visualize the permuted deck.

Here’s the code that made the image above.

spades = list(range(0x1F0A1, 0x1F0AF))
spades.remove(0x1F0AC) # take out the knight
cards = [s + 16*i for s in spades for i in range(4)]

a = Permutation.unrank_lex(52, n)
p = a(cards)

for i in range(4):
    for j in range(13):
        print(chr(p[13*i + j]), end="")
    print()

The code above is plenty fast, but Permutation has methods rank_nonlex and unrank_nonlex that run in O(n) time, which could be useful for n much larger than 52.

I wanted to say more about this, and look at implementing the same code in SymPy. I was curious how well SymPy would do because I’ve noticed that LLMs often generate SymPy code since it’s an open source CAS.

Wilf [1] describes the infinite product above as the exponential generating function (egf) of f(n, k), the number of permutations of n objects that have a kth root. Since egfs have a n! term in the denominator, this is also the ordinary generating function (ogf) of the probability that a randomly chosen permutation on n objects has a kth root.

My first attempt at using Mathematica to probe the generating function was

expq[x_, q_] := MittagLefflerE[q, x^q]	 
p[n_, k_] :=  SeriesCoefficient[	 
    Product[expq[x^m/m, GCD[m, k]], {m, 1, Infinity}], {x, 0, n}]

This hung forever when I tried to use it on a small example. I realized, but apparently Mathematica did not, that Infinity could be replaced by n since terms higher than n do not contribute to the coefficient of xⁿ. With that change, the code ran quickly.

This morning I tried converting the Mathematica code to Sympy; Claude did this in one shot. I also reproduced the table of f(n, k) values on page 150 of [1] to test the code. Since Wilf tabulated f(n, k), not f(n, k)/n!, I multiplied the results by n!.

Here is the output:

k = 2 [1, 1, 3, 12, 60, 270, 1890, 14280, 128520, 1096200]
k = 3 [1, 2, 4, 16, 80, 400, 2800, 22400, 181440, 1814400]
k = 4 [1, 1, 3, 12, 60, 270, 1890, 13020, 117180, 1039500]
k = 5 [1, 2, 6, 24, 96, 576, 4032, 32256, 290304, 2612736]
k = 6 [1, 1, 1, 4, 40, 190, 1330, 8680, 52920, 340200]
k = 7 [1, 2, 6, 24, 120, 720, 4320, 34560, 311040, 3110400]

and here is the SymPy code. I edited the main but the rest is verbatim from Claude.

from sympy import symbols, gcd, factorial, Rational, S

x = symbols('x')

def expq_coeffs(m, q, n):
    """
    Truncated (degree <= n) series coefficients of
        expq(x**m/m, q) = MittagLefflerE(q, (x**m/m)**q)
    Since q is a positive integer:
        E_q(y^q) = sum_j y^(q*j) / (q*j)!
    with y = x**m/m, so the term of degree m*q*j has coefficient
        1 / ( m**(q*j) * (q*j)! ).
    Returns a list c[0..n] of coefficients.
    """
    c = [S.Zero] * (n + 1)
    j = 0
    while m * q * j <= n:
        deg = m * q * j
        c[deg] += Rational(1, m**(q * j) * factorial(q * j))
        j += 1
    return c

def poly_mult_trunc(a, b, n):
    """Multiply two series (lists of coeffs, index = degree) truncated to degree n."""
    c = [S.Zero] * (n + 1)
    for i, ai in enumerate(a):
        if ai == 0:
            continue
        max_j = n - i
        for j2 in range(max_j + 1):
            bj = b[j2]
            if bj != 0:
                c[i + j2] += ai * bj
    return c

def p(n, k):
    """
    SymPy equivalent of:
        expq[x_, q_] := MittagLefflerE[q, x^q]
        p[n_, k_] := SeriesCoefficient[
            Product[expq[x^m/m, GCD[m, k]], {m, 1, n}], {x, 0, n}]
    """
    result = [S.Zero] * (n + 1)
    result[0] = S.One
    for m in range(1, n + 1):
        q = gcd(m, k)
        factor = expq_coeffs(m, q, n)
        result = poly_mult_trunc(result, factor, n)
    return result[n]

# example
if __name__ == "__main__":
    for k in range(2, 8):
        print("k =", k, [factorial(n)*p(n, k) for n in range(1,11)])

[1] Herbert Wilf. Generatingfunctionology. Available online here.

CAFE_hex = 1100 1010 1111 1110_two

I imagine it’s less well known that you can do the same thing with floating point numbers.

I wanted to find the binary representation of a floating point number using Python, and discovered that it has no function to do this. However, there is a method on floats to show a hex representation. For example, here’s the hex representation of π.

>>> import math
>>> (math.pi).hex()
'0x1.921fb54442d18p+1'

Curiously, the p+k part at the end is an exponent of 2, not an exponent of 16. So after we convert 1.921fb54442d18 to binary, we’ll need to multiply by 2, i.e. move the fractional point one space to the right.

So first we convert 1.921fb54442d18_hex to binary by converting 1, 9, 2, etc. each to binary.

1.1001 0010 0001 1111 1011 0101 0100 0100 0100 0010 1101 0001 1000_two

Then after shifting the fraction point to account for the p+1 part we have

π = 11.001001000011111101101010100010001000010110100011000_two

You could use Python’s bin() function to convert the fractional part, interpreted as an integer, to hex, though you may need to pad with 0 bits. For example,

>>> (1.03).hex()
'0x1.07ae147ae147bp+0

>>> bin(0x7ae147ae147)
'0b1111010111000010100011110101110000101000111'

The binary representation of 1.03_ten is

1.000001111010111000010100011110101110000101000111_two

We added a total of five zero bits, four for the 0 after the fractional point and one for converting 7 to 0111_two.

If we let our n elements be the integers 0 through n − 1, then we can represent permutations by what they do to this list of numbers. In Python as a tuple of length n and compose permutations with the following function:

import itertools

def compose(sigma, tau):
    "Return the composition σ ∘ τ (apply τ first, then σ)."
    return tuple(sigma[j] for j in tau)

We can always construct permutations that have square roots by squaring a permutation. If we run the following code

tau = (3, 1, 4, 5, 2, 0)
sigma = compose(tau, tau)

we find σ = (5, 1, 2, 0, 4, 3), and by construction (3, 1, 4, 5, 2, 0) is a square root of &sigma, though it’s not the only one.

The following code shows that σ has four roots.

import itertools

def numroots(sigma):
    n = len(sigma)
    c = 0
    for tau in itertools.permutations(range(n)):
        if sigma == compose(tau, tau):
            c += 1
    return c

print( numroots(sigma) )
print( numroots( (1, 2, 3, 4, 5, 0) ) )

It also shows that the rotation (1, 2, 3, 4, 5. 0) has no roots.

The function numroots has runtime proportional to n! and so it’s not practical for large permutations. There is a theorem that says a permutation σ has a square root if and only if the number of cycles it has of every even length is even. See [1].

We can also define cubes and cube roots of permutations, and higher powers and roots.

How common is it for permutations to have square roots, or cube roots, etc.? If you pick a random permutation on n elements, what is the probability that it has a kth root?

This is a hard question in general, but it is equivalent to finding the coefficient of x^k in the infinite product

This is theorem 4.8.3 in [1]. This theorem was the motivation for writing about exp_q in the previous post.

Although the product is infinite, there’s no need to compute terms in the product that only contribute powers of x higher than you’re interested in. The following Mathematica code will compute the probability that a permutation on n elements has a kth root.

expq[x_, q_] := MittagLefflerE[q, x^q]	 
p[n_, k_] :=  SeriesCoefficient[	 
    Product[expq[x^m/m, GCD[m, k]], {m, 1, n}], {x, 0, n}]

So, for example, the probability that a permutation of 10 elements has a square root is 29/96.

[1] Herbert Wilf. Generatingfunctionology. Available online here.

In general,

The first sum uses Iverson’s bracket notation: a Boolean expression in brackets denotes the function that returns 1 when the expression is true and zero when it is false. Here the bracket equals 1 when q divides n and is zero otherwise.

Closed forms

Let ω = exp(2πi / q). Then

This lets us find closed-form expressions for exp_q(x). For example, when q = 4, ω = i and

Here’s a proof of the identity above:

In the proof we used the identity

which is important in deriving the properties of the discrete Fourier transform.

Differential equations

The first time I saw the function exp_q(x) was in differential equations, though I didn’t know at the time the function had a name.

When a course in differential equations gets to power series solutions, a common example or homework problem is to solve

for k = 3 or 4, i.e. to find a function that equals its third or fourth derivative.

If the initial conditions are

and

the unique solution to

is y(x) = exp_k(x).

Mathematica and Mittag-Leffler

Mathematica does not have a built-in function implementing exp_q(x), but it does have an implementation of the Mittag-Leffler function, and so thanks to a relation between this function and exp_q(x) you can implement the latter as

expq[x_, q_] := MittagLefflerE[q, x^q]

Combinatorics

The first time I saw the notation exp_q(x) was in combinatorics. I had intended to include an application from that book here, but I make that the topic for the next post.

Excel labels columns A through Z, then AA through AZ, then BA through BZ, etc. If this is base 26, then does A correspond to 0? That could work for A through Z, but then what about AA? Then you’d have to say the first A corresponds to 26 but the second A corresponds to 0.

Does Z correspond to 0? If so then the column numbers would be 1 through 25, followed by 0, then 27. And it would mean that columns ZA through ZZ are the same as A through Z.

In fact nothing in Excel column labeling corresponds to 0. The labels cannot be interpreted as a positional number system.

There’s a name for this kind of number system: bijective base 26. The concept extends generally to bijective base b for any positive integer b. The idea is ancient, but the name was coined recently. It has also been called k-adic numbering. For most of history it didn’t have a name.

The motivation behind the name bijective base b is that there is a bijection (a one-to-one correspondence) between these symbols and positive integers; there’s no possibility of leading zeros that would keep the mapping from being a bijection, unlike say 7 and 07 representing the same number.

Excel limits

Before 2007, an Excel file could have a maximum of 2⁸ = 256 columns, and so the largest column label was IV.

Then in 2007 the column limit was increased to 2¹⁴ = 16,384 and the largest column label is XFD.

Conversion code

Converting from column labels to integers is easy; going the other way is a little more complicated.

letter_to_ordinal = lambda c: ord(c) - ord('A') + 1
ordinal_to_letter = lambda n: chr(ord('A') + n - 1)

def label_to_num(label):
    label = label.upper()
    n = 0
    for c in label:
        n = n*26 + letter_to_ordinal(c)
    return n

def num_to_label(n):
    letters = []
    while n > 0:
        n, remainder = divmod(n - 1, 26)
        letters.append(ordinal_to_letter(remainder + 1))
    return ''.join(reversed(letters))

Here’s an online calculator based on the code above.

Tests

The following code verifies the assertions above about the maximum number of Excel columns over time.

assert(num_to_label(256) == "IV")
assert(label_to_num("IV") == 256)

assert(num_to_label(2**14) == "XFD")
assert(label_to_num("XFD") == 2**14)

The conversion routines are not limited to actual Excel labels but work for arbitrarily large integers and bijective base 26 representations. For example, the following code shows that the bijective base 26 representation of Avogadro’s number is MUAEKAUDYDXEWOSDD.

avogadro = 602_214_076_000_000_000_000_000
assert(label_to_num(num_to_label(avogadro)) == avogadro)
print(num_to_label(avogadro))

Related posts

• Radix conversion calculator
• When is floating point radix conversion exact?
• Excel, R, and Unicode

Suppose α and β are two positive numbers such that α/β is irrational. In the previous post, α = 1 and β = π/180. Then the function

f(x) = sin(αx) − sin(βx)

is almost periodic: it is not periodic, but it comes close to being periodic, as close as you’d like provided you’re willing to look over a sufficiently long range of x‘s.

The identity

sin(αx) − sin(βx) = 2 cos((α + β)x/2) sin((α − β)x/2)

shows that f(x) is the product of two periodic functions but is not periodic itself. The periods of the cosine and sine above never coincide because the ratio of their frequencies is irrational.

The zeros of f are not periodic, though they can be divided into two subsequences that are periodic.

But sometimes it doesn’t matter. A trivial example is when x = 0. A more interesting example is

x = 180π/(180 + π) = 3.08770208….

For that value of x,

sin(x) = sin(x°).

In this article I’ll use the common convention of using radians by default and denoting degrees with ° as above.

Note that

x = πx°/180

and so we are interested in solutions to the equation

sin(x) = sin(πx/180)

Now two angles A and B have the same sine if they differ by a multiple of 2π, or if they’re supplementary (i.e. A = π − B), or both. To put it another way, if A and B have the same sine, they are either equal mod 2π or supplementary mod 2π. This means that

sin(x) = sin(πx/180)

if and only if

x = πx/180 + 2πk

or

x = π − πx/180 + 2πk

for some integer k.

Therefore all solutions have the form

x = 360πk/(180 − π)

or

x = 180π(2k + 1)/(180 + π).

Alternative solution

The derivation above is correct, but it occurred to me later that a simpler argument would be to use the identity

sin(A) − sin(B) = 2 cos((A + B)/2) sin((A − B)/2).

Thus A and B have the same sine if

cos((A + B)/2) = 0

or if

sin((A − B)/2) = 0.

These two possibilities correspond to the two families of solutions above.

Density

When reduced modulo 2π, both families are dense in [0, 2π]. This means that for every y in [−1, 1], there is a number x such that

sin(x) = sin(x°) ≈ y

and we can make the approximation as good as we’d like.

Example 1

For example, today is July 22, so let’s set y = 0.722. We’d like to find a value of x such that the sine of x radians and the sine of x degrees both approximately equal 0.722. And let’s say our approximation tolerance is ε = 0.0001.

We can search for a value of x in the first family of solutions by looking for a value of k with

| sin(360πk/(180 − π)) − 0.722 | < 0.0001

and the smallest such k is 96343 and so

x = 360×96343 π/(180 − π) = 616093.78713621…

will do, and sin(x) = 0.72191…

Example 2

Now let’s set y = 0.2026 and look for a solution in the other family of solutions, and this time let’s set ε = 10⁻⁶. The smallest value of k such that

| sin(180π(2k + 1)/(180 + π)) − 0.2026 | < 10⁻⁶

is k = 741141. Then

sin( 4576848.310950611 ) = sin( 4576848.310950611° ) = 0.202600139…

Suppose you have three numbers with uncertain values with a known, or at least purported, sum. The first number could be 31, 41, or 59; the second could be either 26 or 53; the last could be 58, 97, 93, or 23.

The following code enumerates all 3 × 2 × 4 = 24 possibilities and prints their sums.

from itertools import product

# Example input
possibilities = [(31, 41, 59), (26, 53), (58, 97, 93, 23)]

for combo in product(*possibilities):
    total = sum(combo) 
    print(f"Combination {combo} sums to: {total}")

In this example all the sums are unique, though of course that might not happen in practice. If, for example, you know the sum is 187, you know the three numbers were 41, 53, and 93. If the reported sum is 200, you know some assumption has been violated because none of the possible choices add up to 200.

More forensics posts

• Make up your own rules of probability
• Identifiable to man or machine?
• Metadata in photos
• Changing one character in a PDF
• Bits of information in age or birthday

I was curious whether most mathematicians were trying to prove or disprove the conjecture, so I asked Claude.

Before a counterexample to the Jacobian conjecture was found, did most mathematicians believe it was true or false?

Claude’s response was

The premise of this question isn’t quite right — no counterexample to the Jacobian conjecture has been found. It remains an open problem in mathematics: no one has proven it true, and no one has found a counterexample disproving it. … If you encountered a claim that a counterexample was found, do you have a source for that? I’d be happy to look into it, since that would actually be a major result in algebraic geometry if true.

Of course Claude doesn’t know that it solved the conjecture. It didn’t even solve the conjecture. It was an inanimate tool in the hand of a mathematician, just like a piece of chalk or a dry erase marker.

The middle part of Claude’s response was that mathematicians are (were) divided on whether the conjecture is true. So it was not like the Riemann hypothesis, which most people believe to be true, or the P = NP conjecture, which most people believe to be false.

Now what is the Jacobian conjecture? It says that a polynomial function from ℝⁿ to ℝⁿ with constant, non-zero Jacobian determinant has a polynomial inverse. (The conjecture was stated more generally for fields of characteristic 0, in which the derivatives defining the Jacobian would have to be defined algebraically, not in terms of limits.)

Alpöge came up with a counterexample, a polynomial function from ℝ³ to ℝ³ with constant Jacobian determinant −2. The function is

It’s a tedious but simple calculus exercise to show that the determinant equals −2 everywhere. The inverse function theorem says that a function is locally invertible at any point where the Jacobian determinant is non-zero, so Alpöge’s function is locally invertible everywhere.

However, the function takes on some values more than once. For example, (0, 0, −1/4) and (1, −3/2, 13/2) both map to (−1/4, 0, 0). Therefore the function is not invertible globally. So not only does the function not have a polynomial inverse, it doesn’t have an inverse even if you allow non-polynomial functions.

Alpöge’s counterexample disproves the Jacobian conjecture for n = 3. It can trivially be extended to all n > 3 by defining the function to be Alpöge’s function for three variables and the identity for the rest. The conjecture remains open for n = 2.

V = 4πr³ / 3

and the surface area is

A = 4πr²

and so the ratio of volume to area is

V / A = r / 3.

Surprisingly, the same ratio holds for all regular solids if r is the radius of the largest sphere that can be inscribed inside the regular solid.

For example, if the edge of a cube is a, then r = a/2. The volume is 8r³, the area is 24r², and the ratio is r/3.

The relationship between edge length and radius, and between radius and volume, is more complicated for the four other regular solids (tetrahedron, octahedron, dodecahedron, and icosahedron). However, in each case the ratio of volume to area is r/3.

The proof is surprisingly simple. Pick a face and form a pyramid by connecting each face vertex to the center of the inscribed sphere. The pyramid has height r and volume equal to B/3 where B is the area of the base. If the regular solid has f faces, the volume of the solid is fBr / 3 and the area is fB. So the ratio of volume to area is r/3.

The theorem generalizes to n > 3 dimensions. The formula for the volume of a pyramid in n dimensions is Bh/n where B is the (n − 1)-dimensional volume of the base, and so the ratio of n-dimensional volume of a regular solid to (n − 1)-dimensional volume of its boundary is r/n.

Here’s the problem, a variation on the n-queens problem formulated by Martin Gardner.

Can you place five white queens and three black queens on a 5 × 5 board so that no queen of one color attacks a queen of another color?

Here’s the prompt I have Grok 4.5.

Generate SWI Prolog code to solve the following chess puzzle. On a 5 by 5 chessboard, place 5 white queens and 3 black queens so that no queen is attacking a queen of the other color. Have the code enumerate all solutions.

The code ran correctly, giving all eight solutions.

I then asked Grok

Now generate Lean 4 code to solve the same puzzle.

The code did not run on the first or second attempts but ran correctly on the third iteration. That’s better than I’ve seen before. It’s not an apples-to-apples comparison, but it has usually taken many iterations with Claude before Lean code would run.

Solutions

There are eight solutions, each being a flip or a rotation of one fundamental solution. Here’s a representative.

SWI Prolog code

:- use_module(library(clpfd)).

%% Enumerate all ways to place 5 white queens and 3 black queens
%% on a 5x5 board so that no queen attacks a queen of the opposite colour.
%% Same-colour queens are allowed to attack each other.
%%
%% Positions are numbered 1..25 in row-major order:
%%   1  2  3  4  5
%%   6  7  8  9 10
%%  11 12 13 14 15
%%  16 17 18 19 20
%%  21 22 23 24 25

solve :-
    solution(Whites, Blacks),
    print_solution(Whites, Blacks),
    fail.
solve :-
    writeln('No more solutions.').

solution(Whites, Blacks) :-
    length(Whites, 5),
    length(Blacks, 3),
    append(Whites, Blacks, All),
    All ins 1..25,
    all_distinct(All),
    % canonical order so each set of positions is generated once
    Whites = [W1,W2,W3,W4,W5],
    W1 #< W2, W2 #< W3, W3 #< W4, W4 #< W5,
    Blacks = [B1,B2,B3],
    B1 #< B2, B2 #< B3, % every black is safe from every white maplist(safe_from_all(Whites), Blacks), label(All). safe_from_all(Whites, Black) :- maplist(no_attack(Black), Whites). %% Two queens do *not* attack each other. no_attack(P, Q) :- R1 #= (P-1) div 5 + 1, C1 #= (P-1) mod 5 + 1, R2 #= (Q-1) div 5 + 1, C2 #= (Q-1) mod 5 + 1, R1 #\= R2, % different row C1 #\= C2, % different column abs(R1-R2) #\= abs(C1-C2). % different diagonal print_solution(Whites, Blacks) :- format('Whites: ~w~n', [Whites]), format('Blacks: ~w~n', [Blacks]), forall(between(1, 5, R), ( forall(between(1, 5, C), ( Pos is (R-1)*5 + C, ( memberchk(Pos, Whites) -> write('W ')
                      ; memberchk(Pos, Blacks) -> write('B ')
                      ; write('. ')
                      )
                    )),
             nl )),
    nl.

Lean 4 code

-- Lean 4 solution for the 5×5 black-and-white queens puzzle.
-- 5 white queens + 3 black queens; no opposite-colour attacks.

def absDiff (a b : Nat) : Nat :=
  if a ≥ b then a - b else b - a

def attacks (p q : Nat × Nat) : Bool :=
  p.1 == q.1 || p.2 == q.2 || absDiff p.1 q.1 == absDiff p.2 q.2

/-- All 25 squares (0-based, row-major). -/
def allPositions : List (Nat × Nat) :=
  (List.range 5).flatMap fun r =>
    (List.range 5).map fun c => (r, c)

/-- Combinations of length exactly `k` (order-preserving). -/
def combinations {α : Type} : Nat → List α → List (List α)
  | 0,     _      => [[]]
  | _,     []     => []
  | k + 1, x::ys  =>
      (combinations k ys).map (fun zs => x :: zs) ++
      combinations (k + 1) ys

def isAttackedBy (pos : Nat × Nat) (queens : List (Nat × Nat)) : Bool :=
  queens.any (fun q => attacks q pos)

/-- Every solution as a pair (white squares, black squares). -/
def solve : List (List (Nat × Nat) × List (Nat × Nat)) :=
  (combinations 3 allPositions).foldl (fun acc blacks =>
    let forbidden := allPositions.filter fun p =>
      blacks.contains p || isAttackedBy p blacks
    let available := allPositions.filter fun p =>
      !(forbidden.contains p)
    acc ++ (combinations 5 available).map fun whites => (whites, blacks)
  ) []

/-- Pretty-print one board. -/
def printBoard (whites blacks : List (Nat × Nat)) : IO Unit := do
  for r in List.range 5 do
    let mut row := ""
    for c in List.range 5 do
      let p := (r, c)
      row := row ++
        (if whites.contains p then "W "
         else if blacks.contains p then "B "
         else ". ")
    IO.println row
  IO.println ""

/-- Enumerate and display every solution. -/
def main : IO Unit := do
  let sols := solve
  IO.println s!"Found {sols.length} solutions.\n"
  for p in (List.range sols.length).zip sols do
    let i := p.1
    let (whites, blacks) := p.2
    IO.println s!"Solution {i + 1}:"
    printBoard whites blacks

#eval main

grep -w -F -o -f icd10codes.txt notes.txt

Now you might want to combine your list of words into a singular regular expression, for efficiency or possibly for some other reason. Apparently ripgrep does this because when I tried replacing grep with ripgrep in the command above I got an error saying “Compiled regex exceeds size limit of 104857600 bytes.”

Beating brute force

Say you wanted to search for the strings “bluecross”, “blueshield”, and “bluey”. You could simply form the brute force regular expression

bluecross|blueshied|bluey

but that doesn’t take advantage of the fact that all three strings begin with “blue.” A smaller regular expression would be

blue(shield|cross|y)

Finding the shortest regular expression that matches a list of words is a hard problem, but finding a regular expression that’s shorter than brute force is not. The Python package trieregex will do this. According to the documentation,

trieregex creates efficient regular expressions (regexes) by storing a list of words in a trie structure, and translating the trie into a more compact pattern.

Let’s try our blue example with trieregex.

import re
from trieregex import TrieRegEx as TRE

words = ['bluecross', 'blueshield', 'bluey']
tre = TRE(*words) 
print(tre.regex())

This produces the same regular expression as above, except it adds ?: to make the parentheses non-capturing.

blue(?:shield|cross|y)

Prefixes versus suffixes

The library builds a trie data structure using common prefixes. That works well in the example above, but the result is disappointing when we have common suffixes rather than common prefixes. The following code

words = ['javascript', 'typescript']
tre = TRE(*words) 
print(tre.regex())

produces the regular expression

(?:javascript|typescript)

which is no better than brute force, whereas we might have hoped for

(?:java|type)script

HCPCS examples

As mentioned at the top of the post, ripgrep failed to search on a list of ICD-10 codes. The list of HCPCS codes is about 10x smaller, and more compressible. Ripgrep was able to fit all HCPCS codes into a single regex and was able to search the test file much faster than grep. The command

grep -w -F -o -f hcpcs.txt notes.txt

took 73.426 seconds to execute, while the command

rg -w -F -o -f hcpsc.txt notes.txt

took 0.078 seconds, three orders of magnitude faster.

The following code will read a list of HCPCS codes from a file and create a regular expression.

tre = TRE()
with open('hcpcs.txt', 'r') as file:
    for line in file:
        tre.add(line.strip())
print(len(tre.regex()))

This shows that the resulting regular expression has 17,198 characters. The file of codes has 8725 five-character codes, so the regex compresses the code characters by roughly a ratio of 5 to 2.

Let p be an odd prime number. Then half the numbers from 1 through p − 1 are squares and half are not. That is, for half of numbers 1 ≤ k < p, the equation

x² = k mod p

has a solution. The traditional name for these numbers is “quadratic residues” but we can just say “squares” if the context is clear. So, for example, the numbers 1, 2, and 4 are squares mod 7, and the numbers 3, 5, and 6 are not.

If k is a square mod p we will call is a low square if 0 ≤ k < p/2 and a high square if p/2 < k < p.

Signatures

Now let p > 3 be a prime congruent to 3 mod 4. Add up all the low squares mod p and take the remainder mod p. Call this the signature of p. Here’s Python code to make this explicit.

from sympy import isprime, factorint, is_quad_residue

def signature(p):
    assert(p > 3)
    assert(isprime(p))
    assert(p % 4 == 3)
    s = 0
    for k in range(1, 1 + p//2):
        if is_quad_residue(k, p):
            s += k
    return s % p

Inverse signatures

Surprisingly, the signature of each p is unique. Given the signature of p, you can uniquely determine p, and in fact you can do so easily. I ran across this in a paper [1] that presented the results in the form of a parlor trick: have someone pick a prime p such that p = 3 mod 4 and ask them to compute its signature, the sum of the low squares mod p. Then you can quickly tell them what their choice of p was.

Given a signature s, the corresponding prime p is the largest prime factor of 16s + 1.

Not only that,

p = (16s + 1)/m

where m is the smallest of the numbers {3, 7, 11, 15} such that the fraction above is a prime number. In term of Python code, both the following functions should invert the signature of p.

def inverse_signature1(s):
    n = 16*s + 1
    return max(factorint(n).keys())

def inverse_signature2(s):
    n = 16*s + 1
    for m in [3, 7, 11, 15]:
        if n % m == 0 and isprime(n // m):
            return n // m

The following code demonstrates that this is the case for numbers less than 1,000.

for n in range(7, 1000, 4):
    if isprime(n):
        s = signature(n)
        assert(n == inverse_signature1(s))
        assert(n == inverse_signature2(s))        

[1] David M. Bloom. A Quadratic Residues Parlor Trick. Mathematics Magazine, Vol. 71, No. 3 (Jun., 1998), pp. 201–203.

Here’s the treemap for HCPCS codes.

And here’s the treemap for ICD-10 codes.

The sizes of the squares are proportional to the number of codes beginning with that letter. Note that they are not necessarily proportional to how often codes are used.

The HCPCS map omits R and U because these are tiny relative to the rest. The ICD-10 map omits U for the same reason.

Here’s the code that was used to create the HCPCS map.

import matplotlib.pyplot as plt
import squarify

# HCPCS
data = {
    "G": 2010,
    "J": 1232,
    "L": 940,
    "A": 862,
    "E": 671,
    "Q": 639,
    "C": 619,
    "S": 533,
    "M": 506,
    "V": 212,
    "K": 175,
    "T": 114,
    "H": 94,
    "P": 59,
    "B": 51,
  # "U": 5,
  # "R": 3,
}

labels = list(data.keys())
sizes = list(data.values())

# Labels are just the letters (no counts)
display_labels = labels

# Color map — one distinct color per box
colors = plt.cm.tab20.colors[: len(labels)]

fig, ax = plt.subplots(figsize=(12, 8))
squarify.plot(
    sizes=sizes,
    label=display_labels,
    color=colors,
    alpha=0.85,
    ax=ax,
    text_kwargs={"fontsize": 30, "weight": "bold"},
    pad=True,
)
ax.axis("off")

plt.tight_layout()
plt.savefig("treemap.png", dpi=72)
plt.show()

The code to create the ICD-10 map differs only in its data.

# ICD-10
data = {
    "S": 31052, 
    "T": 10090, 
    "M":  6665, 
    "V":  4086, 
    "H":  3330, 
    "O":  2437, 
    "Y":  1590, 
    "I":  1427, 
    "Z":  1411, 
    "W":  1290, 
    "C":  1226, 
    "L":  1000, 
    "E":   971, 
    "Q":   894, 
    "F":   871, 
    "K":   857, 
    "N":   836, 
    "D":   824, 
    "R":   773, 
    "G":   700, 
    "A":   573, 
    "X":   495, 
    "B":   495, 
    "P":   463, 
    "J":   360, 
  # "U":     3,
}  

HCPCS stands for Healthcare Common Procedure Coding System, and is pronounced “hick picks.” When most people say HCPCS, they technically mean HCPCS Level II, and that’s what I mean here.

The format of a HCPCS code is simple: one letter and four digits. In regex terms,

    [A-Z]\d{4}

Not all letters are used, so you can get more specific and say

    [A-CEGHJ-MP-V][0-9]\d{4}

Some sources say no codes begin with U, but there are currently five codes that begin with U.

When I was doing some research on HCPCS codes recently using AI, I was told there is a D code for dentistry, but that was a hallucination.

HCPCS codes can also have modifiers. These consist of a letter and either a letter or digit:

    [A-Z][A-Z0-9]

Not all letters actually appear in modifiers—I, O, W, and Y are missing—so you could be more specific. At the time of writing there are 384 official modifiers.

Modifiers are often stored in a separate column in a database, but in text you’ll see a HCPCS code optionally followed by a dash and a modifier. So a regex to match HCPCS codes with possible modifiers would be

    [A-CEGHJ-MP-V][0-9]\d{4}(-[A-Z][A-Z0-9])?

This regex will have some false positives, but it should not have false negatives: every real HCPCS code should match.

Of course you could search against a complete list of HCPCS codes. This would be more accurate and slower. I did a test similar to the one in the previous post and found a search with the regex above took 20 milliseconds, while a search against the list of HCPCS codes took 46 seconds.

However, the regex searched for possible modifiers and the exhaustive search only looked for unmodified HCPCS codes. A complete list of HCPCS codes with possible modifiers would be tedious to create because some combinations of codes and modifiers make no sense. And I imagine that some combinations that would make sense are not used in practice.

Regular expressions usually do not exactly match what you’re looking for and nothing else. They have error false positives and false negatives. But they also have advantages, and context determines whether the advantages make the error rates tolerable.

The post mentioned above gave the following regular expression for ICD-10 diagnosis codes.

    [A-TV-Z][0-9][0-9AB]\.?[0-9A-TV-Z]{0,4}

As cryptic as this may look at first glance, it’s straight-forward to interpret. It says that an ICD-10 code

1. Begins with a capital letter, excluding U
2. Followed by a digit
3. Followed by a digit or A or B
4. Optionally followed by a period
5. Followed by up to 4 digits or capital letters, excluding U.

Speed

Now suppose you want to scan a text document for ICD-10 codes. One approach would be to use the regex above. Another would be to compare every alphanumeric sequence in the document to a list of ICD-10. Currently this list has 74,719 codes.

I tested both approaches on a 800kb text file. The regex search

    egrep -o '[A-TV-Z][0-9][0-9AB]\.?[0-9A-TV-Z]{0,4}' notes.txt

took 18 milliseconds. Searching against the list of codes

    grep -w -F -o -f icd10codes.txt notes.txt

took 386 seconds, about six and a half minutes or five orders of magnitude longer.

Error rates

The regex

    [A-TV-Z][0-9][0-9AB]\.?[0-9A-TV-Z]{0,4}

had a false negative rate of zero at the time it was written. I tested the regex against the current list of codes with the following command.

    egrep -v '[A-TV-Z][0-9][0-9AB]\.?[0-9A-TV-Z]{0,4}' icd10codes.txt

The -v flag reverses the sense of the search, reporting lines that do not match the regular expression. This returned three matches: U070, U071, and U099. So 3 out of 74,719 valid ICD-10 codes were reported as invalid.

Codes beginning with U are designated for provisional/emergency/special purposes, but these three have become essentially permanent. A change in the application of the ICD-10 standard caused an error in the regular expression.

But the change would also have caused an error in code that did an exhaustive search against the list of ICD-10 codes at the time. In fact, every new code not starting with U would also be reported in error. So the regex is actually more future-proof than an exhaustive search. Presumably the simplified regex

    [A-Z][0-9][0-9AB]\.?[0-9A-Z]{0,4}

will remain valid for the foreseeable future.

We’ve looked at false negatives. What about false positives? That depends on context. The false positive rate when searching medical notes is low: a word matching the regex above in a medical record is most likely an ICD-10 code. But the number of conceivable false positives is enormous. If you were searching a file of randomly generated alphanumeric text, the regex matches would overwhelmingly be false positives [1].

The number of strings matching

    [A-Z][0-9][0-9AB]\.?[0-9A-Z]{0,4}

would be

26 × 10 × 12 × (1 + 36 + 36² + 36³ + 36⁴) = 5,390,127,600.

Out of over five billion strings matching the regular expression, only around 75,000 are valid ICD-10 codes. So a naive theoretical calculation would say the false positive rate is 99.9986%, whereas in practice the false positive rate is very low, though there’s no way to say a priori exactly how low.

Related posts

• HCPCS (“hick pics”) codes
• NPI numbers and checksum
• HIPAA expert determination

[1] You could argue that all positives would be false positives in this context because you’re looking at noise. You couldn’t find an ICD code, though you could find a string of characters that coincides with an ICD code. That may sound like a pedantic distinction, but it matters in the context of evaluating deidentification quality: you want to find instances of PHI, not instances of strings that match the character pattern of PHI.

The ICD-10-CM standard is divided into 21 chapters, which generally correspond to the first letter of a code. However, a chapter may contain blocks beginning with more than one letter, and codes starting with a single letter, namely D, can span two chapters.

Here’s a diagram I made to visualize the relationship between chapters and initial letters of codes.

Notice there’s no letter U on the diagram. That’s because U is reserved for special/provisional codes.

Related posts

• Rare and strange ICD-10 codes
• Regex to match ICD-9 and ICD-10 codes
• New ICD-10 code weirdness

That led the previous post which looked at posterior means for three Bayesian models, showing how the posterior mean is a weighted average of the prior mean and the mean of the new data. The weights are precisions, which means something different for each model.

For the beta-binomial model, variance may increase when seeing unexpected data (details here), but precision always increases.

For the normal-normal model precision is the reciprocal of variance. Every new data point makes precision go up and posterior variance go down.

The Poisson-gamma model may be the most interesting. As stated in the previous post, if data has a Poisson distribution with parameter λ, and λ has a gamma(α₀, β₀) prior distribution, then the posterior distribution on λ after observing k events over time t has a gamma(α₀ + k, β₀ + t) posterior distribution. Therefore the posterior variance is

(α₀ + k) / (β₀ + t)².

Note the posterior variance is an increasing function of k and a decreasing function of t. This means that the posterior variance increases every time an event is observed, and it decreases quadratically between observations.

Here’s an illustration. I simulated data from a Poisson process with λ and used a gamma(1, 1) prior on λ. Here’s a plot of the posterior variance.

How much should new data change your prior beliefs? When prior judgment and new information are in conflict, which one should be given the benefit of the doubt?

Bayesian data models provide a framework for making such decisions quantitative and objective. The choice of a data model is somewhat subjective—whether it’s a Bayesian model or not—but given a Bayesian model, the rules for updating the representation of your beliefs are objective. As some put it, you “turn the Bayesian crank.” A likelihood model and a prior on parameters together specify how new data changes the prior distribution into a posterior distribution.

We will make this more concrete with three examples.

Normal-normal model

Suppose that data X has a normal distribution with unknown mean μ and known variance σ², and we assume that a priori μ has a normal distribution with mean μ₀ and variance σ₀².

After observing x, the posterior distribution on μ also has a normal distribution, but with a different mean and variance. Its mean is somewhere between the prior mean and x. We will ignore the change in the variance for this post.

The posterior mean of μ is

This equation becomes more understandable when we introduce precisions τ = 1/σ² and τ₀ = 1/σ₀².

Then we have

which you can read as saying the posterior mean is the weighted average of the prior mean and x, with the weights given by the precision. Intuitively, you take the weighted mean of your conclusions from previous data and new data, weighting the mean according to how much confidence you have in each.

Beta-binomial model

Now let’s switch over to a different data model. Now assume X is a binary random variable, with probability of success p and probability of failure 1 − p, and we assume p has a beta(a, b) distribution.

After observing s successes and f failures, the posterior mean of the distribution on p becomes

We can rewrite this as

This says that the posterior mean is the weighted average of the prior mean a/(a + b) and the mean of the data s/n. The weights are the prior effective sample size a + b and the sample size of the new data n. In this example (effective) sample size is playing the role that precision played in the normal-normal model above.

Gamma-Poisson model

Suppose data have a Poisson distribution with parameter λ, and λ has a gamma(α₀, β₀) prior distribution [1]. And suppose you observe k events over time t. Then the posterior distribution of λ given the data has a gamma(α₀ + k, β₀ + t) prior distribution and the mean of the posterior distribution is given by

As before, the posterior mean is a weighted average of the prior mean and new data, and the weights are interpretable as some sort of measure of confidence, namely time. The variable t is directly time and the parameter β₀ is sort of an effective time, just as a + b is an effective sample size for the beta distribution.

Common thread

In each example the posterior mean is the weighted average of the prior mean and the mean of the data, with the weights given by a precision. However, precision means something different in each example. In the normal-normal model, precision is the reciprocal of variance, but in the beta-binomial model precision is sample size and in the Poisson-gamma model precision is time.

What all three examples have in common is that they are conjugate models using distributions from the “exponential family” of probability distributions. In technical terms, precision is the multiplicative factor on the sufficient statistic in the exponent of the posterior kernel.

Related posts

• Does additional data always reduce posterior variance?
• Diagram of probability distribution relationships
• Bayesian statistics consulting

[1] There are multiple conventions for parameterizing the gamma distribution. Here we’re using the shape-rate parameterization, where the mean is α/β.

Gilbreath’s conjecture is simple, but it’s also kinda weird. As I wrote before,

Paul Erdős speculated that Gilbreath’s conjecture is true but it would be 200 years before anyone could prove it. I find Erdős’s conjecture more interesting than Gilbreath’s conjecture.

The conjecture is hard in a way that, say, solving a nasty-looking differential equation is not. Over the last three centuries, mathematics has developed quite a toolbox for solving differential equations. But Gilbreath’s conjecture is just odd enough that it’s not at all clear what kind of tool might be useful in approaching it.

Terence Tao has a new blog post announcing a paper that he and two coauthors wrote on a random model intended to mimic Gilbreath’s calculation on primes. This random model is more sophisticated than the little game Gilbreath was playing, but it’s also much more amenable to analysis by established techniques. Tao’s post gives a heuristic explanation for why Gilbreath’s conjecture is plausible, but then adds

However, it seems well beyond current technology to try to make these heuristics rigorous; even the first step … is far out of reach.

Related posts

• Moessner’s Magic
• Closed form solutions to PDEs

The theorem corresponding to the diagram is interesting, but I found reproducing the diagram more interesting.

The segment AB is a diameter and the line CD is perpendicular to the diameter.

Assume the outer circle is a unit circle. I guessed C = (cos(1), sin(1)) and made the following diagram.

I guessed the value of C by eyeballing it, but in retrospect this would have been a convenient value for the creator of the original diagram to have chosen.

Drawing the blue circle inscribed in the triangle was easy using the equations for the center and radius from this post. Drawing the other two circles, the green and orange circles, was harder. They are also inscribed circles, but not inscribed in a triangle. They’re inscribed in a three-sided figure with two perpendicular sides and a circular arc.

The radius r of the green circle is the distance from the center of the circle to each of its tangent lines. Also, the distance from the origin to the center of the circle must be 1 − r. This is enough information to set up a quadratic equation for r. The same reasoning applies to the orange circle.

The original diagram comes from [1] and the theorem it illustrates says the diameter of the blue circle equals the sum of the radii of the green and orange circles.

Python code

In case you’re interested, here’s the code that created the diagram.

#!/usr/bin/env -S uv run --script

# /// script
# dependencies = ["numpy", "matplotlib"]
# ///

import numpy as np
import matplotlib.pyplot as plt

def connect(A, B, color='gray'):
    plt.plot([A[0], B[0]], [A[1], B[1]], color=color, linewidth=2)

def circle(c, r, color='gray'):
    t = np.linspace(0, 2*np.pi)
    plt.plot(c[0] + r*np.cos(t), c[1] + r*np.sin(t), color=color, linewidth=2)

def quadratic(a, b, c):
    det = b**2 - 4*a*c
    return ((-b - det**0.5)/(2*a), (-b + det**0.5)/(2*a))

A = np.array([-1, 0])
B = np.array([ 1, 0])
C = np.array([np.cos(1), np.sin(1)])
a = np.linalg.norm(B - C)
b = np.linalg.norm(A - C)
c = np.linalg.norm(B - A)
s = (a + b + c)/2

circle([0,0], 1)
connect(A, B,)
connect(A, C)
connect(C, B)
connect(C, C*np.array([1, -1]))

center = (a*A + b*B + c*C)/(2*s)
radius = 0.5*a*b/s
circle(center, radius, 'C0')

Ex = C[0]
roots = quadratic(1, 2 + 2*Ex, Ex**2 - 1)
r = roots[1] # Smaller root is negaive
print(roots)
center = (r + Ex, -r)
circle(center, r, 'C1')

roots = quadratic(1, 2 - 2*Ex, Ex**2 - 1)
r = roots[1] # Smaller root is negaive
center = (Ex - r, -r)
circle(center, r, 'C2')

plt.gca().set_aspect("equal")
plt.axis("off")
plt.show()

[1] Leon Bankoff. A Geometrical Coincidence. Mathematics Magazine, Vol. 37, No. 5 (Nov., 1964), p. 324.