Finding Integer Approximations

Question

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.

Answer 1

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]]

Answer 2

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!

Answer 3

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])
        print(f"angle (deg): {np.rad2deg(np.arcsin(err)):.5g}°\n")
    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.