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 of the Fourier coefficients as elements of an infinite vector then the Hilbert transform can be represented as multiplying by an infinite block matrix. I rarely see infinite matrices except in older math books.…

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) where x, y, u, and v are all real-valued. Not only that, but if f is an elementary function, so are u and v. Here “elementary” has a technical meaning, but essentially it means functions that you could evaluate on a scientific…

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 equations here. The equations come from a paper by Henry G. Baker, cited in the linked page. I wrote up Baker’s equations in LaTeX, then used ChatGPT to generate Python code from the LaTeX…

“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 back to some point in the domain. Often this works as expected, but sometimes it doesn’t. In the previous post I said You can relate each trig function “foo” with its hyperbolic…

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 rotate z a quarter turn to the left (i.e. multiply by i) and then take the cosine. For another example, sinh(z) = −i sin(iz). This says that you can calculate the hyperbolic sine of z by rotating z to the left, taking the sine, and then rotating to the right. You can relate each trig function…

How should we define √(z² − 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. When x is a non-negative real number, √x is defined to be the non-negative real number whose square is x. When x is a complex number √x is defined to be a function that extends √x from the real line to the complex plane by analytic continuation. But we can’t extend √x as an analytic function to the entire complex plane…

The previous post derived the identity and said in a footnote that the identity holds at least for x > 1 and y > 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 for x > 1 and y > 1, and more, but not everywhere. If we combine the two square roots and plot again we still get a valid identity for x > 1 and y > 1, but the plot changes. This is because √a √b does not necessarily equal…

Markov numbers are integer solutions to x² + y² + z² = 3xyz. The Wikipedia article on Markov numbers mentions that Don Zagier studied Markov numbers by looking the approximating equation x² + y² + z² = 3xyz + 4/9 which is equivalent to f(x) + f(y) = f(z) where f(t) is defined as arccosh(3t/2). It wasn’t clear to me why the two previous equations are equivalent, so I’m writing this post to show that they are equivalent. Examples Before showing the equivalence of Zagier’s two equations, let’s…

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 = 0xFFFFFFFF def xorshift128() -> int: global a, b, c, d t = d s = a t ^= (t << 11) & MASK t ^= (t >> 8) & MASK s ^= (s >>…

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 compilers, at least of GCC and Clang, is that you cannot initialize a 128-bit integer to a 128-bit integer literal. You can’t directly print a 128-bit integer either, which is why the…

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 conventional random number generator that can pass Big Crush,” a well-respected test for pseudorandom number generators. Implementing lehmer64 The lehmer64 generator can…

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]. One way to see that the convergence is quick is to note that and John Baez pointed out in the comments that p(q, ∞) = φ(1/q) where φ is the Euler function. Euler was extremely prolific, and many things are named…

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 last element h because the right-shift destroyed information about h. A left-shift doesn’t fully recover a right-shift, and yet surely left shift and right shift are in some sense inverses. Yesterday I wrote a…

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 you mean by the inverse. Suppose you have a matrix M filled with 0s and 1s and you’re looking for a matrix N such that MN is the identity matrix. Do you want the entries of N to also be 0s and 1s? And when you…

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 most of basic linear algebra works the same over every field [1]. In particular, we can do linear algebra over a finite field, and we’re interested in the most finite of finite…

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 more common and more efficient, but the linear algebra approach is conceptually simpler. How MT works There are numerous variations on the Mersenne Twister. We’ll focus on…

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 when f has a fairly simple expression, the expression for κ can be complicated. If we define then the level set of f(x, y) = c is an equilateral triangle when c = −4. The level sets are smoothed triangles for −4 < c < 0. The curvature of these level sets at any point is given by Simplification But there is one instance…

The previous post constructed a triangular analog of the squircle, the unit circle in the p-norm where p is typically around 4. The case p = 2 is a Euclidean circle and the limit as p → ∞ is a Euclidean square. The previous post introduced three functions Li(x, y) such the level set of each function forms a side of a triangle. Then it introduced a soft penalty for each L being away 1, and the level sets of that penalty function formed the rounded triangles we were looking for. Another approach…

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 at the corners. See, for example, this post. Suppose the sides of our triangle are given by L1(x, y) = 1 for i = 1, 2, 3. For example, We design a function f(x, y) as a soft penalty for…

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 keys. In particular, I want the key functionality, not the key name, to be the same across operating systems. For example, the Command key on a Mac does what the Control…

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 my conversation and summarized it in a tweet: category theory is just some straightforwardly useful stuff for some purposes in some fields! you can elegantly simplify…

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 couple bytes longer. (See why here.) Maybe adding one non-ASCII character could cause the file to include a glyph it didn’t before. I did a couple experiments. I made a minimal LaTeX file with only the text See pages 9-10.…

I saw a post on X that plotted the function (log x)² + (log y)² = 1. Of course the plot of x² + y² = 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, 3, …, 10. ContourPlot[Log[x]^2 + Log[y]^2, {x, 0, 10}, {y, 0, 10}, Contours -> Range[10]] The dark blue contour near the origin reminded me of a guitar pick, so I decided to take a stab at creating an equation for the shape of a guitar pick. I…

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 part of a more general series The first-order approximation has error O(x4), as shown in the earlier post. Adding the second-order term makes the error O(x6), and adding the…

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 ramp function for obvious reasons. The function is simply r(x) = max(0, x). Pointwise derivative The pointwise derivative would be 0 for x < 0, 1…

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 right multiple you can cancel out the second-order term as well. For an example, I’d like to approximate the Bessel function J0(x) by the better known cosine function.…

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 length b of the chord of half the arc, illustrated below. Here θ is the central angle of the arc. Then the length of the arc, rθ, is approximately a = rθ ≈ b²/(c + 4b). If the arc is moderately small, the approximation is very accurate. This approximation is simple, accurate, and not obvious, much like the one in this post Derivation Let φ = 4θ. Then the angle…

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 equation and then adjusting the period. You can also find an exact solution, but not in terms of elementary functions; you have to use Jacobi elliptic functions. These are…

There’s a nice formula for the nth derivative of a product. It looks a lot like the binomial theorem. There is also a formula for the nth 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. And here’s the general rule in all its glory. Source: V. F. Ivanoff. The nth Derivative of a Fractional Function. The American Mathematical Monthly, Vol. 55, No. 8 (Oct., 1948), p.…

The equation of motion for a pendulum is the differential equation where g 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 rewrite the equation as Questions This raises a lot of questions, or at least it should. Why not leave sin θ alone? What justifies replacing sin θ with just θ? How small does θ have…