Math

Lunar period approximations

Posted on by John

The date of Easter

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.

The gap between Eastern and Western Easter

Posted on by John

Today is Orthodox Easter. Western churches celebrated Easter last week. Why are the Eastern and Western dates of Easter different? Is Eastern Easter always later than Western Easter? How far apart can the two dates be?

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 4. 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

[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.

Distribution of digits in fractions

Posted on by John

There’s a lot of mathematics just off the beaten path. You can spend a career in math and yet not know all there is to know about even the most basic areas of math. For example, this post will demonstrate something you may not have seen about decimal forms of fractions.

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

[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

Random hexagon fractal

Posted on by John

I recently ran across a post on X describing a process for creating a random fractal. First, pick a random point c inside a hexagon.

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

Root prime gap

Posted on by John

I recently found out about Andrica’s conjecture: the square roots of consecutive primes are less than 1 apart.

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

√pn+1 − √pn < 1.

This has been empirically verified for primes up to 2 × 1019.

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

pn+1 < 1 + 2√pn  + pn,

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

pn+1 < 2pn.

 

Kalman and Bayes average grades

Posted on by John

This post will look at the problem of updating an average grade as a very simple special case of Bayesian statistics and of Kalman filtering.

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

m = (x1 + x2 + x3 + … + xn) / n.

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

m′ = (x1 + x2 + x3 + … + xn + xn+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 + xn+1) / (n + 1).

You could split that into

m′ = w1 m + w2 xn+1

where w1 = n/(n + 1) and w2 = 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 xn+1. [2]

You could also rewrite the equation above as

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

where K = 1/(n + 1) and Δ = xn+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

[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.

Hyperbolic version of Napier’s mnemonic

Posted on by John

I was looking through an old geometry book [1] and saw a hyperbolic analog of Napier’s mnemonic for spherical trigonometry. In hindsight of course there’s a hyperbolic analog: there’s a hyperbolic analog of everything. But I was surprised because I’d never thought of this before. I suppose the spherical version is famous because of its practical use in navigational calculations, while the hyperbolic analog is of more theoretical interest.

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.

Pentagonal numbers are truncated triangular numbers

Posted on by John

Pentagonal numbers are truncated triangular numbers. You can take the diagram that illustrates the nth pentagonal number and warp it into the base of the image that illustrates the (2n − 1)st triangular number. If you added a diagram for the (n − 1)st triangular number to the bottom of the image on the right, you’d have a diagram for the (2n − 1)st triangular number.

In short,

Pn = T2n − 1 − Tn.

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

The proof follows immediately from the definition of pentagonal numbers

Pn = (3n² − n)/2

and triangular numbers

Tn = (n² − n)/2.

Related posts

Computing sine and cosine of complex arguments with only real functions

Posted on by John

Suppose you have a calculator or math library that only handles real arguments but you need to evaluate sin(3 + 4i). What do you do?

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.

\begin{align*} \sin(x + iy) &= \sin x \cosh y + i \cos x \sinh y \\ \cos(x + iy) &= \cos x \cosh y - i \sin x \sinh y \\ \end{align*}

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

\begin{align*} \sin iz &= i \sinh z \\ \cos iz &= \cosh z \end{align*}

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

Lebesgue constants

Posted on by John

I alluded to Lebesgue constants in the previous post without giving them a name. There I said that the bound on order n interpolation error has the form

ch^{n+1} + \lambda \delta

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

a = x_0 < x_1 < x_2 < \cdots < x_{n-1} < x_n = b

define

\ell_j(x) := \prod_{\begin{smallmatrix}i=0\\ j\neq i\end{smallmatrix}}^{n} \frac{x-x_i}{x_j-x_i}

Then the Lebesgue function is defined by

\lambda_n(x) = \sum_{j=0}^n |\ell_j(x)|

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

\Lambda_n(T)=\max_{x\in[a,b]} \lambda_n(x)

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

\Lambda_n \sim \frac{2^{n+1}}{n \log n}

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

\Lambda_n \approx \frac{2}{\pi} \log(n + 1) + 1

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.