[{"id":247017,"date":"2026-05-11T10:30:30","date_gmt":"2026-05-11T15:30:30","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247017"},"modified":"2026-05-11T10:31:05","modified_gmt":"2026-05-11T15:31:05","slug":"inverse-shift","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/11\/inverse-shift\/","title":{"rendered":"Inverse shift"},"content":{"rendered":"

What is the inverse of shifting a sequence to the right? Shifting it to the left, obviously.<\/p>\n

But wait a minute. Suppose you have a sequence of eight bits<\/p>\n

abcdefgh<\/em><\/p>\n

and you shift it to the right. You get<\/p>\n

0abcdefg.<\/em><\/p>\n

If you shift this sequence to the left you get<\/p>\n

abcdefg<\/em>0<\/p>\n

You can’t recover the last element\u00a0h<\/em> because the right-shift destroyed information about\u00a0h<\/em>.<\/p>\n

A left-shift doesn’t fully recover a right-shift, and yet surely left shift and right shift are in some sense inverses.<\/p>\n

Yesterday<\/a> I wrote a post about representing bit manipulations, including shifts, as matrix operators. The matrix corresponding to shifting right by\u00a0k<\/em> bits has 1s on the\u00a0k<\/em>th diagonal above the main diagonal and 0s everywhere else. For example, here is the matrix for shifting an 8-bit number right two bits. A black square represents a 1 and a white square represents a 0.<\/p>\n

<\/p>\n

This matrix isn’t invertible. When you’d like to take the inverse of a non-invertible matrix, your kneejerk response should be to compute the pseudoinverse. (Technically the Moore-Penrose pseudoinverse<\/a>. There are other<\/a> pseudoinverses, but Moore-Penrose is the most common.)<\/p>\n

As you might hope\/expect, the pseudoinverse of a right-shift matrix is a left-shift matrix. In this case the pseudoinverse is simply the transpose, though of course that isn’t always the case.<\/p>\n

<\/p>\n

If you’d like to prove that the pseudoinverse of a matrix that shifts right by k<\/em> places is a matrix that shifts left by k<\/em> places, you don’t have to compute the pseudo inverse per se: you can verify your guess. This post<\/a> gives four requirements for a pseudoinverse. You can prove that left shift is the inverse of right shift by showing that it satisfies the four equations.<\/p>\n","protected":false},"excerpt":{"rendered":"

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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[5],"tags":[198],"class_list":["post-247017","post","type-post","status-publish","format-standard","hentry","category-computing","tag-linear-algebra"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247017","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247017"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247017\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247017"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247017"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247017"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247016,"date":"2026-05-11T08:25:13","date_gmt":"2026-05-11T13:25:13","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247016"},"modified":"2026-05-11T08:25:13","modified_gmt":"2026-05-11T13:25:13","slug":"random-binary-matrices","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/11\/random-binary-matrices\/","title":{"rendered":"Probability that a random binary matrix is invertible"},"content":{"rendered":"

The two latest posts have involved invertible matrices with 0 and 1 entries. If you fill an n<\/em> \u00d7 n<\/em> matrix with 0s and 1s randomly, how likely is it to be invertible?<\/p>\n

What kind of inverse?<\/h2>\n

There are a couple ways to find the probability that a binary matrix is invertible, depending on what you mean by the inverse.<\/p>\n

Suppose you have a matrix M<\/em> filled with 0s and 1s and you’re looking for a matrix N<\/em> such that MN<\/em> is the identity matrix. Do you want the entries of N<\/em> to also be 0s and 1s? And when you multiply the matrices, are you doing ordinary integer arithmetic or are you working mod 2?<\/p>\n

In the previous posts we were working over GF(2), the field with two elements, 0 and 1. All the elements of a matrix are either 0 or 1, and arithmetic is carried out mod 2. In that context there’s a nice expression for the probability a square matrix is invertible.<\/p>\n

If you’re working over the real numbers, the probability of binary matrix being invertible is higher. One way to see this is that the inverse of a binary matrix is allowed to be binary but it isn’t required to be.<\/p>\n

Another way to see this is to look at determinants. If you think of a matrix M<\/em> as a real matrix whose entries happen to only be 0 or 1, M<\/em> is invertible if and only its determinant is non-zero. But if you think of M<\/em> as a matrix over GF(2), the entries are either 0 or 1 out of necessity, and M<\/em> is invertible if and only if its determinant, computed in<\/em> GF(2), is non-zero. If the determinant of M<\/em> as a real matrix is a non-zero even number, then M<\/em> is invertible as a real matrix but not as a matrix over GF(2).<\/p>\n

Probability of invertibility in GF(2)<\/h2>\n

Working over GF(2), what is the probability that a random matrix is invertible? Turns out it’s just as easy to answer a more general question: what is the probability that a random n<\/em> \u00d7\u00a0n<\/em> matrix over GF(q<\/em>), a finite field with q<\/em> elements, is invertible? This is<\/p>\n

\"\\prod_{i=1}^n\\left(<\/p>\n

When q<\/em> = 2 and n<\/em> = 8 this probability is 0.289919. The probability is roughly the same for all larger values of n<\/em>, converging to approximately 0.288788 as n<\/em> \u2192 \u221e.<\/p>\n

Probability of invertibility in \u211d<\/h2>\n

What is the probability that an 8 \u00d7 8 matrix with random 0 and 1 entries is invertible as a real matrix? We can estimate this by simulation.<\/p>\n

import numpy as np\r\n\r\ndef simulate_prob_invertible_real(n, numreps=1000):\r\n    s = 0\r\n    for _ in range(numreps):\r\n        M = np.random.randint(0, 2, size=(n, n))\r\n        det = np.linalg.det(M)\r\n        if abs(det) > 1e-9:\r\n            s += 1\r\n    return s\/numreps\r\n<\/pre>\n
When n<\/em> = 8, I got 0.5477 when running the code with 10,000 reps.<\/p>\n
When n<\/em> = 32, I got a probability of 1. Obviously it is possible for a 32 \u00d7 32 binary matrix to be singular, but it’s very unlikely: it didn’t happen in 10,000 random draws.<\/p>\n","protected":false},"excerpt":{"rendered":"
The two latest posts have involved invertible matrices with 0 and 1 entries. If you fill an n \u00d7 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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[9],"tags":[198,105],"class_list":["post-247016","post","type-post","status-publish","format-standard","hentry","category-math","tag-linear-algebra","tag-probability"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247016","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247016"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247016\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247016"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247016"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247016"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247014,"date":"2026-05-10T13:51:40","date_gmt":"2026-05-10T18:51:40","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247014"},"modified":"2026-05-10T20:31:20","modified_gmt":"2026-05-11T01:31:20","slug":"the-linear-algebra-of-bit-twiddling","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/10\/the-linear-algebra-of-bit-twiddling\/","title":{"rendered":"The linear algebra of bit twiddling"},"content":{"rendered":"
The previous post<\/a> 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.<\/p>\n
The theorems of linear algebra generally hold independent of the field of scalars. Typically the field is \u211d or \u2102, 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 fields GF(2), the field with just two elements, 0 and 1.<\/p>\n
In GF(2), addition corresponds to XOR. We will denote this by \u2295 to remind us that although it’s addition, it’s not the usual addition, i.e. 1 \u2295 1 = 0. Similarly, multiplication corresponds to AND. We’ll work with 8-bit numbers to make the visuals easier to see.<\/p>\n
Shifting a number left one bit corresponds to multiplication by a matrix with 1’s below the diagonal main. Shifting left by k<\/em> bits is the same as shifting left by 1 bit k<\/em> times, so the the matrix representation for\u00a0x<\/em> << k<\/em> is the\u00a0k<\/em>th power of the matrix representation of shifting left once. This matrix has 1s on the\u00a0k<\/em>th diagonal below the main diagonal. Below is the matrix for shifting left two bits,\u00a0x<\/em> <<\u00a0k<\/em>.<\/p>\n
<\/p>\n
Right shifts are the mirror image of left shifts. Here’s the matrix for shifting right two bits, x<\/em> >> k<\/em>.<\/p>\n
<\/p>\n
Shifts are not fully invertible because bits either fall off the left or the right end. The steps in the Mersenne Twister are invertible because shifts are always XOR’d with the original argument. For example, although the function that takes x<\/em> to x<\/em> >> 2 is not invertible, the function that takes x<\/em> to x<\/em> \u2295 (x<\/em> >> 2) is invertible. This operation corresponds to the matrix below.<\/p>\n
<\/p>\n
This is an upper triangular matrix, so its determinant is the product of the diagonal elements. These are all 1s, so the determinant is 1, and the matrix is invertible.<\/p>\n
Bitwise AND multiplies each bit of the input by the corresponding bit in another number known as the mask. The bits aligned with a 1 are kept and the bits aligned with a 0 are cleared. This corresponds to multiplying by a diagonal matrix whose diagonal elements correspond to the bits in the mask. For example, here is the matrix that corresponds to taking the bitwise AND with 10100100.<\/p>\n
<\/p>\n
Each of the steps in the Mersenne Twister tempering process are invertible because they all correspond to triangular matrices with all 1’s on the diagonal. For example, the line<\/p>\n
y ^= (y <<  7) & 0x9d2c5680 <\/pre>\n
says to shift the bits of y<\/code> left 7 places, then zero out the elements corresponding to 0s in the mask, then XOR the result with y<\/em>. In matrix terms, we multiply by a lower triangular matrix with zeros on the main diagonal, then multiply by a diagonal matrix that zeros out some of the terms, then add the identity matrix. So the matrix corresponding to the line of code above is lower triangular, with all 1s on the diagonal, so it is invertible. <\/p>\n
[1] Until you get to eigenvalues. Then it matters whether the field is algebraically complete, which no finite field is.<\/p>\n","protected":false},"excerpt":{"rendered":"
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 \u211d or \u2102, but […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[5,9],"tags":[198],"class_list":["post-247014","post","type-post","status-publish","format-standard","hentry","category-computing","category-math","tag-linear-algebra"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247014","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247014"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247014\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247014"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247014"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247014"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247013,"date":"2026-05-10T12:32:02","date_gmt":"2026-05-10T17:32:02","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247013"},"modified":"2026-05-10T13:52:53","modified_gmt":"2026-05-10T18:52:53","slug":"reverse-mersenne-twister","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/10\/reverse-mersenne-twister\/","title":{"rendered":"Reverse engineering Mersenne Twister with Linear Algebra"},"content":{"rendered":"
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.<\/p>\n
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.<\/p>\n
How MT works<\/h2>\n
There are numerous variations on the Mersenne Twister. We’ll focus on the original version that had internal state consists of vector x<\/code> of length 640 filled with 32-bit numbers. The ideas in this post would apply equally to all MT versions.<\/p>\n
The first call to MT returns a “tempered” version of\u00a0x[0]<\/code>. The next call returned a tempered version of x[1]<\/code>, and so on. After every 640 calls, the state is scrambled. This scrambling is where the “twist” in the name Mersenne Twister comes from. (The Mersenne part comes from the fact that the period of an MT generator is a Mersenne prime<\/a>.)<\/p>\n
Tempering<\/h2>\n
Here is Python code for performing the tempering step.<\/p>\n
def temper(y):\r\n    y ^= (y >> 11) \r\n    y ^= (y <<  7) & 0x9d2c5680 \r\n    y ^= (y << 15) & 0xefc60000 \r\n    y ^= (y >> 18)  \r\n    return y\r\n<\/pre>\n
Each step is reversible, and so the temper<\/code> function is reversible.<\/p>\n
Because the tempering step is reversible, the first output can be used to infer the first element of the internal state, the second output to infer the second element, and so on. After 640 calls one can know the entire internal state and predict the rest of the generator’s output from then on.<\/p>\n
Linear algebra<\/h2>\n
The bitwise operations above all correspond to linear operations over GF(2), the field with just two elements, 0 and 1. Addition in this field is XOR and multiplication is AND.<\/p>\n
Each step corresponds to multiplying a vector of 32 bits on the left by a 32 \u00d7 32 matrix with entries that are 0’s and 1’s, with the understanding that the sum of two bits is their XOR and the product of two bits is their AND. Equivalently, arithmetic is carried out mod 2. So you can compute the matrix-vector product as ordinary integers if you then reduce every integer mod 2.<\/p>\n
We will find the matrix M<\/em> that corresponds to the temper operation and prove that it is invertible by finding its inverse. This proves that tempering is invertible, and one could compute the inverse of tempering by multiplying by M<\/em>\u22121<\/sup>, though it would be more efficient to undo temporing by bit twiddling.<\/p>\n
One way to recover a matrix is to multiplying by unit vectors e<\/em>i<\/em><\/sub> where the i<\/em>th component of e<\/em>i<\/em><\/sub> is 1 and the rest of the components are zero. Then<\/p>\n
M<\/em> e<\/em>i<\/em><\/sub><\/p>\n
is the\u00a0i<\/em>th column of M<\/em>.<\/p>\n
So we can find the\u00a0n<\/em>th column of\u00a0M<\/em> by tempering 1 <<\u00a0n<\/em> = 2n<\/em><\/sup>.<\/p>\n
M = np.zeros((32, 32), dtype=int)\r\nfor i in range(32):\r\n    t = temper(1 << (31-i))\r\n    s = f\"{t:032b}\"\r\n    for j in range(32):\r\n        M[j, i] = int(s[j])\r\n<\/pre>\n
Let’s generate a random number and compute its tempered form two ways: directly and matrix multiplication.<\/p>\n
x = random.getrandbits(32)\r\ny = temper(x)\r\nprint(f\"{y:032b}\")\r\n\r\nvx = np.array([int(b) for b in f\"{x:032b}\"]) # vector form of x\r\nvy = M @ vx % 2 # vector form of y\r\nprint(\"\".join(str(b) for b in vy))\r\n<\/pre>\n
Both produce the same bits:<\/p>\n
10100101100101101100110101000110\r\n10100101100101101100110101000110\r\n<\/pre>\n
We can find the matrix representation of the untemper function by inverting the matrix M<\/em>. However, we need to invert it over the field GF(2), not over the integers or reals.<\/p>\n
import galois\r\nGF2 = galois.GF(2)\r\nMinv = np.linalg.inv(GF2(M))\r\n<\/pre>\n
Here are visualizations of M<\/em> and its inverse using a black square for a 1 and a white square for a 0.<\/p>\n
M<\/em>:<\/p>\n
<\/p>\n
M<\/em>\u22121<\/sup>:<\/p>\n
<\/p>\n
The next post<\/a> will back up and look at the linear algebra of the components that comprise the Mersenne Twister.<\/p>\n","protected":false},"excerpt":{"rendered":"
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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[9],"tags":[198,118],"class_list":["post-247013","post","type-post","status-publish","format-standard","hentry","category-math","tag-linear-algebra","tag-rng"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247013","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247013"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247013\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247013"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247013"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247013"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247012,"date":"2026-05-08T08:54:18","date_gmt":"2026-05-08T13:54:18","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247012"},"modified":"2026-05-08T08:54:18","modified_gmt":"2026-05-08T13:54:18","slug":"calculating-curvature","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/08\/calculating-curvature\/","title":{"rendered":"Calculating curvature"},"content":{"rendered":"
Curvature is conceptually simple but usually difficult to calculate. For a level set curve f<\/em>(x<\/em>,\u00a0y<\/em>) = c, such as in the previous couple posts, the equation for curvature is<\/p>\n
\"\\kappa<\/p>\n
Even when\u00a0f<\/em> has a fairly simple expression, the expression for \u03ba can be complicated.<\/p>\n
If we define<\/p>\n
\"f(x,<\/p>\n
then the level set of f<\/em>(x<\/em>,\u00a0y<\/em>) = c is an equilateral triangle when\u00a0c<\/em> = \u22124. The level sets are smoothed triangles for \u22124 <\u00a0c<\/em> < 0.<\/p>\n
<\/p>\n
The curvature of these level sets at any point is given by<\/p>\n
\"\\frac{\\left|2<\/p>\n
Simplification<\/h2>\n
But there is one instance in which curvature is easy to calculate. For the graph of a function\u00a0g<\/em>(x<\/em>) =\u00a0y<\/em>, the curvature is approximately the absolute value of the second derivative of\u00a0g<\/em>, provided the first derivative is small.<\/p>\n
\"<\/p>\n
At a local maximum or local minimum of\u00a0g<\/em>(x<\/em>) the approximation is exact because the first derivative is zero.<\/p>\n
Max and min curvature of smoothed triangles<\/h2>\n
This means that in the example above, we can calculate the maximum and minimum curvature of the level sets. The maximum curvature occurs in the corners and the minimum occurs in the middle of the sides. By symmetry there are three maxima and three minima, but we can take the ones corresponding to\u00a0x<\/em> = 0 for convenience. There we find the curvature is simply<\/p>\n
\"|6<\/p>\n
When\u00a0x<\/em> = 0, we have<\/p>\n
\"f(x,y)<\/p>\n
and so the maximum and minimum curvature are the two roots of the cubic equation for\u00a0y<\/em> that lie in the interval [\u22121, 2]. (There’s another root greater than 2.)<\/p>\n
For example, when c<\/em> = \u22123, the roots are 0.8794, 1.3473, and 2.5321. The first root corresponds to the minimum curvature, the second to the maximum, and the third is outside our region of interest. It follows that the minimum curvature is 0.7237 and the maximum is 14.0838.<\/p>\n
When c<\/em> = \u22121 the minimum and maximum curvature are 2.80747 and 9.91622 respectively,<\/p>\n
Since c<\/em> = \u22124 corresponds to the triangle, the minimum curvature is 0 and the maximum is infinite. As c<\/em> increases, the minimum and maximum curvature come together because the level set is becoming more round.<\/p>\n","protected":false},"excerpt":{"rendered":"
Curvature is conceptually simple but usually difficult to calculate. For a level set curve f(x,\u00a0y) = c, such as in the previous couple posts, the equation for curvature is Even when\u00a0f has a fairly simple expression, the expression for \u03ba can be complicated. If we define then the level set of f(x,\u00a0y) = c is […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[9],"tags":[48],"class_list":["post-247012","post","type-post","status-publish","format-standard","hentry","category-math","tag-differential-geometry"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247012","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247012"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247012\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247012"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247012"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247012"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247009,"date":"2026-05-07T13:19:24","date_gmt":"2026-05-07T18:19:24","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247009"},"modified":"2026-05-08T08:48:09","modified_gmt":"2026-05-08T13:48:09","slug":"smoothed-polygons","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/07\/smoothed-polygons\/","title":{"rendered":"Smoothed polygons"},"content":{"rendered":"
The previous post constructed a triangular analog of the squircle, the unit circle in the\u00a0p<\/em>-norm where\u00a0p<\/em> is typically around 4. The case\u00a0p<\/em> = 2 is a Euclidean circle and the limit as p<\/em> \u2192 \u221e is a Euclidean square.<\/p>\n
The previous post introduced three functions L<\/em>i<\/em><\/sub>(x<\/em>,\u00a0y<\/em>) such that the level set of each function<\/p>\n
\"\\{(x,y)<\/p>\n
forms a side of a triangle. Then it introduced a soft penalty for each L<\/em> being away 1, and the level sets of that penalty function formed the rounded triangles we were looking for.<\/p>\n
Another approach would be to change the\u00a0L<\/em>‘s slightly so that the sides are the levels sets L<\/em>i<\/em><\/sub>(x<\/em>,\u00a0y<\/em>) = 0. The advantage to this formulation is that the product of three numbers is 0 if and only if one of the numbers is zero. That means if we define<\/p>\n
\"f(x,\u00a0y)
\nthen the set of points<\/p>\n
\"\\{(x,<\/p>\n
corresponds to the three lines when c<\/em> = 0 and the level sets\u00a0for small\u00a0c<\/em> > 0 are nearly triangles. The level sets will be smooth if the gradient is non-zero, i.e. c<\/em> is not zero.<\/p>\n
This approach is not unique to triangles<\/strong>. You could create smooth approximations any polygon by multiplying linear functions that define the sides. Or you could do something similar with curved arcs.<\/p>\n
If we define our\u00a0L’<\/em>s by<\/p>\n
\"\\begin{align*}<\/p>\n
then our curves will be the level sets of<\/p>\n
\"f(x,<\/p>\n
A few level sets are shown below. The level set for\u00a0c<\/em> = 0 is the straight lines.<\/p>\n
Note the level sets are not connected. If you’re interested in approximate triangles, you want the components that are inside the triangle.<\/p>\n
<\/p>\n
 <\/p>\n","protected":false},"excerpt":{"rendered":"
The previous post constructed a triangular analog of the squircle, the unit circle in the\u00a0p-norm where\u00a0p is typically around 4. The case\u00a0p = 2 is a Euclidean circle and the limit as p \u2192 \u221e is a Euclidean square. The previous post introduced three functions Li(x,\u00a0y) such that the level set of each function forms […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[9],"tags":[],"class_list":["post-247009","post","type-post","status-publish","format-standard","hentry","category-math"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247009","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247009"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247009\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247009"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247009"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247009"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247007,"date":"2026-05-06T14:53:29","date_gmt":"2026-05-06T19:53:29","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247007"},"modified":"2026-05-06T14:53:29","modified_gmt":"2026-05-06T19:53:29","slug":"triangular-analog-of-the-squircle","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/06\/triangular-analog-of-the-squircle\/","title":{"rendered":"Triangular analog of the squircle"},"content":{"rendered":"
TimF left a comment on my guitar pick post<\/a> 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.<\/p>\n
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<\/a>.<\/p>\n
Suppose the sides of our triangle are given by L<\/em>1<\/sub>(x<\/em>, y<\/em>) = 1 for i<\/em> = 1, 2, 3. For example, <\/p>\n
\\begin{align*}\nL_1(x, y) &= y \\\\\nL_2(x, y) &= \\phantom{-}\\frac{\\sqrt{3}}{2}x - \\frac{1}{2}y \\\\\nL_3(x, y) &= -\\frac{\\sqrt{3}}{2}x - \\frac{1}{2}y \\\\\n\\end{align*}\n<\/p>\n
We design a function f<\/em>(x<\/em>, y<\/em>) as a soft penalty for points not being on one of the sides and look at the set of points f<\/em>(x<\/em>, y<\/em>) = 1.<\/p>\n
f(x, y) = \\left( \\sum_{i=1}^3 \\exp(p (L_i(x, y) - 1)) \\right)^{1\/p}<\/p>\n
You might recognize this as a Lebesgue norm, analogous to the one used to define a squircle.<\/p>\n
The larger p<\/em> is, the heavier the penalty for being far from a side and the closer the level set f<\/em>(x<\/em>, y<\/em>) = 1 comes to being a triangle.<\/p>\n
<\/p>\n","protected":false},"excerpt":{"rendered":"
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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[9],"tags":[],"class_list":["post-247007","post","type-post","status-publish","format-standard","hentry","category-math"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247007","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247007"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247007\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247007"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247007"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247007"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247006,"date":"2026-05-06T11:32:20","date_gmt":"2026-05-06T16:32:20","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247006"},"modified":"2026-05-06T11:32:20","modified_gmt":"2026-05-06T16:32:20","slug":"unified-config-files","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/06\/unified-config-files\/","title":{"rendered":"Unified config files"},"content":{"rendered":"
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.<\/p>\n
Unified keys<\/h2>\n
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<\/a>. In particular, I want the key functionality<\/b>, not the key name<\/b>, to be the same across operating systems. For example, the Command key on a Mac does what the Control key does on Windows and Linux. I have my machines set up so the key in the lower left corner, whatever you call it, does things like copy, paste, and cut.<\/p>\n
Unified config files<\/h2>\n
I use bash everywhere even though zshell is the default on Mac OS. But Linux and MacOS are sufficiently different that I have two .bashrc<\/code> files, one for each OS. However, both .bashrc<\/code> files source a common file .bash_aliases<\/code> for aliases. You can set that up by putting the following code in your .bashrc<\/code> file.<\/p>\n
if [ -f ~\/.bash_aliases ]; then\r\n    . ~\/.bash_aliases\r\nfi\r\n<\/pre>\n
Sometimes you need some kind of branching logic even for two machines running the same OS. For example, if you have two computers, named Castor and Pollux, you might have code like the following in your .bashrc<\/code>.<\/p>\n
HOSTNAME_SHORT=$(hostname -s)\r\nif [ \"$HOSTNAME_SHORT\" = \"Castor\" ]; then\r\n...\r\nelif [ \"$HOSTNAME_SHORT\" = \"Pollux\" ]; then\r\n...\r\n<\/pre>\n
One problem with using hostname is that the OS can change the name on you. At least MacOS will do this; I don’t know whether other operating systems will. If you give your machine an uncommon name then this is unlikely to happen.<\/p>\n
I have one init.el<\/code> file for Emacs. It contains some branching logic for OS:<\/p>\n
(cond\r\n    ((string-equal system-type \"windows-nt\") ; Microsoft Windows\r\n     ...)\r\n    ((string-equal system-type \"gnu\/linux\") ; linux\r\n     ...)\r\n    ((string-equal system-type \"darwin\") ; Mac\r\n     ...)\r\n)\r\n<\/pre>\n
You can also branch on hostname.<\/p>\n
(if (string-equal (system-name) \"Castor\")\r\n...)\r\n(if (string-equal (system-name) \"Pollux\")\r\n...)\r\n<\/pre>\n
This isn’t difficult, but in the short run it’s easier to just make one ad hoc edit to a config file at a time, letting the files drift apart. Along those lines, see bicycle skills<\/a>. <\/p>\n","protected":false},"excerpt":{"rendered":"
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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[5],"tags":[],"class_list":["post-247006","post","type-post","status-publish","format-standard","hentry","category-computing"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247006","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247006"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247006\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247006"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247006"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247006"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247005,"date":"2026-05-06T06:37:42","date_gmt":"2026-05-06T11:37:42","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247005"},"modified":"2026-05-06T06:37:42","modified_gmt":"2026-05-06T11:37:42","slug":"category-mythology","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/06\/category-mythology\/","title":{"rendered":"The mythology of category theory"},"content":{"rendered":"
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.<\/p>\n
Later I ran across the following post<\/a> from Qiaochu Yuan. It felt as if he had overheard my conversation and summarized it in a tweet:<\/p>\n
category theory is just some straightforwardly useful stuff for some purposes in some fields! you can elegantly simplify and streamline some proofs. then there is the mythology of category theory, which is some other thing entirely, mostly wishful thinking and projection afaict<\/p><\/blockquote>\n
His phrase “the mythology of category theory” gives a name to this idea that category theory can deliver specific outputs without specific inputs. It helps to distinguish CT as a scrapbook of patterns from CT as sorcery.<\/p>\n","protected":false},"excerpt":{"rendered":"
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 […]<\/p>\n","protected":false},"author":1,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":{"_acf_changed":false,"footnotes":""},"categories":[9],"tags":[144],"class_list":["post-247005","post","type-post","status-publish","format-standard","hentry","category-math","tag-category-theory"],"acf":[],"aioseo_notices":[],"_links":{"self":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247005","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/comments?post=247005"}],"version-history":[{"count":0,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/posts\/247005\/revisions"}],"wp:attachment":[{"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/media?parent=247005"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/categories?post=247005"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/www.johndcook.com\/blog\/wp-json\/wp\/v2\/tags?post=247005"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}},{"id":247003,"date":"2026-05-05T17:45:26","date_gmt":"2026-05-05T22:45:26","guid":{"rendered":"https:\/\/www.johndcook.com\/blog\/?p=247003"},"modified":"2026-05-05T18:36:26","modified_gmt":"2026-05-05T23:36:26","slug":"changing-one-character-in-a-pdf","status":"publish","type":"post","link":"https:\/\/www.johndcook.com\/blog\/2026\/05\/05\/changing-one-character-in-a-pdf\/","title":{"rendered":"Changing one character in a PDF"},"content":{"rendered":"
I saw a post<\/a> on X saying<\/p>\n
Changing a hyphen to an en-dash increases your PDF file size by ~10 bytes.<\/p>\n
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<\/a>.) Maybe adding one non-ASCII character could cause the file to include a glyph it didn’t before.<\/p>\n
I did a couple experiments. I made a minimal LaTeX file with only the text<\/p>\n
    See pages 9-10.<\/pre>\n
and another with<\/p>\n
    See pages 9--10.<\/pre>\n
(In LaTeX, a hyphen compiles to a hyphen and two hyphens compile to an en-dash.)<\/p>\n
I compiled both files using pdflatex<\/code>. The PDF with the hyphen was 13172 bytes and the one with the en-dash was 13099 bytes. Replacing the hyphen with the en-dash made the file 73 bytes smaller<\/b>.<\/p>\n
I repeated the experiment with a Libre Office ODT document. Changing the hyphen to an en-dash reduced the file size from 13548 bytes to 13514 bytes.<\/p>\n
Then I went back to LaTeX and pasted a paragraph of lorem ipsum text into both files. Now the file with the hyphen produced a 18131 byte PDF and the file with the en-dash produced a 18203 byte PDF. So in this instance changing the hyphen to an en-dash increased<\/b> the size of the file by 72 bytes.<\/p>\n
After adding the lorem ipsum text to the ODT files, the PDF resulting from the file with the hyphen was 23 bytes larger, 17846 bytes versus 17823 bytes.<\/p>\n
PDFs are complicated. Changing one character can make the file bigger or smaller, by an unpredictable amount. It depends, among other things, on what software was used to create the PDF. Even changing one letter in an all-ASCII text can change the size of the PDF, I suppose due to some internal text compression and aesthetic control characters. I don’t pretend to understand what’s going on inside a PDF.<\/p>\n
Related posts<\/h2>\n