John D. Cook

My post from yesterday on permutation roots ends with a Mathematica code for finding the probability that a permutation of n elements has a kth root. This is done by finding the coefficient of xn in the generating function 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…

It’s well known that you can convert the base 16 (hex) representation of an integer to the base 2 (binary) representation by simply converting each digit from hex to binary. For example, CAFEhex = 1100 1010 1111 1110two 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…

Let σ be a permutation on n elements. If there is a permutation τ such that applying τ twice has the same effect on the list of elements as applying σ once, we say σ = τ² and τ is a square root of σ. 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…

The function expq(x) is defined by taking the power series for exp(x) and keeping only the terms whose index is a multiple of q. For example, exp2(x) keeps only the even-numbered terms in the exponential power series and so equals cosh(x). 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…

I was working with a wide spreadsheet from a client the other day and I had to convert between Excel column labels and column numbers. I had never paid attention to how Excel labels columns and implicitly thought it was base 26 using letters rather than digits. But then I realized that’s not right. Excel labels columns A through Z, then AA trough 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…

This post takes a more abstract view of the previous post. That post looked at the concrete question of whether a number ever has the same sine in radians as in degrees. The relation between radians and degrees is irrelevant except that π/180 is an irrational number. 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…

Ordinarily the sine of x radians and the sine of x degrees are very different numbers. Having your calculator in radian mode when it should be in degree mode, or vice versa, results in a major error. 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 vale 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° =…

I recently had a project in which I had to reverse engineer a data analysis. There was some ambiguity regarding which of several possibilities someone chose for several of the variables, something analogous to the following example. 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…

A couple days ago, Levent Alpöge, a mathematician working at Anthropic, discovered a counterexample to the Jacobian conjecture using Claude Fable 5. 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…

The volume of a sphere of radius r is 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…

I’ve written several posts about using Claude or ChatGPT to generate Prolog or Lean code to solve a chess puzzle. I didn’t think Grok would be up to the task, though I didn’t try it. I’ve heard good things about Grok 4.5, so I gave it a shot. It did great. 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…

Suppose you want to search for a list of words. If you’re using grep, you can add the -f flag provide a file of regular expressions, and you can add the -F to tell it that the regular expressions are in fact just words. I did something like this a couple days ago when searching for diagnosis codes. 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…

Squares, high and low 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…

Quick follow up to the previous two posts on ICD-10 codes and HCPCS codes. This post uses Python’s squarify library to create treemaps visualizing how many codes begin with each letter. 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…

Since I revisited my old post on ICD code matching, I thought I’d revisit by post on HCPCS codes too. 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…

Seven years ago I wrote a post about regular expressions to match diagnosis codes. I wanted to revisit that post looking at speed and error rates. 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.…

I’ve been thinking about ICD-10 codes; they come up a lot in my work. 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…

A few days ago I wrote a post entitled Does additional data always reduce posterior variance?. In a nutshell, the answer is no, not always. 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),…

Common sense says that what you believe after seeing new data should be some sort of compromise between what you believed before and what the new data says. You don’t want to ignore previous information or new information. 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…

Years ago I wrote about Gilbreath’s conjecture. It’s a simple conjecture; you could explain it to anyone who understands what prime numbers are. See the linked post for a description of the problem. 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…

I ran across a geometry theorem with the following diagram. 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…

I ran across the approximation e ≈ 2721/1001 recently. What makes this remarkable is its accuracy relative to the size of the denominator. You can create a trivial approximation just by truncating a decimal expansion e ≈ 2718/1000 but this is only good to four significant figures, but 2721/1001 is good to seven, almost eight, significant figures. e = 2.71828182… 2721/1001 = 2.71828171… The comparison is more impressive in binary. $ bc -l >>> obase=2 >>> 2721/1001 10.10110111111000010100… >>>…

A discussion over lunch today brought up the fact that additional data does not always decrease the size of a confidence interval. This post will look at this from a Bayesian perspective. In general, new information reduces your uncertainty regarding whatever you’re estimating. The posterior distribution becomes more concentrated as more data are collected. That’s what happens “in general” but does it necessarily happen every time you get new data? Conceivably if you get surprising data, data…

This morning I wrote a post that included the central Delannoy numbers. The nth central Delannoy number Dn counts the number of ways a king can move from one corner of a chessboard to the diagonally opposite corner without backtracking. The more general Delannoy numbers Dm,n are the analogy for an m × n rectangular board, not necessarily square. Dm,n is also the number of possible sequence alignments for a strand of DNA with m base pairs and a strand with n base pairs [1]. At each step in the…

Imagine the following dialog. Professor: f is a function of a real variable x that takes a real parameter k. Student: What’s a parameter? Professor: It’s a constant that can vary. Student: Then if it can vary, isn’t it a variable? Professor: Sorta, but no not really. This conversation plays out over and over, and unfortunately it often ends as it does above, with the student confused. Here’s how I believe the conversation should continue. Professor: You’re absolutely right that f is a function…

Golden rectangles The defining property of golden rectangle is that if you stick a square on its longer side, you get another golden rectangle. The smaller vertical rectangle is similar to the larger horizontal rectangle. This means φ / 1 = (1 + φ) / φ which tells us φ² = 1 + φ and so the golden ratio φ equals (1 + √5)/2. Silver rectangles A silver rectangle is one that if you stick two squares on its longer side you get another rectangle with the same aspect ratio. This tells us σ / 1 = (1 +…

Here’s kind of a strange problem with an interesting solution: find a function f such that the derivative of f equals the inverse of f for all positive x. f ′(x) = f−1(x) This is a differential equation, but a very unusual one, one that cannot be solved using any of the techniques taught in a class on differential equations. The unique solution is f(x) = φ(x / φ)φ where φ is the golden ratio. What an unexpected appearance of the golden ratio! The problem was proposed by H. L. Nelson and solved…

The calculator utility bc has a minimal math library. For example, there’s no tangent function because you’re expected take the ratio of sine and cosine. (The Gnu version of bc goes have a function for tangent, but the POSIX version does not.) And yet bc includes support for Bessel functions J(x). The bc function j takes two arguments. Is the first argument n or x? Grok said the function arguments are j(n,x). I thought I should run man bc just to make sure, and it said j(x, n) Returns the…

Here’s a crazy bash one-liner I found via an article by Peter Krumins: echo {w,t,}h{e{n{,ce{,forth}},re{,in,fore,with{,al}}},ither,at} This prints 30 English words: when, whence, whenceforth, where, wherein, wherefore, wherewith, wherewithal, whither, what, then, thence, thenceforth, there, therein, therefore, therewith, therewithal, thither, that, hen, hence, henceforth, here, herein, herefore, herewith, herewithal, hither, hat This post will explain how the one-liner works. Bash brace…

The previous post looked at how many digits are in the reduced fraction for the nth harmonic number. I was curious about how long the cycle of digits in a harmonic number might be. I wrote about the period length for the digits of fractions almost a decade ago. This post includes code so I can apply it to harmonic demoninators. from sympy import lcm, factorint, n_order def period(n): factors = factorint(n) exp2 = factors.get(2, 0) exp5 = factors.get(5, 0) r = max(exp2, exp5) d = n // (2**exp2 *…