John D. Cook

The first post in the recent series of posts on Hadamard matrices describes a way of constructing new Hadamard matrices from two other Hadamard matrices by taking their Kronecker product. Starting with a Hadamard matrix H0 and a Hadamard matrix G, you can construct a sequence of Hadamard matrices by Hn+1 = G ⊗ Hn for positive integers n. This is known as the generalized Sylvester method. Let pn be the proportion of 1s in Hn and let q be the proportion of 1s in G. Then you can show that the…

Error correcting codes are most simply described in terms of the errors they can certainly correct. For example, the Hadamard code used for the Mariner 9 probe to Mars encoded each 6-bit pixel to a 32-bit codeword in such a way that the original pixel could be recovered if no more than 7 bits were corrupted in transit. What is the probability that a pixel could be repaired if corrupted? That depends on your probability model. We will assume that the probability of each bit being flipped is p…

Hadamard matrices are in the news following the recent announcement of a newly discovered Hadamard matrix. I’ve written three posts on Hadamard matrices recently, one as a sort of introduction and two on applications: the error correcting code used in the Mariner 9 probe and constructing sphere packings. A Hadamard matrix is an orthogonal matrix with all entries equal to ±1. Jacques Hadamard conjectured that there exist Hadamard matrices of order 4n for all positive integers n. It’s necessary…

Yesterday Levent Alpöge announced that he and his colleagues had discovered a new Hadamard matrix using Claude AI. That motivated a post I wrote this morning on how to construct Hadamard matrices. I mentioned in that post that these matrices arise in applications. This evening I gave an example, describing how NASA used a Hadamard matrix of order 32 to transmit photos from the Mariner 9 spacecraft in 1971. This post will give another application: sphere packing. Conway and Sloane [1] give a…

NASA set Mariner 9 to photograph Mars in 1971. The images had to be encoded for transmission using an error-correcting code, otherwise they would be significantly corrupted when they were received on Earth. The images were encoded for transmission using a code based on Hadamard matrices, specifically a (32, 6, 16) Hadamard code. This means that each 6-bit pixel value was encoded as a 32-bit code word, with all code words differing in at least 16 positions. The previous post explained a way to…

A Hadamard matrix is an orthogonal matrix whose entires are all either 1 or − 1. For example is a Hadamard matrix of order 2. True to Stigler’s law of eponymy, James Joseph Sylvester investigated Hadamard matrices before Jacques Hadamard. Sylvester saw how to bootstrap the example above into more examples. If H is a Hadamard matrix, then the partitioned matrix Sylvester’s construction can be generalized as follows. If Hm is a Hadamard matrix of order m and Hn is a Hadamard matrix of order n,…

Suppose you’ve never worked at the command line and you’re reading a book about the bash shell. You read that !$ is a shortcut to refer to the last word of the previous command. That little fact will almost certainly not stick in your head for a couple reasons. First, you probably see no need for such a shortcut. Second, the syntax seems completely arbitrary. But then after you use the command line a while, you might begin to notice a pattern. You often have to run a pair of commands with the…

I was reading a blog post on boat names because it was on Hacker News this morning. The post contained a link to a data set on dog names in NYC and I poked around the data a little. The top names were not at all what I expected, but then again this is limited to NYC; it’s not a sample across the US. These were the top 10 names: Bella Luna Max Charlie Coco Lola Rocky Milo Teddy Lucy I wondered if the name frequencies might fit a power-law distribution. They do not, but they follow a log-normal…

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…