John D. Cook

The comments to the posts here are generally very helpful. I appreciate your contributions to the site. I wanted to offer a tip for those who leave comments and are frustrated by the way the comments appear, especially those who write nicely formatted snippet of code only to see the formatting lost. There is a […]

Take a positive integer x, square each of its digits, and sum. Now do the same to the result, over and over. What happens? To find out, let’s write a little Python code that sums the squares of the digits. def G(x): return sum(int(d)**2 for d in str(x)) This function turns a number into a […]

Heron’s formula computes the area of a triangle given the length of each side. where If you have a very thin triangle, one where two of the sides approximately equal s and the third side is much shorter, a direct implementation Heron’s formula may not be accurate. The cardinal rule of numerical programming is to […]

I recently stumbled on a paper [1] that looks at a cubic equation that comes out of a problem in orbital mechanics: σx³ = (1 + x)² Much of the paper is about the derivation of the equation, but here I’d like to focus on a small part of the paper where the author looks […]

A couple days ago I wrote about Hamming codes and said that they are so-called perfect codes, i.e. codes for which Hamming’s upper bound on the number of code words with given separation is exact. Not only are Hamming codes perfect codes, they’re practically the only non-trivial perfect codes. Specifically, Tietavainen and van Lint proved […]

The previous post mentioned adding a parity bit to a string of bits as a way of detecting errors. The parity of a binary word is 1 if the word contains an odd number of 1s and 0 if it contains an even number of ones. Codes like the Hamming codes in the previous post […]

The previous post looked at how to choose five or six letters so that their Morse code representations are as distinct as possible. This post will build on the previous one to introduce Hamming codes. The problem of finding Hamming codes is much simpler in some ways, but also more general. Morse code is complicated […]

A century ago the German army used a field cipher that transmitted messages using only six letters: A, D, F, G, V, and X. These letters were chosen because their Morse code representations were distinct, thus reducing transmission error. The ADFGVX cipher was an extension of an earlier ADFGV cipher. The ADFGV cipher was based […]

ChaCha is a CSPRING, a cryptographically secure pseudorandom number generator. When used in cryptography, ChaCha typically carries out 20 rounds of its internal scrambling process. Google’s Adiantum encryption system uses ChaCha with 12 rounds. The runtime for ChaCha is proportional to the number of rounds, so you don’t want to do more rounds than necessary […]

Sometimes you need to count the number of 1’s in a stream of bits. The most direct application would be summarizing yes/no data packed into bits. It’s also useful in writing efficient, low-level bit twiddling code. But there are less direct applications as well. For example, three weeks ago this came up in a post […]

You might hear hysteresis described as a phenomena where the solution to a differential equation depends on its history. But that doesn’t make sense: the solution to a differential equation always depends on its history. The solution at any point in time depends (only) on its immediately preceding state. You can take the state at […]

There are two ways to deidentify data to satisfy HIPAA: Safe Harbor, § 164.514(b)(2), and Expert Determination, § 164.514(b)(1). And for reasons explained here, you may need to be concerned with HIPAA even if you’re not a “covered entity” under the statute. To comply with Safe Harbor, your data may not contain any of eighteen […]

Linear congruence random number generators have the form xn+1 = a xn + b mod p Inverse congruence generators have the form xn+1 = a xn-1 + b mod p were x-1 means the modular inverse of x, i.e. the value y such that xy = 1 mod p. It is possible that x = […]

This is the third post in a series on Runge-Kutta methods. The first post in the series introduces Runge-Kutta methods and Butcher tableau. The next post looked at Fehlberg’s adaptive Runge-Kutta method, first published in 1969. This post looks at a similar method from Dormand and Prince in 1980. Like Fehlberg’s method, the method of […]

This post brings together several themes I’ve been writing about lately: caching function evaluations, error estimation, and Runge-Kutta methods. A few days ago I wrote about how Runge-Kutta methods can all be summarized by a set of numbers called the Butcher tableau. These methods solve by evaluating f at some partial step, then evaluating f […]

This post will look at two numerical integration methods, the trapezoid rule and Romberg’s algorithm, and memoization. This post is a continuation of ideas from the recent posts on Lobatto integration and memoization. Although the trapezoid rule is not typically very accurate, it can be in special instances, and Romberg combined it with extrapolation to […]

I was browsing through SciPy documentation this evening and ran across a function in scipy.misc called electrocardiogram. What?! It’s an actual electrocardiogram, sampled at 360 Hz. Presumably it’s included as convenient example data. Here’s a plot of the first five seconds. I wrote a little code using it to turn the ECG into an audio […]

The HIPAA privacy rule only applies to “covered entities.” This generally means insurance plans, healthcare clearinghouses, and medical providers. If your company is using heath information but isn’t a covered entity per the HIPAA statute, there are a couple reasons you might still need to pay attention to HIPAA [1]. The first is that […]

The previous post explained that Lobatto’s integration method is more efficient than Gaussian quadrature when the end points of the interval need to be included as integration points. It mentioned that this is an advantage when you need to integrate over a sequence of contiguous intervals, say [1, 2] then [2, 3], because the function […]

A basic idea in numerical integration is that if a method integrates polynomials exactly, it should do well on polynomial-like functions [1]. The higher the degree of polynomial it integrates exactly, the more accurate we expect it will be on functions that behave like polynomials. The best known example of this is Gaussian quadrature. However, […]

If you know one numerical method for solving ordinary differential equations, it’s probably Euler’s method. If you know two methods, the second is probably 4th order Runge-Kutta. It’s standard in classes on differential equations or numerical analysis to present Euler’s method as conceptually simple but inefficient introduction, then to present Runge-Kutta as a complicated but […]

I wrote a blog post of sorts, spread over several tweets, about plotting functions of a complex variable. Plot of cosine over a square centered at the origin of the complex plane. Color = phase. Height = magnitude. pic.twitter.com/rOT7tAkfM9 — Analysis Fact (@AnalysisFact) February 11, 2020 Here are a couple of the images from the […]

Richard Arenstorf (1929–2014) discovered a stable periodic orbit between the Earth and the Moon which was used as the basis for the Apollo missions. His orbit is a special case of the three body problem where two bodies are orbiting in a plane, i.e. the Earth and the Moon, along with a third body of […]

Having watched a few episodes of Phineas and Ferb, when I see “Brusselator” I imagine Dr. Doofenschmertz saying “Behold! The Brusselator!” But the Brusselator is considerably older than Phineas and Ferb. It goes back to Belgian scientists René Lefever and Grégoire Nicolis in 1971 [1] who combined “Brussels” and “oscillator” to name the system after […]

In recent posts I’ve written about using RNG test suites on the output of the μRNG entropy extractor. This is probably the last post in the series. I’ve looked at NIST STS, PractRand, and DIEHARDER before. In this post I’ll be looking at TestU01. TestU01 includes three batteries of tests: Small Crush, Crush, and Big […]

The first well-known suite of tests for random number generators was George Marsalia’s DIEHARD battery. The name was a pun on DieHard car batteries. Robert G. Brown took over the DIEHARD test suite and called it DIEHARDER, presumably a pun on the Bruce Willis movie. I’ve written lately about an entropy extractor that creates a […]

Yesterday I wrote about my experience using NIST STS to test an entropy extractor, a filtering procedure that produces unbiased bits from biased sources. This post will look at testing the same entropy extractor using the Practically Random (PractRand) test suite. The results were much worse this time, which speaks to the limitations of both […]

We all learn as children that there are 60 seconds in a minute, 60 minutes in an hour, 24 hours in a day, and 365 days in a year. Then things get more complicated. There are more like 365.25 days in a year, hence leap years. Except that’s quite right either. It’s more like 365.242 […]

Around this time last year I wrote about the entropy extractor used in μRNG. It takes three biased random bit streams and returns an unbiased bit stream, provided each input stream as has least 1/3 of a bit of min-entropy. I’ve had in the back of my mind that I should go back and run […]

Here’s an interesting theorem I ran across recently. The number of odd integers in the nth row of Pascal’s triangle equals 2b where b is the number of 1’s in the binary representation of n. Here are the first few rows of Pascal’s triangle: 1 1 1 1 2 1 1 3 3 1 1 […]

There is no precise definition of what it means for a differential equation to be stiff, but essentially it means that implicit methods will work much better than explicit methods. The first use of the term [1] defined stiff equations as equations where certain implicit methods, in particular BDF, perform better, usually tremendously better, than […]

I keep running into Mittag-Leffler. A couple days ago I wrote about his polynomials. Today I ran across his regularization method for summing a divergent series. Before that I wrote about his generalization of the exponential function, which is closely related to his summation method. The exponential function has power series where we’ve written the […]

I’ve made a couple minor changes to my page that converts between frequency and pitch. (The page also includes Barks, a psychoacoustic unit of measure.) If you convert a frequency in Hertz to musical notation, the page used to simply round to the nearest note in the chromatic scale. Now the page will also tell […]

A while back I wrote about the Mittag-Leffler function which is a sort of generalization of the exponential function. There are also Mittag-Leffler polynomials that are a sort of generalization of the powers of x; more on that shortly. Recursive definition The Mittag-Leffler polynomials can be defined recursively by M0(x) = 1 and for n […]

The Lorenz system is a canonical example of chaos. Small changes in initial conditions eventually lead to huge changes in the solutions. And yet discussions of the Lorenz system don’t simply show this. Instead, they show trajectories of the system, which make beautiful images, but do not demonstrate the effect of small changes to initial […]

A few days ago I wrote about bifurcation for a discrete system. This post looks at a continuous example taken from the book Exploring ODEs. We’ll consider solutions to the differential equation for two different initial conditions: y(0) = 0.02, y‘(0) = 0 and y(0) = 0.05, y‘(0) = 0. We’ll solve the ODE with […]

Goodhart’s law says “When a measure becomes a target, it ceases to be a good measure.” That is, when people are rewarded on the basis of some metric, they’ll learn how to improve that metric, but not necessarily in a way that increases what you’re after. Here are three examples of Goodhart’s law related to […]

There’s a worn-out analogy in software development that you cannot build a skyscraper the same way you build a dog house. The idea is that techniques that will work on a small scale will not work on a larger scale. You need more formality to build large software systems. The analogy is always applied in […]

Cobweb plots are a way of visualizing iterations of a function. For a function f and a starting point x, you plot (x, f(x)) as usual. Then since f(x) will be the next value of x, you convert it to an x by drawing a horizontal line from (x, f(x)) to (f(x), f(x)). In other […]

California’s new CCPA (California Consumer Privacy Act) may become more like HIPAA. In particular, a proposed amendment would apply HIPAA’s standards of expert determination to CCPA. According to this article, The California State Senate’s Health Committee recently approved California AB 713, which would amend the California Consumer Privacy Act (CCPA) to except from CCPA requirements […]

If you enter any number into a calculator and repeatedly press the cosine key, you’ll eventually get 0.73908513, assuming your calculator is in radian mode [1]. And once you get this value, taking the cosine more times won’t change the number. This is the first example of a fixed point of an iteration that many […]

The previous post looked at computing recurrence relations. That post ends with a warning that recursive evaluations may nor may not be numerically stable. This post will give examples that illustrate stability and instability. There are two kinds of Bessel functions, denoted J and Y. These are called Bessel functions of the first and second kinds […]

A video by Raymond Hettinger points out that simultaneous assignment makes it much easier to understand code that evaluates a recurrence relation. His examples were in Python, but the same principle applies to any language supporting simultaneous evaluation. The simplest example of simultaneous evaluation is swapping two variables: a, b = b, a Compare this […]

The other day I ran across the biochemical pathways poster from Roche. This is the first of two posters. Both posters are about 26 square feet in area. Below is a close up of about one square foot in the middle of the first poster. I’d seen this before, I think while I was working […]

This post is a follow-on to the post on how to make the logistic bifurcation diagram below. That post plotting the attractors for iterations of f(x) = r x(1 – x). This post will plot a few trajectories over time for differing starting points and varying values of r. For values of r less than […]

If you look at the Wikipedia article on Weierstrass functions, you’ll find a line that says “the relation between the sigma, zeta, and ℘ functions is analogous to that between the sine, cotangent, and squared cosecant functions.” This post unpacks that sentence. Weierstrass p function First of all, what is ℘? It’s the Weierstrass elliptic […]

List comprehensions in Python let you create a list declaratively, much like the way you would describe the set in English. For example, [x**2 for x in range(10)] creates a list of the squares of the numbers 0 through 9. If you wanted the sum of these numbers, you could of course say sum( [x**2 […]

I looked back to see what posts I had written on this date in years past. January 14 doesn’t have any particular significance for me, so these post are a representative cross section of the blog. Some years I didn’t write anything on January 14, and some years I wrote posts that are not worth […]

The image below is famous. I’ve seen it many times, but rarely with a detailed explanation. If you’ve ever seen this image but weren’t sure exactly what it means, this post is for you. This complex image comes from iterating a simple function f(x) = r x(1 – x) known as the logistic map. The […]

Every condition that all primes satisfy has a corresponding primality test and corresponding class of pseudoprimes. The most famous example is Fermat’s little theorem. That theorem says that for all primes p, a certain congruence holds. The corresponding notion of pseudoprime is a composite number that also satisfies Fermat’s congruence. To make a useful primality […]