John D. Cook

Thenth harmonic number is the sum of the reciprocals of the firstn positive integers. Hn = 1 + 1/2 + 1/3 + 1/4 + ... + 1/n The product of all the denominators isn!, so you could write Hn as a fraction Hn = p/q wherep =n! Hn is an integer and q =n!. While [...]The post Writing down harmonic numbers first appeared on John D. Cook.

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 [...]The post Hart's theorem first appeared on John D. Cook.

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 theincircleand the radii of the [...]The post Incircles and Excircles of Pythagorean triangles first appeared on John D. Cook.

In this post we find all Pythagorean triples that contain consecutive numbers, all Pythagorean triples (a,b,c) such thata + 1 =b orb + 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 [...]The post Consecutive Pythagorean triangle sides first appeared on John D. Cook.

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 [...]The post The Star Trek lemma first appeared on John D. Cook.

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 [...]The post Regular expressions that work everywhere" first appeared on John D. Cook.

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. For a jinc" analog, see this paper. *** Every time I see the name Lobachevsky I [...]The post Lobachevsky's integral formula first appeared on John D. Cook.

Then queens problem is to place on ann *n chessboardn queens so that none attacks any other. This means there is only one queen on every horizontal, vertical, and diagonal line. Whenn is a prime number >= 5, it is sufficient to place the queens on a line that has slope 2, 3, 4, ..., [...]The post Queens on a prime order board first appeared on John D. Cook.

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 [...]The post All pieces on a 6 by 5 board first appeared on John D. Cook.

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 [...]The post Formalizing a ring theorem with Lean 4 and Claude first appeared on John D. Cook.

Nearly everyone who has seen partial fraction decomposition was introduced to it as a way to compute integrals. IfP(x) andQ(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 [...]The post Partial fraction decomposition first appeared on John D. Cook.

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^2 -n)/2 and the triangular numbers by Tn= (n^2 +n)/2 then Pn=T2n- 1-Tn- 1. Here's a visualization [...]The post Three examples suffice first appeared on John D. Cook.

Thenth pentagonal number Pn is the number of dots in diagrams like those below with n concentric pentagons. We have the formula Pn = (3n^2 - n)/2 wheren is a positive integer. Ifn is an integer but not positive, the equation above defines ageneralized pentagonal number. If you're given ann, you can easily compute Pn. [...]The post Testing pentagonal numbers first appeared on John D. Cook.

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 [...]The post Quaternion Rotations, Claude, and Lean first appeared on John D. Cook.

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 [...]The post Writing Prolog with ChatGPT first appeared on John D. Cook.

Back when the US government classified strong encryption as munitions," RSA public key cryptography was illegal to export. In 1995, Adam Back protested this by creating a terse, obfuscated implementation of RSA in Perl code and used it as an email signature. The code was also printed on T-shirts. The shirt was classified as munitions [...]The post RSA munitions T-shirt first appeared on John D. Cook.

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 [...]The post Solving a chess puzzle with Claude and Prolog first appeared on John D. Cook.

I ran an experiment today to see whether Claude [1] 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 [...]The post Formally proving a calculation with Claude and Lean first appeared on John D. Cook.

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^2) (1 + cos(sin(x) + x))/2. Some people said online that that the approximation is simply due to the first [...]The post Pulling on a thread first appeared on John D. Cook.

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 ofE 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 [...]The post Aitken acceleration before Aitken first appeared on John D. Cook.

An earlier 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 expressE as a function ofe andM, but from different perspectives. Bessel thought of his solution as a sum of sines in [...]The post The Laplace limit first appeared on John D. Cook.

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 [...]The post A crank formula for first appeared on John D. Cook.

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 [...]The post From Kepler to Bessel first appeared on John D. Cook.

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, [...]The post Mr. Bessel's eponymous functions first appeared on John D. Cook.

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 [...]The post The Latin of Linux first appeared on John D. Cook.

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. For reasons I wrote about here, this means that the trapezoid rule should find its integral very efficiently. In general, the error in the trapezoid [...]The post Integrating smooth periodic functions first appeared on John D. Cook.

I was thinking more about the cosine approximation to the Gaussian exp(-z^2) (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. Ifz =iy, the right side grows much faster than the left, behaving likeexp(exp(y)). This led to [...]The post Partitions over permutations first appeared on John D. Cook.

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. [...]The post Naively summing an alternating series first appeared on John D. Cook.

There is still active discussion on X about the approximation exp(-x^2) (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 [...]The post It's not just Taylor series first appeared on John D. Cook.

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 [...]The post Subscribe by email first appeared on John D. Cook.

The function (1 + cos(x))/2 gives a fair approximation to the Gaussian density exp(-x^2) 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 [...]The post Another Gaussian approximation first appeared on John D. Cook.

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 [...]The post Spot checking polynomial identities first appeared on John D. Cook.

Canonical example The sample variance of a set of numbers is defined in terms of 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 sample variance can [...]The post Online (one-pass) algorithms first appeared on John D. Cook.

Kullback-Leibler divergence Kullback-Leibler divergence is defined for two random variablesX and Y by K-L divergence is non-negative, and it's zero if and only if X andY 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 [...]The post Turning K-L divergence into a metric first appeared on John D. Cook.

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 parametersa andb to fit the curve? We can rewrite [...]The post The Meta logo and fitting Besace curves first appeared on John D. Cook.

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 ofn samples from aN(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 [...]The post Calculating the expected range of normal samples first appeared on John D. Cook.

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 [...]The post Expected IQ spread on a jury first appeared on John D. Cook.

The previous post linked to a post I wrote a few years ago about the Hilbert transform and Fourier series. That post says that if the Fourier series of a function is then the Fourier series of its Hilbert transform is When I looked back at that post I thought about how if you thought [...]The post Hilbert transform as an infinite matrix first appeared on John D. Cook.

The previous post announced some notes I wrote up based on an article by Henry Baker implementing functions of a complex variable in terms of functions of a real variable. That is, it finds functions u(x, y) and v(x, y) such that f(x + iy) = u(x, y) + i v(x, y) wherex,y,u, and vare [...]The post Real and imaginary parts first appeared on John D. Cook.

A couple months ago I wrote about how to compute the sine and cosine of a complex number using only real functions of real variables using the equations You can do something analogous for all the elementary functions, though some of the equations are quite a bit more complicated than the ones above. See the [...]The post Building complex functions out of real parts first appeared on John D. Cook.

You can't always get what you want. But sometimes you get what you need." - The Rolling Stones Circular functions and hyperbolic functions aren't invertible, but we invert them anyway. These functions map many points in the domain to each point in the range, and we invert them by mapping a point in the range [...]The post Couth and uncouth function pairs first appeared on John D. Cook.

The difference between a circular function and a hyperbolic function is a rotation or two. For example, cosh(z) = cos(iz). You can read that as saying that to find the hyperbolic cosine of z, first you rotatez a quarter turn to the left (i.e. multiply by i) and then take the cosine. For another example, [...]The post Circular and hyperbolic functions differ by rotations first appeared on John D. Cook.

How should we define (z^2 - 1)? Well, you could square z, subtract 1, and take the square root. What else would you do?! The question turns out to be more subtle than it looks. Whenx is a non-negative real number, x is defined to be the non-negative real number whose square isx. Whenx is [...]The post Square root of x^2 - 1 first appeared on John D. Cook.

The previous post derived the identity and said in a footnote that the identity holds at least forx > 1 andy > 1. That's true, but let's see why the footnote is necessary. Let's have Mathematica plot The plot will be 0 where the identity above holds. The plot is indeed flat forx > 1 [...]The post Closer look at an identity first appeared on John D. Cook.

Markov numbers are integer solutions to x^2 +y^2 +z^2 = 3xyz. The Wikipedia article on Markov numbers mentions that Don Zagier studied Markov numbers by looking the approximating equation x^2 +y^2 +z^2 = 3xyz+ 4/9 which is equivalent to f(x) +f(y) = f(z) wheref(t) is defined as arccosh(3t/2). It wasn't clear to me why the [...]The post Approximating Markov's equation first appeared on John D. Cook.

I've written a couple posts lately about reverse engineering the internal state of a random number generator, first Mersenne Twister then lehmer64. This post will look at xorshift128, implemented below. import random # Seed the generator state a: int = random.getrandbits(32) b: int = random.getrandbits(32) c: int = random.getrandbits(32) d: int = random.getrandbits(32) MASK = [...]The post Recovering the state of xorshift128 first appeared on John D. Cook.

If you look very closely at my previous post, you'll notice that I initialize a 128-bit integer with a 64-bit value. The 128-bit unsigned integer represents the internal state of a random number generator. Why not initialize it to a 128-bit value? I was trying to keep the code simple. A surprising feature of C [...]The post Initialize and print 128-bit integers in C first appeared on John D. Cook.

A couple days ago I wrote about hacking the Mersenne Twister. I explained how to recover the random number generator's internal state from a stream of 640 outputs. This post will do something similar with the lehmer64 random number generator. This generator is very simple to implement. Daniel Lemire found it to be the fastest [...]The post Hacking the lehmer64 RNG first appeared on John D. Cook.

This morning I wrote a post about the probability that a random matrix over a finite field is invertible. If the field has q elements and the matrix has dimensions n * n then the probability is In that post I made observation that p(q, n) converges very quickly as a function of n [1]. [...]The post Euler function first appeared on John D. Cook.

What is the inverse of shifting a sequence to the right? Shifting it to the left, obviously. But wait a minute. Suppose you have a sequence of eight bits abcdefgh and you shift it to the right. You get 0abcdefg. If you shift this sequence to the left you get abcdefg0 You can't recover the [...]The post Inverse shift first appeared on John D. Cook.