John D. Cook

Suppose you were able to go back in time, to an era before computers, and give someone contemporary cryptography. Encryption methods that are essentially unbreakable now would certainly be unbreakable then. But there’s a catch: not only do attackers not have computers, neither do users. If you told someone about RSA encryption, for example, you’d lose them right after you said “First you find a couple 1000-digit primes.” But there’s no need for using 1000-digit primes if no attacker has a…

The following extraordinary paragraph comes from Knuth’s TAOCP Volume 4A, right before the last set of exercises. Many of the exercises below ask a modern reader to find and/or to correct errors in the literature of bygone days. The point is not to gloat over how smart we are in the 21st century; the point is rather to understand that even the pioneers of a subject can stumble. One good way to learn that a set of ideas is not really as simple as it might seem to today’s computer scientists and…

In science and engineering classes, you might describe a system using Newton’s laws and end up with a differential equation. You then solve the differential equation, analytically or numerically, to see how the solutions behave. You might also do the opposite, especially in a mathematics class: look at what differential equation a set of functions satisfy in order to understand those functions. Bessel functions came out of solving differential equations from astronomy. But then they turned out…

I was reading a book on the history of the discovery of the structure of DNA [1] and was surprised by a few passing references to Bessel functions. According to Claude, When X-rays are diffracted by a helical structure, the resulting diffraction pattern breaks into a series of horizontal “layer lines.” Cochran, Crick, and Vand showed mathematically that the diffracted amplitude on the n-th layer line is proportional to a Bessel function of the first kind, order n: Jn(2π r R) where r is the…

At the end of my post on how not to calculate cosine I said that the first step in calculating cosine, particularly cosine of a large number, would be to do range reduction. This post will present a simple range reduction method by Cody and Waite that is adequate for moderately large arguments. If you want to compute the sine or cosine of an angle x you could start by reducing x mod 2π since that would not change the result. However, accurately reducing a number mod 2π is not trivial; that’s…

I have a program that shares files between my laptop and my phone. It works well, except for apostrophes. When I type an apostrophe ' on my laptop, it becomes ’ on my phone. And when I type 's on my phone, it becomes 痴 on my laptop. Apparently the phone turns the apostrophe (U+0027) into a right single quote (U+2019), then bungles bytes in the UTF-8 encoding of U+2019 as three Windows-1252 characters. The bytes E28099hex are interpreted as â (E2hex), € (80hex), and ™ (99hex). When I type 's…

Calculus professors with no experience in numerical computing will tell students that computers calculate trig functions with power series. They don’t. I worked on the implementation of trig functions in hardware, and I can assure you we didn’t just use power series. Power series are an excellent way calculate functions near the center of the series, such as computing sine for small angles. But the further you get from the center, the less useful power series are. Let’s suppose you want to…

In a footnote to the previous post, I said that Python’s math library can calculate the logarithm of extremely large numbers but not the cosine. This post will expand on that comment. In this post I’ll use n = 200! as my example rather than 1000! nbecause this value of N is larger than the largest representable floating point number but small enough to be more convenient to work with. Suppose someone calculates 200! for you: 78865786736479050355236321393218506229513597768717326329474253324435\…

The previous post pointed out that the following code such as the following unexpectedly works. >>> from math import log, factorial >>> log(factorial(1000)) 5912.128178488163 If you don’t find this unexpected, note that if you replace math.log with numpy.log the code will fail [1]. Functions like natural logarithm operate on real numbers. Real numbers are represented as floating point numbers in programming languages, and 1000! factorial is too large to represent as a standard floating point…

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

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…

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…

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

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…

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. The sequence c(n) is also the number of constrains 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…

Second order linear differential equations with polynomial coefficients form their own area of study. This seems like a narrow class of equations, but it’s very important in applications. 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 and too…

Combinatorial problems can be interesting for their own sake, but they are more interesting when there is a connection to a problem outside combinatorics, and the more unexpected the connection the better. 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…

The textbook version of the Runge-Kutta method for solving differential equations has 4 stages and has 4th order error. For lower order versions of RK the number of stages s also matches the order of the error p. But in order to achieve error on the order of p ≥ 5, you need more than p stages. This is known as the Butcher barrier. 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…

I was digging into the Runge-Kutta method for solving differential equations and a line from [1] piqued my curiosity. 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 if it had been available in 1900. I go into Runge-Kutta methods in this post.…

A couple years ago I wrote about how to compute the inverse of factorial. I used that code in writing the previous post because the post required solving the equation ⌊log2(n!)⌋ = b given b. That is, given a number of bits b, find the smallest value of n such that n! > 2b. 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…

The previous post looked at the idea of storing a cryptographic key in the order of a deck of cards. A deck of 52 cards can store 225 bits of data because ⌊log2(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…

The latest issue of Paged Out! has an article by Stephen Hewitt “An off-line backup of your cryptographic key using playing cards.” The idea is to use a deck of 52 to store a 128-bit cryptographic key. To erase the key, shuffle the deck. In principle you could store a 225-bit key as a permutation of 52 cards because log2(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…

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…