John D. Cook

The previous post looked at writing the harmonic numbers as reduced fractions and estimating the number of digits in the numerator and denominator based on asymptotics. This is a follow up post with plots. We’ll choose our base b to be 2. And we’ll look at the total number of bits in both the numerator and denominator, which we will use as the height of the fractions. First, let’s look at the actual and estimated heights, using the estimates from the previous post. Next let’s look at the…

The nth harmonic number is the sum of the reciprocals of the first n positive integers. Hn = 1 + 1/2 + 1/3 + 1/4 + … + 1/n The product of all the denominators is n!, so you could write Hn as a fraction Hn = p/q where p = n! Hn is an integer and q = n!. While p/q is a way to write Hn as a fraction, it’s not the most efficient because p and n! will have common factors. If we write Hn as a reduced fraction, the denominator will be the least common multiple of the integers 1 through n. That number…

Hart’s theorem says If a triangle be formed by the arcs of three circles, the inscribed and the three escribed circles are all tangent to a new circle or line. Here “triangle” means a three-sided figure whose sides are portions of a circle. The inscribed circle is the largest circle that can fit inside the three-sided figure. The “escribed” circles are analogous to the excircles in the previous post: you extend two sides and find a circle that is tangent to the triangle side and the extended…

This post will reveal the connection between my two previous posts: one on the Star Trek lemma and one on Pythagorean triples. In the process of writing the latter, I looked at the Wikipedia article on Pythagorean triples and noticed this curious paragraph. In every Pythagorean triangle, the radius of the incircle and the radii of the three excircles are positive integers. Specifically, for a primitive triple the radius of the incircle is r = n(m − n), and the radii of the excircles opposite…

In this post we find all Pythagorean triples that contain consecutive numbers, all Pythagorean triples (a, b, c) such that a + 1 = b or b + 1 = c. a + 1 = b George Osborne wrote a paper [1] addressing the question of when the squares of two consecutive numbers is also a square. Geometrically this is asking for primitive Pythagorean triples for which the legs are consecutive integers. He proved that the sequence shorter legs satisfies the recurrence relation with initial conditions u0 = 0 and u1…

I was reading an article this evening and saw a footnote to a book by Arthur Baragar [1]. This caught my eye because he was my officemate at UT for a year. I found his book on Archive.org and was surprised to see “The Star Trek Lemma” in the table of contents. What could this be? It’s a theorem that goes back to Euclid that applies to an angle formed by connecting a point to two other points on a circle. The theorem says “The measure of an inscribed angle is half the measure of the arc it…

The most frustrating aspect of regular expressions is that implementations vary. Features supported in one tool may not be supported at all in another tool, or they may be supported with slightly different syntax. I learned regular expressions in the context Perl, a maximalist regex environment. This led to frustration when features I expect to work are missing. One way around this is to use Perl analogs of other tools, but this is very non-standard. I want to be able to send colleagues and…

Let f be an even function with period π. Then the following remarkable theorem by Lobachevsky holds. This theorem is useful in Fourier analysis and signal processing. It’s useful to know even in the special case f(x) = 1. Related posts Sinc and jinc integrals Sinc and jinc sums The post Lobachevsky’s integral formula first appeared on 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,…