Pope Gregory, the Calendar,

and Continued Fractions

Friedrich Eisenbrand

Abstract. The success of many activities of modern civilization crucially depends on careful planning. Some activities should be carried out during a certain period of the year. For example: When is the right time of the year to sow, when is the right time to plow? It is thus no surprise that calendars are found in literally every ancient civilization. The earth revolves around the sun in about 365.2422 days. An accurate calendar can thus not provision the same number of days every year if the calendar should be synchronous with the seasons. This article is about the problem of approximating a given number by a rational number with small denominator, continued fractions and their relationship to the Gregorian calendar with its leap-year rule that is still in use today and keeps the calendar synchronized for a very long time.

2010 Mathematics Subject Classification: 11J70, 11Y16, 11A55 Keywords and Phrases: Calendar, Diophantine approximation, continued fractions

The Julian calendar and Gregory’s reform

The number 365.2422 is close to 365 + 1/4. If this was precisely the duration of one year in days, then the following rule would result in an exact calendar.

Each year that is divisible by 4 consists of 366 days and each other year consists of 365 days.

The mean duration of a calendar year is thus 365 + 1/4. In other words, each year that is divisible by 4 will be a leap year. This leap year rule was imposed by Julius Cesar in 45 B.C. Already at this time, astronomers calculated the duration of a year in days fairly accurately and it was clear that the calendar would be behind by one day in roughly 130 years.

In 1582, when the Julian calendar was evidently out of sync by a large extent, pope Gregory the XIII imposed the following calendar reform. As before, every year that is divisible by 4 is a leap-year, except for those divisible by 100 but not by 400. The mean duration of a year of the Gregorian calendar is thus 365 + 97/400.

Best approximations

What is the mathematical challenge behind the design of an accurate leap-year rule? The task is to approximate the number 0.2422 by a rational number p/q with p, q ∈ ℕ_{+} such that q as well as the error E = |.2422 − p/q| is small. The mean duration of a calendar year is then 365 + p/q if the calendar provisions p leap years every q years. The smaller the q, the simpler should be the leap-year rule. In the Julian calendar, p/q = 1/4. The rule “Each year divisible by four is a leap year” is easy to remember. In 1/E years, the calendar will then be ahead by one day or behind by one day depending on whether p/q is smaller or larger than 0.2422.

Finding a convenient and sufficiently accurate leap-year rule is related to approximating a real number α ∈ ℝ_{≥0} by a rational number p/q in a good way. In the following we always assume that p is a natural number or 0 and that q is a positive natural number when we speak about the representation p/q of a rational number. The rational number p/q is a best approximation of α if for any other rational number p′/q′ ≠ p/q one has

|α − p/q| < |α − p′/q′|

if q′ ≤ q. Going back to the calendar problem, this makes sense. If there exists an approximation p′/q′ of 0.2422 with q′ ≤ q that results in a smaller error, then we could hope that we can find a leap year rule that accommodates for p′ leap years in q′ years instead of the one that accommodates for p leap years in q years that is just as easy to remember. Furthermore, the calendar would be more accurate.

Continued fractions

Continued fractions have been used to approximate numbers for a very long time and it seems impossible to attribute their first use to a particular researcher or even to a particular ancient civilization. Keeping the best approximation problem in mind however, the application of continued fractions seems natural.

Suppose our task is to approximate α ∈ ℝ_{≥0} by a rational number with small denominator. If α is not a natural number then we can re-write

α = ⌊α⌋ + (α − ⌊α⌋)

= ⌊α⌋ + \frac{1}{1/(α − ⌊α⌋)}.

Pope Gregory, the Calendar, and Continued Fractions 89

The number β = 1/(α − ⌊α⌋) is larger than one. If β is not a natural number, one continues to expand the number β and obtains

α = ⌊α⌋ + \frac{1}{⌊β⌋ + \frac{1}{1/(β − ⌊β⌋)}}.

The continued fraction expansion of α is inductively defined as the sequence α if α ∈ ℕ and ⌊α⌋, a_{1}, a_{2}, … otherwise, where a_{1}, a_{2}, … is the continued fraction expansion of 1/(α − ⌊α⌋). On the other hand, a finite sequence of integers b_{0}, …, b_{n}, all positive, except perhaps b_{0} gives rise to the continued fraction

⟨b_{0}, …, b_{n}⟩ = b_{0} + \frac{1}{b_{1} + \frac{1}{... + \frac{1}{b_{n}}}}.

If the sequence a_{0}, a_{1}, … is the continued fraction expansion of α ∈ ℝ_{≥0} and if its length is at least k + 1, then the k-th convergent of α is the continued fraction

⟨a_{0}, …, a_{k}⟩ = a_{0} + \frac{1}{a_{1} + \frac{1}{... + \frac{1}{a_{k}}}}

Let us compute the first convergents of the number α = 365.2422. Clearly, a_{0} is 365. To continue, it is convenient to represent α as a rational number α = 1826211/5000. Clearly α − ⌊α⌋ is the remainder of the division of 1826211 by 5000 divided by 5000. One has

1826211 = 5000 · 365 + 1211.

Thus we continue to expand 5000/1211 and obtain a_{1} = 4. The remainder of the division of 5000 by 1211 is 156 which means that we next expand 1211/156 which results in a_{2} = 7. The remainder of this division is 119 and we next expand 156/119 resulting in a_{3} = 1, then 119/37 yielding a_{4} = 3 and 37/8 yields a_{5} = 4.

At this point we can record an important observation. If α = p/q is a rational number, then its continued fraction expansion is precisely the sequence of quotients of the division-with-remainder steps that are carried out by the Euclidean algorithm on input p and q. Also, for arbitrary real α ∈ ℝ_{≥0}, the function f_{k}(x) = ⟨a_{0}, …, a_{k−1}, x⟩ defined for x > 0 is strictly increasing in x if k is even and decreasing if k is odd. Furthermore, if k is even, then a_{k} is the largest integer with ⟨a_{0}, …, a_{k}⟩ ≤ α and if k is odd then a_{k} is the largest integer such that ⟨a_{0}, …, a_{k}⟩ ≥ α.

The quality of the Gregorian calendar

The third convergent of 365.2422 is

365 + \frac{1}{4 + \frac{1}{7 + \frac{1}{1}}} = 365 + 8/33.

According to Rickey [6], the Persian mathematician, philosopher and poet Omar Khayyam (1048 - 1131) suggested a 33-year cycle where the years 4, 8, 12, 16, 20, 24, 28 and 33 should be leap years. Thus the mean-duration of a year according to his suggestion would be exactly the value of the third convergent. How does this compare to the mean duration of a year of the Gregorian calendar. We calculate both error terms

E_{1} = |365.2422 − 365 + 8/33| = 0.000224242424242432
E_{2} = |365.2422 − 365 + 97/400| = 0.000299999999999995

and surprisingly, one finds that Omar Khayyam’s leap-year rule is more accurate. Using the third convergent, his calendar will be imprecise by one day in roughly 4459.45 years, whereas Gregory’s calendar will be off by one day in “only” 3333.33 years. Still the leap-year rule of the Gregorian calendar is convenient, as it relates nicely with our decimal number system and is simple to remember. However, why is it a good idea to approximate a number by its convergent? What is the relation of the convergents of a number with its best approximations?

Best approximations and convergents

We now explain the relationship of convergents of α ∈ ℝ_{≥0} and best approximations. The subject is nicely treated in [2]. Let a_{0}, a_{1}, … be a sequence of natural numbers where again all are positive except perhaps a_{0} and consider the two sequences g_{k} and h_{k} that are inductively defined as

(\begin{matrix}g_{−1} & g_{−2} \\ h_{−1} & h_{−2} \end{matrix}) = (\begin{matrix}1 & 0 \\ 0 & 1 \end{matrix}) , (\begin{matrix}g_{k} & g_{k−1} \\ h_{k} & h_{k−1} \end{matrix}) = (\begin{matrix}g_{k−1} & g_{k−2} \\ h_{k−1} & h_{k−2} \end{matrix}) (\begin{matrix}a_{k} & 1 \\ 1 & 0 \end{matrix}) , k ≥ 0.

(1)

It follows from a simple inductive argument that, if β_{k} is the number β_{k} = g_{k}/h_{k}, then one has ⟨a_{0}, …, a_{k}⟩ = β_{k} for k ≥ 0.

Now the process of forming convergents admits a nice geometric interpretation. Notice that, since the a_{i} are integers and since the determinant of

(\begin{matrix}g_{k} & g_{k−1} \\ h_{k} & h_{k−1} \end{matrix})  (2)

\begin{matrix}v_{k−1} \\ v_{k} \\ h = α · g \\ g \\ h \\ 0 \\ h_{k−1} \end{matrix}

Figure 1: An illustration of the geometric interpretation of convergents

is 1, such a matrix (2) is a basis of the standard lattice ℤ^{2}. This means that each vector in ℤ^{2} can be obtained by multiplying the matrix (2) with an integral 2-dimensional vector and conversely, the result of such a multiplication is always an integral 2-dimensional vector. If v_{k} = ( \begin{matrix}g_{k} \\ h_{k} \end{matrix} ) then the line with slope α through 0 is sandwiched between the vectors v_{k} and v_{k−1} in the positive orthant, see Figure 1. In Figure 1, the rational number g_{k−1}/h_{k−1} is larger than α. Since there is no integer point in the shaded region, any other rational number p/q ≥ α with p/q − α ≤ g_{k−1}/h_{k−1} − α must have a denominator that is larger than h_{k−1}. One says that g_{k−1}/h_{k−1} is a best approximation from above. Similarly, g_{k}/h_{k} is a best approximation from below. At this point it is already clear that one of the convergents is a best approximation.

Next we show that the following best approximation problem can be solved in polynomial time.

Given a rational number α ∈ ℚ_{>0} and a positive integer M, compute the best approximation of α with denominator bounded by M, i.e., compute a rational number p/q with p ≤ M such that |α−p/q| is minimum.

The algorithm is described in [2], see also [1], and is as follows. One computes the convergents α as long as the denominator (h-component) of the latest convergent is bounded by M. Since the denominators double every second round, the number of steps is bounded by the encoding length of M. Suppose that this is the k-th convergent and we denote the columns of the matrix (2) again by v_{k} and v_{k−1}. In the next round, the new first column would be v_{k−1} + a_{k+1} · v_{k} but the h-component of this vector exceeds M. Instead, one computes now the largest μ ∈ ℕ_{0} such that the h-component of v_{k−1} + μ · v_{k} does not exceed M. If we denote the resulting vector by u then still u, v_{k} is a basis of ℤ^{2} but the second component of u + v_{k} exceeds M. The situation is depicted in Figure 2. Any rational number p/q that approximates α better

\begin{matrix}u \\ v_{k} \\ h ≤ M \\ h = α · g \\ 0 \\ g \\ h \end{matrix}

Figure 2: An illustration of the algorithm solving the best approximation problem

than u and v_{k} is in the cone C spanned by u and v_{k}

C = {λ_{1}u + λ_{2}v_{k}: λ_{1}, λ_{2} ≥ 0}.

But if this rational number is different from the one represented by u and v, then λ_{1} and λ_{2} must be strictly positive. However, since u and v_{k} form a lattice-basis, λ_{1} and λ_{2} are positive integers and thus the h-component q of the corresponding vector exceeds M. Thus u or v_{k} is a solution to the bestapproximation problem.

Further historical remarks

Continued fractions are a true classic in mathematics and it is impossible to give a thorough historical account. In this final section I content myself with a very brief discussion of computational issues related to best approximations and continued fractions and some recent results. The simultaneous best approximation problem is the high-dimensional counterpart to the best approximation problem that we discussed. Here, one is given a rational vector and a denominator bound and the task is to find another rational vector where each component has the same denominator that is bounded by the prescribed denominator bound. The objective is to minimize the error in the ℓ_{∞}-norm. Lagarias [3] has shown that this problem is NP-hard and applied the LLL-algorithm [4] to approximate this optimization problem. Variants of this simultaneous best approximation problem are also shown to be hard to approximate [7]. Schönhage [8] showed how to compute convergents in a quasilinear amount of bit-operations. Recently Novocin, Stehlé and Villard [5] have shown that a variant of LLL-reduction depends on the bit-size of the largest input coefficient in a similar way.

References

[1] M. Grötschel, L. Lovász, and A. Schrijver. Geometric Algorithms and Combinatorial Optimization, volume 2 of Algorithms and Combinatorics. Springer, 1988.

[2] A. Ya. Khinchin. Continued fractions. Dover Publications Inc., Mineola, NY, russian edition, 1997. Reprint of the 1964 translation. The first Russion edition was published in 1935.

[3] J. C. Lagarias. The computational complexity of simultaneous diophantine approximation problems. SIAM J. Computing, 14(1):196–209, 1985.

[4] A. K. Lenstra, H. W. Lenstra, and L. Lovász. Factoring polynomials with rational coefficients. Math. Annalen, 261:515–534, 1982.

[5] Andrew Novocin, Damien Stehlé, and Gilles Villard. An lll-reduction algorithm with quasi-linear time complexity: extended abstract. In Lance Fortnow and Salil P. Vadhan, editors, STOC, pages 403–412. ACM, 2011.

[6] V. Frederick Rickey. Mathematics of the Gregorian calendar. The Mathematical Intelligencer, 7(1):53–56, 1985.

[7] Carsten Rössner and Jean-Pierre Seifert. Approximating good simultaneous Diophantine approximations is almost NP-hard. In Mathematical foundations of computer science 1996 (Cracow), volume 1113 of Lecture Notes in Comput. Sci., pages 494–505. Springer, Berlin, 1996.

[8] A. Schönhage. Schnelle Berechnung von Kettenbruchentwicklungen. Acta Informatica, 1:139–144, 1971.

Friedrich Eisenbrand

EPFL SB IMA MA C1 573 1015 Lausanne Switzerland

                  friedrich.eisenbrand@epfl.ch

94

Documenta Mathematica · Extra Volume ISMP (2012)