Draft. Aimed at High Primes and Misdemeanors.

Doubly focused enumeration of locally square polynomial values

Daniel J. Bernstein

Department of Mathematics, Statistics, and Computer Science (M/C 249) The University of Illinois at Chicago Chicago, IL 60607–7045

                                     djb@cr.yp.to

Abstract. Let f be a nonconstant squarefree polynomial. Which of the values f(c + 1), f(c + 2), …, f(c + H) are locally square at all small primes? This paper presents an algorithm that answers this question in time H/M^{2+o(1)} for an average small c as H → ∞, where M = H^{1/log_{2} log H}. In contrast, the usual method takes time H/M^{1+o(1)}. This paper also presents the results of two record-setting computations: an enumeration of locally square integers up to 24·2^{64}, and an enumeration of locally square values of x^{3} + y^{7} for small x and y.

1 Introduction

A rational number is locally square at a prime p if it is a square in the p-adic field Q_{p}. In particular:

• A nonzero integer is locally square at 2 if and only if it is a square modulo 2, 4, 8, 16, …: equivalently, it has the form 2^{2e}s where s is an odd square modulo 8, i.e., s mod 8 = 1.

• A nonzero integer is locally square at an odd prime p if and only if it is a square modulo p, p^{2}, p^{3}, p^{4}, …: equivalently, it has the form p^{2e}s where s is a nonzero square modulo p, i.e., s^{(p−1)/2} mod p = 1.

A nonnegative rational number is a square if and only if it is locally square at every prime p.

Consider the problem of enumerating all integers x in a specified interval such that a specified polynomial in x is locally square at all primes below a specified bound. Consider, for example, the problem of enumerating all integers x between 0 and 10^{17} − 1 such that x^{3} + 1 is locally square at all primes in {2, 3, …, 251}.

One can simply check, for each of the 10^{17} values of x, the local squareness of x^{3} + 1. This is an example of what I call “unfocused enumeration.”

2000 Mathematics Subject Classification. Primary 11Y16.

A standard improvement is to first compute the possibilities for x modulo some small integer m, then focus on the corresponding arithmetic progressions. For example, x mod 5 must be in {0, 2, 4}, so one can focus attention on the 10^{17}/5 values of x with x mod 5 = 0, then the 10^{17}/5 values of x with x mod 5 = 2, and finally the 10^{17}/5 values of x with x mod 5 = 4, avoiding 40% of the work. Even better, there are only 18 possibilities for x mod 55, so one can focus attention on 18 arithmetic progressions modulo 55, avoiding more than 67% of the work. These are examples of what I call “focused enumeration.”

The number of x’s drops exponentially with the number of prime divisors of m: each additional prime reduces the number of x’s by a factor of approximately 2. But the number of arithmetic progressions, and thus the overhead of considering each arithmetic progression, grows more than exponentially with the number of prime divisors. The minimum computation time is achieved for m somewhere around 10^{14}; the best choice of m depends on the relative costs of considering an arithmetic progression and considering an x.

The point of this paper is that another technique, which I call “doubly focused enumeration,” makes the overhead much smaller, allowing m to be chosen much larger. For example, one can reasonably take m ≈ 3.1 · 10^{25} as 8 times the product of all primes between 3 and 67, reducing the number of x’s by a factor of about 66123.

Section 2 of this paper explains doubly focused enumeration in a more general setting. Section 3 returns to locally square polynomial values:

• It presents an algorithm that uses doubly focused enumeration to figure out which of f(c+1), f(c+2), …, f(c+H) are locally square at all primes p ≤ h, given a nonconstant squarefree polynomial f, an integer c, a positive integer H, and an integer h ≥ 2 log H.

• It proves that this algorithm typically takes time H/M^{2+o(1)} where M = H^{1/lg log H}; here lg = log_{2} as usual. More precisely: The average time for the algorithm, when c is a uniform random element of an interval of length about H^{2}, is H/M^{2+o(1)} for fixed f as H → ∞, provided that c has at most M^{o(1)} digits and h is in M^{o(1)}. I do not know how to prove the same bound for a fixed c.

For comparison: Focused enumeration takes time H/M^{1+o(1)}.

Section 4 and Section 5 report, as examples of doubly focused enumeration, two record-setting computations. The first computation showed that every non-square positive integer below 24 · 2^{64} is locally non-square at some prime in {2, 3, …, 283}; this is numerical evidence for a standard conjecture related to fast deterministic primality proving. The second computation, which is in progress right now, is scanning 10^{20} small pairs (x, y) to find locally square values of x^{3} + y^{7} with x, y coprime. I should have the second computation done by mid-May; perhaps I’ll also extend the first computation.

2 Doubly focused enumeration

Consider the general problem of finding all integers x ∈ [1, H] such that x mod m_{1} ∈ S_{1} and x mod m_{2} ∈ S_{2}. Here H is a positive integer; m_{1} and m_{2} are coprime positive integers; S_{1} is a subset of Z/m_{1}; and S_{2} is a subset of Z/m_{2}.

This section presents three solutions to this problem. The solutions do not have standard names; I call them “unfocused enumeration,” “focused enumeration,” and “doubly focused enumeration,” as in Section 1.

In common applications, focused enumeration is asymptotically faster than unfocused enumeration, and doubly focused enumeration is asymptotically faster than focused enumeration.

The following sizes are typical for applications: H ≈ 10^{20}; m_{1} ≈ m_{2} ≈ 10^{14}; and #S_{1} ≈ #S_{2} ≈ 10^{11}, so that H(#S_{1}/m_{1})(#S_{2}/m_{2}) ≈ 10^{14}. In many situations, one can prove that the number of outputs is approximately H(#S_{1}/m_{1})(#S_{2}/m_{2}); see, for example, Section 3.

Unfocused enumeration. The first method is to consider the possibilities x = 1, x = 2, x = 3, and so on, checking for each x in turn whether x mod m_{1} ∈ S_{1} and x mod m_{2} ∈ S_{2}.

The sets S_{1} and S_{2} can be represented in the obvious way as circular arrays of m_{1} and m_{2} bits respectively. There is very little work for each new x: check the next bit in each array, and record x on the rare occasions that both bits are 1. Common general-purpose computers can check 32 or 64 values of x simultaneously.

Focused enumeration. The second method is to generate, for each r ∈ S_{1}, the arithmetic progression of x ∈ [1, H] such that x mod m_{1} = r, and then check for each x successively whether x mod m_{2} ∈ S_{2}.

The advantage of focused enumeration over unfocused enumeration is that the number of operations drops from H to about H(#S_{1}/m_{1}) + m_{1}. The disadvantage is that S_{2} is no longer checked sequentially.

Doubly focused enumeration. The third method uses the following special case of the explicit Chinese remainder theorem: every integer x ∈ [1, H] may be written as a difference between a reasonably small multiple of m_{1} and a reasonably small multiple of m_{2}. More precisely: x may be written in the form a_{1} − a_{2}, where a_{1} is a multiple of m_{1} in [m_{1}, H + (m_{1} − 1)m_{2}] and a_{2} is a multiple of m_{2} in [0, (m_{1} − 1)m_{2}]. Notice that x mod m_{1} ∈ S_{1} if and only if −a_{2} mod m_{1} ∈ S_{1}, and x mod m_{2} ∈ S_{2} if and only if a_{1} mod m_{2} ∈ S_{2}.

Here is the algorithm. Enumerate, in increasing order, the multiples a_{1} of m_{1} in [m_{1}, H + (m_{1} − 1)m_{2}] such that a_{1} mod m_{2} ∈ S_{2}. Simultaneously enumerate, in increasing order, the multiples a_{2} of m_{2} in [0, (m_{1} −1)m_{2}] such that −a_{2} mod m_{1} ∈ S_{1}. Merge these two lists to see all differences a_{1} − a_{2} in [1, H].

The advantage of doubly focused enumeration over focused enumeration is that the number of operations drops from about H(#S_{1}/m_{1}) + m_{1} to, typically, about H(#S_{1}/m_{1})(#S_{2}/m_{2})+m_{1} +m_{2}. The disadvantage is that each operation is fairly complicated: for example, a multidigit comparison.

This idea is so simple that it must have been written down before. However, I have not been able to locate it in the literature, and it is certainly not widely known in the context of enumerating locally square polynomial values.

Further factorization. In most applications, one can factor m_{1} and m_{2} into much smaller pieces, and correspondingly factor S_{1} and S_{2}. Consider, for example, coprime positive integers m_{11} and m_{12} and sets S_{11} and S_{12} such that m = m_{11}m_{12} and S_{1} = {r : r mod m_{11} ∈ S_{11}, r mod m_{12} ∈ S_{12}}.

For unfocused enumeration: One can store S_{11} and S_{12} instead of S_{1}, using m_{11} + m_{12} bits of memory instead of m_{1} bits of memory. One then checks bits of S_{1} by checking the corresponding bits of S_{11} and S_{12}.

For focused enumeration: One can enumerate values r ∈ S_{1} more quickly than trying each r ∈ Z/m_{1}, by applying the explicit Chinese remainder theorem to each pair in S_{11} × S_{12}. This well-known technique again avoids the need to store S_{1}, and reduces the number of operations from about H(#S_{1}/m_{1}) + m_{1} to about H(#S_{1}/m_{1}) + #S_{1} + m_{11} + m_{12}.

For doubly focused enumeration: The only extra difficulty is that values r ∈ S_{1} need to be enumerated in increasing order. This can be done without much memory; see, e.g., [4].

3 Enumeration of locally square polynomial values

Fix a nonconstant squarefree polynomial f in one variable over Z. Consider the problem of finding all integers x ∈ [1, H] such that f(c + x) is locally square at all primes p ≤ h, given a positive integer H, an integer c, and an integer h ≥ 2 log H.

Here is an algorithm for solving this problem. Select coprime positive integers m_{1} and m_{2} such that p ≤ h for every prime p dividing m_{1}m_{2}. Define S_{1} as the set of r ∈ Z/m_{1} such that f(c + r) is a square in Z/m_{1}, and define S_{2} as the set of r ∈ Z/m_{2} such that f(c + r) is a square in Z/m_{2}. Enumerate all integers x ∈ [1, H] such that x mod m_{1} ∈ S_{1} and x mod m_{2} ∈ S_{2}, as explained in Section 2. Check, for each such x, whether f(c + x) is locally square at all primes p ≤ h.

The rest of this section shows that this algorithm takes time H/M^{2+o(1)} for an average c under mild assumptions, if m_{1} and m_{2} are selected properly. Here M = H^{1/lg log H}, as in Section 1, and o(1) is as H → ∞ for f fixed.

Assume for simplicity that m_{1} and m_{2} are squarefree; that they are each in H/M^{2+o(1)}; that each prime p dividing m_{1}m_{2} is smaller than 2 log H; and that the number of p’s is in (2+o(1))(log H)/log log H = (2+o(1)) lg M. Theorem 3.1 below explains one way to construct m_{1} and m_{2} satisfying these conditions.

Assume also that c has at most M^{o(1)} digits and that h is in M^{o(1)}. The basic operations in the algorithm—checking whether f(c + x) is a square modulo various primes, or is locally square at various primes—then take time M^{o(1)} with negligible memory. (One can speed up the algorithm by using more memory, as discussed in Section 2, and by choosing m_{1} and m_{2} somewhat larger. However, the speedup is only M^{o(1)}, which is not visible at the level of precision of this analysis.)

There are three bottlenecks in the algorithm:

• Checking whether f(c + x) is locally square at all primes p ≤ h, for each x ∈ [1, H] such that x mod m_{1} ∈ S_{1} and x mod m_{2} ∈ S_{2}, i.e., for each x ∈ [1, H] such that f(c + x) is a square modulo m_{1}m_{2}. If c is a uniform random element of an interval of length m_{1}m_{2} then the average number of candidates x is at most H/2^{(2+o(1)) lg M} = H/M^{2+o(1)} by Theorem 3.3 below.

• Checking whether a_{1} mod m_{2} is in S_{2}, for each multiple a_{1} of m_{1} between m_{1} and H + (m_{1} − 1)m_{2}. There are H/M^{2+o(1)} multiples to check.

• Checking whether −a_{2} mod m_{1} is in S_{1}, for each multiple a_{2} of m_{2} between

0 and (m_{1} − 1)m_{2}. There are H/M^{2+o(1)} multiples to check. The total time is H/M^{2+o(1)} if c is a uniform random element of an interval of length m_{1}m_{2} ≈ H^{2}.

Theorem 3.1 Let H be a positive integer. Define u = log H; assume that 2u ≥ 41. Define M = exp(u/lg u). Let v be the largest integer such that the product m of the primes in [1, v] satisfies m ≤ H^{2}/M^{4}. Let v_{1} be the largest integer such that the product m_{1} of the primes in [1, v_{1}] satisfies m_{1} ≤ H/M^{2}. Then there are (2 + o(1))u/log u prime divisors of m; all prime divisors of m are smaller than 2u; and both m_{1} and m/m_{1} are in H/M^{2+o(1)}.

Proof The prime number theorem implies that v is in (1+o(1)) log(H^{2}/M^{4}) = (2 + o(1))u. Thus there are (2 + o(1))u/log 2u = (2 + o(1))u/log u prime divisors of m. Furthermore, m(v + 1) > H^{2}/M^{4} by definition of v, so m is in H^{2}/M^{4+o(1)}; by the same argument, m_{1} is in H/M^{2+o(1)}; so m/m_{1} is also in H/M^{2+o(1)}.

Apply one of the Rosser-Schoenfeld theorems from [15]: the product of the primes in [1, 2u] is larger than exp(2u(1 − 1/ log 2u)) since 2u ≥ 41. But log 2u > log u > (1/2) lg u, so exp(2u/log 2u) < exp(4u/ lg u) = M^{4}; hence the product of the primes in [1, 2u] is larger than H^{2}/M^{4}. Consequently v < 2u.

Theorem 3.2 Let d be a positive integer. Let f be a polynomial of degree d over Z. Let p be an odd prime number that does not divide the leading coefficient of f and does not divide the discriminant of f. Then there are at most (p + (d − 1)\sqrt{p} + d)/2 elements r of Z/p such that f(r) is a square in Z/p.

Proof Define χ : Z/p → {−1, 0, 1} as the Legendre symbol modulo p. Define X_{i} = {r ∈ Z/p : χ(f(r)) = i}. Define a as the leading coefficient of f. Define g as the polynomial a^{−1}f over the field Z/p.

Check the conditions of Weil’s theorem, as stated in [11, Theorem 5.41]: χ is a multiplicative character of Z/p of order 2; g is a monic polynomial over Z/p of positive degree (namely, degree d); g is not a square, because otherwise p would divide the discriminant of f; and g has at most d (in fact, exactly d) distinct roots in its splitting field over Z/p.

Conclusion: X_{1} − X_{−1} = ∑_{r} χ(f(r)) = ∑_{r} χ(ag(r)) is at most (d − 1)\sqrt{p}. But X_{1} + X_{0} + X_{−1} is exactly p; and X_{0} is exactly the number of roots of f in Z/p, which is at most d. Add: 2X_{1} + 2X_{0} ≤ p + (d − 1)\sqrt{p} + d.

Theorem 3.3 Let f be a nonconstant squarefree polynomial over Z. Then there is a function ϵ : N → R, with ϵ ∈ o(1), such that f(r) is a square in Z/m with probability at most 2^{−(1+ϵ(k))k} if m is a squarefree positive integer, k is the number of prime divisors of m, and r is a uniform random element of Z/m.

Proof Write d = deg f ≥ 1. Note that the discriminant of f is nonzero. Define α as the maximum of the following quantities: d; \sqrt{2}; \sqrt{p} for the primes p dividing the discriminant of f; and \sqrt{p} for the primes p dividing the leading coefficient of f.

I claim that if p is prime and r mod p is uniform then f(r) is a square modulo p with probability at most (1 + α/\sqrt{p})/2. Indeed, if p > d^{2} is an odd prime that divides neither the discriminant of f nor the leading coefficient of f, then the probability is at most

\frac{1}{2} (1 + \frac{d − 1}{\sqrt{p}} + \frac{d}{p}) < \frac{1}{2} (1 + \frac{d}{\sqrt{p}}) ≤ \frac{1}{2} (1 + \frac{α}{\sqrt{p}})

by Theorem 3.2. Otherwise \sqrt{p} ≤ α by definition of α; the probability is at most 1 ≤ (1 + α/\sqrt{p})/2.

The probability that f(r) is a square modulo m is the product, over the primes p dividing m, of the probability that f(r) is a square modulo p; which is, in turn, at most the the product of (1 + α/\sqrt{p})/2 for these k primes p; which is, in turn, at most (1 + α/\sqrt{2})(1 + α/\sqrt{3})(1 + α/\sqrt{4}) ⋯ (1 + α/\sqrt{k + 1})/2^{k}; in other words, at most 2^{−(1+ϵ(k))k}, where ϵ(0) = 0 and

ϵ(k) = \frac{−1}{k} (lg (1 + \frac{α}{\sqrt{2}}) + lg (1 + \frac{α}{\sqrt{3}}) + ⋯ + lg (1 + \frac{α}{\sqrt{k + 1}}))

for k ≥ 1. The sum log(1+α/\sqrt{2})+⋯+log(1+α/\sqrt{k + 1}) ≤ α/\sqrt{2}+⋯+α/\sqrt{k + 1} is bounded by a multiple of \sqrt{k + 1}, so ϵ(k) → 0 as k → ∞.

4 Example: locally square integers

Let x be a positive non-square integer in 1 + 8Z. What is the smallest odd prime p such that x^{(p−1)/2} mod p ≠ 1? In other words, what is the smallest odd prime p such that x is divisible by p or locally non-square at p or both?

It is widely conjectured that p/log p is at most (1 + o(1)) lg x, where lg = log_{2}. In fact, no examples are known in which p/log p is larger than lg x. A proof of an explicit bound such as p/log p ≤ 2 lg x would imply, among other things, that there is a deterministic primality-proving algorithm taking essentially cubic time. See [7], [2], [13, pages 130–136], and [14, Section 2].

I have verified that p ≤ 283 for all x < 24 · 2^{64} ≈ 4.4 · 10^{20}. Here are the cutoffs for each p between 149 and 281 inclusive:

p ≤ 149, p/log p < 29.777 if x <  26250887023729 ≈ 1.492 · 2^{44}
p ≤ 157, p/log p < 31.051 if x <  112434732901969 ≈ 1.598 · 2^{46}
p ≤ 173, p/log p < 33.571 if x <  178936222537081 ≈ 1.271 · 2^{47}
p ≤ 181, p/log p < 34.818 if x <  696161110209049 ≈ 1.237 · 2^{49}
p ≤ 193, p/log p < 36.674 if x <  2854909648103881 ≈ 1.268 · 2^{51}
p ≤ 197, p/log p < 37.288 if x <  6450045516630769 ≈ 1.432 · 2^{52}
p ≤ 211, p/log p < 39.426 if x <  11641399247947921 ≈ 1.292 · 2^{53}
p ≤ 227, p/log p < 41.844 if x < 190621428905186449 ≈ 1.323 · 2^{57}
p ≤ 229, p/log p < 42.145 if x < 196640248121928601 ≈ 1.364 · 2^{57}
p ≤ 233, p/log p < 42.745 if x < 712624335095093521 ≈ 1.236 · 2^{59}
p ≤ 239, p/log p < 43.642 if x < 1773855791877850321 ≈ 1.539 · 2^{60}
p ≤ 241, p/log p < 43.940 if x < 2327687064124474441 ≈ 1.009 · 2^{61}
p ≤ 251, p/log p < 45.427 if x < 6384991873059836689 ≈ 1.385 · 2^{62}
p ≤ 257, p/log p < 46.315 if x < 8019204661305419761 ≈ 1.739 · 2^{62}
p ≤ 263, p/log p < 47.199 if x < 10198100582046287689 ≈ 1.106 · 2^{63}
p ≤ 277, p/log p < 49.254 if x < 69848288320900186969 ≈ 1.893 · 2^{65}
p ≤ 281, p/log p < 49.838 if x < 208936365799044975961 ≈ 1.416 · 2^{67}
p ≤ 283, p/log p < 50.129 if x < 24 · 2^{64} = 1.5 · 2^{68}  (not maximal)

This computation took ten days, about 1.2 · 10^{15} clock cycles, on a Pentium 4 running at 1406MHz.

Another way to phrase the same result: Every non-square positive integer below 24 · 2^{64} is locally non-square at some prime in {2, 3, …, 283}. Indeed, write the integer in the form xy^{2} where x is squarefree. If x is a square then the original integer is a square, contradiction. If x ∉ 1 + 8Z then xy^{2} is locally non-square at 2. If x is a non-square in 1 + 8Z then, by this computation, there is an odd prime p ≤ 283 for which x^{(p−1)/2} mod p ≠ 1. If x^{(p−1)/2} mod p = p − 1 then x is locally non-square at p. If x^{(p−1)/2} mod p = 0 then x is divisible by p but, being squarefree, not by p^{2}, so it is locally non-square at p.

A series of previous computations, initiated by Kraitchik in 1924 and continued by Lehmer, Lehmer, Shanks, Patterson, Williams, Stephens, and Lukes, showed with considerably more effort that p ≤ 281 for all x up to about 7 · 10^{19} ≈ 2^{66}. See [8], [9], [10], [16], [13, page 134], and [14]. For example, the computation of Lukes, Patterson, and Williams in [14] was a focused enumeration of all small y such that 1 + 24y is a non-unit square modulo m_{1} = 5 · 7 · 11 · 13; there are about H/27 such values of y in [1, H].

My computation was a doubly focused enumeration, as explained in Section 2 and Section 3, of all small y such that 1 + 24y is a non-unit square modulo both m_{1} = 41 · 43 · 47 · 53 · 59 · 61 · 71 · 73 and m_{2} = 5 · 7 · 11 · 13 · 17 · 19 · 23 · 29 · 31 · 37 · 67. There are about H/1386487 such values of y in [1, H].

Further speed improvements are possible. I could use somewhat larger moduli m_{1} and m_{2}, for example, especially if I balance the primes.

A note on terminology: pseudosquares. Lehmer in [9] defined, for each prime q ≥ 3, the corresponding “pseudo-square” L_{q} as the smallest positive nonsquare integer x ∈ 1 + 8Z for which p > q. The same terminology is used in [17, page 522], [13], and [14].

On the other hand, Lehmer, Lehmer, and Shanks in [10] defined a “pseudosquare” for q as any non-square integer x ∈ 1 + 8Z for which p > q. The same terminology is used in [5]. The “Table of Pseudo-Squares” in [10, page 434] includes one “least solution” column and one “least prime solution” column. A “prime pseudo-square” in this terminology does not have a short name in the previous terminology.

“Pseudo-squares” have an unrelated definition in [1] and [3]: a sequence of “pseudo-squares” is a sequence of integers in which the nth integer is close to n^{2} in the usual metric.

5 Example: locally square values of x^{3} + y^{7}

The problem here is to find all coprime integer pairs (x, y) with y ∈ [1, 1000] and x ∈ [−y^{3}, 10^{17} −1−y^{3}] such that x^{3} +y^{7} is locally square at all primes p ≤ 251. This is inspired by the well-known problem of finding all coprime integer pairs (x, y) for which x^{3} + y^{7} is a square.

Define m = 3 · 5 · 7 · 8 · 11 · 13 · 17 · 19 · 23 · 29 · 31 · 37 · 41 · 43 · 47 · 53 · 59 · 61 · 67 ≈ 3.1433 · 10^{25}. For each y separately, I am doing a doubly focused enumeration of all x ∈ [−y^{3}, 10^{17} − 1 − y^{3}] such that x is coprime to gcd {y, m} and x^{3} + y^{7} is a square modulo m. This enumeration is producing only a small fraction of the original 10^{20} pairs (x, y). I am then checking, for each enumerated x, whether x^{3} + y^{7} is locally square at all primes p ≤ 251, and whether x is coprime to y.

A subsequent draft of this paper will report the results.

Previously known results: 1^{3} + 0^{7} = 1^{2}; (−1)^{3} + 1^{7} = 0^{2}; 0^{3} + 1^{7} = 1^{2}; 2^{3} +1^{7} = 3^{2}; 17^{3} +2^{7} = 71^{2}; 76271^{3} +17^{7} = 21063928^{2}; (−1414)^{3} +65^{7} = 2213459^{2}; (−9262)^{3} + 113^{7} = 15312283^{2}. The last three were discovered in 1993 by Beukers and Zagier.

I should look at a wider range of y’s, including some negative y’s.

References

[1] A. O. L. Atkin, On pseudo-squares, Proceedings of the London Mathematical Society, Third Series 14a (1965), 22–27. ISSN 0024–6115. MR 34:2547.

[2] Eric Bach, Lorenz Huelsbergen, Statistical evidence for small generating sets, Mathematics of Computation 61 (1993), 69–82. ISSN 0025–5718. MR 93k:11089.

[3] R. Balasubramanian, D. S. Ramana, Atkin’s theorem on pseudo-squares, Institut Mathématique, Publications, Nouvelle Série 63 (1998), 21–25. ISSN 0350–1302. MR 99e:11012.

[4] Daniel J. Bernstein, Enumerating solutions to p(a) + q(b) = r(c) + s(d), Mathematics of Computation 70 (2001), 389–394. ISSN 0025–5718. Available from http://cr.yp.to/papers.

html.

[5] Nathan D. Bronson, Duncan A. Buell, Congruential sieves on FPGA computers, in [6] (1994), 547–551. MR 95k:11165.

[6] Walter Gautschi (editor), Mathematics of Computation 1943–1993: a half-century of computational mathematics, American Mathematical Society, Providence, 1994. ISBN 0–8218– 0291–7. MR 95j:00014.

[7] Marshall Hall, Quadratic residues in factorization, Bulletin of the American Mathematical Society 39 (1933), 758–763. ISSN 0273–0979. Available from http://cr.yp.to/bib/entries.

html#1933/hall.

[8] Derrick H. Lehmer, The mechanical combination of linear forms, American Mathematical Monthly 35 (1928), 114–121. ISSN 0002–9890. Available from http://links.jstor.

org/sici?sici=0002-9890%28192803%2935%3A3%3C114%3ATMCOLF%3E2.0.CO%3B2-Z.

[9] Derrick H. Lehmer, A sieve problem on “pseudo-squares”, Mathematical Tables and Other Aids to Computation 8 (1954), 241–242. ISSN 0891–6837. MR 16,113e.

[10] Derrick H. Lehmer, Emma Lehmer, Daniel Shanks, Integer sequence having prescribed quadratic character, Mathematics of Computation 24 (1970), 433–451. ISSN 0025–5718. MR 42:5889.

[11] Rudolf Lidl, Harald Niederreiter, Finite fields, 2nd edition, Encyclopedia of Mathematics and its Applications, 20, Cambridge University Press, Cambridge, 1997. ISBN 0–521–39231– 4. MR 97i:11115.

[12] John H. Loxton (editor), Number theory and cryptography, London Mathematical Society Lecture Note Series, 154, Cambridge University Press, Cambridge, 1990. ISBN 0–521–39877– 0. MR 90m:11003.

[13] Richard F. Lukes, C. D. Patterson, Hugh C. Williams, Numerical sieving devices: their history and some applications, Nieuw Archief voor Wiskunde Series 4 13 (1995), 113–139. ISSN 0028–9825. MR 96m:11082. Available from http://cr.yp.to/bib/entries.html#1995/

lukes.

[14] Richard F. Lukes, C. D. Patterson, Hugh C. Williams, Some results on pseudosquares, Mathematics of Computation 65 (1996), 361–372. ISSN 0025–5718. MR 96e:11010.

[15] J. Barkley Rosser, Lowell Schoenfeld, Approximate formulas for some functions of prime numbers, Illinois Journal of Mathematics 6 (1962), 64–94. ISSN 0019–2082. MR 25:1139.

[16] A. J. Stephens, Hugh C. Williams, An open architecture number sieve, in [12] (1990), 38–75. MR 1 055 399.

[17] Hugh C. Williams, Jeffrey O. Shallit, Factoring integers before computers, in [6] (1994), 481–531. MR 95m:11143.