# Prime divisors of a Babylonian sequence

Consider the odd natural number sequence (see OEIS A001601): $$S_0 = 3 \\ S_{k+1} = 2S_k^2 - 1$$

These are numerators of reduced rational approximations to $$\sqrt 2$$ generated by the usual Newton iterations starting at $$1$$. Knuth's paper "Ancient Babylonian Algorithms" (Comm. of the ACM, 1972) sets out a case that cuneiform tablets from "about 1800 - 1600 B.C." depict a sexagesimal floating-point implementation of this method for extracting square roots.

$$1, \frac{1+\frac{2}{1}}{2} = \frac{3}{2}, \frac{\frac{3}{2}+\frac{2}{3/2}}{2} = \frac{17}{12}, \; \ldots$$

Two conjectures about the prime divisors of these numerators arise from this previous Question:
1. Each $$S_k$$ is square-free.
2. For $$k\ge 1$$ all prime divisors of $$S_k$$ are congruent to $$1 \bmod 4$$.

A proof of these conjectures will immediately imply two similar conjectures posed in that previous Question, which references a "modified Lucas-Lehmer sequence" of even numbers. (See also OEIS A003423.) [NB: Somos suggested that Conjecture 2 is implied by Conjecture 1. In discussing the circle of ideas to verify this, an outright proof of Conjecture 2 was found. It is briefly sketched in the CW Answer appended below.]

Low hanging fruit

Because a great deal is known about prime divisors of Fibonacci numbers, my initial expectation was that there would be considerable research relevant to the conjectures posed above. For the first conjecture, square-freeness of the entries, the only related material I found was this previous MSE Question, involving a different initial condition:

Search square factors in Lucas–Lehmer sequence

There David Speyer supplies a heuristic estimate that the number of terms that are not square-free is finite and "pretty small".

Regarding the second conjecture, there are two easy partial observations. Certainly $$S_k$$ itself is $$1 \bmod 4$$ for each $$k\ge 1$$ as the induction argument shows (basis $$S_1 = 17$$):

$$S_{k+1} \equiv 2 S_k^2 - 1 \equiv 2\cdot 1 - 1 \equiv 1 \bmod 4$$

A similar argument shows that the terms $$S_k$$ are pairwise coprime. Let $$p$$ be a prime that divides $$S_k$$, i.e. $$S_k \equiv 0 \bmod p$$. Then:

$$S_{k+1} \equiv 2 S_k^2 - 1 \equiv -1 \bmod p$$

and thereafter (by induction for $$j \gt 1$$):

$$S_{k+j} \equiv 2 - 1 \equiv 1 \bmod p$$

Evidently $$p$$ will not divide any entry subsequent to $$S_k$$, so those entries are coprime to $$S_k$$. (We didn't use primality of $$p$$, only that $$p\gt 1$$, but it was enough to establish coprimality of the $$S_k$$.)

Computational evidence

Faced with a paucity of theoretical results, we can always hope that just factoring the entries $$S_k$$ for $$k\ge 1$$ may quickly give a negative result, either by turning up a square factor or a prime factor congruent to $$3 \bmod 4$$.

The initial entry $$S_0=3$$ and the next three are primes congruent to $$1 \bmod 4$$:

\begin{aligned} S_1 &= 17 \\ S_2 &= 577 \\ S_3 &= 665857 \end{aligned}

By the earlier proof (or inspection!) these entries are congruent to $$1 \bmod 4$$.

The next four entries (spaced into six digit blocks for easier reading) are composite and have been fully factored:

\begin{aligned} S_4 &= 886731\,088897 \\ &= 257 \cdot 1409 \cdot 2\,448769 \\ S_5 &= 1\,572584\,048032\,918633\,353217 \\ &= 11777 \cdot 2\,393857 \cdot 55\,780318\,173953 \\ S_6 &= 4\,946041\,176255\,201878\,775086\,487573\,351061\,418968\,498177 \\ &= 7681 \cdot 1\,492993 \cdot 431\,302713\,980890\,947612\,633357\,964569\,696769 \\ S_7 &= \mathbf{C}98 = 7460\,967982\,148609 \cdot \mathbf{P}82 \end{aligned}

The digit lengths of entries $$S_k$$ increase rapidly, at least doubling with each step. For the last entry above we've used an abbreviation scheme common in factoring large numbers, indicating a composite number with $$\mathbf C$$ and a prime number with $$\mathbf P$$, followed respectively by the number of digits needed. Readers can verify the "small" prime factor shown is congruent to $$1 \bmod 4$$, and so the other prime factor must be also.

For the next entry we have as yet only a partial factorization, with two prime divisors congruent to $$1 \bmod 4$$ and a large composite factor:

$$S_8 = \mathbf{C}196 = 79873 \cdot 719073\,884161 \cdot \mathbf{C}179$$

This is a "hard row to hoe," computationally speaking, but there is an alternative approach that finds all the prime divisors of this Babylonian sequence up to a moderate limit, say less than a million, with considerably less effort. This and other partial results will be collected in a CW Answer below.

• @Somos: It would be nice to connect the two conjectures, but we only need an even number of factors of the form $4k+3$ to get a product of the form $4k+1$, e.g. $7\cdot 11=4 (19)+1$. – hardmath Apr 5 at 20:22
• Ouch! Of course. – Somos Apr 5 at 20:25
• @Somos: But if we knew something more about the entries, your suggestion would be apt. If not only congruent to $1 \bmod 4$, the $S_k$ were a sum of two squares, then any prime factors congruent to $3 \bmod 4$ would have to appear as even powers. Thus Conjecture 1 would imply Conjecture 2, as you wished. – hardmath Apr 6 at 4:39
• See recent OEIS sequence 307372 for sum of squares. – Somos Apr 6 at 4:42
• @Somos: Amazingly it seems your suggestion has busted the problem open. A little scratch work shows that indeed $S_k = T_k^2 + 1$, where $T_1 = 4$ and thereafter, $T_k = 2S_{k-2}T_{k-1}$. So each $S_k$ is a sum of two squares. – hardmath Apr 6 at 5:41

A closed form solution

The terms $$S_k$$ are iterates of the second degree Tchebychev polynomial (of the first kind):

$$f(x) = 2x^2 - 1$$

with initial value $$3$$. It follows that:

$$S_k = \cosh (2^k \cosh^{-1} 3)$$

as outlined in answer to this previous Question or this other Question. Sadly this explicit solution seems to be useless for our purposes.

A product formula for entries

For $$k\ge 1$$ we show by induction:

$$S_{k+1} = \left[ 4^{k+2} \prod_{i=0}^{k-1} S_i^2 \right] + 1$$

Indeed if we adopt the usual convention that an empty product is $$1$$, then this hold for $$k=0$$ as well. The first few terms, $$k=0,1,2$$ , serve as our basis step:

\begin{aligned} S_1 &= 4^2\cdot 1 + 1 = 17 \\ S_2 &= 4^3 \cdot S_0^2 + 1 = 577 \\ S_3 &= 4^4 (S_0 \cdot S_1)^2 + 1 = 665857 \end{aligned}

The induction step can be taken as follows:

\begin{aligned} S_{k+2} &= 2S_{k+1}^2 -1 \\ &= 2\left( 2S_k^2 - 1 \right) S_{k+1} - 1 \\ &= 4 S_k^2 S_{k+1} - 2 S_{k+1} - 1 \\ &= 4 S_k^2 \left(4^{k+2} \prod_{i=0}^{k-1} S_i^2 + 1 \right) - 2\left(2S_k^2 - 1\right) - 1 \\ &= \left[4^{k+3} \prod_{i=0}^k S_i^2\right] + 4S_k^2 - 4S_k^2 + 2 - 1 \\ &= \left[4^{k+3} \prod_{i=0}^k S_i^2\right] + 1 \end{aligned}

Note that the induction hypothesis is used here only once, applied to $$k$$ on our way to establishing the case $$k+1$$. The other identities used are the sequence's defining relations $$S_{k+1} = 2S_k^2 - 1$$.

Coprimality of terms revisited

Using the product formula above gives more precise information about the coprimality of terms in the sequence:

\begin{aligned} S_{k+1} &\equiv -1 \bmod S_k^2 \\ S_{k+1} &\equiv +1 \bmod S_j^2 \text{ for } 0\le j \lt k \end{aligned}

The first of these follows from the defining relation of the sequence, and the second follows from the product formula established above.

Proof of Conjecture 2

The product formula above says that for $$k \ge 1$$, $$S_k$$ is of the form $$n^2 + 1$$. It follows that modulo any prime divisor $$p$$ of such $$S_k$$, $$-1$$ is a quadratic residue. Thus any odd prime divisor $$p$$ will be congruent to $$1 \bmod 4$$. See also this previous Question.

Strengthening Conjecture 2

For $$k\ge 1$$, the prime divisors of $$S_k$$ are actually congruent to $$1\bmod 8$$, a slight strengthening of Conjecture 2. Given that prime $$p$$ divides $$S_k = 2S_k^2 - 1$$, it must be that $$2$$ is a quadratic residue mod $$p$$. Then the second supplement to the law of quadratic reciprocity implies $$p$$ is either congruent to $$\pm 1 \bmod 8$$. Only $$p \equiv 1 \bmod 8$$ is consistent with $$p \equiv 1 \bmod 4$$.

Empirically for prime divisors $$p$$ of $$S_k$$ the top power of $$2$$ that divides $$p-1$$ grows with increasing $$k$$, as the following table summarizes:

$$\begin{array}{|c|c|c|} k & \text{prime?} & S_k \text{ prime factors } p \\ \hline 1 & \text{yes} & 2^4 + 1 \\ \hline 2 & \text{yes} & 2^6 3^2 + 1 \\ \hline 3 & \text{yes} & 2^8 3^2 17^2 + 1 \\ \hline 4 & \text{no} & (2^8+1)(2^7\cdot 11 + 1)(2^7\cdot 3\cdot 7\cdot 911 + 1) \\ \hline 5 & \text{no} & (2^9\cdot 23+1)(2^8 3^2\cdot 1039 + 1)(2^8\cdot 3041\cdot 5147\cdot 13921 + 1) \\ \hline 6 & \text{no} & \begin{array}{c}(2^9\cdot 3\cdot 5 +1)(2^{11}3^6 +1)\\(2^9\cdot 3\cdot 7^3\cdot 13\cdot 263\cdot 2147503\cdot 111497203601498327476613 +1)\end{array} \\ \hline 7 & \text{no} & \begin{array}{c} (2^{10}\cdot 7286101545067 + 1)\\(2^{10}\cdot 11\cdot 197\cdot 14268312047910853080882581\cdot 207118459764708621027403333949253228941663842143239 +1)\end{array} \\ \hline \end{array}$$

This suggests more can be said about congruence of such $$p$$ mod higher powers of two.