Square-free values of reducible polynomials

We calculate admissible values of r such that a square-free polynomial with integer coefficients, no fixed prime divisor and irreducible factors of degree at most 3 takes infinitely many values that are a product of at most r distinct primes.


Introduction
Let H ∈ Z[x] be a non-constant square-free polynomial of degree h with κ irreducible factors. Unless there is good reason to suppose otherwise, one expects that H(n) is squarefree infinitely often. In fact, assuming that H has no fixed prime divisor, one expects that H(n) is a product of precisely κ distinct primes for infinitely many values of n ∈ Z. We shall assume that H has no fixed prime divisor in all that follows, by which we mean that there is no prime p such that p | H(n) for every n ∈ Z.
When κ = 1 and h = 3 it follows from work of Erdős [11] that H(n) is square-free infinitely often. This fact has not yet been extended to a single irreducible quartic polynomial, although work of Granville [13] handles irreducible polynomials of arbitrary degree under the abc conjecture. On the other hand, thanks to work of Hooley [15], in the setting of irreducible cubic polynomials we have an asymptotic formula for the number of n x for which H(n) is square-free. Recent work of Reuss [19] has substantially improved this asymptotic formula, with the outcome that we now have a power saving in the error term.
The topic of almost-prime values of polynomials has been studied extensively by several authors. The DHR weighted sieve of Diamond, Halberstam and Richert (see [5,6]) gives an efficient means of producing values of r such that Ω(H(n)) r for infinitely many n ∈ Z, where Ω(m) is the total number of prime factors of m counted with multiplicity. The following result gives sufficient conditions on polynomials H, not necessarily irreducible, under which H(n) takes infinitely many values that are simultaneously square-free and have at most r prime factors, for relatively small values of r. Theorem 1.1. Assume that H has no fixed prime divisor and that each irreducible factor of H has degree at most 3. Then there exists r ∈ N such that h \ κ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16  prime numbers. The set of all possible sequences obtained from a given starting prime p has the natural structure of a directed graph, say G p , which one can then ask questions about; e.g., is G p a tree? Applying Theorem 1.1 to the polynomial it was shown in [2] that G p is not a tree for a positive proportion of p, and in fact there are infinitely many coprime pairs (a i , q i ) ∈ N 2 , with (q i , q j ) = 1 for i = j, such that G p has a loop of height at most 13 for any p ≡ a i (mod q i ).
For any z 1 let P (z) = p z p and let δ > 0 be a small parameter to be chosen in due course. Our starting point in the proof of Theorem 1.1 is the observation that where Σ 1 = # n x : Ω(H(n)) r, (H(n), P (x δ )) = 1 and Σ 2 = # n x : ∃ p > x δ s.t. p 2 | H(n) . Our strategy is now evident. First we will seek a lower bound for Σ 1 , for which we shall invoke the DHR weighted sieve [5] in §3. Next, in §4, we shall prove that there exists η > 0 such that provided that each irreducible factor of H has degree at most 3. The key input here will come from the work of Reuss [19] that we have already mentioned and which is based on the approximate determinant method. Finally, in §5, we shall turn to the question of optimising the value of r for which the lower bound in Theorem 1.1 holds. A well-known feature of the weighted sieve is that any admissible value of r can only be derived through a precise analysis of certain integrals involving complicated functions σ, f, F that arise as solutions to certain differential-delay problems. These are all described in Theorem 3.1 below. Wheeler, in his thesis [21], gave algorithms for solving the differential-delay problem via Gaussian quadrature with rigorous error bounds, and some of the theoretical work concerning the solution of this problem (see [6] for a comprehensive treatment) relies on his computations. However, [21] is rather difficult to obtain and there appears to have been no published account of Wheeler's algorithms or source code. With this in mind, in §5 we describe an independent implementation based on polynomial approximations. Thanks to the exponential gains in computer power since Wheeler's thesis was written, we can afford to leave more of the work to the computer, to the point that our error estimates require nothing more complicated than the alternating series test or the Lagrange form of the remainder term in Taylor's theorem. We rely heavily throughout on Fredrik Johansson's excellent library Arb [17] for arbitrary-precision interval arithmetic. Thus, we make no assumptions about round-off error, so that, modulo bugs in the code or computer hardware, our results are rigorous. The end result of our computations is the online table [3] of minimum values of r for every pair κ, h with κ h 50, from which Table 1 is derived. The reader wishing to compute r for numbers outside of the tables is welcome to make use of our program, also available at [3].
We conclude the introduction by proving one more result used in [2]. Namely, if we drop the condition Ω(H(n)) r then we can promote the lower bound in Theorem 1.1 to an asymptotic formula by adapting work of Hooley [16,Chapter 4] and invoking Reuss' work. To do so we write n x where Σ 2 is as above and N(x) (resp. E(x)) denotes the number of n x for which H(n) is not divisible by p 2 for any prime p

Polynomial congruences
An important function in our work is the multiplicative arithmetic function ̺(q) that was defined in (1.2). Our assumption that H has no fixed prime divisor is equivalent to the statement that ̺(p) < p for every prime p. We henceforth assume that It will pay dividends to consider the general function ̺(G; q) = #{ν mod q : G(ν) ≡ 0 mod q}, associated to any irreducible polynomial G ∈ Z[x] of degree d and any q ∈ N. Note that ̺(G; q) is a multiplicative function of q and ̺(q) = ̺(H; q). Assuming that p does not divide the content of G, we have the basic estimate for any prime p and any k ∈ N.

The DHR weighted sieve
We begin by recalling the set-up for a version of the DHR weighted sieve [5] that we shall use. Let A be a finite sequence of integers. Under suitable hypotheses, the DHR weighted sieve produces a lower bound for the cardinality #{a ∈ A : Ω(a) r, (a, P (z)) = 1}, for suitable r and z. In fact, in the literature, it is usually recorded as a lower bound for #{a ∈ A : Ω(a) r}, but the results are actually valid for the restricted cardinality in which any a ∈ A is forbidden to have prime factors less than z. (See the proof of [6, Thm. 11.1] for elucidation of this point.) In order to get a lower bound for Σ 1 we take A = {H(n) : n x} and z = x 1/v , for a parameter v > 1.
When κ = 1, one can solve the differential-delay problem explicitly, and this leads to the conclusion that r = h + 1 is admissible for every h 1 (see [20]). Moreover, for κ = 1 and h 2 we may take r = h by the prime number theorem for arithmetic progressions and [18]. These values are reflected in Table 1.
For κ 2 we refer to [5]. One sees that we have a sieve of dimension κ, that all of the hypotheses of [5, Thm 1] are met with X = x, ω(d) = ̺(d), τ = 1 and that the choice µ = h is acceptable. For any v, w ∈ R >0 satisfying v > w + 1, it therefore follows that and f, F is a pair of solutions to the following differential-delay problem: Theorem 3.1 (Diamond, Halberstam and Richert [8,9,10]). Let κ > 1 be a real number, and let σ : R >0 → R be the continuous solution of the system Then there are real numbers α > β > 2 such that the system increases monotonically, and Now it is clear that if δ < 1/v. We therefore conclude that Σ 1 ≫ x/(log x) κ provided that r ⌊R(v, w)⌋ and δ < 1/v, which is plainly satisfactory for Theorem 1.1.

Large square divisors
The goal of this section is to prove the upper bound (1.1) for Σ 2 . Our treatment is modelled on an argument of Hooley [16,Chapter 4]. We partition the interval (x δ , ∞) into We clearly have where Σ 2,i is the number of n x such that p 2 | H(n) for some prime p ∈ I i . We now make the further assumption that δ < 1/11. Once taken together, the following three results suffice to establish (1.1), as required to complete the proof of Theorem 1.1.
Proof. Breaking the sum into residue classes modulo p 2 , we find that  Proof. The contribution to Σ 2,3 from n x such that H(n) = 0 is O(1), which is satisfactory. Hence we may focus on n for which H(n) = 0. We assume that x is large so that p ∤ ∆, whenever p > x 1+δ , where ∆ is given by (2.2). But then p 2 | H(n) implies that p 2 | H i (n) for some i ∈ {1, . . . , κ}, since p doesn't divide any of the resultants Res(H i , H j ). In particular only the irreducible factors of degree 3 occur, since Breaking into residue classes modulo m the inner cardinality is seen to be We make the change of variables n = ν + mu for |u| x/m + 1. 6 At this point we call upon work of Heath-Brown [14,Thm. 15] to estimate the inner cardinality. Given ε > 0 and an absolutely irreducible polynomial F ∈ Z[u, v] of degree d, this shows that there are at most Here T is defined to be the maximum of U e 1 V e 2 , taken over all monomials u e 1 v e 2 which appear in F (u, v) with non-zero coefficient. Moreover, the implied constant is only allowed to depend on d and ε. We shall apply this bound with U = x/m + 1 and F (u, v) = H i (ν + mu) − mv 2 . In particular we may take T V 2 . Thus it follows that , for any ε > 0. We now factorise m = st where s is cube-free and t is cube-full, with (s, t) = 1. Then We complete the proof of the lemma on redefining the choice of ε. Proof. As previously, if x is large enough we have p ∤ ∆ whenever p ∈ I 2 . Thus say. Suppose first that i ∈ {1, . . . , κ} is such that d i = deg(H i ) 2. Then we can assume that H i (n) = p 2 m with m ≪ x 2 /x 2(1−δ) = x 2δ . Arguing as in the proof of the previous lemma we find that Breaking m into its cube-free and cube-full part, we are easily led to the satisfactory contribution Σ in the notation of [19,Lemma 3]. Appealing to [19, page 283] with d = 3 one finds that N (X; A, B) ≪ M 2/3 X 2/3+ε , where M is any choice of parameter satisfying the condition Lemma 7]. A moment's thought shows that M = X 4/11+δ is admissible, which leads to the bound N (X; A, B) ≪ X 10/11+2δ/3+ε . We conclude the proof of the lemma by inserting this into our bound for Σ (i) 2,2 and taking ε sufficiently small in terms of δ.

Numerical estimates
In this section, we describe how to compute rigorous numerical approximations of (3.1). Since our main application requires only integer values of κ 2, we make this restriction for convenience. However, our approach should work just as well for arbitrary real κ 1, and in fact one could likely treat κ as a variable. This might be useful, for instance, for extending the results of [7] to small κ. (On the other hand, recent work of Franze [12] and Blight [1] has shown that lower-bound sieves that are superior to the DHR sieve are possible once κ 3, so such an extension is likely to be of limited interest. where h n (z) = i+j=n f i (z)g j (z) for n < N and We obtain a bounding interval for h N by evaluating the above in interval arithmetic with z replaced by Θ ([0, 1]). • If f is known to be bounded away from 0 then we may compute a representation for its reciprocal via where k n is defined for n < N by the relation To compute a bounding interval for k N we again use interval arithmetic, replacing f (z) in the denominator by a known bounding interval.
• We have Using the Lagrange form of the remainder term in Taylor's theorem, we have Given a bounding interval for f n , this yields one for R n,N via interval arithmetic. • Let I n be a bounding interval for f n . Then we have Similarly, 5.2. Solving the differential-delay system. Next we describe how to compute the solution to Theorem 3.1 for integral κ 2, beginning with α and β. These are found via a number of auxiliary functions. First, we have p, q : R >0 → R defined by the differential-delay Then for α, β we may take any solution to the pair of equations (5.5) (Conjecturally, the solution is unique.) We describe the computation of each piece in turn.

5.2.1.
Computing p. By [6, (12.11)], we have the following integral representation for p: where We split the integral (5.6) as Using the series in (5.7), for any positive integer N and any z ∈ [0, 1], we have On the other hand, differentiating the integral representation in (5.7), we have Comparing series term by term, we see that for any fixed x > 0, the sequence Ein (k) (x)/k! is alternating and decreasing in magnitude. Hence,  For small u 0 , the worst case error in (5.10) occurs when m ≈ N/u 0 , and is roughly of size u −N 0 . If we aim for B bits of precision for all u 0 2, then a sensible choice of parameters is N = B, M = ⌈ B log 2 u 0 ⌉. 5.2.2. Computing q. Since we have restricted to positive integer values of κ, it turns out that q is a polynomial of degree 2κ − 1. If we put q(u) = 2κ−1 n=0 a n u n then, solving (5.3), we find that the coefficients a n are given recursively by n r a n for r < 2κ − 1.
Once the a n have been computed, we work out the coefficients of q(u 0 + u) for integers u 0 1 by polynomial arithmetic. In order to limit the precision loss, we perform these computations in rational arithmetic first before converting to floating point intervals. where A = κ!(2e γ ) κ , which is an expansion of the above type provided that N κ. For u 0 = 1 we similarly have Next, suppose that an expansion of the desired type is known for u 0 −2 and u 0 −1. Solving (3.2), we have Using the expansion around u 0 − 1 to estimate σ(u 0 ) and applying the formulas and algorithms from §5.1 in a straightforward way, we obtain an expansion for σ(u 0 + u).

5.2.4.
Computing 1/σ. For u 0 3, we use (5.2), together with the fact that σ is monotonic [6, (14.6)] to obtain a bounding interval, to compute an expansion for 1/σ(u 0 + u). For u 0 ∈ {1, 2}, this method gives a poor approximation when u is close to 1 since, in each case, the analytic function that agrees with σ(u 0 + u) on [0, 1] has a zero at or just to the left of u = −1. This can be mitigated in various ways, e.g. by using expansions around u 0 = 3/2 or 5/2 instead. In order to keep our later algorithms as simple and uniform as possible, we use the following ad hoc methods to obtain better polynomial approximations of 1/σ(u 0 + u) for u 0 ∈ {1, 2}. First, for u 0 = 1 and u ∈ [0, 1], we have We truncate the series at N and express it in the form N n=0 c n u n by polynomial arithmetic. Since the series has positive terms, the greatest error occurs at u = 0. Hence, to account for the error, it suffices to replace the constant term c 0 by Θ([c 0 , 1]).
Turning to u 0 = 2, we solve (3.2) to find, for u ∈ [0, 1], (5.11) We compute this series using the arithmetic of §5.1 and invert it to obtain a degree-(N + 1) approximation for 1/σ (3 − u). Since the error in (5.11) is concentrated in the highestdegree term, the same is true of our approximation to the inverse. Hence, in the degree-N polynomial part we may replace u by 1 − u without excessive precision loss. (There is some loss due to cancellation between terms, but that can be compensated for by increasing the working precision.) In the final term, we simply replace u N +1 by Θ([0, 1]) and incorporate it into the constant term. Finally, for u 0 = 0, 1/σ(u 0 + u) = Au −κ has an unavoidable singularity at u = 0. However, the computation of 1/σ in this region is only needed in order to compute the integrals v−1 w F (u) du and v−1 w F (u) v−u du, and we can evaluate the contribution from u ∈ [w, 1] by elementary means, as detailed below. 5.2.5. Computing Π and Ξ. For integers u 0 2, we compute expansions for Π(u 0 + u) and Ξ(u 0 + u) via (5.4), using our expansions for p, q, 1/σ and a straightforward application of the operations from §5.1.

5.2.6.
Computing α and β. With these preliminaries in place, we can now proceed with the computation of α and β. First, eliminating the integral terms in (5.5), we see that α is a root of the function l(u) = Π(u) − 2 q(u − 1) + Ξ(u)p(u − 1).
Our computed values of α and β for κ 50 are displayed to 20 decimal place accuracy at [3].

5.2.7.
Computing F and f . We compute expansions for F (u 0 + u) and f (u 0 + u) by a similar strategy to that for σ. Note that for u 0 ⌊α⌋, F (u 0 + u) changes behaviour at u = {α} = α − ⌊α⌋, and similarly for f (u 0 + u) when u 0 ⌊α⌋ + 1. Hence, for each u 0 , we compute two approximations for F (u 0 + u) and f (u 0 + u), one valid for u ∈ [0, {α}] and one for u ∈ [{α}, 1]. Solving (3.3), we have where c is chosen to ensure continuity at u = 0 or u = {α}, and similarly with the roles of F and f reversed. Applying this iteratively yields the desired expansions.

5.3.
Optimising R(v, w). Differentiating (3.1), we find that Since F (u) > 1 for all u, (5.12) is a strictly increasing function of w > 0. It tends to −∞ as w → 0 + and, by the above, is positive at w = β − 1; hence, for each fixed v > β, there is a unique w = w(v) ∈ (0, β − 1) at which R(v, w) is minimal. Given a value of v, we compute w(v) as follows. First we compute c(v) = hvf (v) − κ v−1 2 F (u) du. If c(v) 0 then w(v) 2, and we have which we then solve for w(v). Otherwise, w(v) ∈ (2, β − 1), and we find the root of (5.12) by bisection.
Replacing w by w(v) in (3.1) and applying (5.12), we have To compute the integral, we split it over the intervals [u 0 , u 0 + {α}] and [u 0 + {α}, u 0 + 1] for each relevant integer u 0 , taking the appropriate subinterval around the endpoints w(v) and v − 1. For each interval [u 1 , u 2 ] ⊆ [u 0 , u 0 + 1] with u 0 2 we use the approximation Note that we always have u 2 −u 0 v−u 0 1 2 , so this is accurate to at least N bits. When w(v) < 2, we compute the contribution from [w(v), 2] directly, via .
Thus, we have reduced the problem to one of minimising the single-variable function R(v) = R(v, w(v)). Presumably, R(v) has a unique local (and global) minimum on (β, ∞); although we are not aware of a proof of this, we found no violations of it in practice. Empirically, the optimal v is always larger than α + 3. A good upper bound that is uniform in all parameters is not as easy to describe, but by trial and error we found that v < 200 for every κ h 50. Hence, to compute the table in [3], we first worked out the expansions for F and f for all u 0 < 200. We then used a simple linear search through the points v = ⌈α⌉ + n for n = 1, 2, 3 . . . to find one at which R(v) < min(R(v − 1), R(v + 1)). Finally, we zoomed in on the minimum by a simple bisection algorithm. 5.4. Working precision. We conclude this section with a few words on numerical precision. Since our method for computing f and F is iterative, the total precision loss accumulates as u 0 increases. This mainly affects large values of κ, for which the optimal v can be large. We counteracted this effect by increasing the working precision with κ. Empirically we found that using 12(κ+10) bits of precision was enough to determine the minimum value of R(v, w) to at least 20 significant (decimal) digits, for all κ h 50. The demand for high precision is the main obstacle to extending our computations to larger values of κ, though we expect that κ in the hundreds would still be feasible.