John D. Cook

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 […]

Fermat’s little theorem says that if p is a prime number, then for any integer b, bp-1 = 1 (mod p). This gives a necessary but not sufficient test for a number to be prime. A number that satisfies the equation above but is not prime is called a pseudoprime [1] for base b. For example, […]

I ran across an article this morning about the booming market for low-tech tractors. Farmers are buying up old tractors at auction because the older tractors are cheaper to buy, cheaper to operate, and cheaper to repair. The article opens with the story of Kris Folland, a Minnesota farmer who bought a 1979 John Deere […]

A number is said to be “smooth” if all its prime factors are small. To make this precise, a number is said to be y-smooth if it only has prime factors less than or equal to y. So, for example, 1000 is 5-smooth. The de Bruijn function φ(x, y) counts the number of y-smooth positive […]

A sufficient statistic summarizes a set of data. If the data come from a known distribution with an unknown parameter, then the sufficient statistic carries as much information as the full data [0]. For example, given a set of n coin flips, the number of heads and the number of flips are sufficient statistics. More […]

The number of elements in a finite field must be a prime power, and for every prime power there is only one finite field up to isomorphism. The finite field with 256 elements, GF(28), is important in applications. From one perspective, there is only one such field. But there are a lot of different isomorphic […]

Someone left a comment this morning on my blog post on sinc and jinc integrals regarding the area of the lobes. It would be nice to have the values of integrals of each lobe, i.e. integrals between 0 and multiples of pi. Anyone knows of such a table? This post will include Python code to […]

It’s easy to manipulate CSV files with basic command line tools until you need to do a join. When your data is spread over two different files, like two tables in a normalized database, joining the files is more difficult unless the two files have the same keys in the same order. Fortunately, the xsv […]

This post shows how to export an Excel file to a CSV file using in2csv from the csvkit package. You could always use Excel itself to export an Excel file to CSV but there are several reasons you might not want to. First and foremost, you might not have Excel. Another reason is that you […]

Sometimes you’re in the flow using the command line and you’d like to briefly switch over to Python without too much interruption. Or it could be the other way around: you’re in the Python REPL and need to issue a quick shell command. One solution would be to run your shell and Python session in […]

Vitalii Tymchyshyn and Andrii Khlevniuk posted a new paper here entitled “On the average number of birthdays in birthday paradox setup.” This paper gives a different perspective on the birthday problem and a new proof. In the usual formation of the birthday problem, one asks the probability that everyone in a group of size n […]

I’ve playing around with driven Van der Pol oscillators. The cosine term is the driving function. You can see a variety of behavior depending on the interaction between the driving frequency and the natural frequency, along with the other parameters. Sometimes you get periodic solutions and sometimes you get chaotic behavior. Here’s the Python code […]

A few days ago I wrote about how to solve differential equations with SciPy’s ivp_solve function using Van der Pol’s equation as the example. Van der Pol’s equation is The parameter μ controls the amount of nonlinear damping. For any initial condition, the solution approach a periodic solution. The limiting periodic function does not depend […]