John D. Cook

The n queens problem is to place on an n × n chessboard n queens so that none attacks any other. This means there is only one queen on every horizontal, vertical, and diagonal line. When n is a prime number ≥ 5, it is sufficient to place the queens on a line that has slope 2, 3, 4, …, n − 2. (The slope cannot be 1 because that’s a diagonal. And it cannot be n − 1 because n − 1 = −1 mod n is also a diagonal.) [1] Here we imagine the top and bottom edge being identified. Geometrically, this makes…

I’ve written a couple posts lately on getting an LLM to generate code to solve chess problems. The first used Claude to generate Prolog and the second used ChatGPT to generate Prolog. This post will use Claude to generate Z3/Python code. The puzzle is one I’ve written about before: Place all the pieces—king, queen, two bishops, two knights, and two rooks—on a 6 × 5 chessboard, with the requirement that the two bishops be on opposite colored squares and no piece is attacking another.…

I’ve been testing Claude’s ability to generate Lean 4 code to prove theorems. I’ve written about a couple experiments that verified calculations. I did not write about my failed attempt to get Claude to formalize a proof of the pqr theorem for seminorms. This time I asked Claude to formally prove the theorem from the previous post on partial fraction decomposition in an abstract setting. Unlike the previous experiments, I did not give it a proof or calculation to verify; I asked it to generate…

Nearly everyone who as seen partial fraction decomposition was introduced to it as a way to compute integrals. If P(x) and Q(x) are polynomials, then you can break their ratio P(x)/Q(x) into a sum of terms that can each be integrated in closed form. As with most topics in a calculus class, partial fractions go by in a blur. This post will look at partial fractions more generally. Computation Every polynomial with real coefficients can be factored into a product of linear and irreducible…

You can’t prove a theorem by just checking a few examples. Except sometimes you can. A few weeks ago I wrote Pentagonal numbers are truncated triangular numbers. In a nutshell, if the pentagonal numbers are defined by Pn = (3n² − n)/2 and the triangular numbers by Tn = (n² + n)/2 then Pn = T2n − 1 − Tn − 1. Here’s a visualization of the equation. Note that the equation asserts that two quadratic polynomials are equal. If the two polynomials are equal at three points, then they’re equal…

The nth pentagonal number Pn is the number of dots in diagrams like those below with n concentric pentagons. We have the formula Pn = (3n² − n)/2 where n is a positive integer. If n is an integer but not positive, the equation above defines a generalized pentagonal number. If you’re given an n, you can easily compute Pn. But suppose you’re given a large number x. How would you determine if it is a pentagonal number? And if it is a pentagonal number, how would you find n such that x = Pn? If x =…

I got an email message this afternoon reporting a typo in a blog post from about a year ago on converting between quaternions and rotation matrices [1]. The email said exactly where the typo was, but I decided to see whether Claude would find it. Specifically, I prompted Sonnet 4.6 Medium with the following. Write Lean code to verify the two theorems at the top of this post: https://www.johndcook.com/blog/2025/05/07/quaternions-and-rotation-matrices/ That is, prove that the expressions given in…

A few days ago I wrote about using Claude to solve a chess puzzle by writing Prolog code. This morning I tried a similar chess puzzle with ChatGPT. The task is to place a queen, king, rook, bishop, and knight on a 4 by 4 chessboard so no piece attacks another. Of course there’s not a lot of market demand for solving chess puzzles, but there is demand for solving real-world problems that amount to logic puzzles. Experiments with chess puzzles suggest writing Prolog with LLM help might be…

When the US government classified strong encryption as “munitions,” RSA public key cryptography was technically illegal. In 1995, Adam Back protested this by creating a terse, obfuscated implementation of RSA in two lines of Perl code and used it as an email signature. The code was also printed on T-shirts. This was the code: #!/bin/perl -s-- -export-a-crypto-system-sig -RSA-3-lines-PERL $m=unpack(H.$w,$m."\0"x$w),$_=`echo "16do$w 2+4Oi0$d*-^1[d2%Sa…

Prolog is the original logic programming language. The name comes from programming in logic. More specifically, the name comes from programmation en logique because the inventor of the language, Philippe Roussel, is French. Prolog has its advantages and disadvantages. One of the advantages is that the language represents logical problems directly. One of the disadvantages is that the syntax can be quirky. But if an LLM is writing the code, or at least helping to write the code, the syntax…

I ran an experiment today to see whether Claude could generate Lean code to prove a calculation at the bottom of this post, six lines of calculus. I started with this prompt This page contains a mathematical proof that a Fourier coefficient, a_n, is given in terms of a Bessel function. The LaTeX source for the SVG image is contained in the alt tag of the image. Generate a formal proof of the result using Lean. and give it the URL of the post. It produced a proof without any further input, but…

Often there’s a thread running through a sequence of my posts. Sometimes I make this explicit and sometimes I don’t. The latest thread started with this post commenting on a tweet that observed that exp(−x²) ≈ (1 + cos(sin(x + x))/2. Some people said online that that the approximation is simply due to the first few terms of the Taylor series on both sides matching up, so I wrote a follow up post explaining that it’s not that simple. The series for the left hand side alternates and converges…

Kepler solved his eponymous equation M = E − e sin(E) by finding a fixed point of E = M + e sin(E). So guess a value of E and stick it into the right hand side. Then plug that value into the right hand side again. Kepler said a couple iterations should be enough. And a couple iterations are enough if the eccentricity e is small and you don’t need much accuracy. The rate of convergence is determined by e. Kepler implicitly had in mind small values of e because he wasn’t aware of anything…

The previous post discussed how to solve Kepler’s equation M = E − e sin(E) using a sine series. You could also solve Kepler’s equation using a power series, which Lagrange did in 1771. Both approaches express E as a function of e and M, but from different perspectives. Bessel though of his solution as a sum of sines in M, with coefficients that depend on e. Lagrange thought of his solution as a power series in e whose coefficients involve sines in M. You can rearrange the terms of either…

I ran across a cranky formula for π based on physical constants here and decided to play around with it. The source describes λ as “wavelength (chosen in the microwave region)” and I thought perhaps you could chose a value of λ to make the equation work. But as a comment pointed out, the bracketed expression is simply 2hc², independent of λ, due to Planck’s blackbody law. That means we can simplify the expression above to Now the values of h and c are known. In fact, they’re now exactly known…

The previous post very briefly said that the integral representation for Bessel functions was motived by solving Kepler’s equation. This post will go into more detail. Kepler’s equation There are multiple ways to describe the position of a planet in an elliptical orbit around a star. For historical reasons, these descriptions have arcane names such as mean anomaly, true anomaly, and eccentric anomaly. This post explains how these three are related. For this post, it is enough to say that often…

Yesterday I wrote a post showing that the trapezoid rule evaluates the integral very efficiently. But how do we know what the exact integral is for comparison? If you ask Mathematica, it will tell you the integral equals −2π J1(1) where J1 is a Bessel function. This may seem like rabbit out of a hat, but it’s actually a simple calculation given the integral definition of Bessel functions: Since cosine is even, we can write our integral over [−π, π] as twice the integral over [0, π]. Then a…

One reason people study Latin is that it is the ancestor of many modern languages. English derives from West Germanic languages, not from Latin, but much of English vocabulary, perhaps as much as 60%, derives from Latin, either directly or indirectly through French. Knowing a bit of Latin makes sense of many things that would otherwise seem completely arbitrary, such as why the symbols for gold, silver, and lead are Au, Ag, and Pb respectively. Similarly, ed(1) is the Latin of Linux [1]. Many…

Several posts lately have looked at the function f(x) = cos(sin(x) + x). This post will look at the function from a different angle. It’s a smooth function with period 2π, and it’s very flat at odd multiples of π, i.e. the first five derivatives are zero. For reasons I wrote about here, this means that the trapezoid rule should be very efficient for integrating this function. In general, the error in the trapezoid rule is on the order of 1/N² where N is the number of integration points. To be…

I was thinking more about the cosine approximation to the Gaussian exp(−z²) ≈ (1 + cos(sin(z) + z))/2 that I wrote about last week. The two expressions above are close along the real axis but not along the imaginary axis. If z = iy, the right side grows much faster than the left, behaving like exp(exp(y)). This led to me looking up the power series for the double exponential function exp(exp(y)). This is an interesting series because the coefficient of xn is e Bn / n! where Bn is the nth Bell…

Suppose you run across the power series for the exponential function and decide to code it up. Good idea: you’ll probably learn something, though maybe not what you expect. Maybe you decide a tolerance of 10−12 is good enough, and so you sum the terms until the next term to add is below the tolerance. from math import factorial, exp def naive_exp(x): tolerance = 1e-12 s = 0 n = 0 while True: delta = x**n / factorial(n) s += delta if abs(delta) < tolerance: return s n += 1 You want to try your…

There is still active discussion on X about the approximation exp(−x²) ≈ (1 + cos(sin(x) + x))/2 and some are saying this can just be explained by Taylor series: the series for the two sides differ for the first time at the x6 term, so that’s why you get a good approximation. As I wrote yesterday, that’s only part of it. If it were just about Taylor series you could use exp(−x²) ≈ 1 − x² + x4/2 which also has error O(x6). But this approximation is only good for fairly small x, say in [−0.5,…

Readers have subscribed to this blog via email almost from its beginning in 2008, but how they have subscribed has changed several times. I’ve used several services to provide email subscription that have come and gone. For the past two years I’ve been using Substack to send out emails announcing new blog posts. That has worked out well. Substack delivers email reliably, and that’s all I wanted. I’m not active on Substack other than using it to an email. I give a brief introduction to the…

The function (1 + cos(x))/2 gives a fair approximation to the Gaussian density exp(−x²) You can make the approximation much better by raising it to a power. The function ((1 + cos(x))/2)4 gives a good lower bound and ((1 + cos(x))/2)3.5597 gives a good upper bound. More on that here. There are other ways of improving the cosine approximation to the Gaussian. Yesterday I came across one I hadn’t seen before, adding a sin(x) term to x. (1 + cos(sin(x) + x))/2 This function matches the first few…

If a polynomial identity holds at a few random points, it’s very like true. We’ll make this statement more precise, but first let’s look at some applications. You may want to test an identity that naturally presents itself as a statement that two polynomials are equal. Or you might use something like the binomial coefficient trick to reframe a problem isn’t obviously an identity about polynomials. And with algebraic circuits, you can reformulate a wide range of computations as polynomial…

Canonical example The variance of a set of numbers is defined as the sum of the squared distances from each point to the mean. So it would seem that you first need to calculate the mean, then go back and compute the squared differences from the mean. And yet the variance can be computed in one pass through the data. You’ll find two equivalent equations in statistics books: the one described above and another based on the sum of the data points and the sum of the data points squared. While this…

Kullback-Leibler divergence Kullback-Leibler divergence is defined for two random variables X and Y by K-L divergence is non-negative, and it’s zero if and only if X and Y have the same distribution. But it is not a metric, for reasons explained here. For one thing, it’s not symmetric. Jeffreys divergence We can fix the symmetry problem by defining The J above stands for Jeffreys, for Harold Jeffreys. J is called either the symmetrized K-L divergence or Jeffreys’ divergence. It’s still a…

I saw a post yesterday saying that the Meta logo is a Besace curve. A Besace curve has the implicit form and the parametric form where t ranges over [0, 2π]. So given a Besace curve, such as the Meta logo, how do you find the parameters a and b to fit the curve? We can rewrite the parametric expression for x as a sine with a phase shift (see notes here) where Also, we can rewrite the parametric expression for y as Now the extreme values of x and y are easier to see. The maximum value of x is A…

The previous post looked at the expected IQ range in a jury of 12. This post will look more generally at computing the expected range of n samples from a N(0, 1) random variable. This will give the expected range in units of σ, i.e. multiply the results by σ if your σ isn’t 1. As mentioned in the previous post, the expected range is given by where φ and Φ are the PDF and CDF of a standard normal. The integral can be calculated in closed form for n ≤ 5, but in general it requires numerical…

There’s been some discussion online lately about how a large difference in IQ makes it difficult for two people to communicate. There have been studies that confirm this effect. The difficulty is not insurmountable, but it takes deliberate effort to overcome. Someone dismissed this communication difficulty by pointing out that the expected difference in IQ between two individuals is around 17, suggesting that most communication is between people who differ by more than one standard deviation in…