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 kneejerk 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 temporing 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.

Even when f has a fairly simple expression, the expression for κ can be complicated.

If we define

then the level set of f(x, y) = c is an equilateral triangle when c = −4. The level sets are smoothed triangles for −4 < c < 0.

The curvature of these level sets at any point is given by

Simplification

But there is one instance in which curvature is easy to calculate. For the graph of a function g(x) = y, the curvature is approximately the absolute value of the second derivative of g, provided the first derivative is small.

At a local maximum or local minimum of g(x) the approximation is exact because the first derivative is zero.

Max and min curvature of smoothed triangles

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 x = 0 for convenience. There we find the curvature is simply

When x = 0, we have

and so the maximum and minimum curvature are the two roots of the cubic equation for y that lie in the interval [−1, 2]. (There’s another root greater than 2.)

For example, when c = −3, 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.

When c = −1 the minimum and maximum curvature are 2.80747 and 9.91622 respectively,

Since c = −4 corresponds to the triangle, the minimum curvature is 0 and the maximum is infinite. As c increases, the minimum and maximum curvature come together because the level set is becoming more round.

The previous post introduced three functions L_i(x, y) such that the level set of each function

forms a side of a triangle. Then it introduced a soft penalty for each L being away 1, and the level sets of that penalty function formed the rounded triangles we were looking for.

Another approach would be to change the L‘s slightly so that the sides are the levels sets L_i(x, y) = 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

then the set of points

corresponds to the three lines when c = 0 and the level sets for small c > 0 are nearly triangles. The level sets will be smooth if the gradient is non-zero, i.e. c is not zero.

This approach is not unique to triangles. You could create smooth approximations any polygon by multiplying linear functions that define the sides. Or you could do something similar with curved arcs.

If we define our L’s by

then our curves will be the level sets of

A few level sets are shown below. The level set for c = 0 is the straight lines.

Note the level sets are not connected. If you’re interested in approximate triangles, you want the components that are inside the triangle.

A squircle is not exactly a square with rounded corners. The sides are continuously curved, but curved most at the corners. See, for example, this post.

Suppose the sides of our triangle are given by L₁(x, y) = 1 for i = 1, 2, 3. For example,

We design a function f(x, y) as a soft penalty for points not being on one of the sides and look at the set of points f(x, y) = 1.

You might recognize this as a Lebesgue norm, analogous to the one used to define a squircle.

The larger p is, the heavier the penalty for being far from a side and the closer the level set f(x, y) = 1 comes to being a triangle.

Unified keys

One thing I do is remap keys so that the same key does the same thing everywhere, to the extent that’s practical. This requires remapping keys. In particular, I want the key functionality, not the key name, to be the same across operating systems. For example, the Command key on a Mac does what the Control 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.

Unified config files

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 files, one for each OS. However, both .bashrc files source a common file .bash_aliases for aliases. You can set that up by putting the following code in your .bashrc file.

if [ -f ~/.bash_aliases ]; then
    . ~/.bash_aliases
fi

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.

HOSTNAME_SHORT=$(hostname -s)
if [ "$HOSTNAME_SHORT" = "Castor" ]; then
...
elif [ "$HOSTNAME_SHORT" = "Pollux" ]; then
...

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.

I have one init.el file for Emacs. It contains some branching logic for OS:

(cond
    ((string-equal system-type "windows-nt") ; Microsoft Windows
     ...)
    ((string-equal system-type "gnu/linux") ; linux
     ...)
    ((string-equal system-type "darwin") ; Mac
     ...)
)

You can also branch on hostname.

(if (string-equal (system-name) "Castor")
...)
(if (string-equal (system-name) "Pollux")
...)

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.

Later I ran across the following post from Qiaochu Yuan. It felt as if he had overheard my conversation and summarized it in a tweet:

category theory is just some straightforwardly useful stuff for some purposes in some fields! you can elegantly simplify and streamline some proofs. then there is the mythology of category theory, which is some other thing entirely, mostly wishful thinking and projection afaict

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.

Changing a hyphen to an en-dash increases your PDF file size by ~10 bytes.

My first thought was that it had something to do with hyphen being an ASCII character and an en-dash not. Changing a hyphen to an en-dash would make a UTF-8 encoded text file a couple bytes longer. (See why here.) Maybe adding one non-ASCII character could cause the file to include a glyph it didn’t before.

I did a couple experiments. I made a minimal LaTeX file with only the text

    See pages 9-10.

and another with

    See pages 9--10.

(In LaTeX, a hyphen compiles to a hyphen and two hyphens compile to an en-dash.)

I compiled both files using pdflatex. 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.

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.

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 the size of the file by 72 bytes.

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.

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.

Related posts

• How text case can change QR code size
• Handyman lorem ipsum
• PDF metadata

(log x)² + (log y)² = 1.

Of course the plot of

x² + y² = 1

is a circle, but I never thought what taking logs would do to the shape.

Here’s what the contours look like setting the right hand side equal to 1, 2, 3, …, 10.

ContourPlot[Log[x]^2 + Log[y]^2, {x, 0, 10}, {y, 0, 10}, 
    Contours -> Range[10]]

The dark blue contour near the origin reminded me of a guitar pick, so I decided to take a stab at creating an equation for the shape of a guitar pick.

I wanted to rotate the image so the axis of symmetry for the pick is vertical, so I replaced x and y with x + y and x − y.

The aspect ratio was too wide, so I experimented with

log(y + kx)² + log(y − kx)² = r²

where increasing k increases the height-to-width ratio. After a little experimentation I settled on k = 1.5 and r = 1.

This has an aspect ratio of roughly 5:4, which is about what I measured from a photo of a guitar pick.

Updating: refining the fit

After posting this article on X, Paul Graham replied with a photo of a Fender guitar pick with the image above overlaid. The fit was fairly good, but the aspect ratio wasn’t quite right.

So then I did a little research. The shape referred to in this post is known as the “351,” but even for the 351 shape the aspect ratio varies slightly between picks.

Setting k = 1.6 gives a better fit to Paul Graham’s pick.

The blue line represents my fit using k = 1.5 and the red line represents my fit using k = 1.6.

The first-order approximation has error O(x⁴), as shown in the earlier post. Adding the second-order term makes the error O(x⁶), and adding the third-order term makes it O(x⁸).

I’ve written a few times about cosine approximations to the normal probability density. For example, see this post. We could use the same idea as the series above to approximate the normal density with a series of powers of cosine. This gives us

and as before, the first, second, and third order truncated series have error O(x⁴), O(x⁶), and O(x⁸).

The general theory behind what’s going on here is an extension of Bürmann’s theorem. The original version of the theorem relies on a series inversion theorem that in turn relies on the approximating function, in our case cos(x) − 1, not having zero derivative at the center of the series. But there is a more general form of Bürmann’s theorem based on a more general form of series inversion. We will always need a more general version of the theorem when working with even functions because even functions have zero derivative at zero.

Here’s another example, this time using the Bessel function J₁, an odd function, which does use the original version of Bürmann’s theorem to approximate J₁ by powers of sine.

In this case truncating the series after sin^k(x) gives an error O(x^{k + 2}).

You can find more on Bürmann’s theorem in Whittker and Watson.

The function is simply r(x) = max(0, x).

Pointwise derivative

The pointwise derivative would be 0 for x < 0, 1 for x > 0, and undefined at x = 0. So except at 0, the pointwise derivative of the ramp function is the Heaviside function.

In a real analysis course, you’d simply say r′(x) =H(x) because functions are only defined up to equivalent modulo sets of measure zero, i.e. the definition at x = 0 doesn’t matter.

Distributional derivative

In distribution theory you’d identify the function r(x) with the distribution whose action on a test function φ is

Then the derivative of r would be the distribution r′ satisfying

for all smooth functions φ with compact support. You can prove using integration by parts that the above equals the integral of φ from 0 to ∞, which is the same as the action of H(x) on φ.

In this case the distributional derivative of r is the same as the pointwise derivative of r interpreted as a distribution. This does not happen in general [1]. For example, the pointwise derivative of H is zero but the distributional derivative of H is δ, the Dirac delta distribution.

For more on distributional derivatives, see How to differentiate a non-differentiable function.

Subgradient

The subgradient of a function f at a point x, written ∂f(x), is the set of slopes of tangent lines to the graph of f at x. If f is differentiable at x, then there is only one slope, namely f′(x), and we typically say the subgradient of f at x is simply f′(x) when strictly speaking we should say it is the one-element set {f′(x)}.

A line tangent to the graph of the ReLU function at a negative value of x has slope 0, and a tangent line at a positive x has slope 1. But because there’s a sharp corner at x = 0, a tangent at this point could have any slope between 0 and 1.

My dissertation was full of subgradients of convex functions. This made me uneasy because subgradients are not real-valued functions; they’re set-valued functions. Most of the time you can blithely ignore this distinction, but there’s always a nagging suspicion that it’s going to bite you unexpectedly.

[1] When is the pointwise derivative of f as a function equal to the derivative of f as a distribution? It’s not enough for f to be continuous, but it is sufficient for f to be absolutely continuous.

I wanted to see if I could get anything interesting by turning the trick in the previous post into a technique. The trick created a high-order approximation by subtracting a multiple one even function from another. Even functions only have even-order terms, and by using the right multiple you can cancel out the second-order term as well.

For an example, I’d like to approximate the Bessel function J₀(x) by the better known cosine function. Both are even functions.

J₀(x) = 1 − x²/4 + x⁴/64 + …
cos(x) = 1 − x²/2 + x⁴/24 + …

and so

2 J₀(x) − cos(x) = 1 − x⁴/96 + …

which means

J₀(x) ≈ (1 + cos(x))/2

is an excellent approximation for small x.

Let’s try this for a couple examples.

J₀(0.2) = 0.990025 and (1 + cos(0.2))/2 = 0.990033.

J₀(0.05) = 0.99937510 and (1 + cos(0.05))/2 = 0.99937513.

Here θ is the central angle of the arc. Then the length of the arc, rθ, is approximately

a = rθ ≈ 12 b²/(c + 4b).

If the arc is moderately small, the approximation is very accurate.

This approximation is simple, accurate, and not obvious, much like the one in this post

Derivation

Let φ = θ/4. Then the angle between the chords b and c is φ. This follows from the inscribed angle theorem, illustrated below.

There are two right triangles in the diagram above that have an angle φ: a smaller triangle with hypotenuse b and a larger triangle with hypotenuse 2r. From the smaller triangle we learn

cos(φ) = c / 2b

and from the larger triangle we learn

sin(φ) = b / 2r.

Now expand in power series.

c / 2b = cos(φ) = 1 − φ²/2! + φ⁴/4! − …
2b / a = sin(φ) / φ = 1 − φ²/3! + φ⁴/5! − …

If we multiply 2b / a by 3 and subtract c / 2b then the φ² terms cancel out and we get

6b / a − c / 2b = 2 − φ⁴/60 + …

and so

6b / a − c / 2b ≈ 2

to a very high degree of accuracy when φ is small. The approximation follows by solving for a.

Example

Let θ = π/3 and so φ = 0.26…, not a particularly small value of φ, but small enough for the approximation to work well.

Set r = 1 so a = θ. Then

b = 2 sin(π/12) = 0.51764

and

c = 2b cos(π/12) = 1.

Now in application, we know b and c, not θ, and so pretend we measured b = 0.51764 and c = 1. Then we would approximate a by

12b²/(c + 4b) = 1.04718

while the exact value is 1.04720. Unless you can measure lengths to more than four significant figures, the approximation may has well be exact because approximation error would be less than measurement error.

[1] J. M. Bruce. Approximation to a Circular Arc. The American Mathematical Monthly. Vol. 49, No. 3 (March 1942), pp. 184–185

If the initial displacement is small enough, you can simply replace sin θ with θ. If the initial displacement is larger, you can improve the accuracy quite a bit by solving the linearized equation and then adjusting the period.

You can also find an exact solution, but not in terms of elementary functions; you have to use Jacobi elliptic functions. These are functions somewhat analogous to trig functions, though it’s not helpful to try to pin down the analogies. For example, the Jacobi function sn is like the sine function in some ways but very different in others, depending on the range of arguments.

We start with the differential equation

θ″(t) + c² sin( θ(t) ) = 0

where c² = g/L, i.e. the gravitational constant divided by pendulum length, and initial conditions θ(0) = θ₀ and θ′(0) = 0. We assume −π < θ₀ < π.

Then the solution is

θ(t) = 2 arcsin( a cd(ct | m ) )

where a = sin(θ₀/2), m = a², and cd is one of the 12 Jacobi elliptic functions. Note that cd, like all the Jacobi functions, has an argument and a parameter. In the equation above the argument is ct and the parameter is m.

The last plot in the previous post was misleading, showing roughly equal parts genuine difference and error from solving the differential equation numerically. Here’s the code that was used to solve the nonlinear equation.

from scipy.special import ellipj, ellipk
from numpy import sin, cos, pi, linspace, arcsin
from scipy.integrate import solve_ivp

def exact_period(θ):
    return 2*ellipk(sin(θ/2)**2)/pi

def nonlinear_ode(t, z):
    x, y = z
    return [y, -sin(x)]    

theta0 = pi/3
b = 2*pi*exact_period(theta0)
t = linspace(0, 2*b, 2000)

sol = solve_ivp(nonlinear_ode, [0, 2*b], [theta0, 0], t_eval=t)

The solution is contained in sol.y[0].

Let’s compare the numerical solution to the exact solution.

def f(t, c, theta0):
    a = sin(theta0/2)
    m = a**2
    sn, cn, dn, ph = ellipj(c*t, m)
    return 2*arcsin(a*cn/dn)

There are a couple things to note about the code. First,SciPy doesn’t implement the cd function, but it can be computed as cn/dn. Second, the function ellipj returns four functions at once because it takes about as much time to calculate all four as it does to compute one of them.

Here is a plot of the error in solving the differential equation.

And here is the difference between the exact solution to the nonlinear pendulum equation and the stretched solution to the linear equation.

There is also a formula for the nth derivative of a quotient, but it’s more complicated and less known.

We start by writing the quotient rule in an unusual way.

Applying the quotient rule twice gives the following.

And here’s the general rule in all its glory.

Source: V. F. Ivanoff. The nth Derivative of a Fractional Function. The American Mathematical Monthly, Vol. 55, No. 8 (Oct., 1948), p. 491

where g is the acceleration due to gravity and ℓ is the length of the pendulum. When this is presented in an introductory physics class, the instructor will immediately say something like “we’re only interested in the case where θ is small, so we can rewrite the equation as

Questions

This raises a lot of questions, or at least it should.

1. Why not leave sin θ alone?
2. What justifies replacing sin θ with just θ?
3. How small does θ have to be for this to be OK?
4. How do the solutions to the exact and approximate equations differ?

First, sine is a nonlinear function, making the differential equation nonlinear. The nonlinear pendulum equation cannot be solved using mathematics that students in an introductory physics class have seen. There is a closed-form solution, but only if you extend “closed-form” to mean more than the elementary functions a student would see in a calculus class.

Second, the approximation is justified because sin θ ≈ θ when θ is small. That’s true, but it’s kinda subtle. Here’s a post unpacking that.

The third question doesn’t have a simple answer, though simple answers are often given. An instructor could make up an answer on the spot and say “less than 10 degrees” or something like that. A more thorough answer requires answering the fourth question.

I address how nonlinear affects the solutions in a post a couple years ago. This post will expand a bit on that post.

Longer period

The primary difference between the nonlinear and linear pendulum equations is that the solutions to the nonlinear equation have longer periods. The solution to the linear equation is a cosine. Solving the equation determines the frequency, amplitude, and phase shift of the cosine, but qualitatively it’s just a cosine. The solution to the corresponding nonlinear equation, with sin θ rather than θ, is not exactly a cosine, but it looks a lot like a cosine, only the period is a little longer [1].

OK, the nonlinear pendulum has a longer period, but how much longer? The period is increased by a factor f(θ₀) where θ₀ is the initial displacement.

You can find the exact answer in my earlier post. The exact answer depends on a special function called the “complete elliptic integral of the first kind,” but to a good approximation

The earlier post compares this approximation to the exact function.

Linear solution with adjusted period

Since the nonlinear pendulum equation is roughly the same as the linear equation with a longer period, you can approximate the solution to the nonlinear equation by solving the linear equation but increasing the period. How good is that approximation?

Let’s do an example with θ₀ = 60° = π/3 radians. Then sin θ₀ = 0.866 but θ₀, in radians, is 1.047, so we definitely can’t say sin θ₀ is approximately θ₀. To make things simple, let’s set ℓ = g. Also, assume the pendulum starts from rest, i.e. θ'(0) = 0.

Here’s a plot of the solutions to the nonlinear and linear equations.

Obviously the solution to the nonlinear equation has a longer period. In fact it’s 7.32% longer. (The approximation above would have estimated 7.46%.)

Here’s a plot comparing the solution of the nonlinear equation and the solution to the linear equations with period stretched by 7.32%.

The solutions differ by less than the width of the plotting line, so it’s too small to see. But we can see there’s a difference when we subtract the two solutions.

Here’s a plot of the solutions to the nonlinear and linear equations.

Update: The plot above is misleading. Part of what it shows is numerical error from solving the pendulum equation. When we redo the plot using the exact solution the error is about half as large. And the error is periodic, as we’d expect. See this post for more on the exact solution using Jacobi functions and the differential equation solver that was used to make the original plot.

Related posts

• Mechanical vibrations
• Bowie integrator and the nonlinear pendulum

[1] The period of a pendulum depends on its length ℓ, and so we can think of the nonlinear term effectively replacing ℓ by a longer effective length ℓ_eff.

in radians. The previous post worked in degrees, but here we’ll use radians.

If the triangle is oblique rather than a right triangle, there an approximation for the angle A that doesn’t require inverse trig functions, though it does require square roots. The approximation is derived in [1] using the same series that is the basis of the approximation in the earlier post, the power series for 2 csc(x) + cot(x).

For an oblique triangle, the approximation is

where s is the semiperimeter.

For comparison, we can find the exact value of A using the law of cosines.

and so

Here’s a little Python script to see how accurate the approximation is.

from math import sqrt, acos

def approx(a, b, c):
    "approximate the angle opposite a"
    s = (a + b + c)/2
    return 6*sqrt((s - b)*(s - c)) / (2*sqrt(b*c) + sqrt(s*(s - a)))

def exact(a, b, c):
    "exact value of the angle opposite a"    
    return acos((b**2 + c**2 - a**2)/(2*b*c))

a, b, c = 6, 7, 12
print( approx(a, b, c) )
print( exact(a, b, c) )

This prints

0.36387538476776243
0.36387760856668505

showing that in our example the approximation is good to five decimal places.

[1] H. E. Stelson. Note on the approximate solution of an oblique triangle without tables. American Mathematical Monthly. Vol 56, No. 2 (February, 1949), pp. 84–95.

A ≈ a 172° / (b + 2c).

The approximation above only involves simple arithmetic. No trig functions. Not even a square root. It could be carried out with pencil and paper or even mentally. And yet it is surprisingly accurate.

If we use the 3, 4, 5 triangle as an example, the exact value of the smallest angle is

A = arctan(3/4) × 180°/π ≈ 36.8699°

and the approximate value is

A ≈  3 × 172° / (4 + 2×5) = 258°/7 ≈ 36.8571°,

a difference of 0.0128°. When the angle is more acute the approximation is even better.

Derivation

Where does this magical approximation come from? It boils down to the series

2 csc(x) + cot(x) = 3/x + x³/60 + O(x⁴)

where x is in radians. When x is small,  x³/60 is extremely small and so we have

2 csc(x) + cot(x) ≈ 3/x.

Apply this approximation with csc(x) = c/a and cot(x) = b/a. and you have

x ≈ 3a/(b + 2c)

in radians. Multiply by 180°/π to convert to degrees, and note that 540/π ≈ 172.

Discovery

It’s unmotivated to say “just expand 2 csc(x) + cot(x) in a series.” Where did that come from?

There’s a line in [1] that says “It can been seen, either from tables or from a consideration of power series that the radian measure of a small angle lies approximately one-third of the way from the sine to the tangent.” In other words

3x ≈ 2 sin(x) + tan(x).

You can verify that by adding the power series and noting that the cubic terms cancel out.

But that’s just the beginning. The author then makes the leap to conjecturing that if the weighted arithmetic mean gives a good approximation, maybe the weighted harmonic mean gives an even better approximation, and that leads to considering

2 csc(x) + cot(x) ≈ 3/x.

Extension

See the next post for an extension to oblique triangles. Not as simple, but based on the same trick.

[1] J. S. Frame. Solving a right triangle without tables. The American Mathematical Monthly, Vol. 50, No. 10 (Dec., 1943), pp. 622-626

The models feel subjectively much smarter. They can accomplish a much broader range of tasks. They seem to have a larger, more comprehensive in-depth view of the code base and what you are trying to accomplish. They can locate more details in obscure parts of the code related to the specific task at hand. By rough personal estimate, I would say that they helped me with 20% of my coding effort last August and about 60% now. Both of these figures may be below what is possible, I may not be using them to fullest capability.

This being said, they are not a panacea. Sometimes they need help and direction to where to look for a problem. Sometimes they myopically look at the trees and don’t see the forest, don’t step back to see the big picture to understand something is wrong at a high level. Also they can overoptimize to the test harness. They can also generate code that is not conceptually consistent with existing code.

They can also generate much larger amounts of code than necessary. Some warn that that this will lead to explosion of coding debt. On the other hand, you can also guide the coding agent to refactor and improve its own code—quickly. In one case I’ve gotten it to reduce the size of one piece of code to less than half its original size with no change in behavior.

I use OpenAI Codex rather than Claude Code.  I’m glad to hear some technically credible people think this a good choice. Though maybe I should try both.

My work is a research project for which the code itself is the research product, so I can’t give specifications of everything in advance; writing the code itself is a process of discovery. Also, I want a code base that remains human-readable. So I am deeply involved in discussion with the coding agent that will sometimes go off to do a task for some period of time. It is not my desire to treat the agent as what some have called a dark software factory.

Some say they fear that using a coding agent will result in forgetting how to write code without one. I have felt this, but from having exercised this muscle for such a very long time I don’t think it is a skill I would easily forget. The flip side of this argument is, you might even learn new coding idioms from observing what the coding agent writes, that is a good thing.

Some say they haven’t written a line of code in weeks because the coding agent does it for them. I don’t think I’ll ever stop writing code, any more than I will stop scribbling random ideas on the back of an envelope, or typing out some number of lines of code by which I discover some new idea in real time while I am typing. Learning is multisensory.

My hat’s off to the developers who are able to keep many different agents in flight at the same time. Personally I have difficulty handling the cognitive load of thinking deeply about each one and doing a mental context switch across many agents. Though I am experimenting with running more than one agent so I can be doing something while another agent is working.

I continue to be astounded by productivity gains in some situations. I recently added a new capability, which normally I think would’ve taken me two months to learn a whole new algorithmic method and library to use. With the coding agent, it took me four days to get this done, on the order of 10X productivity increase. Though admittedly things are not always that rosy.

In short, the tools are getting better and better. I’m looking forward to what they will be like a few months or a year from now.

The previous post contained an illustration using a cubic polynomial and three evenly spaced parallel lines. This post uses a fifth degree polynomial, and shows that the parallel lines need not be evenly spaced.

In this post

f(x, y) = y³ + y − x (x + 1) (x + 2) (x − 3) (x − 2).

Here’s an example of three lines that each cross the curve five times.

The lines are y = x + k where k = 0.5, −0.5, and −3. The coordinates of the centroids are (0.4, 0.9), (0.4, -0.1), and (0.4, -2.6).

And to show that the requirement that the lines cross five times is necessary, here’s a plot where one of the parallel lines only crosses three times. The top line is now y = x + 2 and the centroid on the top line moved to (0.0550019, 2.055).

Why NF4

LLM parameters have a roughly Gaussian distribution, and so evenly spaced numeric values are not ideal for parameters. Instead, you’d like numbers that are closer together near 0.

The FP4 floating point numbers, described in the previous post, are spaced 0.5 apart for small values, and the larger values are spaced 1 or 2 apart. That’s hardly a Gaussian distribution, but it’s closer to Gaussian than a uniform distribution would be. NF4 deliberately follows more of a Gaussian distribution.

QLoRA

The QLoRA formats [1], unlike FP4, are not analogs of IEEE numbers. The bits are not interpreted as sign, exponent, and mantissa, but rather as integers to be used as indexes. An NFn number is an index into a list of 2ⁿ real numbers with Gaussian spacing. To put it another way, the numbers represented by NFn have uniformly distributed z-scores.

That makes sense at a high level, but the paper [1] is hard to follow in detail. It says

More formally, we estimate the 2^k values q_i of the data type as follows:

where Q_X(·) is the quantile function of the standard normal distribution N(0, 1).

The paper doesn’t give the range of i but it says there are 2^k values, implying that i runs from 0 to 2^k −1 or from 1 to 2^k. Either way runs into infinite values since Q(0) = −∞ and Q(1) = ∞. We could avoid infinities by letting i run from 1 to 2ⁿ − 1.

The next sentence is puzzling.

A problem for a symmetric k-bit quantization is that this approach does not have an exact representation of zero, which is an important property to quantize padding and other zero-valued elements with no error.

I understand the desire to represent 0 exactly, but the equation above has an exact representation of 0 when i = 2^{n − 1}. Perhaps the authors had in mind that i takes on the values ½, 1 + ½, 2 + ½, …, 2ⁿ − ½. This would be reasonable, but a highly unusual use of notation. It seems that the real problem is not the lack of a representation of 0 but an unused index, with i running from 1 to 2ⁿ − 1.

To be fair, the first sentence quoted above says “we estimate the 2^k values …” and so the equation above may not be intended as a definition but as motivation for the actual definition.

Reproducing NF4

The authors give a procedure for using 2ⁿ values of i and obtaining an exact representation of 0, and they give a list of NF4 values in Appendix E. I was not able to get the two to match. I implemented a few possible interpretations of the procedure described in the paper, and each approximates the list of values in the appendix, but not closely.

The following code, written with the help of ChatGPT, reverse engineers the NF4 values to 8 decimal places, i.e. to the precision of a 32-bit floating point number.

from scipy.stats import norm

Q = norm.ppf

α  = 0.9677083
Z  = Q(α)
δ1 = (α - 0.5)/7
δ2 = (α - 0.5)/8

q = [0]*16
for i in range(7):
    q[i] = -Q(α - i*δ1)/Z
for i in range(8):
    q[i+8] = Q(0.5 + (i+1)*δ2)/Z
    
# Values given in Appendix E
NF4 = [
    -1.0,
    -0.6961928009986877,
    -0.5250730514526367,
    -0.39491748809814453,
    -0.28444138169288635,
    -0.18477343022823334,
    -0.09105003625154495,
    0.0,
    0.07958029955625534,
    0.16093020141124725,
    0.24611230194568634,
    0.33791524171829224,
    0.44070982933044434,
    0.5626170039176941,
    0.7229568362236023,
    1.0
]

# Compare 
for i in range(16):
    print(i, NF4[i] - q[i])

The magic number α = 0.9677083 is a mystery. I asked ChatGPT to look into this further, and it said that bitsandbytes uses α = 929/960 = 0.9677083333333333. When I use this value for α the precision is about the same, which is fine. However, the values in the paper were given to 16 decimal places, so I thought it might be able to match the values to more precision.

Quibbles over the exact values of NF4 aside, the NF4 format works well in practice. Models. quantized to 4 bits using NF4 perform better than models quantized to other 4-bit formats on some benchmarks.

Related posts

• Anatomy of a floating point number
• Anatomy of a posit number
• Half precision numbers

[1] QLoRA: Efficient Finetuning of Quantized LLMs by Tim Dettmers, Artidoro Pagnoni, Ari Holtzman, and Luke Zettlemoyer. https://arxiv.org/abs/2305.14314.

Programmers were grateful for the move from 32-bit floats to 64-bit floats. It doesn’t hurt to have more precision, and some numerical problems go away when you go from 32-bits to 64-bits. (Though not all. Something I’ve written about numerous times.)

Neural networks brought about something extraordinary: demand for floating point numbers with less precision. These networks have an enormous number of parameters, and its more important to fit more parameters into memory than it is to have higher precision parameters. Instead of double precision (64 bits) developers wanted half precision (16 bits), or even less, such as FP8 (8 bits) or FP4 (4 bits). This post will look at 4-bit floating point numbers.

Why even bother with floating point numbers when you don’t need much precision? Why not use integers? For example, with four bits you could represent the integers 0, 1, 2, …, 15. You could introduce a bias, subtracting say 7 from each value, so your four bits represent −7 through 8. Turns out it’s useful to have a more dynamic range.

FP4

Signed 4-bit floating point numbers in FP4 format use the first bit to represent the sign. The question is what to do with the remaining three bits. The notation ExMy denotes a format with x exponent bits and y mantissa bits. In the context of signed 4-bit numbers,

x + y = 3

but in other contexts the sum could be larger. For example, for an 8-bit signed float, x + y = 7.

For 4-bit signed floats we have four possibilities: E3M0, E2M1, E1M2, and E0M3. All are used somewhere, but E2M1 is the most common and is supported in Nvidia hardware.

A number with sign bit s, exponent e, and mantissa m has the value

(−1)^s 2^e−b (1 + m/2)

where b is the bias. The purpose of the bias is to allow positive and negative exponents without using signed numbers for e. So, for example, if b = 1 and e = 1, 2, or 3 then the exponent part 2^e−b can represent 1, 2, or 4.

The bias impacts the range of possible numbers but not their relative spacing. For any value of bias b, the E3M0 format is all exponent, no mantissa, and so its possible values are uniformly distributed on a log scale. The E0M3 format is all mantissa, so its values are uniformly distributed on a linear scale. The E1M2 and E2M1 formats are unevenly spaced on both log and linear scales.

There is an exception to the expression above for converting (s, e, m) into a real number when e = 0. In that case, m = 0 represents 0 and m = 1 represents ½.

Table of values

Since there are only 16 possible FP4 numbers, it’s possible to list them all. Here is a table for the E2M1 format.

Bits s exp m  Value
-------------------
0000 0  00 0     +0
0001 0  00 1   +0.5
0010 0  01 0     +1
0011 0  01 1   +1.5
0100 0  10 0     +2
0101 0  10 1     +3
0110 0  11 0     +4
0111 0  11 1     +6
1000 1  00 0     -0
1001 1  00 1   -0.5
1010 1  01 0     -1
1011 1  01 1   -1.5
1100 1  10 0     -2
1101 1  10 1     -3
1110 1  11 0     -4
1111 1  11 1     -6

Note that even in this tiny floating point format, there are two zeros, +0 and −0, just like full precision floats. More on that here.

Pychop library

The Python library Pychop emulates a wide variety of reduced-precision floating point formats. Here is the code that used Pychop to create the table above.

import pychop

# Pull the format metadata from Pychop.
spec = pychop.MX_FORMATS["mxfp4_e2m1"]
assert (spec.exp_bits, spec.sig_bits) == (2, 1)

def e2m1_value(s: int, e: int, m: int) -> float:
    sign = -1.0 if s else 1.0

    # Subnormal / zero
    if e == 0:
        return sign * (m / 2.0)

    # Normal
    return sign * (2.0 ** (e - 1)) * (1.0 + m / 2.0)

def display_value(bits: int, x: float) -> str:
    if bits == 0b0000:
        return "+0"
    if bits == 0b1000:
        return "-0"
    return f"{x:+g}"

rows = []
for bits in range(16):
    s = (bits >> 3) & 0b1
    e = (bits >> 1) & 0b11
    m = bits & 0b1
    x = e2m1_value(s, e, m)

    rows.append(
        {
            "Bits": f"{bits:04b}",
            "s": s,
            "exp_bits": f"{e:02b}",
            "m": m,
            "Value": display_value(bits, x),
        }
    )

# Pretty-print the table.
header = f"{'Bits':<4} {'s':>1} {'exp':>3} {'m':>1} {'Value':>6}"
print(header)
print("-" * len(header))
for row in rows:
    print(
        f"{row['Bits']:<4} " f"{row['s']:>1} "
        f"{row['exp_bits']:>3} "
        f"{row['m']:>1} "
        f"{row['Value']:>6}"
    )

Other formats

FP4 isn’t the only 4-bit floating point format. There’s a surprisingly large number of formats in use. I intend to address another format my next post.

Update: See the next post for a discussion of NF4, a format whose representable numbers more closely matches the distribution of LLM weights.

Newton defined a diameter to be any line that crosses the zero set of f exactly n times. If

f(x, y) = x² + y² − 1

then the zero set of f is a circle and diameters of the circle in the usual sense are diameters in Newton’s sense. But Newton’s notion of diameter is more general, including lines the cross the circle without going through the center.

Newton’s theorem of diameters says that if you take several parallel diameters (in his sense of the word), the centroids of the intersections of each diameter with the curve f(x, y) = 0 all line on a line.

To illustrate this theorem, let’s look at the elliptic curve

y² = x³ − 2x + 1,

i.e. the zeros of f(x, y) = y² − (x³ − 2x + 1). This is a third degree curve, and so in general a straight line will cross the curve three times [2].

The orange, green, and red lines are parallel, each intersecting the blue elliptic curve three times. The dot on each line is the centroid of the intersection points, the center of mass if you imagine each intersection to be a unit point mass. The centroids all lie on a line, a vertical line in this example though in general the line could have any slope.

I hadn’t seen this theorem until I ran across it recently when skimming [3]. Search results suggest the theorem isn’t widely known, which is surprising for a result that goes back to Newton.

Update: See this post for another illustration using a fifth degree polynomial.

Related posts

• What is an elliptic curve?
• The fundamental theorem of algebra
• Finding where two quadratic curves intersect

[1] Bézout’s theorem says a curve of degree m and a curve of degree n will always intersect in mn points. But that includes complex roots, adds a line at infinity, and counts intersections with multiplicity. So a line, a curve of degree 1, will intersect a curve of degree n at n points in this extended sense.

[2] See the description of Bézout’s theorem in the previous footnote. In the elliptic curve example, the parallel lines meet at a point at infinity. A line that misses the closed component of the elliptic curve and only passes through the second component has 1 real point of intersection but there would be 2 more if we were working in ℂ² rather than ℝ².

In algebraic terms, the system of equations

y² = x³ − 2x + 1
3y = 2x + k

has three real solutions for small values of k, but for sufficiently large values of |k| two of the solutions will be complex.

[3] Mathematics: Its Content, Methods, and Meaning. Edited by A. D. Aleksandrov, A. N. Kolmogorov, and M. A. Lavrent’ev. Volume 1.

If you know the distance to a second satellite as well, then you can find two circles that contain your location. The two circles intersect at two points, and you know that you’re at one of two possible positions. If you know your approximate position, you may be able to rule out one of the intersection points.

If you know the distance to three different satellites, now you know three circles that you’re standing on, and the third circle will only pass through one of the two points determined by the first two satellites. Now you know exactly where you are.

Knowing the distance to more satellites is even better. In theory additional observations are redundant but harmless. In practice, they let you partially cancel out inevitable measurement errors.

If you’re not on the earth’s surface, you’re still at the intersection of n spheres if you know the distance to n satellites. If you’re in an airplane, or on route to the moon, the same principles apply.

Errors and corrections

How do you know the distance to a satellite? The satellite can announce what time it is by its clock, then when you receive the announcement you compare it to the time by your clock. The difference between the two times tells you how long the radio signal traveled. Multiply by the speed of light and you have the distance.

However, your clock will probably not be exactly synchronized with the satellite clock. Observing a fourth satellite can fix the problem of your clock not being synchronized with the satellite clocks. But it doesn’t fix the more subtle problems of special relativity and general relativity. See this post by Shri Khalpada for an accessible discussion of the physics.

Numerical computation

Each distance measurement gives you an equation:

|| x – s_i || = d_i

where s_i is the location of the ith satellite and d_i is your distance to that satellite. If you square both sides of the equation, you have a quadratic equation. You have to solve a system of nonlinear equations, and yet there is a way to transform the problem into solving linear equations, i.e. using linear algebra. See this article for details.

Related posts

• The navigation triangle
• Dutton’s Navigation and Piloting

A quick search didn’t turn up anything about Artzt parabolas [1], but apparently the parabolas go through pairs of vertices with tangents parallel to the sides.

The general form of a conic section is

ax² + bxy + cy² + dx + ey + f = 0

and the constraint b² = 4ac means the conic will be a parabola.

We have 6 parameters, each determined only up to a scaling factor; you can multiply both sides by any non-zero constant and still have the same conic. So a general conic has 5 degrees of freedom, and the parabola condition b² = 4ac takes us down to 4. Specifying two points that the parabola passes through takes up 2 more degrees of freedom, and specifying the slopes takes up the last two. So it’s plausible that there is a unique solution to the problem.

There is indeed a solution, unique up to scaling the parameters. The following code finds parameters of a parabola that passes through (x_i, y_i) with slope m_i for i = 1, 2.

def solve(x1, y1, m1, x2, y2, m2):
    
    Δx = x2 - x1
    Δy = y2 - y1
    λ = 4*(Δx*m1 - Δy)*(Δx*m2 - Δy)/(m1 - m2)**2
    k = x2*y1 - x1*y2

    a = Δy**2 + λ*m1*m2
    b = -2*Δx*Δy - λ*(m1 + m2)
    c = Δx**2 + λ
    d =  2*k*Δy + λ*(m1*y2 + m2*y1 - m1*m2*(x1 + x2))
    e = -2*k*Δx + λ*(m1*x1 + m2*x2 - y1 - y2)
    f = k**2 + λ*(m1*x1 - y1)*(m2*x2 - y2)

    return (a, b, c, d, e, f)

[1] The page said “Artz” when I first looked at it, but it has since been corrected to “Artzt”. Maybe I didn’t find anything because I was looking for the wrong spelling.

and the constant 1. The following equations, taken from the paper’s supplement, show how to bootstrap addition, subtraction, multiplication, and division from the eml function.

See the paper and supplement for how to obtain constants like π and functions like square and square root, as well as the standard circular and hyperbolic functions.

Related posts

• Bootstrapping a small math library
• Toffoli gates are all you need

The church fixed Easter to be the first Sunday after the first full moon after the Spring equinox. They were choosing a date in the Roman (Julian) calendar to commemorate an event whose date was known according to the Jewish lunisolar calendar, hence the reference to equinoxes and full moons.

The previous post explained why the Eastern and Western dates of Easter differ. The primary reason is that both churches use March 21 as the first day of Spring, but the Eastern church uses March 21 on the Julian calendar and the Western church uses March 21 on the Gregorian calendar.

But that’s not the only difference. The churches chose different algorithms for calculating when the first full moon would be. The date of Easter doesn’t depend on the date of the full moon per se, but the methods used to predict full moons.

This post will show why determining the date of the full moon is messy.

Lunation length

The moon takes between 29 and 30 days between full moons (or between new moons, which are easier to objectively measure). This period is called a lunation. The average length of a lunation is L = 29.530588853 days. This is not a convenient number to work with, and so there’s no simple way of reconciling the orbital period of the moon with the rotation period of the earth [1]. Lunar calendars alternate months with 29 and 30 days, but they can’t be very accurate, so they have to have some fudge factor analogous to leap years.

The value of L was known from ancient times. Meton of Athens calculated in 432 BC that 235 lunar cycles equaled 19 tropical years or 6940 days. This corresponds to L ≈ 29.5319. Around a century later the Greek scholar Callippus refined this to 940 cycles in 76 years or 27,759 days. This corresponds to L ≈ 29.53085.

The problem wasn’t knowing L but devising a convenient way of working with L. There is no way to work with lunations that is as easy as the way the Julian (or even the more complicated Gregorian) calendar reconciles days with years.

Approximations

Let’s look at the accuracy of several approximations for L. We’d like an approximation that is not only accurate in an absolute sense, but also accurate relative to its complexity. The complexity of a fraction is measured by a height function. We’ll use what’s called the “classic” height function: log( max(n, d) ) where n and d are the numerator and denominator of a fraction. Since we’re approximating a number bigger than 1, this will be simply log(n).

We will compare the first five convergents, approximations that come from the continued fraction form of L, and the approximations of Meton and Callippus. Here’s a plot.

And here’s the code that produced the plot, showing the fractions used.

from numpy import log
import matplotlib.pyplot as plt

fracs = [
    (30, 1), 
    (59, 2),
    (443, 15),
    (502, 17),
    (1447, 49),
    (6940, 235),
    (27759, 940)
]

def error(n, d):
    L = 29.530588853    
    return abs(n/d - L)

for f in fracs:
    plt.plot(log(f[0]), log(error(*f)), 'o')
plt.xlabel("log numerator")
plt.ylabel("log error")
plt.show()

The approximation 1447/49 is the best by far, both in absolute terms and relative to the size of the numerator. But it’s not very useful for calendar design because 1447 is not nicely related to the number of days in a year.

[1] The time between full moons is a synodic month, the time it takes for the moon to return to the same position relative to the sun. This is longer than a sidereal month, the time it takes the moon to complete one orbit relative to the fixed stars.