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π J₁(1) where J₁ is a Bessel function. This may seem like rabbit out of a hat, but it’s actually a simple calculation given the integral definition of Bessel functions:

Since cosine is even, we can write our integral over [−π, π] as twice the integral over [0, π]. Then a change of variables turns this into the definition of J_n(z) with n = 1 and z = 1.

A deeper question is what have we accomplished by just giving a new name to essentially the same problem we started with. Another question is why in the world are Bessel functions defined as above.

As for what we’ve accomplished, we’ve related out integration problem to a very well-studied function. Bessel functions have been studied for two centuries and it’s easy to find software to evaluate them. Even the usually minimalist command line calculator bc has a function j(x, n) for evaluating J_n(x) for integer values of n. We could calculate our integral to 50 decimal places as follows.

~$ bc -l
>>> scale = 50
>>>  -8*a(1)*j(1,1)
-2.76491937476833705153256665538788207487495025542883

Note that bc doesn’t have a value of π built in, but a(x) evaluates the arctangent function, and π = 4 arctan(1).

There are multiple ways of defining Bessel functions. The three main ways would be in terms of their power series, in terms of the differential equations they solve, and in terms of their integral representation. Friedrich Bessel defined what we now call Bessel functions of the first kind, the J_n functions, in terms of their integral representations.

Why did Bessel care about these integrals? They came out of his calculations in celestial mechanics. One example of this is solving Kepler’s equation with Fourier series; the Fourier coefficients are given by Bessel functions. Bessel functions had occurred in applications before Mr. Bessel drew attention to them and studied them methodically.

Mathematics is developed inductively but taught deductively. It’s common for things to be kicked around for years before someone decides they deserve a name and systematic study. See this post on the central limit theorem for another example. The CLT is older than the Gaussian distribution, even older than Gauss.

Related posts

• FM radio and Bessel functions
• Vibrating circular membranes
• Fourier-Bessel series

Knowing a bit of Latin makes sense of many things that would otherwise seem completely arbitrary, such as why the symbols for gold, silver, and lead are Au, Ag, and Pb respectively.

Similarly, ed(1) is the Latin of Linux [1]. Many conventions in command line utilities follow conventions that go back to the ed(1) line editor. They may go back even further. Just as Latin didn’t come out of nowhere, neither did ed(1), but you can’t go back indefinitely. It’s convenient to start history somewhere, and this post will start with ed(1) just as much discussion of Western linguistics starts with Latin.

The following are features of ed(1) that live on in sed, awk, grep, vi, perl, bash, etc.

1. Using slashes to delimit regular expressions
2. Using $ to indicate the end of a line or the end of a file
3. The pattern of specifying address + action or address range + action
4. Using regular expressions as address ranges
5. Using \1, \2, etc to refer to regex captures
6. Using & to refer to the entire matched text
7. The g/regexp/command pattern
8. Using p for printing lines, as in g/re/p
9. The commands a, c, d, i, j, l, p, q, r, and w in vi
10. ! for shell escape

[1] Because the name “ed” is so short, and looks so much like the name Ed, it’s convenient to use its full Unix name ed(1). The parenthesized number is used to disambiguate different things that have the same name, such as the user command kill(1) and the system call kill(2). There is no ed(2) or any other higher-numbered ed. The number is there to make the name stand out, not to disambiguate anything.

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 rule is on the order of 1/N² where N is the number of integration points. To be more specific, the error in integrating a function f over [a, b] with N points is bounded by

(b − a)³ M / 12N²

where M is the maximum absolute value of the second derivative of f. So in our case we should expect the error to be less than 82.67/N². In fact we do much better than that. The error does not decrease quadratically, as it does in general, but exponentially.

With just 16 integration points, we’ve reached the limit of floating point representation.

exp(−z²) ≈ (1 + cos(sin(z) + z))/2

that I wrote about last week. The two expressions above are close along the real axis but not along the imaginary axis.

If z = iy, the right side grows much faster than the left, behaving like exp(exp(y)).

This led to me looking up the power series for the double exponential function exp(exp(y)). This is an interesting series because the coefficient of xⁿ is

e B_n / n!

where B_n is the nth Bell number, which equals the number of ways to partition a set of n labeled items [1]. And of course n! is the number of ways to permute a set of n labeled items. So the nth coefficient in the power series for exp(exp(y)) is the ratio of the number of partitions to permutations for a set of n labeled things, multiplied by e.

The number of ways to partition a set of n things grows quickly as n increases, almost as quickly as the number of permutations, and so the series for the double exponential function converges very slowly.

Computing

SymPy has a function bell for computing Bell numbers, so you could compute the ratio of partitions to permutations as follows.

from sympy import bell, factorial
f = lambda n: bell(n)/factorial(n)

This returns a number of type sympy.core.numbers.Rational and so the result is exact. You can cast it to float for convenience.

Asymptotics

If we look at only the terms in the asymptotic series for log B_n and log n! that grow with n we have

log B_n ~ n log n − n log log n
log n! ~ n log n − ½ log n

and so

log( B_n / n! ) ~ ½ log n − n log log n

There’s also an asymptotic series for log( B_n / n! ) involving the Lambert W function:

log( B_n / n! ) ~ n/r − 1 − n log r

where r = W(n).

Related posts

• Mr. Bell and Bell numbers
• Estimating partition numbers
• Richard Stanley’s twelvefold way (pdf)

[1] It’s important that the items are labeled. Partition numbers are the number of partitions of an unlabeled set. Partition numbers are much smaller than Bell numbers.

Maybe you decide a tolerance of 10⁻¹² is good enough, and so you sum the terms until the next term to add is below the tolerance.

from math import factorial, exp

def naive_exp(x):
    tolerance = 1e-12
    s = 0
    n = 0
    while True:
        delta = x**n / factorial(n)
        s += delta
        if abs(delta) < tolerance:
            return s
        n += 1

You want to try your program out, so you compute e by calling the function at 1. If you compare this to calling exp(1) you find that you got all the digits correct.

Now you try computing exp(-20). Calling naive_exp(-20) gives

    5.47893091802112e-10

but calling exp(-20) gives

    2.061153622438558e-09

Don’t brush things like this as flukes or compiler bugs [1]. This is your golden opportunity to learn something.

Maybe you add a print statement to see the intermediate values of the sum stored in the variable s. If you do, you’ll see that the partial sums oscillate wildly before settling down.

Maybe that seems wrong, but then you look more carefully at the series. The nth term is xⁿ/n!. Since x is negative, the terms alternate in sign. And the absolute values of the term get bigger before they get smaller. When x = −20, each numerator is 20 times larger than the previous, and each denominator is n times larger than the previous. So the terms will get bigger until n > 20. So the wild oscillations are real, not a bug.

The largest partial sum is 21822593.77927747 in absolute value. You know that exp(−20) is a very small number, so there’s going to have to be a lot of cancellation before the partial sums settle down to a small number. Maybe you’ve heard that cancellation is where numerical calculations lose precision. If not, now you know!

Look again at the largest partial sum. There are eight figures to the right of the decimal point. The code is printing out results to as much precision as it has, so the error at this point is on the order of 10⁻⁸. We’re trying to compute a number on the order of 10⁻⁹, and if any digits in our result are correct, it would be a coincidence.

If you go back and try your code on x = −22, the result is even worse, giving a negative result for a quantity that for theoretical reasons cannot be negative. But you can see why: you’re asking the code to compute a number that is closer to zero than the accuracy of the code.

Computers don’t represent numbers in base 10 internally, but the argument above is sufficient in this case. If you want to dig deeper, look into the anatomy of a floating point number.

There is a simple way around the problem above, but discovering it sooner would short-circuit the learning process. You could calculate exp(−20) as 1/exp(20) and avoid all the cancellation because the series for exp(20) does not alternate.

[1] Compilers do have bugs occasionally, but it’s orders of magnitude more likely that something is wrong with your code.

exp(−x²) ≈ (1 + cos(sin(x) + x))/2

and some are saying this can just be explained by Taylor series: the series for the two sides differ for the first time at the x⁶ term, so that’s why you get a good approximation. As I wrote yesterday, that’s only part of it.

If it were just about Taylor series you could use

exp(−x²) ≈ 1 − x² + x⁴/2

which also has error O(x⁶). But this approximation is only good for fairly small x, say in [−0.5, 0.5], whereas the approximation at the top of the post is good over [−4, 4]. When x = 4, the error in the cosine approximation is 0.002579 but the error in the Taylor approximation is 113, five orders of magnitude larger.

If the accuracy of the cosine approximation were due to Taylor series, then we’d expect the function

exp(−x²) − (1 + cos(sin(x) + x))/2

to be small not just over the interval [−4, 4] but over a disk of radius 4 in the complex plane. But it’s not. When x = 4i the function is on the order of 10¹³.

Both the cosine approximation and the Taylor approximation are good for small disks, say of radius 0.5. They’re both bad for much larger disks, and in fact the cosine approximation is worse.

For the past two years I’ve been using Substack to send out emails announcing new blog posts. That has worked out well. Substack delivers email reliably, and that’s all I wanted. I’m not active on Substack other than using it to an email. I give a brief introduction to the latest two or three blog posts in each email, and sometimes I include additional ideas that occurred to me later.

If you’d like to get blog post announcements and a little extra commentary via email, sign up for my free Substack newsletter here.

If you’d like to learn about new posts sooner, you can subscribe via RSS or follow me on X or Mastodon.

(1 + cos(x))/2

gives a fair approximation to the Gaussian density

exp(−x²)

You can make the approximation much better by raising it to a power. The function

((1 + cos(x))/2)⁴

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 improving the cosine approximation to the Gaussian. Yesterday I came across one I hadn’t seen before, adding a sin(x) term to x.

(1 + cos(sin(x) + x))/2

This function matches the first few terms of the power series for exp(−x²) and has an error on the order of x⁶/240. You can’t see the difference between the two functions in a plot for −4 ≤ x ≤ 4.

***

There’s a tension between the previous two statements. If the error in on the order of x⁶/240 then we’d expect the error to be huge at x = 4. We have

4⁶/240 = 17.07

and yet

exp(−4²) − ((1 + cos(4 + sin(4)))/2) = −0.002579,

i.e. the error is between 3 and 4 orders of magnitude smaller than we might expect.

We have an alternating series, so the truncation error should be roughly equal to the first term after the truncation, right? No, the alternating series theorem doesn’t apply because the absolute values of the terms in the series are not decreasing yet for x = 4. The terms have to decrease eventually because the series has infinite radius of convergence, but they’re not decreasing at the 6th term; the terms will get much larger in absolute value before they get smaller.

The basic alternating series theorem gives only an upper bound on truncation error, but there are extensions that also give a lower bound. I wrote about these extensions a few weeks ago. But they don’t apply here because the terms have not started decreasing in absolute value.

Update: See further discussion in the post It’s not just Taylor series.

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 trick to reframe a problem that isn’t obviously an identity about polynomials. And with algebraic circuits, you can reformulate a wide range of computations as polynomial identities; this is widely used in zero-knowledge proofs.

The theorem alluded to at the top of the post is the Schwartz-Zippel lemma. It is formulated in terms of the probability of a non-zero polynomial P evaluating to zero. To prove that two polynomials Q₁ and Q₂ are equal, you can show that

P = Q₁(x) − Q₂(x) = 0.

Schwartz-Zippel lemma

Let F be a (typically large) finite field and let P be a non-zero polynomial in n variables

P(x₁, x₂, x₃, …, x_n)

of total degree d. If we choose the x‘s randomly from F then the probability that P evaluates to zero is no more than d/|F|. [1]

If the total degree d is small relative to the size of the field, then the probability of P evaluating to zero is small. As long as d is less than |F|, you can evaluate the polynomial k times to make

(d / |F|)^k

as small as you’d like. If d isn’t too large, and F is large, like the integers mod p = 2²⁵⁵ − 19 used in cryptography, one polynomial evaluation might be enough to give convincing evidence that the polynomial is zero.

[1] The Schwartz-Zippel lemma in its full generality applies to polynomials over an integral domain R with variables drawn from S, a finite subset of R. Here we’re setting R = S = F.

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 be computed in one pass through the data.

You’ll find two equivalent equations in statistics books: the one described above and another based on the sum of the data points and the sum of the data points squared.

While this equation is theoretically correct, it is numerically unstable. Code that directly implements this equation can return a negative value for a quantity that is theoretically positive. I’ve seen this happen with real data, causing a program to crash when taking the square root of the variance to get the standard deviation.

However, there is an algorithm that computes mean and variance in one pass that is accurate and numerically stable. This algorithm was developed by B. P. Welford in 1962. I discuss Welford’s algorithm and give code for implementing it here.

Online algorithms

Welford’s algorithm is known in computer science as an “online” algorithm. This term was coined well before the Internet. For example, see the paper [1] from 1965.

But of course now “online” means something else, and so the technical and colloquial uses of “online algorithm” have split. Technical literature uses the phrase to describe the kinds of algorithms in this post. Most people would take “online algorithm” to mean code that runs on a remote server. You may see “streaming algorithm” as a contemporary technical term, but I’d still search on “online algorithm” to find papers.

Computing higher moments online

Welford’s algorithm computes the first two moments, mean and variance, of a data set online. It is also possible to compute skewness and kurtosis online, as well as higher moments.

Online regression

Simple linear regression is closely related to calculating mean and variance, and it is possible to compute simple regression coefficients online. I have some old notes on this here.

This post was motivated by an email asking me about multiple regression. It is also possible to compute multiple regression coefficients online, but I haven’t done this. I found a couple references, [2] and [3], but I have not read them. There is a simple procedure for two predictor variables but I believe things get a little more complicated with three or more predictors, requiring a recursive least squares algorithm.

Related posts

The notion of online algorithms is closely related to the notion of a fold in functional programming. Here are several posts on computing things with folds.

• Insertion sort as a fold
• Kalman filters as folds
• Computing higher moments with a fold

[1] One-Tape, Off-Line Turing Machine Computations by F. C. Hennie. Information and Control. 8, 553-578 (1965). Available here. In this paper Hennie writes “In an on-line computation the input data are supplied to the machine, one symbol at a time, at a special input terminal. … In an off-line computation all of the input symbols are written on one of the machine’s tapes prior to the start of the computation.

[2] Arthur Albert and Robert W. Sittler, “A Method for Computing Least Squares Estimators that Keep Up with the Data,” Journal of the Society for Industrial and Applied Mathematics, Series A: Control, 3(3), 384–417, 1965. DOI: 10.1137/0303026.

[3] Petre Stoica and Per Ashgren. Exact initialization of the recursive least-squares algorithm. Int. J. Adapt. Control Signal Process. 2002; 16:219&ndashh;230.

Kullback-Leibler divergence is defined for two random variables X and Y by

K-L divergence is non-negative, and it’s zero if and only if X and Y 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 defining

The J above stands for Jeffreys, for Harold Jeffreys. J is called either the symmetrized K-L divergence or Jeffreys’ divergence. It’s still a divergence, not a distance.

A distance (metric) d has to have four properties:

1. d(x, x) = 0
2. d(x, y) > 0 if x ≠ y
3. d(x, y) = d(y, x)
4. d(x, z) ≤ d(x, y) + d(y, z)

K-L divergence satisfies the first two properties. Jeffreys’ divergence satisfies the first three, but not the last one, the triangle inequality.

To show that J doesn’t satisfy the triangle inequality, let X, Y, and Z be Bernoulli random variables with p equal to 0.1, 0.2, and 0.3 respectively. Then the following Python code shows that the divergence from X to Y, plus the divergence from Y to Z, is less than the divergence from X to Z. This would be like saying you could get from LA to NYC faster by having a layover in Denver rather than taking a direct flight.

from math import log

kl = lambda p, q: p*log(p/q) + (1-p)*log((1-p)/(1-q))
j  = lambda p, q: kl(p, q) + kl(q, p)

a = j(0.1, 0.2)
b = j(0.2, 0.3)
c = j(0.1, 0.3)
print(a + b, c)

This prints 0.135 and 0.270.

Jensen-Shannon distance

Jensen-Shannon distance turns K-L divergence into a metric as follows. First, define the random variable M to be the average of X and Y. Then average the K-L divergence from M to each of X and Y. This defines the Jensen-Shannon divergence. It’s still not a metric, but its square root is, which defines the Jensen-Shannon distance.

The following code gives an example of Jensen-Shannon distance satisfying the triangle inequality.

def d(p, q):
    m = 0.5*(p + q)
    jsd = 0.5*kl(p, m) + 0.5*kl(q, m) 
    return jsd**0.5

a = d(0.1, 0.2)
b = d(0.2, 0.3)
c = d(0.1, 0.3)
print(a + b, c)

This prints 0.1817 and 0.1801. Now a layover makes the trip longer.

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 parameters a and b to fit the curve?

We can rewrite the parametric expression for x as a sine with a phase shift (see notes here)

where

Also, we can rewrite the parametric expression for y as

Now the extreme values of x and y are easier to see. The maximum value of x is A and the minimum value is −A. The maximum value of y is A(cos(φ) + 1)/2 and the minimum value is A(cos(φ) − 1)/2.

We can simplify the cosine of an arctangent (see here) to find the height, i.e. the difference between the maximum and minimum y value, in terms of a and b.

Then the height is given by

The width is given by

and so

and

Now the Meta logo is drawn with a thick line, and the line width isn’t constant. It’s a little fuzzy what the height and width of the middle of the curve are, but I estimated h = 120 and w = 200 from one image. This leads to b = 20 and a = 97.98.

The Mathematica code

ParametricPlot[{a Cos[t] + 
   b Sin[t], -Sin[t] ( a Cos[t] + b Sin[t])}, {t, 0, 2 Pi}, 
 PlotStyle -> Thickness[0.05]]

produces the following image.

This is reminiscent of the Meta logo, but not a great match. I suspect the logo is not exactly a Besace curve. You could tinker with the a and b parameters and the aspect ratio to get a closer match. The logo may have been inspired by a Besace curve and then drawn by hand.

As mentioned in the previous post, the expected range is given by

where φ and Φ are the PDF and CDF of a standard normal. The integral can be calculated in closed form for n ≤ 5, but in general it requires numerical integration [1].

The following Python code can compute d_n.

from scipy.stats import norm
from scipy.integrate import quad
import numpy as np

def d(n):
    integrand = lambda x: x*norm.pdf(x)*norm.cdf(x)**(n-1)
    res, info = quad(integrand, -np.inf, np.inf)
    return 2*n*res

For large n we have the asymptotic approximation

which we could implement in Python by

def approx(n):
    return 2*norm.ppf((n - 0.375)/(n + 0.25))

For very large n the asymptotic expression may be more accurate than the integral due to numerical integration error.

Here are a few example values.

|-----+-------|
|   n |   d_n |
|-----+-------|
|   2 | 1.128 |
|   3 | 1.693 |
|   5 | 2.326 |
|  10 | 3.078 |
|  12 | 3.258 |
|  23 | 3.858 |
|  50 | 4.498 |
| 100 | 5.015 |
|-----+-------|

[1] Order Statistics by H. A. David. John Wiley & Sons. 1970.

Someone dismissed this communication difficulty by pointing out that the expected difference in IQ between two individuals is around 17, suggesting that most communication is between people who differ by more than one standard deviation in IQ. But this calculation assumes people are chosen at random, which they usually are not. People tend to live around and work around others of similar intelligence.

However, a jury is a random sample. It’s not a perfect random sample. For one thing, it starts with a random sample of people who are registered to vote, or who have a driver’s license, not all individuals. Furthermore, the pool of potential jurors is reduced to a jury through the process of voir dire, which is not random.

For this post I will make the simplifying assumption that a jury is a random sample from a population with normally distributed IQ with standard deviation σ = 15. The mean doesn’t matter here, but you could assume it’s 100 if you’d like.

By symmetry, the expected range of n samples from a normal random variable is twice the maximum. For n = 12 the range is about 3.26σ, which corresponds to nearly 50 IQ points.

This suggests there’s usually a big spread of IQ on a jury. Even if IQ doesn’t measure intelligence, it measures something, and that something varies a lot over 12 people chosen at random [1].

Related posts

• Explaining probability to a jury
• Expert testimony
• Seven questions a statistician can answer for a lawyer

[1] In case you’re interested in the technical details, the expected range of n samples from a standard normal random variable is given by

where φ and Φ are the PDF and CDF of a standard normal. Multiply this by σ to get the range of a normal random variable with standard deviation σ. As for how to calculate d_n, see the next post.

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. Apparently they were more fashionable a few decades ago than they are now. I suppose the notation falls between two stools, too concrete for some tastes and not concrete enough for others. The former folks would prefer something like H and the latter would prefer the sum above.

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 calculator. A couple of the functions might be unfamiliar, such as sgn and atan2, but there are no functions like the gamma function that are defined in terms of integrals.

One application of Baker’s equations would be to bootstrap a math library that doesn’t support complex numbers into one that does. But the equations could be useful in pure math when you’d like to have a convenient expression for the real or imaginary part of a function.

The real and imaginary parts of a complex analytic function are harmonic functions. So the functions on the right hand side of Baker’s equations satisfy Laplace’s equation:

u_xx + u_yy = 0

and

v_xx + v_yy = 0.

Furthermore, the functions u and v form harmonic conjugate pairs, meaning each is the Hilbert transform of the other.

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 to numerically verify the equations and my typesetting of them. This caught a few typos on my part.

The test code evaluated the equations at points from each quadrant. All matched NumPy, implying that Baker and NumPy use the same branch cuts on inverse functions.

***

This post is part of a thread that has gone on for a few days. Maybe it’s the last post in the thread; we’ll see.

It all started with a post on Markov’s equation

x² + y² + z² = 3 xyz

and an approximation to the equation that has a closed-form solution. That led to the identity

cosh( arccosh(a) + arccosh(b) ) = ab + √(a² − 1) √(b² − 1).

The approximation to Markov’s equation only needed the identity to be valid for real a and b greater than 1. But when I looked closer at the identity I found several complications with branch cuts. The identity doesn’t hold everywhere using the principle branch of the square root function. But if you define √(z² − 1) to have a branch cut along [−1, 1] then the equation holds everywhere in the complex plane. And that led to my writing up some notes on how to define all the elementary inverse functions in terms of log.

Someone reading these posts suggested I look at a paper that mentioned “couth” and “uncouth” function pairs, which led to this post and its warm up.

I find all this interesting because it’s an advanced perspective on a questions that are latent in an intro calculus class. What exactly do functions like arccos mean and why where they defined as they were? These are fairly deep and interesting questions that are swept under the rug, and swept there for good reason. A calculus class has to cover an enormous amount of material and there’s no time to dwell on fine points. Some of my favorite posts look back leisurely on things that go by in a blur when you’re a student.

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 counterpart “fooh” by applying one of the functions to iz then multiplying by a constant c that depends on foo: c = i for sin and tan, c = 1 for cos and sec, and c = −i for csc and cot.

In symbols,

c foo(z) = fooh(iz).

Let’s suppose foo and fooh are invertible, ignoring any complications, and solve foo(z) = w for z. We get

i foo⁻¹(w) = fooh⁻¹(cw)

or using “arc” nomenclature for inverse functions

i arcfoo(w) = arcfooh(cw).

When the naive calculation above holds, except possibly at a finite number of points, we say the pair (foo, fooh) is couth. Otherwise we say the pair is uncouth. These term were coined by Robert Corless and his coauthors in their paper [1].

Whether the pair (foo, fooh) is couth depends not only on foo and fooh, but also on the details of how arcfoo and arcfooh are defined.

In Python’s NumPy library, the pairs (sin, sinh) and (tan, tanh) are couth, but the pair (cos, cosh) is uncouth.

Numpy doesn’t define the reciprocal functions sec, sech, csc, csch, cot, and coth. I used to find that annoying, but I’m beginning to think that was wise. These functions cause problems. For example, there may be two reasonable ways to define these functions, one of which forms a couth pair and one of which forms an uncouth pair.

For example, how should you define cot and coth? There would be no disagreement over the definition

cot = lambda x: 1/tan(x)

but there are at least two definitions of inverse coth that you’ll find in practice:

arccot = lambda z: 0.5*pi - arctan(z)
arccot = lambda z: arctan(1/z).

Both definitions have their advantages, but the former is uncouth while the latter is couth. You can verify that both definitions are the same at z = 1 but not at z = −1.

With the following definitions, the pairs (cos, cosh) and (sec, sech) are uncouth and the rest are couth.

from numpy import *

csc     = lambda x: 1/sin(x)
sec     = lambda x: 1/cos(x)
cot     = lambda x: 1/tan(x)
csch    = lambda x: 1/sinh(x)
sech    = lambda x: 1/cosh(x)
coth    = lambda x: 1/tanh(x)

arccot  = lambda z: arctan(1/z)
arcsec  = lambda z: arccos(1/z)
arccsc  = lambda z: arcsin(1/z)
arccoth = lambda z: arctanh(1/z)
arcsech = lambda z: arccosh(1/z)
arccsch = lambda z: arcsinh(1/z)

[1] “According to Abramowitz and Stegun” or arccoth needn’t be uncouth. Robert M. Corless et al. ACM SIGSAM Bulletin, Volume 34, Issue 2, pp 58 – 65 https://doi.org/10.1145/362001.362023

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 “foo” with its hyperbolic counterpart “fooh” by applying one of the functions to iz then multiplying by a constant c that depends on foo: c = i for sin and tan, c = 1 for cos and sec, and c = −i for csc and cot.

Note that if the constant for foo is c, the constant for 1/foo is 1/c. So, for example, the constant for tan is i and the constant for cot is 1/i = −i.

We have four groups of equations, depending on whether the left side of the equation is foo(iz), fooh(iz), foo(z), or fooh(z).

This post was written as a warm-up for the next post on couth and uncouth function pairs.

foo(iz)

fooh(iz)

foo(z)

fooh(z)

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 ℂ. We have to choose to make a “cut” somewhere, and the conventional choice is to make a cut along the negative real axis.

Using the principal branch

The “principal branch” of the square root function is defined to be the unique function that analytically extends √x from the positive reals to ℂ \ (−∞, 0].

Assume for now that by √x we mean the principal branch of the square root function. Now what does √(z² − 1) mean? It could mean just what we said at the top of the post: we square z, subtract 1, and apply the (principal branch of the) square root function. If we do that, we must exclude those values of z such that (z² − 1) is negative. This means we have to cut out the imaginary axis and the interval [−1, 1].

This is what Mathematica will do when asked to evaluate Sqrt[z^2 - 1]. The command

ComplexPlot[Sqrt[z^2 - 1], {z, -2 - 2 I, 2 + 2 I}]

makes the branch cuts clear by abrupt changes in color.

A different approach

Now let’s take a different approach. Consider the function √(z² − 1) as a whole. Do not think of it procedurally as above, first squaring z etc. Instead, think of a it as a black box that takes in z and returns a complex number whose square is z² − 1.

This function has an obvious definition for z > 1. And we can extend this function, via analytic continuation, to more of the complex plane. We can do this directly, not by extending the square root function. But as before, we cannot extend the function analytically to all of ℂ. We have to cut something out. A common choice is to cut out [−1, 1]. This eliminates the need for a branch cut along the imaginary axis.

The function

extends √(z² − 1) the way we want [1].

The Mathematica code

ComplexPlot[Exp[(1/2) (Log[z - 1 ] + Log[z + 1])], {z, -2 - 2 I, 2 + 2 I}]

shows that our function is now continuous across the imaginary axis, though there’s still a discontinuity as you cross [−1, 1].

We used this analytic extension of √(z² − 1) in the previous post to eliminate branch cuts in an identity.

Related posts

• Trig composition table
• sinh( arccosh(x) )
• Duplicating a hand-drawn plot

[1] The principal branch of the logarithm has a cut along the negative real axis. Why does our square root function, defined using log, not have a branch cut along the negative axis?

First of all, the log function, and Mathematica’s implementation of it Log[], isn’t undefined on (−∞, 1), it just isn’t continuous there. The function still has a value. By convention, the value is taken to be the limit of log(z) approaching z from above, i.e. from the 2nd quadrant.

Second, the value of (log(z – 1) + log(z + 1))/2 differs by a factor of 2πi when approaching a value z < −1 from above versus from below. This factor goes away when taking the exponential. So our function f(z) is analytic across (−∞, 1) even though the log functions in its definition are not.

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 √(ab) when the arguments may be negative.

The square root function and the arccosh function are not naturally single-valued functions. They require branch cuts to force them to be single-valued, and the two functions require different branch cuts. I go into this in some detail here.

There is a way to reformulate our identity so that it holds everywhere. If we replace

with

which is equivalent for z > 1, the corresponding identity holds everywhere.

We can verify this with the following Mathematica code.

f[z_] := Exp[(1/2) (Log[z - 1 ] + Log[z + 1])]
FullSimplify[Cosh[ArcCosh[x] + ArcCosh[y]] - x y - f[x] f[y]]

This returns 0.

By contrast, the code

FullSimplify[
 Cosh[ArcCosh[x] + ArcCosh[y]] - x y - Sqrt[x^2 - 1] Sqrt[y^2 - 1]]

simply returns its input with no simplification, unless we add restrictions on x and y. The code

FullSimplify[
 Cosh[ArcCosh[x] + ArcCosh[y]] - x y - Sqrt[x^2 - 1] Sqrt[y^2 - 1], 
 Assumptions -> {x > -1 && y > -1}]

does return 0.

Related posts

• Inverse cosine
• Branch cuts and Common Lisp
• Duplicating a plot from A & S

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 look at an example that shows solutions to his second equation approximate solutions to Markov’s equation.

The following code verifies that (5, 13, 194) is a solution to Markov’s equation.

x, y, z = 5, 13, 194
assert(x**2 + y**2 + z**2 == 3*x*y*z)

With the same x and y above, let’s show that the z in Zagier’s second equation is close to the z above.

from math import cosh, acosh

f = lambda t: acosh(3*t/2)
g = lambda t: cosh(t)*2/3
z = g(f(x) + f(y))
print(z)

This gives z = 194.0023, which is close to the value of z in the Markov triple above.

Applying Osborn’s rule

Now suppose

f(x) + f(y) = f(z)

which expands to

arccosh(3x/2) + arccosh(3y/2)  = arccosh(3z/2).

It seems sensible to apply cosh to both sides. Is there some identity for cosh of a sum? Maybe you recall the equation for cosine of a sum:

cos(a + b) = cos(a) cos(b) − sin(a) sin(b).

Then Osborn’s rule says the corresponding hyperbolic identity is

cosh(a + b) = cosh(a) cosh(b) − sinh(a) sinh(b).

Osborn’s rule also says that the analog of the familiar identity

sin²(a) + cos²(b) = 1

is

sinh²(a) = cosh²(b) − 1.

From these two hyperbolic identities we can show that [1]

cosh( arccosh(a) + arccosh(b) ) = ab + √(a² − 1) √(b² − 1).

Slug it out

The identity derived above is the tool we need to reduce our task to routine algebra.

If

arccosh(3x/2) + arccosh(3y/2)  = arccosh(3z/2)

then

(3x/2)  (3y/2)  + √((3x/2)² − 1) √((3y/2)² − 1) = 3z / 2

which simplifies to Zagier’s equation

x² + y² + z² = 3xyz + 4/9.

Related posts

• Primes and weeds
• sinh( arccosh(x) )
• Three paths converge

[1] The equation holds at least for x > 1 and y > 1, which is enough for this post. More general arguments run into complications due to branch cuts.

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 >> 19) & MASK

    a, b, c, d = (t ^ s) & MASK, a, b, c

    return a

Recovering state

Recovering the internal state of the generator is simple: it’s the four latest outputs in reverse order. This is illustrated by the following chart.

This means that once we’ve seen four outputs, we can predict the rest of the outputs. The following code demonstrates this.

Let’s generate five random values.

out = [xorshift128() for _ in range(5)]

Running

print(hex(out[4]))

shows the output 0xc3f4795d.

If we reset the state of the generator using the first four outputs

d, c, b, a, _ = out
print(hex(xorshift128()))

we get the same result.

Good stats, bad security

Mersenne Twister and lehmer64 have good statistical properties, despite being predictable. The xorshift128 generator is even easier to predict, but it also has good statistical properties. These generators would be fine for many applications, such as Monte Carlo simulation, but disastrous for use in cryptography.

The post on lehmer64 mentioned at the end that the internal state of PCG64 can also be recovered from its output. However, doing so requires far more sophisticated math and thousands of hours of compute time. Still, it’s not adequate for cryptography. For that you’d need a random number generator designed to be secure, such as ChaCha.

So why not just use a cryptographically secure random number generator (CSPRNG) for everything? You could, but the other generators mentioned in this post use less memory and are much faster. PCG64 occupies an interesting middle ground: simple and fast, but not easily reversible.

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 previous post introduces a function print_u128.

The code

__uint128_t x = 0x00112233445566778899aabbccddeeff;

Produces the following error message.

error: integer literal is too large to be represented in any integer type

The problem isn’t initializing a 128-bit number to a 128-bit value; the problem is that the compiler cannot parse the literal expression

0x00112233445566778899aabbccddeeff

One solution to the problem is to introduce the macro

#define U128(hi, lo) (((__uint128_t)(hi) << 64) | (lo))

and use it to initialize the variable.

__uint128_t x = U128(0x0011223344556677, 0x8899aabbccddeeff);

You can verify that x has the intended state by calling print_u128 from the previous post.

void print_u128(__uint128_t n)
{
    printf("0x%016lx%016lx\n",
           (uint64_t)(n >> 64),      // upper 64 bits
           (uint64_t)n);             // lower 64 bits
}

Then

print_u128(x);

prints

0x00112233445566778899aabbccddeeff

Update. The code for print_u128 above compiles cleanly with gcc but clang gives the following warning.

warning: format specifies type 'unsigned long' but the argument has type 'uint64_t' (aka 'unsigned long long') [-Wformat]

You can suppress the warning by including the inttypes header and modifying the print_u128 function.

Here’s the final code. It compiles cleanly under gcc and clang.

#include <stdio.h>
#include <stdint.h>
#include <inttypes.h>
#define U128(hi, lo) (((__uint128_t)(hi) << 64) | (lo))

void print_u128(__uint128_t n)
{
    printf("0x%016" PRIx64 "%016" PRIx64 "\n",
           (uint64_t)(n >> 64),
           (uint64_t)n);
}

int main(void)
{
    __uint128_t x = U128(0x0011223344556677, 0x8899aabbccddeeff);
    print_u128(x);
    return 0;
}

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 be implemented in C by

__uint128_t g_lehmer64_state;

uint64_t lehmer64() {
    g_lehmer64_state *= 0xda942042e4dd58b5ULL;  
  return g_lehmer64_state >> 64;
}

The analogous code in Python would have to simulate the overflow behavior of a 128-bit integer by reducing the state mod 2¹²⁸ after the multiplication.

Reverse engineering lehmer64 is easier than reverse engineering the Mersenne Twister because only three outputs are needed. However, the theory behind the exploit is more sophisticated. See [1].

The following code sets the state to an arbitrary initial seed value and generates three values.

#include <stdio.h>
#include <stdint.h>

int main(void)
{
    g_lehmer64_state = 0x4789499d78770934; // seed
    for (int i = 0; i < 3; i++) {
        printf("0x%016lx\n", lehmer64());
    }

    return 0;
}

The code prints the following.

0x3d144d12822bcc2e
0x85a67226191a568d
0x53e803dffc88e8f8

Exploiting lehmer64

Here is Python code for recovering the state of the lehmer64 generator given in [1].

def reconstruct(X):
    a = 0xda942042e4dd58b5
    r = round(2.64929081169728e-7 * X[0] + 3.51729342107376e-7 * X[1] + 3.89110109147656e-8 * X[2])
    s = round(3.12752538137199e-7 * X[0] - 1.00664345453760e-7 * X[1] - 2.16685184476959e-7 * X[2])
    t = round(3.54263598631140e-8 * X[0] - 2.05535734808162e-7 * X[1] + 2.73269247090513e-7 * X[2])
    u = r * 1556524 + s * 2249380 + t * 1561981
    v = r * 8429177212358078682 + s * 4111469003616164778 + t * 3562247178301810180
    state = (a*u + v) % (2**128)
    return state

Let’s call reconstruct with the output of our C code.

X = [0x3d144d12822bcc2e, 0x85a67226191a568d, 0x53e803dffc88e8f8]
print( hex( reconstruct(X) ) )

This prints

0x3d144d12822bcc2e1b81101c593761c4

Now for the confusing part: at what point is the number above the state of the generator? Is it the state before or after generating the three values? Neither! It is the state after generating the first value.

We can verify this by modifying the C code as follows and rerunning it.

void print_u128(__uint128_t n)
{
    printf("0x%016lx%016lx\n",
           (uint64_t)(n >> 64),      // upper 64 bits
           (uint64_t)n);             // lower 64 bits
}

int main(void)
{
    g_lehmer64_state = 0x4789499d78770934; // seed
    for (int i = 0; i < 3; i++) {
        printf("0x%016lx\n", lehmer64());
        printf("state: ");
        print_u128(g_lehmer64_state);
    }
 
    return 0;
}

The main goal of [1] is to recover the state of the PCG generator, not the lehmer64 generator. The latter was a side quest. Recovering the state of PCG64 is much harder; the authors estimate it takes about 20,000 CPU-hours. The paper shows that a technique used as part of pursuing their main goal can quickly recover the lehmer64 state.

Related posts

• Recovering the Mersenne Twister internal state with linear algebra
• Testing the PCG random number generator
• Random number generator testing

[1] Charles Bouillaguet, Florette Martinez, and Julia Sauvage. Practical seed-recovery for the PCG Pseudo-Random Number Generator. IACR Transactions on Symmetric Cryptology. ISSN 2519-173X, Vol. 2020, No. 3, pp. 175–196.

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 after him. Several functions are known as Euler’s function, the most common being his totient function in number theory. The Euler function we’re interested in here is

for −1 < x < 1. Usually the argument of φ is denoted “q” but that would be confusing in our context because our q, the number of elements in a field, is the reciprocal of Euler’s q, i.e. x = 1/q.

Euler’s identity [2] (in this context, not to be confused with other Euler identities!) says

This function is easy to calculate because the series converges very quickly. From the alternating series theorem we have

When q = 2 and so x = 1/2, N = 6 is enough to compute φ(x) with an error less than 2⁻⁷⁰, beyond the precision of a floating point number. When q is larger, even fewer terms are needed.

To illustrate this, we have the following Python script.

def phi(x, N):
    s = 0
    for n in range(-N, N+1):
        s += (-1)**n * x**((3*n**2 - n)/2)
    return s

print(phi(0.5, 6))

Every digit in the output is correct.

Related posts

• q-analogs of rising powers
• Another pentagonal number theorem

[1] I didn’t say that explicitly, but I pointed out that p(2, 8) was close to p(2, ∞).

[2] This identity is also known as the pentagonal number theorem because of its connection to pentagonal numbers.

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 post about representing bit manipulations, including shifts, as matrix operators. The matrix corresponding to shifting right by k bits has 1s on the kth 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.

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

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.

If you’d like to prove that the pseudoinverse of a matrix that shifts right by k places is a matrix that shifts left by k places, you don’t have to compute the pseudo inverse per se: you can verify your guess. This post 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.

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 multiply the matrices, are you doing ordinary integer arithmetic or are you working mod 2?

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.

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.

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

Probability of invertibility in GF(2)

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 × n matrix over GF(q), a finite field with q elements, is invertible? This is

When q = 2 and n = 8 this probability is 0.289919. The probability is roughly the same for all larger values of n, converging to approximately 0.288788 as n → ∞.

Probability of invertibility in ℝ

What is the probability that an 8 × 8 matrix with random 0 and 1 entries is invertible as a real matrix? We can estimate this by simulation.

import numpy as np

def simulate_prob_invertible_real(n, numreps=1000):
    s = 0
    for _ in range(numreps):
        M = np.random.randint(0, 2, size=(n, n))
        det = np.linalg.det(M)
        if abs(det) > 1e-9:
            s += 1
    return s/numreps

When n = 8, I got 0.5477 when running the code with 10,000 reps.

When n = 32, I got a probability of 1. Obviously it is possible for a 32 × 32 binary matrix to be singular, but it’s very unlikely: it didn’t happen in 10,000 random draws.

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 fields GF(2), the field with just two elements, 0 and 1.

In GF(2), addition corresponds to XOR. We will denote this by ⊕ to remind us that although it’s addition, it’s not the usual addition, i.e. 1 ⊕ 1 = 0. Similarly, multiplication corresponds to AND. We’ll work with 8-bit numbers to make the visuals easier to see.

Shifting a number left one bit corresponds to multiplication by a matrix with 1’s below the diagonal main. Shifting left by k bits is the same as shifting left by 1 bit k times, so the the matrix representation for x << k is the kth power of the matrix representation of shifting left once. This matrix has 1s on the kth diagonal below the main diagonal. Below is the matrix for shifting left two bits, x << k.

Right shifts are the mirror image of left shifts. Here’s the matrix for shifting right two bits, x >> k.

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 to x >> 2 is not invertible, the function that takes x to x ⊕ (x >> 2) is invertible. This operation corresponds to the matrix below.

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.

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.

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

y ^= (y <<  7) & 0x9d2c5680 

says to shift the bits of y left 7 places, then zero out the elements corresponding to 0s in the mask, then XOR the result with y. 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.

[1] Until you get to eigenvalues. Then it matters whether the field is algebraically complete, which no finite field is.

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 the original version that had internal state consists of vector x of length 640 filled with 32-bit numbers. The ideas in this post would apply equally to all MT versions.

The first call to MT returns a “tempered” version of x[0]. The next call returned a tempered version of x[1], 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.)

Tempering

Here is Python code for performing the tempering step.

def temper(y):
    y ^= (y >> 11) 
    y ^= (y <<  7) & 0x9d2c5680 
    y ^= (y << 15) & 0xefc60000 
    y ^= (y >> 18)  
    return y

Each step is reversible, and so the temper function is reversible.

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.

Linear algebra

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.

Each step corresponds to multiplying a vector of 32 bits on the left by a 32 × 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.

We will find the matrix M 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⁻¹, though it would be more efficient to undo tempering by bit twiddling.

One way to recover a matrix is to multiplying by unit vectors e_i where the ith component of e_i is 1 and the rest of the components are zero. Then

M e_i

is the ith column of M.

So we can find the nth column of M by tempering 1 << n = 2ⁿ.

M = np.zeros((32, 32), dtype=int)
for i in range(32):
    t = temper(1 << (31-i))
    s = f"{t:032b}"
    for j in range(32):
        M[j, i] = int(s[j])

Let’s generate a random number and compute its tempered form two ways: directly and matrix multiplication.

x = random.getrandbits(32)
y = temper(x)
print(f"{y:032b}")

vx = np.array([int(b) for b in f"{x:032b}"]) # vector form of x
vy = M @ vx % 2 # vector form of y
print("".join(str(b) for b in vy))

Both produce the same bits:

10100101100101101100110101000110
10100101100101101100110101000110

We can find the matrix representation of the untemper function by inverting the matrix M. However, we need to invert it over the field GF(2), not over the integers or reals.

import galois
GF2 = galois.GF(2)
Minv = np.linalg.inv(GF2(M))

Here are visualizations of M and its inverse using a black square for a 1 and a white square for a 0.

M:

M⁻¹:

The next post will back up and look at the linear algebra of the components that comprise the Mersenne Twister.