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.

The two latest posts have involved invertible matrices with 0 and 1 entries. If you fill an n * n matrix with 0s and 1s randomly, how likely is it to be invertible? What kind of inverse? There are a couple ways to find the probability that a binary matrix is invertible, depending on what [...]The post Probability that a random binary matrix is invertible first appeared on John D. Cook.

The previous post looked at the tempering step of the Mersenne Twister, formulating a sequence of bit operations as multiplication by a matrix mod 2. This post will look at the components more closely. The theorems of linear algebra generally hold independent of the field of scalars. Typically the field is or , but [...]The post The linear algebra of bit twiddling first appeared on John D. Cook.

The Mersenne Twister (MT) is a random number generator with good statistical properties but bad cryptographic properties. In buzzwords, it's a PRNG but not a CSPRNG. This post will show how the internal state of a MT generator can be recovered from its output. We'll do this using linear algebra. The bit twiddling approach is [...]The post Reverse engineering Mersenne Twister with Linear Algebra first appeared on John D. Cook.

Curvature is conceptually simple but usually difficult to calculate. For a level set curve f(x,y) = c, such as in the previous couple posts, the equation for curvature is Even whenf has a fairly simple expression, the expression for can be complicated. If we define then the level set of f(x,y) = c is [...]The post Calculating curvature first appeared on John D. Cook.

The previous post constructed a triangular analog of the squircle, the unit circle in thep-norm wherep is typically around 4. The casep = 2 is a Euclidean circle and the limit as p is a Euclidean square. The previous post introduced three functions Li(x,y) such that the level set of each function forms [...]The post Smoothed polygons first appeared on John D. Cook.

TimF left a comment on my guitar pick post saying the image was a squircle-ish analog for an isosceles triangle." That made me wonder what a more direct analog of the squircle might be for a triangle. A squircle is not exactly a square with rounded corners. The sides are continuously curved, but curved most [...]The post Triangular analog of the squircle first appeared on John D. Cook.

I try to maintain a consistent work environment across computers that I use. The computers differ for important reasons, but I'd rather they not differ for unimportant reasons. Unified keys One thing I do is remap keys so that the same key does the same thing everywhere, to the extent that's practical. This requires remapping [...]The post Unified config files first appeared on John D. Cook.

Yesterday a friend and I had a conversation about category theory, how it can be a useful pattern description language, but also about how people have unrealistic expectations for it, believing category theory can deliver something for nothing. Later I ran across the following post from Qiaochu Yuan. It felt as if he had overheard [...]The post The mythology of category theory first appeared on John D. Cook.

I saw a post on X saying Changing a hyphen to an en-dash increases your PDF file size by ~10 bytes. My first thought was that it had something to do with hyphen being an ASCII character and an en-dash not. Changing a hyphen to an en-dash would make a UTF-8 encoded text file a [...]The post Changing one character in a PDF first appeared on John D. Cook.

I saw a post on X that plotted the function (logx)^2 + (logy)^2 = 1. Of course the plot of x^2 + y^2 = 1 is a circle, but I never thought what taking logs would do to the shape. Here's what the contours look like setting the right hand side equal to 1, 2, [...]The post The shape of a guitar pick first appeared on John D. Cook.

A couple days ago I wrote a post about turning a trick into a technique, finding another use for a clever way to construct simple, accurate approximations. I used as my example approximating the Bessel function J(x) with (1 + cos(x))/2. I learned via a helpful comment on Mathstodon that my approximation was the first-order [...]The post Approximating even functions by powers of cosine first appeared on John D. Cook.

When a function is not differentiable in the classical sense there are multiple ways to compute a generalized derivative. This post will look at three generalizations of the classical derivative, each applied to the ReLU (rectified linear unit) function. The ReLU function is a commonly used activation function for neural networks. It's also called the [...]The post Three ways to differentiate ReLU first appeared on John D. Cook.

Someone said a technique is a trick that works twice. I wanted to see if I could get anything interesting by turning the trick in the previous post into a technique. The trick created a high-order approximation by subtracting a multiple one even function from another. Even functions only have even-order terms, and by using [...]The post Turning a trick into a technique first appeared on John D. Cook.

Suppose you have an arc a, a portion of a circle of radius r, and you know two things: the length c of the chord of the arc, and the lengthbof the chord of half the arc, illustrated below. Here is the central angle of the arc. Then the length of the arc, r, [...]The post Circular arc approximation first appeared on John D. Cook.

The previous post looks at the nonlinear pendulum equation and what difference it makes to the solutions if you linearize the equation. If the initial displacement is small enough, you can simply replace sin with . If the initial displacement is larger, you can improve the accuracy quite a bit by solving the linearized [...]The post Closed-form solution to the nonlinear pendulum equation first appeared on John D. Cook.

There's a nice formula for thenth derivative of a product. It looks a lot like the binomial theorem. There is also a formula for thenth derivative of a quotient, but it's more complicated and less known. We start by writing the quotient rule in an unusual way. Applying the quotient rule twice gives the following. [...]The post nth derivative of a quotient first appeared on John D. Cook.

The equation of motion for a pendulum is the differential equation whereg is the acceleration due to gravity and is the length of the pendulum. When this is presented in an introductory physics class, the instructor will immediately say something like we're only interested in the case where is small, so we can [...]The post How nonlinearity affects a pendulum first appeared on John D. Cook.

The previous post gave a simple and accurate approximation for the smaller angle of a right triangle. Given aright triangle with sides a,b, andc, wherea is the shortest side andc is the hypotenuse, the angle opposite side a is approximately in radians. The previous post worked in degrees, but here we'll use radians. If the [...]The post Approximation to solve an oblique triangle first appeared on John D. Cook.

Suppose you have a right triangle with sidesa,b, andc, wherea is the shortest side andc is the hypotenuse. Then the following approximation from [1] for the angle Aopposite side a seems too simple and too accurate to be true. In degrees, A a 172 / (b + 2c). The approximation above only involves simple [...]The post Simple approximation for solving a right triangle first appeared on John D. Cook.

AI coding agents improved greatly last summer, and again last December-January. Here are my experiences since my last post on the subject. The models feel subjectively much smarter. They can accomplish a much broader range of tasks. They seem to have a larger, more comprehensive in-depth view of the code base and what you are [...]The post An AI Odyssey, Part 4: Astounding Coding Agents first appeared on John D. Cook.

A few days ago I wrote a post on Newton's diameter theorem. The theorem says to plot the curve formed by the solutions to f(x, y) = 0 where f is a polynomial in x and y of degree n. Next plot several parallel lines that cross the curve at n points and find the [...]The post More on Newton's diameter theorem first appeared on John D. Cook.

The previous post looked at the FP4 4-bit floating point format. This post will look at another 4-bit floating point format, NF4, and higher precision analogs. NF4 and FP4 are common bitsandbytes 4-bit data types. If you download LLM weights from Hugging Face quantized to four bits, the weights might be in NF4 or FP4 [...]The post Gaussian distributed weights for LLMs first appeared on John D. Cook.

In ancient times, floating point numbers were stored in 32 bits. Then somewhere along the way 64 bits became standard. The C programming language retains the ancient lore, using float to refer to a 32-bit floating point number and double to refer to a floating point number with double the number of bits. Python simply [...]The post 4-bit floating point FP4 first appeared on John D. Cook.