MoreRSS

site iconJohn D. CookModify

I have decades of consulting experience helping companies solve complex problems involving applied math, statistics, and data privacy.
Please copy the RSS to your reader, or quickly subscribe to:

Inoreader Feedly Follow Feedbin Local Reader

Rss preview of Blog of John D. Cook

The code that didn’t break

2026-08-06 02:26:25

Last week I wrote a post on hiding cryptographic keys in decks of cards. I wrote some code for that post that shouldn’t work, but before fixing I noticed that it in fact did work.

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

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

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

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

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

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

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

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

math.log(1e309)

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

That explains the following behavior:

>>> 1e309 == 1e310
True

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

The post The code that didn’t break first appeared on John D. Cook.

Enumerating trees and circles

2026-08-05 22:49:35

A few days ago I wrote a post on counting rooted trees. That post looked at the sequence c(n) which counts the number of rooted trees with n nodes. Here one node is distinguished as the root, but the nodes below the root are not distinguished from each other; all that matters is how the nodes are connected.

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

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

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

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

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

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

Here are the nine rooted trees with five nodes:

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

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

The post Enumerating trees and circles first appeared on John D. Cook.

Mathematical alchemy

2026-08-04 21:18:58

After writing the previous post about metallic ratios, I thought about the analogy to alchemy and the attempt to make precious metals out of base metals.

When can you make one metallic ratio out of another? Can you make the golden ratio out of the lead ratio?

Before we can make gold out of lead, we have to say what lead is.

Defining metallic ratios

The metallic ratios M(n) can be defined several ways. The most interesting definition is the number whose continued fraction representation contains all ns. A more prosaic but more convenient definition is the larger number that equals its reciprocal plus n, which can be found using the quadratic formula.

M(n) = n + \cfrac{1}{n+\cfrac{1}{n+\cfrac{1}{n+\cdots}}} = \frac{n + \sqrt{n^2 + 4}}{2}

The golden ratio is M(1), the silver ratio is M(2), and the bronze ratio is M(3).

Gold from silver and bronze?

Can you make the golden ratio out of the silver and bronze ratios? Not by integer arithmetic. The golden ratio involves √5, the silver ratio √2 and the bronze ratio √13. No integer operations on the latter two radicals will produce the former, though you can come arbitrarily close.

Gold from lead

The metallic ratios for n > 3 don’t have standard names, but let’s call M(4) the lead ratio. Can you make the golden ratio out of the lead ratio? Yes you can:

M(1) = (M(4) − 1)/2.

General solution

In general, when can you make M(n) out of M(m)? In abstract terms the question is when the fields

ℚ(√(n² + 4))

and

ℚ(√(m² + 4))

are the same, i.e. when adjoining √(n² + 4) to the rational numbers gives the same field as adjoining √(m² + 4) to the rational numbers. This occurs if and only if

(n² + 4)/( + 4)

is the square of a rational number.

Bronze from copper and tin

Can you make bronze out of copper and tin? Yes, if you define M(36) to be the copper ratio and M(393) to be the tin ratio, because

(3² + 4)/(36² + 4) = (1/10)²

and

(3² + 4)/(292² + 4) = (1/109)².

The post Mathematical alchemy first appeared on John D. Cook.

Ratio of metallic ratios

2026-08-04 19:41:49

The golden ratio is the first and best known of the metallic ratios. I’ve written about the silver ratio a few times, most recently here. And I’ve mentioned the bronze ratio a couple times. The metallic ratios after bronze don’t have standard names.

The nth metallic ratio M(n) is the number whose continued fraction representation contains all ns.

n + \cfrac{1}{n+\cfrac{1}{n+\cfrac{1}{n+\cdots}}} = \frac{n + \sqrt{n^2 + 4}}{2}

When n = 1, 2, and 3 we get the gold, silver, and bronze ratios.

You can approximate any positive real number as a ratio of metallic ratios. To see this, note that for large n, M(n) is approximately n. For any positive rational number a/b,

\lim_{n\to\infty} \frac{M(na)}{M(nb)} = \frac{a}{b}

and so you can make M(na) / M(nb) as close to a/b as you like by taking n large enough. And since the rationals are dense in the reals, you can approximate any positive real number as close as you’d like.

Let’s look for metallic ratios whose ratios approximate π to within 0.001 with the following Python code.

from math import pi, sqrt

M = lambda n: 0.5*(n + sqrt(n**2 + 4))

for n in range(1, 100):
    a = round(pi*n)
    b = n
    r = M(a)/M(b)
    if abs(r - pi) < 0.001:
        print(a, b, r)

This shows

π ≈ M(132) / M(42) = 3.1412…

Could we find smaller numbers that work? The following code shows the answer is no.

k = 132 + 42
# loop over numbers whose sum is less than k
for n in range(1, k):
    for a in range(1, n):
        b = n - a
        r = M(a)/M(b)
        if abs(r - pi) < 0.001:
            print(a, b, r)
            exit()

Related posts

The post Ratio of metallic ratios first appeared on John D. Cook.

Holonomic functions

2026-08-03 04:47:53

Yesterday I wrote that a lot of the special functions that pop up in mathematical physics are solutions to second order linear differential equations with polynomial coefficients. More generally, holonomic functions are defined to be those functions that are the solutions to linear differential equations, of any order, with polynomial coefficients.

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, though 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 post Holonomic functions first appeared on John D. Cook.

Estimating a cumulative sum

2026-08-03 01:57:31

In this post I mentioned two series which I denoted t(n) and c(n). The former is the number of unlabeled rooted trees with n nodes. The latter is the cumulative sum of the former, i.e.

c(n) = t(1) + t(2) + t(3) + \cdots + t(n)

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

t(n) \sim C \frac{a^n}{n^{3/2}}

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.

\begin{align*} c(n) &= \sum_{n=1}^N t(n) \\ &\sim \sum_{n=1}^N C \frac{a^n}{n^{3/2}}\\ &= C \frac{a^N}{N^{3/2}} \sum_{k=0}^{N-1} a^{-k}\left(1 - \frac{k}{N} \right)^{-3/2} \\ &\sim C \frac{a^N}{N^{3/2}} \sum_{k=0}^\infty a^{-k} \\ &= C \frac{a^N}{N^{3/2}} \frac{a}{a-1} \\ &= C \frac{a^{N+1}}{(a-1)N^{3/2}} \end{align*}

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:

The post Estimating a cumulative sum first appeared on John D. Cook.