Strong Pseudoprimes to Twelve Prime Bases

For the purposes of this paper, it suffices to check square-free odd integers. A computation by Dorais and Klyve [3] shows that a pseudoprime to the bases 2 and 3 must be square-free if it is less than $4.4\cdot 10^{31}$. Since the values Zhang conjectured as candidates for $\psi_{12}$ and $\psi_{13}$ are less than $4.4\cdot 10^{31}$, we restrict our search to square-free numbers.

In general, to rule out integers divisible by squares, it suffices to find all Wieferich primes $\leq B^{1/2}$, and for each such prime, perform appropriate tests. This can easily be done in well under $O(B^{2/3})$ time. See §3.6.1 in [2].

To each prime $p$ we associate a vector called its signature. Let $\nu=(a_{1},a_{2},\ldots,a_{m})$ for $a_{i}$ distinct positive integers. Then the signature of $p$ is

where $v$ denotes valuation with respect to $2$ and $\mbox{ord}_{p}(a)$ is the multiplicative order of $a$ modulo $p$.

If $t>2$, then this theorem greatly limits the number of candidate $k$ values we need to check. In the case $t=2$ the initial sub-product $k$ is just a single prime, so this theorem does not, at first, appear to help. However, it does play a role in sieving for $p_{t}$. For the rest of this paper, we use $\nu$ as the vector containing the first $m=11$ or $12$ prime numbers.

As above, let $k=p_{1}\cdots p_{t-1}$. We explain in this subsection how to find all possible choices for $p_{t}$ to form $n=kp_{t}$ as a strong pseudoprime to the first $m$ prime bases using a simple GCD computation.

Let $b$ be one of the $a_{i}$ from $\nu$. The pseudoprimality condition $b^{kp_{t}-1}\equiv 1\pmod{kp_{t}}$ implies that

i.e., that $p_{t}$ divides $b^{k-1}-1$. Considering multiple bases, we know that $p_{t}$ must divide

However, since we are concerned with strong pseudoprimes, we can consider only the relevant algebraic factor of $b^{k-1}-1$. Following Bleichenbacher’s notation [2], let $k-1=u2^{d}$ for $u$ odd and

This theorem allows us to construct the product of all possible $p_{t}$ values for a given $k$ by computing a GCD. To carry this out, we figure out which of the $h(a_{i},k)$ values would be the smallest, compute that, and then compute the next smallest $h$ value modulo the smallest, perform a GCD, and repeat until we get an answer $\leq k$, or use up all the $a_{i}$ base values.

These computations start off easy as $k$ is small and progressively get harder (in an essentially linear fashion) as $k$ grows. As a practical matter, these computations only need to be done once. As an example, [8] found

to be a strong pseudoprime to the first eight prime bases. When considering $p_{1}=206135341$, they sieved (which we explain in the next section) for $p_{2}=412270681$. It is reasonable to ask whether there is a larger prime that can be paired with $p_{1}$ that would form a strong pseudoprime to eight bases. In our computation we check by computing a GCD to answer “no.”

• Input:

  Integers $t$ and $k$, where $t$ is the number of prime divisors and $k$ is the product of $t-1$ primes with matching signature, and the base vector $\nu$ containing the $m$ smallest primes.

• 1:

  Let $b\in\nu$ give the smallest estimated value for $h(b,k)$ as defined above.

• 2:

  If $b=a_{1}$, let $i=2$ else let $i=1$.

• 3:

  Compute $h(b,k)$.

• 4:

  Compute $x=h(a_{i},k)\bmod h(b,k)$ using modular exponentiation.

• 5:

  Set $x=\gcd(h(b,k),x)$;

• 6:

  while  $x>k$ and $i<m$  do

  • 7:

    Set $i=i+1$; if $a_{i}=b$ set $i=i+1$ again.

  • 8:

    Compute $y=h(a_{i},k)\bmod x$.

  • 9:

    Set $x=\gcd(x,y)$.

• 10:

  end while

• 11:

  if $x<k$ then

  • 12:

    We can rule out $k$.

• 13:

  else

  • 14:

    We factor $x$ and check each prime $p_{t}\mid x$ with $p_{t}>k$ to see if $kp_{t}$ is a strong pseudoprime.

• 15:

  end if

Computing $h(b,k)$ is the bottleneck of this algorithm. When using $m=11$ or higher, we never found anything to factor in the last step.

As above, let $k=p_{1}\cdots p_{t-1}$ and we try to construct $n=kp_{t}$ by finding $p_{t}$. In order to search up to a bound $B$, we need to search primes $p_{t}$ in the interval $[p_{t-1},B/k]$. We combine two types of sieving to search this interval quickly. The first we call $\lambda$-sieving, and the second we call signature sieving. For each prime $p$ define

By the r-rank Artin Conjecture [11], we expect with very high probability that $\lambda_{p}=p-1$. Let $\lambda=\operatorname{lcm}\left\{\lambda_{p_{1}},\ldots,\lambda_{p_{t-1}}\right\}$. Since $a^{kp_{t}-1}\equiv 1\pmod{n}$ for any $a\in\nu$ this means $kp_{t}-1$ is a multiple of $\lambda$. So, $p_{t}\equiv k^{-1}\pmod{\lambda}$. This tells us that $p_{t}$ lies in a specific residue class modulo $\lambda$, and we call searching for $p_{t}$ using this value $\lambda$-sieving.

We can further narrow down the residue classes that $p_{t}$ may lie in with respect to primes in $\nu$. The following two propositions appeared originally in [7].

One may consider higher reciprocity laws, and Jaeschke did so [7]. We found quadratic reciprocity sufficient.

The authors of [8] consider searching for $p_{t}$ modulo $\operatorname{lcm}\{\lambda,9240\}$. They considered $30$ residue classes modulo $9240=8\cdot 3\cdot 5\cdot 7\cdot 11$ that arose from the propositions above. With the inclusion of each additional prime from $\nu$, we effectively rule out half of the cases to consider, thereby doubling the speed of the sieving process.

Ideally, we would like to include every prime in $\nu$ for the best performance, but with normal sieving, the sieve modulus becomes quite large and requires we keep track of too many residue classes to be practical. Instead, we adapted the space-saving wheel datastructure [16, 15], which had been used successfully to sieve for pseudosquares. Sieving for primes $p_{t}$ with specified quadratic character modulo a list of small primes is the same algorithmic problem. This wheel sieve uses a data structure that takes space proportional to $\sum_{a\in\nu}a$ instead of space proportional to $\prod_{a\in\nu}a$ as used in traditional sieving. It does, however, produce candidate primes $p_{t}$ out of order. This is not an issue, so long as we make sure the sieve modulus does not exceed $B$. In practice, we dynamically include primes from $\nu$ so that $1000\cdot k\cdot\operatorname{lcm}\{\lambda,8,3,5,7,11,13,17,19,23,29,31\}<\psi_{13}$, and we favor the inclusion of smaller primes.

• Input:

  $k$ in factored form as $k=p_{1}p_{2}\cdots p_{t-1}$, $\lambda_{p_{i}}$ for $i<t$, search bound $B$, and base vector $\nu$ of length $m$.

• 1:

  Compute $\lambda_{k}=\operatorname{lcm}\{\lambda_{p_{i}}\}$.

• 2:

  Set $w=1$. Compute wheel modulus $w$ as follows:

• 3:

  for  $i=2,\ldots,m$  do

  • 4:

    if  $k\lambda_{k}a_{i}w<B/1000$  then

    • 5:

      if  $a_{i}$ does not divide $\lambda_{k}$  then

      • 6:

        Set $w=w\cdot a_{i}$.

    • 7:

      end if

  • 8:

    end if

• 9:

  end for

• 10:

  if  $\lambda_{k}$ divisible by 4  then

  • 11:

    Set $\lambda=\lambda_{k}$

• 12:

  else

  • 13:

    Set $\lambda=\lambda_{k}/2$ and

  • 14:

    Set $w=8w$

• 15:

  end if

• 16:

  Let $\sigma$ be the signature of $p_{1}$.

• 17:

  Build a wheel with modulus $w$ so that for each prime power $q\mid w$ we have primes $p$ generated by the wheel with $(p/q)$ consistent with $\sigma$ as in Propositions 2.5 and 2.6.

• 18:

  for  each residue $r$ modulo $w$ generated by the wheel  do

  • 19:

    Use the Chinese Remainder Theorem to compute $x\bmod w\lambda$ such that $x\equiv r\bmod w$ and $x\equiv k^{-1}\bmod\lambda$.

  • 20:

    Sieve for primes in the interval $[k+1,B/k]$ that are $\equiv x\bmod w\lambda$.

  • 21:

    For each such probable prime $p_{t}$ found, perform a strong pseudoprime test on $kp_{t}$.

• 22:

  end for

Note that if $k$ consists entirely of primes that are $\equiv 3\bmod 4$, then all primes $p_{t}$ found by the sieve, when using all of $\nu$, will have the correct signature.

In order to compute GCDs or sieve as described above, it is necessary to construct integers $k=p_{1}p_{2}\cdots p_{t-1}$ in an efficient way, where all the primes $p_{i}$ have matching signatures. We do this by storing all small primes in a hash table datastructure that supports the following operations:

• •

  Insert the prime $p$ with its signature $\sigma_{p}$. We assume any prime being inserted is larger than all primes already stored in the datastructure. (That is, insertions are monotone increasing.) We also store $\lambda_{p}$ with $p$ and its signature for use later, thereby avoiding the need to factor $p-1$.

• •

  Fetch a list of all primes from the datastructure, in the form of a sorted array $s[\,]$, whose signatures match $\sigma$. The fetch operation brings the $\lambda_{p}$ values along as well.

We then use this algorithm to generate all candidate $k$ values for a given $t$ and $\nu$:

• Input:

  $t>2$, search bound $B$, cutoff $X<B$, base vector $\nu$ of length $m$.

• 1:

  Let $T$ denote our hash table datastructure

• 2:

  Let $a_{m+1}$ denote the smallest prime not in $\nu$.

• 3:

  for  each prime $p\leq\sqrt{B/a_{m+1}^{t-2}}$ (with $p-1$ in factored form)  do

  • 4:

    Compute $\lambda_{p}$ from the factorization of $p-1$;

  • 5:

    $s[\,]:=T$.fetch($\sigma_{p}$); List of primes with matching signature

  • 6:

    for  $0\leq i_{1}<\cdots<i_{t-2}\leq s$.length  do

    • 7:

      Form $k=s[i_{1}]\cdots s[i_{t-2}]\cdot p$;

    • 8:

      if  $k\leq X$  then

      • 9:

        Use $k$ in a GCD computation to find matching $p_{t}$ values;

    • 10:

      else

      • 11:

        Use $k$ for a sieving computation to find matching $p_{t}$ values;

    • 12:

      end if

  • 13:

    end for

  • 14:

    if  $p\leq(B/a_{m+1}^{t-3})^{1/3}$  then

    • 15:

      $T$.insert($p$,$\sigma_{p}$,$\lambda_{p}$);

  • 16:

    end if

• 17:

  end for

Note that the inner for-loop above is coded using $t-2$ nested for-loops in practice. This is partly why we wrote separate programs for each value of $t$. To make this work as a single program that handles multiple $t$ values, one could use recursion on $t$.

Also note that we split off the GCD and sieving computations into separate programs, so that the inner for-loop was coded to make sure we always had either $k\leq X$ or $X<k<B/p$ as appropriate.

We implemented this datastructure as an array or table of $2^{m}$ linked lists, and used a hash function on the prime signature to compute which linked list to use. The hash function $h:\sigma\rightarrow 0..(2^{m}-1)$ computes its hash value from the signature as follows:

• •

  If $\sigma=(0,0,0,\ldots,0)$, the all-zero vector, then $h(\sigma)=0$, and we are done.

• •

  Otherwise, compute a modified copy of the signature, $\sigma^{\prime}$. If $\sigma$ contains only 0s and 1s, then $\sigma^{\prime}=\sigma$. If not, find the largest integer entry in $\sigma$, and entries in $\sigma^{\prime}$ are 1 if the corresponding entry in $\sigma$ matches the maximum, and 0 otherwise.

  For example, if $\sigma=(3,4,0,4,2,1,2,4)$ (our example from earlier), then $\sigma^{\prime}=(0,1,0,1,0,0,0,1)$.

• •

  $h(\sigma)$ is the value of $\sigma^{\prime}$ viewed as an integer in binary.

  So for our example, $h((3,4,0,4,2,1,2,4))=01010001_{2}=64+16+1=81$.

Computing a hash value takes $O(m)$ time.

Each linked list is maintained in sorted order; since the insertions are monotone, we can simply append insertions to the end of the list, and so insertion time is dominated by the time to compute the hash value, $O(m)$ time. Note that the signature and $\lambda_{p}$ is stored with the prime for use by the fetch operation.

The fetch operation computes the hash value, and then scans the linked list, extracting primes (in order) with matching signature. If the signature to match is binary (which occurs roughly half the time), we expect about half of the primes in the list to match, so constructing the resulting array $s[\,]$ takes time linear in $m$ multiplied by the number of primes fetched. Note that $\lambda$ values are brought along as well. Amortized over the full computation, fetch will average time that is linear in the size of the data returned. Also note that the average cost to compare signatures that don’t match is constant.

The total space used by the data structure is an array of length $2^{m}$, plus one linked list node for each prime $\leq B^{1/3}$. Each linked list node holds the prime, its $\lambda$ value, and its signature, for a total of $O(2^{m}\log B+mB^{1/3})$ bits, since each prime is $O(\log B)$ bits.

Now that all the pieces are in place, we can describe the algorithm as a whole.

Given an input bound $B$, we compute the GCD/sieve cutoff $X$ (this is discussed in §3 below).

We then compute GCDs as discussed above using candidate $k$ values $\leq X$ for $t=2,3,\ldots$. For some value of $t$, we’ll discover that there are no candidate $k$ values with $t-1$ prime factors, at which point we are done with GCD computations.

Next we sieve as discussed above, using candidate $k$ values $>X$ for $t=2,3,\ldots$. Again, for some $t$ value, we’ll discover there are no candidate $k$ values with $t-1$ prime factors, at which point we are done.

We wrote a total of eight programs; GCD programs for $t=2,3,4$ and sieving programs for $t=2,3,4,5,6$. Also, we did not bother to use the wheel and signature sieving, relying only on $\lambda$ sieving, for $t>3$.

We assume a RAM model, with potentially infinite memory. Arithmetic operations on integers of $O(\log B)$ bits (single-precision) and basic array indexing and control operations all are unit cost. We assume fast, FFT-based algorithms for multiplication and division of large integers. For $n$-bit inputs, multiplication and division cost $M(n)=O(n\log n\log\log n)$ operations [14]. From this, an FFT-based GCD algorithm takes $O(M(n)\log n)=O(n(\log n)^{2}\log\log n)$ operations; see, for example, [17].

Let our modulus $q=\prod_{a_{i}\in\nu}a_{i}$ so that $\phi(q)=\prod_{a_{i}\in\nu}\phi(a_{i})$; each base $a_{i}$ is an odd prime with $\phi(a_{i})=a_{1}-1$, except for $a_{1}=8$ with $\phi(a_{1})=4$. We use our two propositions from §2.3, together with quadratic reciprocity to obtain a list of $(a_{i}-1)/2$ residue classes for each odd prime $a_{i}$ and two residue classes modulo $8$. The total number of residue classes is $2^{-m}\phi(q)$. A prime we are searching for with a matching signature must lie in one of these residue classes.

Let us clarify that when we state in sums and elsewhere that a prime $p$ has signature equivalent to $\sigma$, we mean that the quadratic character $(p/a_{i})$ is consistent with $\sigma$ as stated in Propositions 2.5 and 2.6. We write $\sigma\equiv\sigma_{p}$ to denote this consistency under quadratic character, which is weaker than $\sigma=\sigma_{p}$.

We will choose $m$ proportional to $\log B/\log\log B$ so that $q=B^{c}$ for some constant $0<c<1$.

This follows easily from Lemma 3.1 by partial summation. We could make this tighter by making use of Theorem 4 from [9], and observe that the constants $C(q,a)$ (in their notation) cancel with one another in a nice fashion when summing over different residue classes modulo $q$. See also section 6 of that paper.

Let $\pi_{t,\sigma}(x)$ denote the number of integers $\leq x$ with exactly $t$ distinct prime divisors, each of which has signature equivalent to $\sigma$ of length $m$.

Let $\lambda_{p,m}$ denote the value of $\lambda_{p}$ for the prime $p$ using a vector $\nu$ with the first $m$ prime bases. We need the following lemma from [10, Cor. 2.4]:

Using this, we obtain the following.

In order to prove our $B^{2/3+o(1)}$ running time bound for $t>2$, we make use of the conjecture below. Let ${\mathcal{K}}={\mathcal{K}}(t,\sigma)$ be the set of integers with exactly $t$ distinct prime divisors each of which has signature matching $\sigma$. In other words, $\pi_{t,\sigma}(x)$ counts integers in $\mathcal{K}$ bounded by $x$.

We now have the pieces we need to prove our running time bound.

As above, each pseudoprime candidate we construct will have the form $n=kp_{t}$, where $k$ is the product of $t-1$ distinct primes with matching signature. Again, $\nu=(2,3,5,\ldots)$ is our base vector of small primes of length $m$.