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.

Why the dates differ

Easter is on the first Sunday after the first full moon in Spring [1]. East and West agree on this. What they disagree on is the details of “full moon” and “Spring.” The dates are not based on precise astronomical measurements but rather on astronomical approximations codified long ago.

Spring begins on March 21 for the purposes of calculating Easter. But the Western church uses March 21 on the Gregorian calendar and the Eastern church uses March 21 on the Julian calendar. This mostly accounts for the difference between Eastern and Western dates for Easter. East and West also use slightly different methods of approximating when the moon will be full. More on that in the next post.

Pascha never comes before Easter

The Eastern name for Easter is Pascha. Eastern Pascha and Western Easter can occur on the same day, but otherwise Pascha is always later, never earlier. This is because the Julian year is longer than the Gregorian year, causing fixed dates on the former calendar to occur after the later. Also, the Eastern method of approximating the date of the Paschal full moon gives a later date than the Western method.

The Julian calendar has exactly 365 1/4 days. The Gregorian calendar has 365 97/400 days; centuries are not leap years unless they’re divisible by 400. This complication in the Gregorian calendar was necessary to match the solar year. The date March 21 on the Julian calendar is drifting later in the year from the perspective of the Gregorian calendar, moving further past the astronomical equinox [2].

Size of the gap

Eastern and Western dates of Easter can coincide. The were the same last year, and will be the same again in 2028. The gap is always a whole number of weeks because Easter is always on a Sunday.

The gap is usually 1 week. It can be 0, 4, or 5 weeks, but never 2 or 3 weeks.

This is the pattern for now. Sometime in the distant future the Julian and Gregorian calendars will diverge further than the gaps will increase. Presumably Orthodox churches will make some sort of adjustment before the Julian date March 21 drifts into summer or fall.

Related posts

• Cycle of New Year’s Days
• When was Newton born?
• Martian leap years

[1] The reason for this definition is that Christ was crucified at the time of the Passover. Due to the lunisolar design of the Jewish calendar, this would have been during the first full moon after the Spring equinox. Christ rose from the dead the Sunday following the crucifixion, so Easter is on the first Sunday after the first full moon of Spring.

[2] The Julian and Gregorian calendars currently differ by 13 days, and they’re drifting apart at the rate of 3 days every 400 years. Somewhere around 47,000 years from now the two calendars will agree again, sorta, because the Julian calendar will be a full year behind the Gregorian calendar.

Let p > 5 be a prime number and 0 < k < p. Then the digits in k/p might be the same for all k, varying only by cyclic permutations. This is the case, for example, when p = 7 or p = 17. More on these kinds of fractions here.

The digits in k/p repeat for every k, but different values of k might have sequences of digits that vary by more than cyclic permutations. For example, let’s look at the values of k/13.

>>> for i in range(1, 13):
...   print(i/13)
...
 1 0.0769230769230769
 2 0.1538461538461538
 3 0.2307692307692307
 4 0.3076923076923077
 5 0.3846153846153846
 6 0.4615384615384615
 7 0.5384615384615384
 8 0.6153846153846154
 9 0.6923076923076923
10 0.7692307692307693
11 0.8461538461538461
12 0.9230769230769231

One cycle goes through the digits 076923. You’ll see this when k = 1, 3, 4, 9, 10, or 11. The other cycle goes through 153846 for the rest of the values of k. The cycles 076923 and 153846 are called the distinct repeating sets of 13 in [1].

If we look at fractions with denominator 41, thee are six distinct repeating sets.

02439
04878
07317
09756
12195
14634
26829
36585

You could find these by modifying the Python code above. However, in general you’ll need more than default precision to see the full periods. You might want to shift over to bc, for example.

When you look at all the distinct repeating sets of a prime number, all digits appear almost the same number of times. Some digits may appear one more time than others, but that’s as uneven as you can get. A corollary in [1] states that if p = 10q + r, with 0 < r < 10, then 11 − r digits appear q times, and r − 1 digits appear q + 1 times.

Looking back at the example with p = 13, we have q = 1 and r = 3. The corollary says we should expect 8 digits to appear once and 2 digits to appear twice. And that’s what we see: in the sets 076923 and 153846 we have 3 and 6 repeated twice and the remaining 8 digits appear once.

In the example with p = 41, we have q = 4 and r = 1. So we expect all 10 digits to appear 4 times, which is the case.

Related posts

• Cyclic fractions
• Periods of fractions
• Calendars and continued fractions

[1] James K. Schiller. A Theorem in the Decimal Representation of Rationals. The American Mathematical Monthly
Vol. 66, No. 9 (Nov., 1959), pp. 797-798

I looked into this, expecting it to be approximately true. It’s exactly true in the sense that the speed of light in vacuum is 299,792,458 m/s and the line of latitude 29.9792458° N passes through the pyramid. The exact center of the pyramid is at 29.97917° N, 31.13417° E.

Of course this is a coincidence. Even if you believe that somehow the ancient Egyptians knew the speed of light, the meter was defined four millennia after the pyramid was built.

Then at each subsequent step, pick a random side of the hexagon and create the triangle formed by that side and c. Update c to be the center of the new triangle and plot c.

Note that you only choose a random point inside the hexagon once. After that you randomly choose sides.

Now there are many ways to define the center of a triangle. I assumed the original meant barycenter (centroid) when it said “center”, and apparently that was correct. I was able to create a similar figure.

But if you define center differently, you get a different image. For example, here’s what you get when you use the incenter, the center of the largest circle inside the triangle.

Related posts

• Randomly selecting points in a triangle
• Subdividing a triangle with various centers
• Randomly generated dragon fractal
• The chaos game and the Sierpinski triangle

In symbols, Andrica’s conjecture says that if p_n and p_n+1 are consecutive prime numbers, then

√p_n+1 − √p_n < 1.

This has been empirically verified for primes up to 2 × 10¹⁹.

If the conjecture is true, it puts an upper bound on how long you’d have to search to find the next prime:

p_n+1 < 1 + 2√p_n  + p_n,

which would be an improvement on the Bertrand-Chebyshev theorem that says

p_n+1 < 2p_n.

Because Artemis I was unmanned, it could spend a lot more time in orbit. The Artemis I mission took 25 days while Artemis II will take 10 days. Artemis I took an unusual path, orbiting the moon the opposite direction of the moon’s orbit around earth. This video by Primal Space demonstrates the orbit both from the perspective of earth and from the perspective of the Moon.

Another video from Primal Space describes the orbit of the third stage of Apollo 12. This stage was supposed to orbit around the sun in 1971, but an error sent it on a complicated unstable orbit of the earth, moon, and sun. It returned briefly to earth in 2002 and expected to return sometime in the 2040s.

Related posts

• Finding Lagrange points
• When are Lagrange points stable?
• Bad takes on chaos theory

E ≥ log(2) k_B T

where k_B is the Boltzmann constant and T is the ambient temperature in Kelvin. The lower bound applies no matter how the bit is physically stored. There is no theoretical lower limit on the energy required to carry out a reversible calculation.

In practice the energy required to erase a bit is around a billion times greater than Landauer’s lower bound. You might reasonably conclude that reversible computing isn’t practical since we’re nowhere near the Landauer limit. And yet in practice reversible circuits have been demonstrated to use less energy than conventional circuits. We’re far from the ultimate physical limit, but reversibility still provides practical efficiency gains today.

A Toffoli gate is a building block of reversible circuits. A Toffoli gate takes three bits as input and returns three bits as output:

T(a, b, c) = (a, b, c XOR (a AND b)).

In words, a Toffoli gate flips its third bit if and only if the first two bits are ones.

A Toffoli gate is its own inverse, and so it is reversible. This is easy to prove. If a = b = 1, then the third bit is flipped. Apply the Toffoli gate again flips the bit back to what it was. If ab = 0, i.e. at least one of the first two bits is zero, then the Toffoli gate doesn’t change anything.

There is a theorem that any Boolean function can be computed by a circuit made of only NAND gates. We’ll show that you can construct a NAND gate out of Toffoli gates, which shows any Boolean function can be computed by a circuit made of Toffoli gates, which shows any Boolean function can be computed reversibly.

To compute NAND, i.e. ¬ (a ∧ b), send (a, b, 1) to the Toffoli gate. The third bit of the output will contain the NAND of a and b.

T(a, b, 1) = (a, b, ¬ (a ∧ b))

A drawback of reversible computing is that you may have to send in more input than you’d like and get back more output than you’d like, as we can already see from the example above. NAND takes two input bits and returns one output bit. But the Toffoli gate simulating NAND takes three input bits and returns three output bits.

Related posts

• Fredkin automata
• Machine learning by satisfiability solving
• Minimizing Boolean expressions

There are HIPAA-compliant cloud options, but they’re both restrictive and expensive. Even enterprise options are not “HIPAA compliant” out of the box. Instead, they are “HIPAA eligible” or that they “support HIPAA compliance,” because you still need the right Business Associate Agreement (BAA), configuration, logging, access controls, and internal process around it, and the end product often ends up far less capable than a frontier model. The least expensive and therefore most accessible services do not even allow this as an option.

Specific examples:

• Only sales-managed ChatGPT Enterprise or Edu customers are eligible for a BAA, and OpenAI explicitly says it does not offer a BAA for ChatGPT Business. The consumer ChatGPT Health product says HIPAA and BAAs do not apply. ChatGPT for Healthcare pricing is based on ChatGPT Enterprise, depends on organization size and deployment needs, and requires contacting sales. Even within Enterprise, OpenAI’s Regulated Workspace spec lists Codex and the multi-step Agent feature as “Non-Included Functionality,” i.e. off limits for PHI.
• Google says Gemini can support HIPAA workloads in Workspace, but NotebookLM is not covered by Google’s BAA, and Gemini in Chrome is automatically blocked for BAA customers. If a work or school account does not have enterprise-grade data protections, chats in the Gemini app may be reviewed by humans and used to improve Google’s products.
• GitHub Copilot, despite being a Microsoft product, is not under Microsoft’s BAA. Azure OpenAI Service is, but only for text endpoints. While Microsoft is working on their own models, it is unlikely that they will deviate significantly here.
• Anthropic says its BAA covers only certain “HIPAA-ready” services, namely the first-party API and a HIPAA-ready Enterprise plan, and does not cover Claude Free, Pro, Max, Team, Workbench, Console, Cowork, or Claude for Office. The HIPAA-ready Enterprise offering is sales-assisted only. Bundled Claude Code seats are not covered. AWS Bedrock API calls can work, but this requires extensive configuration and carries its own complexities and restrictions.

Running AI locally is already practical as of early 2026. Open-weight models that approach the quality of commercial coding assistants run on consumer hardware. A single high-end GPU or a recent Mac with enough unified memory can run a 70B-parameter model at a reasonable token speed.

There’s an interesting interplay between economies of scale and diseconomies of scale. Cloud providers can run a data center at a lower cost per server than a small company can. That’s the economies of scale. But running HIPAA-compliant computing in the cloud, particularly with AI providers, incurs a large direct costs and indirect bureaucratic costs. That’s the diseconomies of scale. Smaller companies may benefit more from local AI than larger companies if they need to be HIPAA-compliant.

Related posts

• HIPAA expert determination
• An AI Odyssey
• Queueing and economies of scale

[1] This post is not legal advice. My clients are often lawyers, but I’m not a lawyer.

Suppose you’re keeping up with your average grade in a class, and you know your average after n tests, all weighted equally.

m = (x₁ + x₂ + x₃ + … + x_n) / n.

Then you get another test grade back and your new average is

m′ = (x₁ + x₂ + x₃ + … + x_n + x_n+1) / (n + 1).

You don’t need the individual test grades once you’ve computed the average; you can instead remember the average m and the number of grades n [1]. Then you know the sum of the first n grades is nm and so

m′ = (nm + x_n+1) / (n + 1).

You could split that into

m′ = w₁ m + w₂ x_n+1

where w₁ = n/(n + 1) and w₂ = 1/(n + 1). In other words, the new mean is the weighted average of the previous mean and the new score.

A Bayesian perspective would say that your posterior expected grade m′ is a compromise between your prior expected grade m and the new data x_n+1. [2]

You could also rewrite the equation above as

m′ = m + (x_n+1 − m)/(n + 1) = m + KΔ

where K = 1/(n + 1) and Δ = x_n+1 − m. In Kalman filter terms, K is the gain, the proportionality constant for how the change in your state is proportional to the difference between what you saw and what you expected.

Related posts

• A Bayesian view of Amazon Resellers
• Kalman filters and functional programming
• Kalman filters and bottom-up learning

[1] In statistical terms, the mean is a sufficient statistic.

[2] You could flesh this out by using a normal likelihood and a flat improper prior.

Perilune is sometimes called periselene. The two terms come from two goddesses associated with the moon, the Roman Luna and the Greek Selene. Since the peri- prefix is Greek, perhaps periselene would be preferable. But we’re far more familiar with words associated with the moon being based on Luna than Selene.

The neutral terms for closest and furthest points in an orbit are periapsis and apoapsis. but there are more colorful terms that are specific to orbiting particular celestial objects. The terms perigee and apogee for orbiting earth (from the Greek Gaia) are most familiar, and the terms perihelion and aphelion (not apohelion) for orbiting the sun (from the Greek Helios) are the next most familiar.

The terms perijove and apojove are unfamiliar, but you can imagine what they mean. Others like periareion and apoareion, especially the latter, are truly arcane.

Napier’s mnemonic is a clever way to remember 10 equations in spherical trig. See the linked post for the meanings of the variables.

sin a = sin A sin c = tan b cot B
sin b = sin B sin c = tan a cot A
cos A = cos a sin B = tan b cot c
cos B = cos b sin A = tan a cot c
cos c = cot A cot B = cos a cos b

The hyperbolic analog replaces every circular function of a, b, or c with its hyperbolic counterpart.

sinh a = sin A sinh c = tanh b cot B
sinh b = sin B sinh c = tanh a cot A
cos A = cosh a sin B = tanh b coth c
cos B = cosh b sin A = tanh a coth c
cosh c = cot A cot B = cosh a cosh b

[1] D. M. Y. Sommerville. The Elements of Non-Euclidean Geometry. 1919.

Apollo 8 spent 2 hours and 44 minutes in low earth orbit (LEO) before performing trans-lunar injection (TLI) and heading toward the moon. Artemis II made one low earth orbit before moving to high earth orbit (HEO) where it will stay for around 24 hours before TLI. The Apollo 8 LEO was essentially circular at an altitude of around 100 nautical miles. The Artemis II HEO is highly eccentric with an apogee of around 40,000 nautical miles.

Apollo 8 spent roughly three days traveling to the moon, measured as the time between TLI and lunar insertion orbit. Artemis II will not orbit the moon but instead swing past the moon on a “lunar free-return trajectory” like Apollo 13. The time between Artemis’ TLI and perilune (the closest approach to the moon, on the far side) is expected to be about four days. For Apollo 13, this period was three days.

The furthest any human has been from earth was the Apollo 13 perilune at about 60 nautical miles above the far side of the moon. Artemis is expected to break this record with a perilune of between 3,500 and 5,200 nautical miles.

Related posts

• Sphere of influence
• Arenstorf’s orbit
• Math developed for going to the moon

In short,

P_n = T_{2n − 1} − T_n.

This is trivial to prove algebraically, though the visual proof above is more interesting.

The proof follows immediately from the definition of pentagonal numbers

P_n = (3n² − n)/2

and triangular numbers

T_n = (n² − n)/2.

Related posts

• Partitions and pentagons
• Tetrahedral numbers
• A pentagon, hexagon, and decagon walk into a bar …
  

Quantum computers exist now, but the question is when and if a cryptographically relevant quantum computer (CRQC) is coming. At the moment, a quantum computer cannot factor 21 without cheating (i.e. not implementing circuits that you know a priori won’t be needed). But that could change suddenly. And some believe quantum computers could quickly go from being able to factor numbers with two digits to being able to factor numbers with thousands of digits (i.e. breaking RSA encryption) without much incremental transition between.

The move to post-quantum encryption may be a lot like Y2K, fixing vast amounts of 20th century software that represented years with two digits. Y2K turned out to be a big nothingburger, but only because the world spent half a trillion dollars on preparation to make sure it would be a big nothingburger.

Programmers in the 1970s obviously knew that the year 2000 was coming, but they also knew that they needed to conserve bytes at the time. And they assumed, reasonably but incorrectly, that their software would all be replaced before two-digit dates became a problem.

Programmers still need to conserve bytes, though this is less obvious today. Quantum-resistant signatures and encryption keys are one or two orders of magnitude bigger. This takes up bandwidth and storage space, which may or may not be a significant problem, depending on context. Programmers may conclude that it’s not (yet) worth the extra overhead to use post-quantum encryption. Like their counterparts 50 years ago, they may assume, rightly or wrongly, that their software will be replaced by the time it needs to be.

Moving to post-quantum cryptography ASAP is not a great idea if you can afford to be more strategic. It takes many years to gain confidence that new encryption algorithms are secure. The SIKE algorithm, for example, was a semi-finalist the NIST post-quantum encryption competition, but someone found a way to break it using an hour of computing on a laptop.

Another reason to not be in a hurry is that it may be possible to be more clever than simply replacing pre-quantum algorithms with post-quantum analogs. For example, some blockchains are exploring zero-knowledge proofs as a way to aggregate signatures. Simply moving to post-quantum signatures could make every transaction block 100 times bigger. But replacing a set of signatures by a (post-quantum) zero-knowledge proof of the existence of the signatures, transaction blocks could be smaller than now.

As with Y2K, the move to post-quantum cryptography will be gradual. Some things have already moved, and some are in transition now. You may have seen the following warning when connecting to a remote server.

** WARNING: connection is not using a post-quantum key exchange algorithm.
** This session may be vulnerable to "store now, decrypt later" attacks.
** The server may need to be upgraded. See https://openssh.com/pq.html

Key sizes don’t matter as much to sftp connections as they do to blockchains. And the maturity of post-quantum algorithms is mitigated by OpenSSH using hybrid encryption: well-established encryption (like ECDH) wrapped by newer quantum-resistant encryption (like MK-KEM). If the newer algorithm isn’t as secure as expected, you’re no worse off than if you had only used the older algorithm.

When clocks rolled over from 1999 to 2000 without incident, many people felt the concern about Y2K had been overblown. Maybe something similar will happen with quantum computing. Let’s hope so.

Related posts

• Mixing error-correcting codes and cryptography
• Regression ond PQC
• Blockchains and cryptocurrency

The receive side of the coin is a decision tree for decoding Morse code. The shape is what makes this one interesting.

Decision trees are typically not very compact. Each branch is usually on its own horizontal level, with diagonal lines going down from each node to its children. But by making the lines either horizontal or vertical, the tree fits nicely into a circle.

I thought for a second that the designer had made the choices of horizontal or vertical segments in order to make the tree compact, but that’s not so. The direction of the path through the tree changes when and only when the Morse code switches from dot to dash or dash to dot.

It would be fun to play around with this, using the same design idea for other binary trees.

Related posts

• Family tree numbering
• Q code tree
• Morse code and psychological limits

Not finding the answer immediately on the product page, I thought I’d try clicking the AI shopping assistant helper tool to ask this question.

I waited with anticipation for an answer I would expect be more informative and useful than a standard search result. But it was not to be. The AI tool had nothing worthwhile.

Then I decided on an old-fashioned keyword search across all the product reviews. And, lo and behold, I immediately found several credible reviews addressing my question.

It is not good usability when multiple search mechanisms exist but only one of them is reliable. And it is surprising that a retrieval-based approach (e.g., RAG) could not at least match the effectiveness of a simple keyword search over reviews.

Models are capable, but effective integration can be lacking. Without improvements for cases like this, customers will not be satisfied users of these new AI tools.

Related posts

• An AI Odyssey, Part 1: Correctness Conundrum
• An AI Odyssey, Part 2: Prompting Peril

If you’re using Python, for example, and you don’t have NumPy installed, you can use the built-in math library, but it will not accept complex inputs.

>>> import math
>>> math.sin(3 + 4j)
Traceback (most recent call last):
File "<stdin>", line 1, in <module>
TypeError: must be real number, not complex

You can use the following identities to calculate sine and cosine for complex arguments using only real functions.

The proof is very simple: just use the addition formulas for sine and cosine, and the following identities.

The following code implements sine and cosine for complex arguments using only the built-in Python functions that accept real arguments. It then tests these against the NumPy versions that accept complex arguments.

from math import *
import numpy as np

def complex_sin(z):
    x, y = z.real, z.imag
    return sin(x)*cosh(y) + 1j*cos(x)*sinh(y)

def complex_cos(z):
    x, y = z.real, z.imag
    return cos(x)*cosh(y) - 1j*sin(x)*sinh(y)

z = 3 + 4j
mysin = complex_sin(z)
mycos = complex_cos(z)
npsin = np.sin(z)
npcos = np.cos(z)
assert(abs(mysin - npsin) < 1e-14)
assert(abs(mycos - npcos) < 1e-14)

Related posts

• Why j for imaginary unit?
• Bootstrapping a minimal math library
• Osborn's rule

where h is the spacing between interpolation points and δ is the error in the tabulated values. The constant c depends on the function f being interpolated, and to a lesser extent on n. The constant λ is independent of f but depends on n and on the relative spacing between the interpolation nodes. This post will look closer at λ.

Given a set of n + 1 nodes T

define

Then the Lebesgue function is defined by

and the Lebesgue constant for the grid is the maximum value of the Lebesgue function

The values of Λ are difficult to compute, but there are nice asymptotic expressions for Λ when the grid is evenly spaced:

When the grid points are at the roots of a Chebyshev polynomial then

The previous post mentioned the cases n = 11 and n = 29 for evenly spaced grids. The corresponding values of Λ are approximately 155 and 10995642. So 11th order interpolation is amplifying the rounding error in the tabulated points by a factor of 155, which might be acceptable. But 29th order interpolation is amplifying the rounding error by a factor of over 10 million.

The corresponding values of Λ for Chebyshev-spaced nodes are 2.58 and 3.17. Chebyshev spacing is clearly better for high-order interpolation, when you have that option.

If you want to know the value of a tabulated function between values of x given in the table, you have to use interpolation. Linear interpolation is often adequate, but you could get more accurate results using higher-order interpolation.

Suppose you have a function f(x) tabulated at x = 3.00, 3.01, 3.02, …, 3.99, 4.00 and you want to approximate the value of the function at π. You could approximate f(π) using the values of f(3.14) and f(3.15) with linear interpolation, but you could also take advantage of more points in the table. For example, you could use cubic interpolation to calculate f(π) using f(3.13), f(3.14), f(3.15), and f(3.16). Or you could use 29th degree interpolation with the values of f at 3.00, 3.01, 3.02, …, 3.29.

The Lagrange interpolation theorem lets you compute an upper bound on your interpolation error. However, the theorem assumes the values at each of the tabulated points are exact. And for ordinary use, you can assume the tabulated values are exact. The biggest source of error is typically the size of the gap between tabulated x values, not the precision of the tabulated values. Tables were designed so this is true [1].

The bound on order n interpolation error has the form

c h^{n + 1} + λ δ

where h is the spacing between interpolation points and δ is the error in the tabulated values. The value of c depends on the derivatives of the function you’re interpolating [2]. The value of λ is at least 1 since λδ is the “interpolation” error at the tabulated points.

The accuracy of an interpolated value cannot be better than δ in general, and so you pick the value of n that makes c h^{n + 1} less than δ. Any higher value of n is not helpful. And in fact higher values of n are harmful since λ grows exponentially as a function of n [3].

See the next post for mathematical details regarding the λs.

Examples

Let’s look at a specific example. Here’s a piece of a table for natural logarithms from A&S.

Here h = 10⁻³, and so linear interpolation would give you an error on the order of h² = 10⁻⁶. You’re never going to get error less than 10⁻¹⁵ since that’s the error in the tabulated values, so 4th order interpolation gives you about as much precision as you’re going to get. Carefully bounding the error would require using the values of c and λ above that are specific to this context. In fact, the interpolation error is on the order of 10⁻⁸ using 5th order interpolation, and that’s the best you can do.

I’ll briefly mention a couple more examples from A&S. The book includes a table of sine values, tabulated to 23 decimal places, in increments of h = 0.001 radians. A rough estimate would suggest 7th order interpolation is as high as you should go, and in fact the book indicates that 7th order interpolation will give you 9 figures of accuracy,

Another table from A&S gives values of the Bessel function J₀ in with 15 digit values in increments of h = 0.1. It says that 11th order interpolation will give you four decimal places of precision. In this case, fairly high-order interpolation is useful and even necessary. A large number of decimal places are needed in the tabulated values relative to the output precision because the spacing between points is so wide.

Related posts

• Using a table of logarithms
• Bessel, Everett, and Lagrange interpolation

[1] I say were because of course people rarely look up function values in tables anymore. Tables and interpolation are still widely used, just not directly by people; computers do the lookup and interpolation on their behalf.

[2] For functions like sine, the value of c doesn’t grow with n, and in fact decreases slowly as n increases. But for other functions, c can grow with n, which can cause problems like Runge phenomena.

[2] The constant λ grows exponentially with n for evenly spaced interpolation points, and values in a table are evenly spaced. The constant λ grows only logarithmically for Chebyshev spacing, but this isn’t practical for a general purpose table.

Theorem (Markov): If P(x) is a real polynomial of degree n, and |P(x)| ≤ 1 on [−1, 1] then |P′(x)| ≤ n² on [−1, 1].

If P(x) is a trigonometric polynomial then Bernstein proved that the bound decreases from n² to n.

Theorem (Bernstein): If P(x) is a trigonometric polynomial of degree n, and |P(z)| ≤ 1 on [−π, π] then |P′(x)| ≤ n on [−π, π].

Now a trigonometric polynomial is a truncated Fourier series

and so the max norm of the T′ is no more than n times the max norm of T.

This post and the previous one were motivated by Terence Tao’s latest post on Bernstein theory.

Related posts

• Zeros of trig polynomials
• Convert between real and complex Fourier series
• Systematically solving trig equations

I ran across Mendeleev’s theorem via a paper by Boas [1]. The opening paragraph describes what Mendeleev was working on.

Some years after the chemist Mendeleev invented the periodic table of the elements he made a study of the specific gravity of a solution as a function of the percentage of the dissolved substance. This function is of some practical importance: for example, it is used in testing beer and wine for alcoholic content, and in testing the cooling system of an automobile for concentration of anti-freeze; but present-day physical chemists do not seem to find it as interesting as Mendeleev did.

Mendeleev fit his data by patching together quadratic polynomials, i.e. he used quadratic splines. A question about the slopes of these splines lead to the following.

Theorem (Mendeleev): Let P(x) be a quadratic polynomial on [−1, 1] such that |P(x)| ≤ 1. Then |P′(x)| ≤ 4.

Mendeleev showed his result to mathematician Andrey Markov who generalized it to the following.

Theorem (Markov): If P(x) is a real polynomial of degree n, and |P(x)| ≤ 1 on [−1, 1] then |P′(x)| ≤ n² on [−1, 1].

Both inequalities are sharp with equality if and only if P(x) = ±T_n(x), the nth Chebyshev polynomial. In the special case of Mendeleev’s inequality, equality holds for

T₂(x) = 2x² − 1.

Andrey Markov’s brother Vladimir proved an extension of Andrey’s theorem to higher derivatives,

Related posts

• Periods in the periodic table
• How to memorize the periodic table
• The Brothers Markov

[1] R. P. Boas, Jr. Inequalities for the Derivatives of Polynomials. Mathematics Magazine, Vol. 42, No. 4 (Sep., 1969), pp. 165–174

If A and B are two sorted lists,

    comm A B

prints A − B, B − A, and A ∩ B. You usually don’t want all three, and so comm lets you filter the output. It’s a little quirky in that you specify what you don’t want instead of what you do. And you have to remember that 1, 2, and 3 correspond to A − B, B − A, and A ∩ B respectively.

A couple little scripts can hide the quirks. I have a script intersect

    comm -12 <(sort "$1") <(sort "$2")

and another script setminus

    comm -23 <(sort "$1") <(sort "$2")

that sort the input files on the fly and eliminate the need to remember comm‘s filtering syntax.

The setminus script computes A − B. To find B − A call the script with the arguments reversed.

Embedded regular expression modifiers address one of the environmental complications by putting the modifier in the regular expression itself.

For example, if you want to make a grep search case-insensitive, you pass it the -i flag. But if you want to make a regex case-insensitive inside a Python program, you pass a function the argument re.IGNORECASE. But if you put (?i) at the beginning of your regular expression, then the intention to make the match case-insensitive is embedded directly into the regex. You could use the regex in any environment that supports (?i) without having to know how to specify modifiers in that environment.

I was debugging a Python script this morning that worked under one version of Python and not under another. The root of the problem was that it was using re.findall() with several huge regular expression that had embedded modifiers. That was OK up to Python 3.5, then it was a warning between versions 3.6 and 3.10, and it’s an error in versions 3.11 and later.

The problem isn’t with all embedded modifiers, only global modifiers that don’t appear at the beginning of the regex. Older versions of Python, following Perl’s lead, would let you put a modifier like (?i) in the middle of a regex, and apply the modifier from that point to the end of the expression. In the latest versions of Python, you can either place the modifier at the beginning of the regex, or use a scoped modifier like (?:…) in the middle of the expression.

I didn’t want to edit the regular expressions in my code—some had over a thousand characters—so I changed re.findall() to regex.findall(). The third-party regex module is generally more Perl-compatible than Python’s standard re module.

The most common answer is the Bohr-Mollerup theorem. This theorem says that if f: (0, ∞) → (0, ∞) satisfies

1. f(x + 1) = x f(x)
2. f(1) = 1
3. log f is convex

then f(x) = Γ(x). The theorem applies on the positive real axis, and there is a unique holomorphic continuation of this function to the complex plane.

But the Bohr-Mollerup theorem is not the only theorem characterizing the gamma function. Another theorem was by Helmut Wielandt. His theorem says that if f is holomorphic in the right half-plane and

1. f(z + 1) = z f(z)
2. f(1) = 1
3. f(z) is bounded for {z: 1 ≤ Re z ≤ 2}

then f(x) = Γ(x). In short, Wielandt replaces the log-convexity for positive reals with the requirement that f is bounded in a strip of the complex plane.

You might wonder what is the bound alluded to in Wielandt’s theorem. You can show from the integral definition of Γ(z) that

|Γ(z)| ≤ |Γ(Re z)|

for z in the right half-plane. So the bound on the complex strip {z: 1 ≤ Re z ≤ 2} equals the bound on the real interval [1, 2], which is 1.

To be more precise, assume we have a series of the form

where the a_i are positive and monotonically converge to zero. Then the tail of the series is bounded by its first term:

The more we can say about the behavior of the a_i the more we can say about the remainder. So far we’ve assumed that these terms go monotonically to zero. If their differences

also go monotonically to zero, then we have an upper and lower bound on the truncation error:

If the differences of the differences,

also converge monotonically to zero, we can get a larger lower bound and a smaller upper bound on the remainder. In general, if the differences up to order k of the a_i go to zero monotonically, then the remainder term can be bounded as follows.

Source: Mark B. Villarino. The Error in an Alternating Series. American Mathematical Monthly, April 2018, pp. 360–364.

Related posts

• Euler’s method for accelerating an alternating series
• Cohen’s method for accelerating an alternating series

Theorem: Suppose r is a real number that is not an integer, and the decimal part of r terminates. Then r^k is not an integer for any positive integer k.

Proof: The number r can be written as a reduced fraction a / 10^m for some positive m. If s = r^k were an integer, then

10^mk s = a^k.

Now the left side of this equation is divisible by 10 but the right side is not, and so s = r^k must not be an integer. QED.

The only thing special about base 10 is that we most easily think in terms of base 10, but you could replace 10 with any other base.

What about repeating decimals, like 1/7 = 0.142857142857…? They’re only repeating decimals in our chosen base. Pick the right base, i.e. 7 in this case, and they terminate. So the theorem above extends to repeating decimals.

[1] Eli Leher. √2 is Not 1.41421356237 or Anything of the Sort. The American Mathematical Monthly, Vol. 125, No. 4 (APRIL 2018), page 346.