Finding all S-Diophantine quadruples for a fixed set of primes S

Given a finite set of primes S and an m-tuple \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(a_1,\ldots ,a_m)$$\end{document}(a1,…,am) of positive, distinct integers we call the m-tuple S-Diophantine, if for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1\le i < j\le m$$\end{document}1≤i<j≤m the quantity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a_ia_j+1$$\end{document}aiaj+1 has prime divisors coming only from the set S. For a given set S we give a practical algorithm to find all S-Diophantine quadruples, provided that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|S|=3$$\end{document}|S|=3.


Introduction
Given an irreducible polynomial f ∈ Z[X , Y ] and a set A of positive integers we consider the product Π = a,b∈A a =b f (a, b).
It is an interesting question what is the largest prime divisor P(Π ) of Π or alternatively what is the number of prime divisors ω(Π) of Π in terms of the size of |A|. In case of f (x, y) = x + y this question was intensively studied by several authors starting with Erdős and Turán [4] who established that in this case ω(Π) log |A|. It was shown by Erdős, Stewart and Tijdeman [3] that this is essentially best possible.  In the case of f (x, y) = x y + 1 the investigation was started by Sárközy [8]. It was shown by Győry et.al. that ω (Π) log |A|. However it seems that this lower bound is far from being best possible and the situation is not so clear as in the case f (x, y) = x + y.
One recent approach to this problem was picked up by Szalay and Ziegler [11] who studied so-called S-Diophantine m-tuples. Let S be a finite set of primes and (a 1 , . . . , a m ) a m-tuple of positive, distinct integers with a 1 < · · · < a m . Then we call (a 1 , . . . , a m ) S-Diophantine, if for each 1 ≤ i < j ≤ m the quantity a i a j + 1 has prime divisors coming only from the set S. In other words if A = {a 1 , . . . , a m } and Π has only prime divisors coming form S, then (a 1 , . . . , a m ) is called a S-Diophantine m-tuple. Szalay and Ziegler [11] conjectured that if |S| = 2, then no S-Diophantine quadruple exists. This conjecture has been confirmed for several special cases by Luca, Szalay and Ziegler [5,[11][12][13]16]. Even a rather efficient algorithm has been described by Szalay and Ziegler [13] that finds for a given set S = {p, q} of two primes all S-Diophantine quadruples, if there exist any. However, the algorithm is strictly limited to the case that |S| = 2 and does not work (nor any slightly modified version of the algorithm) for |S| ≥ 3.
In the present paper we establish an algorithm that finds for any set S = {p, q, r } of three primes all S-Diophantine quadruples. Unfortunately the algorithm is not so efficient as its counter part for |S| = 2. However, when implemented in Sage [14] it takes on a usual PC (Intel i7-8700) -using a single core -about six and a half hours to verify: However the case S = {2, 3, 5} is in some sense the slowest instance for small primes p, q and r . Thus within four days of runtime we could establish the following result: Theorem 2 Let S = {p, q, r } with 2 ≤ p < q < r ≤ 100. Then all S-Diophantine quadruples are listed in Table 1.
Thus Theorem 2 raises the following problem: Table 1 all S-Diophantine quadruples with |S| = 3?

Problem 1 Are the S-Diophantine quadruples listed in
Although we present the method only in the case that |S| = 3 obvious modifications to the algorithm would provide an algorithm that would work for any set S of primes with |S| ≥ 3. However, already in the case that |S| = 4 the algorithm would take instead of several hours several months of computation time. Thus a systematic search for S-Diophantine quadruples in the case |S| = 4 seems to be not feasible. Therefore we refer from discussing the case |S| > 3 in this paper.
In the next section we remind some of the auxiliary results we need to deal with S-Diophantine quadruples. In Sect. 3 we establish several results that will allow us to find small upper bounds for the relevant unknowns and will also allow us to apply LLL-reduction to reduce these bounds. The LLL-reduction will be discussed in Sect. 4. In Sect. 5 we will give details how to apply the LLL-reduction method to find upper bounds small enough to enumerate all possible S-Diophantine quadruples for a given set S = {p, q, r } of three primes.

Auxiliary results
Assume that (a, b, c, d) is a S-Diophantine quadruple, with S = {p, q, r }. In particular, we assume that a < b < c < d. Moreover, we will assume that p < q < r . Let us write The almost trivial observation that yields the (non-linear) S-unit equation Similar computations lead us to the following system of S-unit equations , , , , , , .
Proof A proof can be found in [11,Lemma 3].
Proof A proof can be found in [12, Lemma 2.2] except for the case that p = 2. However, in the case that p = 2 we slightly change the proof. That is we assume that min{α 1 , α 2 , α 4 } ≥ 2 and obtain which is an obvious contradiction.
In order to obtain a first upper bound in the next section we will apply results on lower bounds for linear forms in (complex and p-adic) logarithms. Therefore let η = 0 be an algebraic number of degree δ and let be the minimal polynomial of η. Then the absolute logarithmic Weil height is defined by We will mainly need the notation of height for rational numbers. Therefore let us remark that if p/q ∈ Q with p, q integers such that gcd( p, q) = 1, then h( p/q) = max{log | p|, log |q|}. With this basic notation we have the following result on lower bounds for linear forms in (complex) logarithms due to Matveev [6].
In our case of main interest where |S| = 3 we can apply results on linear forms in two p-adic logarithms. For a prime p we denote by Q p the field of p-adic numbers with the standard p-adic valuation ν p . As above let η 1 , η 2 be algebraic numbers over Q and we regard them as elements of the field K p = Q p (η 1 , η 2 ). In our case η 1 and η 2 will be rational integers, thus K p = Q p . Similar as in the case of Matveev's theorem above we have to use a modified height. In particular, we write With these notations at hand, let us state the result due to Bugeaud and Laurent [1, Corollary 1]): Lemma 5 Let b 1 , b 2 be positive integers and suppose that η 1 and η 2 are multiplicatively independent algebraic numbers such that ν p (η 1 ) = ν p (η 2 ) = 0. Put . and E := max log E + log log p + 0.4, 10, 10 log p .
Then we have The next lemma is an elementary, but rather useful result due to Pethő and de Weger [7]. For a proof of Lemma 6 we refer to [9,Appendix B].

A first upper bound
As already mentioned in the previous section we consider the case that |S| = 3, i.e. we have S = {p, q, r } with P = max{ p, q, r } and write First, we prove an upper bound for α i , β i and γ i with i = 1, 2, 4.
Proof Following Stewart and Tijdeman [10] we consider the quantity Note that we have |A | ≤ 2 A, |B | ≤ 2B and |C | ≤ 2C. Let us note that log(1 + x) ≤ 2x provided that |x| ≤ 1/2. Since ab + 1 ≥ 2 we obtain by taking logarithms We apply Matveev's theorem (Lemma 4) with Therefore we have E ≤ max{A , B , C } ≤ 2M provided M > 2 which we clearly may assume. Thus we obtain log s 1 − log 2 < 2.375 · 10 10 log p log q log r log(2M), which yields the upper bound for log s 1 . Since log s 1 = α 1 log p + β 1 log q + γ 1 log r this also yields the upper bounds for α 1 , β 1 and γ 1 .
If we consider instead of T 1 the quantity We end up with the linear form and again an application of Matveev's result yields the same upper bound for log s 2 and also the same upper bounds for α 2 , β 2 and γ 2 .
Finally let us note that which yields after some easy computations the stated upper bounds for log s 4 , α 4 , β 4 and γ 4 .
Let us denote by M 0 and M 4 upper bounds for max{log s 1 , log s 2 } and log s 4 respectively. Then the previous lemma provides upper bounds for M 0 and M 4 .
Proof Following again the arguments of Stewart and Tijdeman [10] we consider the following inequality where A = α 3 − α 5 , B = β 3 − β 5 and C = γ 3 − γ 5 . We apply Matveev's result with On the other hand we have Combining these two inequalities yields the content of the lemma. Now the proof of Proposition 1 is a combination of Lemmas 7 and 8.

Proof of Proposition 1
By inserting the bound for M 0 obtained by Lemma 7 into the bound provided by Lemma 8 we obtain the following inequality M < 9.03 · 10 22 (log P) 5 (log 2M) 2 .
If we put P = 5 and P = 100 we obtain the other bounds stated in the proposition.
Although we have an upper bound for M in Proposition 1 this upper bound provides only rather huge bounds for all the exponents. Moreover, the linear form of logarithms (6) is not suitable for applying standard reduction schemes like LLL-reduction. In particular, there are far too many possible pairs (a, b) left to apply for each possible pair the LLL-reduction to the linear form of logarithms (6) to get in each case a lower bound. Therefore we give a more detailed account on how to find smaller upper bounds for all the exponents.
First, we will prove that all exponents with one possible exception can be bounded in terms of M 0 and M 4 . In particular, Proposition 2 will be useful in the bound reduction process which we will describe in Sect. 5. Before we can state and prove Proposition 2 we prove another helpful lemma: (6) . and such that one of the following two relations holds: (6) and no two indices out of σ (1), (6) and σ (1) and σ (2) are complementary.
Since Lemma 1 we know that α σ (1) = α σ (2) . Assume that σ (1) and σ (2) are not complementary. Then there is exactly one S-unit equation out of the system (1) such that this unit equation contains the index σ (1) but not σ (2). Since α σ (1) is minimal it is equal to one exponent of this unit equation, i.e. α σ (1) If σ (1) and σ (3) are complementary then there is exactly one S-unit equation out of the system (1) that does not contain the indices σ (1) and σ (3) but contains the index σ (2). But since α σ (2) is minimal it is equal to α σ (4) and we have α σ (1) (4) . By exchanging the values of σ (2) and σ (3) the numbers σ (1) and σ (2) are complementary and we have found a permutation that satisfies the second case described in the lemma. A similar argument holds if σ (2) and σ (3) are complementary. Thus we have proved: If σ (1) and σ (2) are not complementary, then either the first case holds or the second case holds after rearranging the order of the α's, i.e. we have found a suitable permutation σ .
Therefore let us assume that σ (1) and σ (2) are complementary. Then there is a unique S-unit equation out of the system (1) such that this unit equation does not contain the indices σ (1) and σ (2). Thus by Lemma 1 we obtain that α σ (3) = α σ (4) and we are in the second case of the lemma.
Similar as in the case of the exponents of p we let τ and ρ be permutations of {1, . . . , 6} such that (6) and (6) respectively. By exchanging the roles of the α's, β's and γ 's we can show that permutations τ and ρ with analogous properties as stated in Lemma 9 exist. Let us fix the permutations σ, ρ and τ . We are now in a position to state the next proposition that will prove to be useful in reducing the huge upper bounds we get from Proposition 1.
Therefore we obtain max{ p α σ (6) , q β τ (6) , r γ ρ(6) } ≤ s 6 = cd + 1 ≤ exp(M 0 + 3M 5 ).  = (a x , b x , c x , a y , b y , c y , a x < a y , and 1 ≤ p a x q b x r c x , p a y q b y r c y ≤ exp(M 4 ), a x = a y , and We define  Proof Let us assume that Case I holds. Then we consider the S-unit equation First, we note that ν p (s y − s x ) = α y = α z and since α x > 0 we have Note that with these notations and assumptions we have after canceling a common factor p α y = p α z the equation with non negative exponents α x , α y , β x , β y , γ x and γ y satisfying α y < α x and 0 ≤ α x log p + β x log q + γ x log r , α y log p + β y log q + γ y log r ≤ M 4 .
Moreover due to Lemma 10 we have either 0 ≤ β z ≤ M 5 log q and 0 ≤ γ z ≤ M or 0 ≤ β z ≤ M and 0 ≤ γ z ≤ M 5 log r . That is (α x , β x , γ x , α y , β y , γ y , β z , γ z ) ∈ M x,y,z . Thus we conclude that Which implies in view of α y ≤ M 4 / log p and (7) the upper bound for α σ (6) . Almost the same arguments apply in the case that α x = α z < α y holds, i.e. that Case II holds. We only have to exchange the roles of x and y.
In Case III, that is we have α x = α y , the proof runs along similar arguments. That is we consider the equation However this time we note that This implies that and therefore the upper bound for α σ (6) . Finally let us note that a non vanishing factor p of s y −s x − 1 ≤ 0, thus we may eliminate all possible powers of p in the formula for C p , which accounts in the usage of p-free parts in the formulas for Case II.
Let us note that we can prove similar results for β τ (6) and γ ρ (6) by exchanging the roles of p, q and r respectively. For the sake of completeness we state the results for β τ (6) and γ ρ (6) in the "Appendix" but do not provide a proof since the proofs are identical after an appropriate permutation of p, q and r .
If we have no explicit bounds for M 4 and M 5 this proposition is hard to apply. However, if we have small, explicit upper bounds for M 4 and M 5 we can estimate C p by applying lower bounds for linear forms in two p-adic logarithms (for details see Sect. 5) which will yield reasonable small bound for M.

The LLL-reduction
In this section we gather some basic facts on LLL-reduced bases, approximation lattices and their applications to Diophantine problems as they can be found in [9,Chapters 5 and 6].
Let L ⊆ R k be a k-dimensional lattice with LLL-reduced basis b 1 , . . . , b k and denote by B be the matrix with columns b 1 , . . . , b k . Moreover, we denote by b * 1 , . . . , b * k the orthogonal basis of R k which we obtain by applying the Gram-Schmidt process to the basis b 1 , . . . , b k . In particular, we have that Further, let us define where · denotes the euclidean norm on R k . It is well known, that by applying the LLL-algorithm it is possible to give in polynomial time a lower bound for l(L, y) ≥c 1 (see e.g. [9, Section 5.4]):

Lemma 11
Let y ∈ R k , z = B −1 y. Furthermore we define -If y / ∈ L let i 0 be the largest index such that z i 0 = 0 and put σ = {z i 0 }, where {·} denotes the distance to the nearest integer.

Finally letc
Then we have In our applications we suppose that we are given real numbers η 0 , η 1 , . . . , η k linearly independent over Q and two positive constantsc 3 ,c 4 such that where the integers x i are bounded by |x i | ≤ X i with X i given upper bounds for 1 ≤ i ≤ k. We write X 0 = max 1≤i≤s {X i }. The basic idea in such a situation, due to de Weger [2], is to approximate the linear form (8) by an approximation lattice. Namely, we consider the lattice L generated by the columns of the matrix whereC is a large constant usually of the size of about X k 0 . Let us assume that we have an LLL-reduced basis b 1 , . . . , b k of L and that we have a lower bound l(L, y) >c 1 with y = (0, 0, . . . , − Cη 0 ). Note thatc 1 can be computed by using the results of Lemma 11. Then we have with these notations the following lemma (e.g. see [9, Lemma VI.1]): . Ifc 2 1 ≥ T 2 + S, then inequality (8) implies that we have either x 1 = x 2 = · · · = x k−1 = 0 and x k = −

Reduction of the bounds
In this section we describe how we can reduce the rather huge bounds for the exponents α i , β i , γ i for i = 1, . . . , 6 for a given set of primes S = {p, q, r } and how it is in practice possible to enumerate for the given set S all S-Diophantine quadruples. We give the full details for the proof of Theorem 1, i.e. the case that S = {2, 3, 5} and will outline the strategy how to find all S-Diophantine quadruples for a general set of three primes S. The reduction process follows in eight steps: Step I We compute an upper bound for M by using Proposition 1. In the case that S = {2, 3, 5} we obtain M < 1.26 · 10 28 . Note that in the case that P < 100 which is relevant for the proof of Theorem 2 we obtain the bound M < 2.88 · 10 30 .
Step II We apply the LLL-reduction method explained in Sect. 4 to the linear forms of logarithms (3) and (5) in order to obtain a small bound for M 0 .
Step III By a direct computation we are able to compute B p , B q and B r by a brute force algorithm. The search implemented in Sage [14] takes on a usual PC a few seconds. In particular we obtain in the case that S = {2, 3, 5} the bounds By Proposition 2 this yields In particular, we have M 5 = M 4 ≤ 272.8.
Step IV In this step we use the theorem of Bugeaud and Laurent (Lemma 5) to estimate the quantities C p and obtain a new upper bound for α σ (6) due to Proposition 3 (and the upper bounds for C q and C r by Propositions 5 and 6 in the "Appendix"). Since the computations for C q and C r are similar we give the details only for C p for Case I note that the upper bounds for Cases II and III can be deduced similarly.
To find an upper bound for C p we split up the set M x,y,z into two subsets M x,y,z , M x,y,z ⊆ M x,y,z with and also consider the two quantities Of course we have C p = max{C p , C p } and we have to find upper bounds for C p (see Step IV (a)) and C p (see Step IV (b)). Similarly we can define C q , C q , C r and C r which will yield bounds for C q and C r respectively.
Step IV (a): In this step we compute an upper bound for C p . We put and b 1 = c z and b 2 = 1. We find that h(η 2 ) ≤ max{b z log q + log( p a y q b y r c y − 1), 5 log q and 0 ≤ a x log p + b x log q + c x log r , a y log p + b y log q + c y log r ≤ M 4 .
Applying Lemma 5 with this data we obtain a new hopefully smaller upper bound for A = α σ (6) .
Step IV (b): In this step we compute an upper bound for C p . In this case we put η 1 = q, η 2 = r c z ( p a y q b y r c y − 1) q b y r c y − p a x −a y q b x r c x and b 1 = b z and b 2 = 1. We find that Again we apply Lemma 5 with this data and we obtain an upper bound for A = α σ (6) in this case.
Similarly we find upper bounds in Cases II and III. We also find upper bounds for B = β τ (6) and C = γ ρ(6) by using Propositions 5 and 6 instead of Propositions 3.
Step V We repeat the process of the Steps II-IV iteratively with the new found bounds for M 0 and M four more times and will find significantly smaller bounds. In the case that S = {2, 3, 5} we obtain after five LLL-reductions the upper bounds M 0 < 33.624, M 4 = M 5 < 67.248 and A < 4.51 · 10 6 , B < 1.3 · 10 6 , C < 1.01 · 10 6 .
We implemented the Steps I-V in Sage. It took a usual desktop PC only a few seconds to obtain the bounds stated in Step V. Before we proceed with Step VI let us summarize what we proved so far in the case that S = {2, 3, 5}: that (a, b, c, d) is a {2, 3, 5}-Diophantine quadruple with a < b < c < d. Then log s 1 , log s 2 < 33.624 and log s 4 < 67.25.
That is we enumerate all S-units s 1 such that the content of Lemma 3 is not violated for the S-Diophantine triple (a, b, c). 5. For all these triples (s 1 , s 2 , s 4 ) we compute a = (s 1 −1)(s 2 −1) We discard all triples (s 1 , s 2 , s 4 ) for which a is not the square of an integer. For the remaining triples we compute a = √ a . 6. We discard all triples (s 1 , s 2 , s 4 ) for which b = s 1 −1 a and c = s 2 −1 a are not integers. 7. We store all triples (a, b, c) in a list Triples.

Let us note that
Step VI is the bottleneck of our algorithm. It took more than six and a half hours to find all these triples using a usual PC (Intel i7-8700) on a single core.
Step VII Let us fix one S-Diophantine triple (a, b, c) from the list Triples computed in Step VI. Let us assume that this triple can be extended to a S-Diophantine quadruple. In this step we reduce the previously found huge upper bound for s 3 = ad + 1 by using the LLL-reduction (Lemma 12). Therefore we consider the inequality and taking logarithms we obtain We apply Lemma 12 to this Diophantine inequality that is we have k = 3, η 1 = log r , η 2 = log q, η 4 = log p, η 0 = log(b/a).
Moreover we put X 1 = A, X 2 = B and X 3 = C, where A, B and C are the bounds for α σ (6) , β τ (6) and γ ρ (6) found in Step V. We distinguish now between two cases.
Case II Let us assume that p, q, r and a/b are multiplicatively dependent and that b/a = p x p q x q r x r . Then Inequality (10) turns into Thus we compute l(L, y) by using Lemma 11 with A = ⎛ ⎝ 1 0 0 0 1 0 C log r C log q C log p ⎞ ⎠ and y = (0, 0, 0) T . We put X 1 = A + |x p |, X 2 = B + |x q | and X 3 = C + |x r | and C sufficiently large and obtain an upper bound for log s 3 . For example in the case of the {2, 3, 5}-Diophantine triple (1,3,5) we have x 2 = x 5 = 0 and x 3 = 1 and obtain with this strategy log s 3 < 39.4 if we chooseC = 10 25 .
Let us note that in the case that S = {2, 3, 5} we obtain in all cases the upper bound log s 3 < 40.
Step VIII Let M 3 denote the upper bound for log s 3 found in Step VII. We enumerate all ad + 1 = s 3 = p α 3 q β 3 r γ 3 with log s 3 < M 3 and for each triple (a, b, c) found in Step VI we compute d = s 3 −1 a . We discard all quadruples (a, b, c, d) for which d is not an integer and for which bd + 1 and cd + 1 are not S-units.
Note that in the case that S = {2, 3, 5} we have M 3 = 40. However a rather quick computer search yields that no triple can be extended to a quadruple and the proof of Theorem 1 is therefore complete.

Remark 1
Finally let us note that the algorithm presented in this paper would also work for arbitrary S with |S| ≥ 3 after some modifications. For instance in Step IV we would have to apply the results of Yu on linear forms in p-adic logarithms (e.g. the results from [15]) instead of the results due to Bugeaud and Laurent [1] on linear forms in two p-adic logarithms. This will lead to much larger bounds for M 0 and M 5 than what we get in the case |S| = 3. Moreover, Step VI in which we search for extendable triples seems to get unfeasible in the case |S| ≥ 4. Thus all together a systematic search for S-Diophantine quadruples or even quintuples in the case |S| ≥ 4 seems to be unfeasible with the method presented in this paper.
Funding Open access funding provided by Paris Lodron University of Salzburg.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. and if α x = α y (Case III) we put M x,y,z = (a x , b x , c x , a y , b y , c y , a z , b z ) ∈ N 8 : 1 ≤ p a x q b x r c x , p a y q b y r c y ≤ exp(M 4 ), c x = c y , and We define C r = max M x,y,z ν r p a z q b z p a y q b y r c y − 1 p a y q b y − p a x q b x r c x −c y − 1 Case I,