Prime divisors of a Babylonian sequence

Question

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. At present the computational evidence supports the strengthened but unproven:

2'. For $k\ge 1$ all prime divisors of $S_k$ are congruent to $1 \bmod 2^{k+3}$.

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$. But this has not been the case.

The initial entry $S_0=3$ and the next three are primes:

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

Indeed those latter prime entries $S_k$ are congruent to $1 \bmod 2^{k+3}$.

The next five 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 \\ S_8 &= \mathbf{C}196 = 79873 \cdot 719073\,884161 \cdot \mathbf{P}66 \cdot \mathbf{P}114 \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.

So far the prime factors of these $S_k$ all satisfy the strengthened Conjecture $2'$ of congruence to $1 \bmod 2^{k+3}$.

This is increasingly 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.

Answer 1

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 8 & \text{no} & \begin{array}{c} (2^{11}\cdot 3\cdot 13 + 1)(2^{11}\cdot 3^4\cdot 5\cdot 23\cdot 37693 + 1) \\ (2^{11}\cdot 3\cdot 5\cdot 31 \\ \cdot 146895790987891709799980308411 454666755530685063771873467473 + 1) \\ (2^{11}\cdot 3\cdot 53 \cdot 97 \cdot 197 \\ \cdot 438921501503783941728573362042998675569676958777 \\ \cdot 700769597973643172674517472532110991945276879083466747+ 1)\end{array} \\ \hline \end{array} $$

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