Finding Integer Approximations

The Saros is the time period for the draconic month ($$T_d$$ = 27.212220815 days), synodic month ($$T_s$$ = 29.530588861 days) and anomalistic month ($$T_a$$ = 27.554549886 days) to approximately match. More specifically, it states that $$242 T_d \approx 223T_s \approx 239 T_a$$ My question is: Is there a non brute force way to derive such integer approximations?

Consider a simpler case, in which we want to get integer approximations for only 2 periods. Take the synodic month and the tropical year ($$T_y$$ = 365.24219 days). Therefore, $$aT_s=bT_y$$ $$\frac{a}{b}=\frac{T_y}{T_s}$$ To derive increasingly accurate approximations, one can use the continued fraction expansion for $$\frac{T_y}{T_s}=12.368266400619...$$, which is (in the WolframAlpha notation) $$12.368266400619=[12; 2, 1, 2, 1, 1, 17, 3, 196, 1, 4, 1, 1, 2, 2,...]$$ From which it can be seen that truncating just before the 17 yields a good approximation $$12.368266400619\approx[12; 2, 1, 2, 1, 1]=\frac{235}{19}$$ Therefore, $$235T_s\approx19T_y$$ Which is known as the metonic cycle.

However, I cannot see a good way to extend this method for 3 periods and beyond.

• BTW, questions related to Saros series are welcome on Astronomy.SE Apr 3 at 4:08

You could try the Lenstra–Lenstra–Lovász (LLL) lattice basis reduction algorithm to calculate a short vector in a lattice. I use Pari GP (p. 382, section 3.10.62 qflll) for the calculation.

We define a lattice with a basis consisting of the columns of the Matrix $$M$$

$$M:=\begin{pmatrix} \frac 1 m& 0& 0 \\0&\frac 1 m & 0 \\0 & 0& \frac 1 m \\T_d &T_s & 0\\ T_d&0&T_a\end{pmatrix}$$ $$m$$ is a scaling factor. A linear combination of these basis vectors

$$\begin{pmatrix} \frac {\lambda_1} m \\ \frac {\lambda_2} m \\ \frac {\lambda_3} m \\{\lambda_1} T_d+ {\lambda_2} T_s \\{\lambda_1} T_d+ {\lambda_3} T_a \end{pmatrix}$$

is of small length, if $${\lambda_1} T_d+ {\lambda_2} T_s$$ and $${\lambda_1} T_d+ {\lambda_3} T_a$$ is small and $$\frac {\lambda_1} m, \frac {\lambda_2} m , \frac {\lambda_3} m$$ are small, too. In this case $${\lambda_2} T_s - {\lambda_3} T_a= ({\lambda_1} T_d+ {\lambda_2} T_s )-({\lambda_1} T_d+ {\lambda_3} T_a)$$ is small, too. Note that $${\lambda_2}$$ and $${\lambda_3}$$ have the same sign.

With Pari/GP online I get the following output

\p20
Td= 27.212220815;
Ts = 29.530588861;
Ta = 27.554549886;
m=1000
M=[1/m,0,0;0,1/m,0;0,0,1/m;Td,Ts,0;Td,0,Ta]
s=qflll(M)
[s[,1],[s[1,1]*Td+s[2,1]*Ts, s[1,1]*Td+s[3,1]*Ta, s[2,1]*Ts-s[3,1]*Ta]]

%2 = [[242, -223, -239]~, [0.036121227000000000000, -0.17998552400000000000, 0.21610675100000000000]]


The coefficients are $$242, 223, 239$$, the differences are $$0.036, 0.180, 0.216$$.

I added the length of the tropical month ($$T_t$$) and the sidereal month ($$Ti$$), from Wiki and got, for m=100000

? \p20
Td= 27.212220815;
Ts = 29.530588861;
Ta = 27.554549886;
Tt=27.321582252
Ti=27.321661554
m=100000
M=[1/m,0,0,0,0;0,1/m,0,0,0;0,0,1/m,0,0;0,0,0,1/m,0;0,0,0,0,1/m;Td,Ts,0,0,0;Td,0,Ta,0,0;Td,0,0,Tt,0;Td,0,0,0,Ti]
s=qflll(M)
[s[,1],[s[1,1]*Td+s[2,1]*Ts, s[1,1]*Td+s[3,1]*Ta, s[1,1]*Td+s[4,1]*Tt,  s[1,1]*Td+s[5,1]*Ti]]

%27 = [[4993, -4601, -4931, -4973, -4973]~, [0.37917983400000000000, -0.86695857100000000000, 0.38999009900000000000, -0.0043787470000000000000]]

• @PM2Ring yes, both months differ only by some seconds. Mar 31 at 18:45
• Here's a live version, running on the SageMathCell server. Apr 1 at 7:53
• @ordptt Here's a Sage / Python version. sagecell.sagemath.org/… Apr 1 at 20:43
• @ordptt Sure, powers of ten are fine. There's no point going too high, since the various month lengths are changing over the long term. Also, the anomalistic month is less important than the other two, it governs the relative sizes of the Sun & Moon, but you need the synodic & draconitic cycles to line up to actually get an eclipse. BTW, in my code, you can type Sage / Python expressions, including stuff like sqrt(2) or pi, in the data box. Apr 1 at 20:53
• @ordptt I think powers of 10 will make sense. Increasing m only by adding 1 will not change anything most of the time. The usage of $T_d$ differs from the usage of $T_s$ and $T_a$ in the algorithm, because $T_d$ is compared to both other numbers directly. $T_a$ and $T_s$ are not compared directly. So you can change the role of $T_d$ and $T_a$ or the role of $T_d$ and $T_s$ and may get different results. Apr 2 at 9:35

Invert the problem, with $$F_a=1/T_a$$ etc. Then we are looking for integers $$a,b,c$$ that satisfy $$\frac{F_a}{a}\approx\frac{F_b}{b}\approx\frac{F_c}{c}\approx \delta.$$ That is, we want find some grid spacing $$δ$$ so that all $$F_{a,b,c}$$ are very close to grid points. The new interpretation of the task is thus to $$\text{ minimize }\|\vec F-δ\vec a\|_2^2\text{ where }\vec a\in\Bbb Z^3,$$ with $$\vec a$$ minimal, that is, co-prime as vector.

This can be done by some kind of Euclidean algorithm, similar to how the continued fraction computation is related to it. In the field of motivations is also the L3 algorithm.

Iterate

• Take the two largest elements and reduce the largest by the second-largest using integer quotients.
• Keep track of the transformation matrix $$M$$, here with the convention that if $$\vec F'$$ is the reduced vector at some stage, then $$M^T\vec F'=\vec F$$ is equal to the original vector.
• Then whenever the (new) maximal element of $$\vec F'$$ at position $$k$$ is decidedly larger than the remaining elements, the vector $$\vec a=M_k$$, the $$k$$th row of $$M$$ is a candidate for the integer coefficient vector.

As error measure use the residual of the minimization task, or closely related, a measure of how parallel $$\vec F$$ and $$\vec a$$ are.

In the following python script these ideas are realized. A report is printed for the case that the residual is substantially improved.

T = [ 27.212220815, 29.530588861, 27.554549886 ]

def largest(F):
k0=0; k1=1
if F[0]<F[1]: k0=1; k1=0
for k in range(2,len(F)):
if F[k]>F[k0]: k1=k0; k0=k; continue
if F[k]>F[k1]: k1=k
return k0,k1

def defect(a):
n2_F=n2_a=0
s_Fa=0
for ak, Tk in zip(a,T): n2_F+=1/Tk**2; n2_a+=ak**2; s_Fa+=ak/Tk
delta = s_Fa/n2_a
return (sum((1/Tk-delta*ak)**2 for ak, Tk in zip(a,T))/n2_F)**0.5

F = [ 1/Ta for Ta in T ]
tol=1e-1
max_F = max(F)
M = np.eye(len(F), dtype=int)
while True:
k0,k1 = largest(F)
err = defect(M[k0])
if err < tol:
tol=0.2*err
print("  a: ",list(M[k0]))
print("a*T: ",[a*Ta for a,Ta in zip(M[k0],T)])
print(" ^F: ",[Fk/F[k0] for Fk in F])
if err < 5e-10: break
q = int(round(F[k0]/F[k1]))
F[k0] -= q*F[k1]
M[k1] += q*M[k0]
if F[k0]<0: F[k0]=-F[k0]; M[k0]=-M[k0]


This gives the log

  a:  [1, 1, 1]
a*T:  [27.212220815, 29.530588861, 27.554549886]
^F:  [0.013482132197497906, 1.0, 0.07171370910340992]
angle (deg): 2.035°

a:  [15, 14, 15]
a*T:  [408.18331222499995, 413.42824405399995, 413.31824829]
^F:  [0.18799937091605323, 0.055664774527038254, 1.0]
angle (deg): 0.34535°

a:  [76, 70, 75]
a*T:  [2068.12878194, 2067.14122027, 2066.59124145]
^F:  [1.0, 0.29609021698213056, 0.31916673511916593]
angle (deg): 0.018052°

a:  [242, 223, 239]
a*T:  [6585.35743723, 6585.321316003, 6585.537422754]
^F:  [0.14353663918949092, 1.0, 0.0779374555912061]
angle (deg): 0.00082299°

a:  [3783, 3486, 3736]
a*T:  [102943.831343145, 102943.63276944599, 102943.798374096]
^F:  [0.15830991529461017, 0.06102937657324801, 1.0]
angle (deg): 4.7152e-05°

a:  [20928, 19285, 20668]
a*T:  [569497.35721632, 569497.406184385, 569497.437043848]
^F:  [1.0, 0.3855057117532664, 0.3167237386175655]
angle (deg): 3.3867e-06°

a:  [66325, 61118, 65501]
a*T:  [1804850.545554875, 1804850.530006598, 1804850.572082886]
^F:  [0.49417557377594934, 0.21716709153510322, 1.0]
angle (deg): 5.4577e-07°

a:  [285986, 263534, 282433]
a*T:  [7782314.18199859, 7782314.204894774, 7782314.187952638]
^F:  [0.2755545984556782, 1.0, 0.05364004447339751]
angle (deg): 6.9819e-08°

a:  [1255666, 1157087, 1240066]
a*T:  [34169460.46188779, 34169460.4734079, 34169460.458932474]
^F:  [1.0, 0.3709551370058305, 0.19466212784696205]
angle (deg): 1.0192e-08°


where indeed the reference result is in the 4th record.

For $$3$$ periods

Divide $$aT_1 = bT_2 = cT_3$$ into $$2$$ simpler cases: $$aT_1 = bT_2$$ and $$bT_2 = cT_3$$.

It follows that $$\frac{a}{b} = \frac{T_2}{T_1} \; \text{and} \; \frac{b}{c} = \frac{T_3}{T_2}$$

The appropriate continuation can be best demonstrated with an example.

Consider the case provided where $$T_1 = 27.212220815, T_2 = 29.530588861$$ and $$T_3 = 27.554549886$$.

We find appropriate approximations for $$\frac{a}{b}$$ and $$\frac{b}{c}$$ using the continued fraction method described in the question, getting

$$\frac{a}{b} \approx [1; 11, 1, 2, 1, 4, 3, 5, 1, 31, 1, 1, 6] = \frac{2381497}{2194532}$$ $$\frac{b}{c} \approx [0; 1, 13, 1, 16, 1, 27, 3, 5, 1, 1] = \frac{879525}{942599}$$

We can now multiply $$\frac{a}{b}$$ appropriately to make the denominator in $$\frac{a}{b}$$ equal to the numerator in $$\frac{b}{c}$$ as follows: $$\frac{a}{b} \times \frac{2194532}{879525} \approx \frac{2381497}{2194532} \times \frac{2194532}{879525} = \frac{2381497}{879525}$$ $$\implies \frac{\frac{2194532}{879525}a}{b} \approx \frac{2381497}{879525}$$

We can now equate $$\frac{2194532}{879525}a = 2381497$$ $$\implies a = \frac{2094586148925}{2194532}$$

Hence $$a : b : c : : \frac{2094586148925}{2194532} : 879525 : 942599$$ $$\implies a : b : c : : 2094586148925 : 1930145757300 : 2068563668668$$

Therefore, $$2094586148925T_1 \approx 1930145757300T_2 \approx 2068563668668T_3$$.

Obviously, these numbers are larger than $$242, 223$$ and $$239$$ but that is because I chose to truncate the continued fraction beyond a certain point. If you like, you can also take the values found for $$a, b, c$$ and find smaller integer approximations for $$a : b : c$$.

For $$n$$ periods

The above argument can be easily generalised as follows:

Consider $$a_1T_1 = a_2T_2 = a_3T_3 = a_4T_4 = ... = a_nT_n$$.

Simply take $$a_1T_1 = a_2T_2, \, a_2T_2 = a_3T_3, \, a_3T_3 = a_4T_4, ... , a_{n-1}T_{n-1} = a_{n}T_{n}$$, getting

$$\frac{a_1}{a_2} = \frac{T_2}{T_1}, \, \frac{a_2}{a_3} = \frac{T_3}{T_2}, ..., \frac{a_{n-1}}{a_n} = \frac{T_n}{T_{n-1}}$$

Then multiply each of the fractions appropriately so that the necessary numerators and denominators match. After that, find the value for $$a_1$$ by equating the numerators. Next, consider the ratio $$a_{1} : a_{2} : a_{3} : ... : a_{n}$$, multiply to make everything an integer and voila, your answer is ready!

• That is a good way to solve the problem. However, it is very far from being the best. The 3 periods case, for example, is equivalent to let $b$ and $c$ equal the numerator and denominator of the continued fraction approximation of $b/c$ and solving for $a$. However, that only guarantees that $bT_2=cT_3$ is a good approximation, but not overall. In fact, by using your method with approximations close to the size of what I presented in the question results in $a,b,c=272,251,269$, which is around the 2000th place for the first 2508 best approximations (done by brute force for $a=1,2,...,2508$) Mar 27 at 23:43
• Just to detail a bit more: There should be a balance. For some approximated $a,b,c$, if $bT_2=cT_3$ is a really good approximation, but $aT_1=bT_2$ or $aT_1=cT_3$ are horrible approximations, the overall relation $aT_1=bT_2=cT_3$ is a bad approximation. Another way to see this is by noting that the equation is equivalent to the following two: $bT_2-aT_1=0$ and $cT_3-aT_1=0$, such that the expression $(bT_2-aT_1)^2+(cT_3-aT_1)^2$ is minimized for good approximations. Mar 27 at 23:50
• @ordptt Yes, you are correct. Unfortunately, this is the best I could come up with. I'll be looking forward to seeing better approaches by more experienced mathematicians.
– user905694
Mar 28 at 4:57
• No problem! It sure is a difficult problem that doesn't seem to have a straightforward answer. Mar 28 at 5:27