\documentclass[12pt,reqno]{article} \usepackage[usenames]{color} \usepackage{amssymb} \usepackage{graphicx} \usepackage{amscd} \usepackage[colorlinks=true, linkcolor=webgreen, filecolor=webbrown, citecolor=webgreen]{hyperref} \definecolor{webgreen}{rgb}{0,.5,0} \definecolor{webbrown}{rgb}{.6,0,0} \usepackage{color} \usepackage{fullpage} \usepackage{float} \usepackage{psfig} \usepackage{graphics,amsmath,amssymb} \usepackage{amsthm} \usepackage{amsfonts} \usepackage{latexsym} \usepackage{epsf} \setlength{\textwidth}{6.5in} \setlength{\oddsidemargin}{.1in} \setlength{\evensidemargin}{.1in} \setlength{\topmargin}{-.5in} \setlength{\textheight}{8.9in} \newcommand{\seqnum}[1]{\href{http://oeis.org/#1}{\underline{#1}}} \begin{document} \begin{center} \epsfxsize=4in \leavevmode\epsffile{logo129.eps} \end{center} \theoremstyle{plain} \newtheorem{theorem}{Theorem} \newtheorem{corollary}[theorem]{Corollary} \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \theoremstyle{definition} \newtheorem{definition}[theorem]{Definition} \newtheorem{example}[theorem]{Example} \newtheorem{conjecture}[theorem]{Conjecture} \theoremstyle{remark} \newtheorem{remark}[theorem]{Remark} \begin{center} \vskip 1cm{\LARGE\bf A Wieferich Prime Search up to $6.7 \times 10^{15} $ } \vskip 1cm \large Fran{\c c}ois G.\ Dorais \\ Department of Mathematics\\ University of Michigan\\ 530 Church Street\\ Ann Arbor, MI 48109\\ USA\\ \href{mailto:dorais@umich.edu}{\tt dorais@umich.edu}\\ \ \\ Dominic Klyve\\ Department of Mathematics\\ Central Washington University\\ 400 E. University Way\\ Ellensburg, WA 98926 \\ USA\\ \href{mailto:klyved@cwu.EDU}{\tt klyved@cwu.EDU} \\ \end{center} \vskip .2 in \newcommand{\Z}{\mathbb{Z}} \newcommand{\M}{\mathbb{M}} \def\coloneqq{:=} \begin{abstract} A \emph{Wieferich prime} is a prime $p$ such that $2^{p-1} \equiv 1 \pmod {p^2}$. Despite several intensive searches, only two Wieferich primes are known: $p = 1093$ and $p=3511$. This paper describes a new search algorithm for Wieferich primes using double-precision Montgomery arithmetic and a memoryless sieve, which runs significantly faster than previously published algorithms, allowing us to report that there are no other Wieferich primes $p < 6.7 \times 10^{15}$. Furthermore, our method allowed for the efficent collection of statistical data on Fermat quotients, leading to a strong empirical confirmation of a conjecture of Crandall, Dilcher, and Pomerance. Our methods proved flexible enough to search for new solutions of $a^{p-1} \equiv 1 \pmod {p^2}$ for other small values of $a$, and to extend the search for Fibonacci-Wieferich primes. We conclude, among other things, that there are no Fibonacci-Wieferich primes less than $p < 9.7 \times 10^{14}$. \end{abstract} \begin{section}{Introduction} During the 19th and 20th centuries, several different classes of prime numbers were identified and studied because of their relationship to Fermat's Last Theorem (FLT). Most notable among these are Fibonacci-Wieferich primes, Wilson primes, and Wieferich primes. For example, if the first case of Fermat's Last Theorem holds for a prime $p$, then $p$ must be a Fibonacci-Wieferich prime~\cite{SS92}. Similarly, any exponent $p$ which permits a solution to the FLT equation $x^p+y^p = z^p$, while being coprime to $x, y$, and $z$, is a Wieferich prime. The relationship between Wilson primes and FLT is more complex; see \cite{L38}. Although it no longer makes sense to search for a solution to the equation in Fermat's Last Theorem, the questions inspired by these classes of primes still remain, and we can now turn our attention to the study of these primes for their own sake. Of the three above-mentioned classes of primes, Wieferich primes are perhaps the simplest to define. Define the \emph{Fermat quotient} of $2 \bmod p$ to be \begin{equation*} q_2(p) = \frac{a^{p-1}-1}{p}. \end{equation*} It is known from Fermat's little theorem that for any prime $p$, $2^{p-1} -1$ is always divisible by $p$, and therefore $q_2(p)$ is always an integer. If $q_2(p)$ vanishes modulo $p$ (that is, if $2^{p-1} -1$ is divisible by $p^2$), $p$ is said to be \emph{Wieferich}. These primes were first studied by Arthur Wieferich, who in 1909 related them to Fermat's Last Theorem~\cite{W09}. Although Wieferich himself found no example of such a prime, W. Meissner~\cite{M13} in 1913 found that $1093$ was Wieferich, and in 1922 N.~G.~W.~H. Beeger~\cite{B22} showed that 3511 was Wieferich, also. Since 1922, however, no new examples have been found. \end{section} \begin{section}{Previous Searches} Exhaustive searches for new Wieferich primes began with Beeger, and continue today. The last eighty years have seen first computers, then new algorithmic techniques, and finally distributed computing applied to the search for Wieferich primes. Because we have been unable to find a comprehensive summary of the history of these searches in the literature, one has been compiled in Table~\ref{tbl:PreviousSearches}. \begin{table}[h] \begin{center} \caption{Previous Wieferich prime searches\newline} \label{tbl:PreviousSearches} \begin{tabular}{r@{}llr} \multicolumn{2}{c}{Search bound} & Author & Year \\ \hline \multicolumn{2}{c}{ $16000$} & Beeger \cite{B39} & 1940 \\ \multicolumn{2}{c}{ $50000$} & Fr{\"o}berg \cite{F58} & unknown \\ \multicolumn{2}{c}{$100000$} & Kravitz \cite{K60} & 1960\\ \multicolumn{2}{c}{$200183$} & Pearson \cite{P64} & 1964\\ \multicolumn{2}{c}{$500000$} & Riesel \cite{R64} & 1964\\ $ 3$ & $\times 10^7$ & Fr{\"o}berg \cite{F68c} & 1968 \\ $ 3$ & $\times 10^9$ & Brillhart, Tonascia, and Weinberger \cite{BTW69} & 1971 \\ $ 6$ & $\times 10^9$ & Lehmer \cite{L81} & 1981 \\ $ 6.1$ & $\times 10^{10} $ & Clark & c.\ 1996 \\ $ 4$ & $\times 10^{12}$ & Crandall, Dilcher, and Pomerance \cite{CDP97} & 1997\\ $ 4.6$ & $\times 10^{13}$ & Brown and McIntosh \cite{BM2001} & 2001\\ $ 2$ & $\times 10^{14}$ & Crump \cite{C2002} & 2002 \\ $ 1.25$ & $\times 10^{15}$ & Knauer and Richstein \cite{KR2005} & 2005 \end{tabular} \end{center} \end{table} The most recent of these searches, that of Knauer and Richstein, used a distributed approach, and incorporated more than 250 client computers during the course of their search, but in order to include as many computers as possible they were unable to use many of the standard optimizations sometimes used in a search for Wieferich primes. Notably, their code assumed only a 32-bit processor on client machines. \end{section} \begin{section}{Improved Search Methods} Our search employed a number of new algorithmic enhancements not used in previous searches. For computations modulo $p^2$, we used a new ``double-precision'' variant Montgomery arithmetic. Finally, we used a new type of ``memoryless'' sieve to quickly eliminate composites. \subsection{Faster Arithmetic} Our first task was to reduce the time for computing Fermat quotients to about $1\mu\mathrm{s}$ (for the machines we had at hand). This was accomplished by using \emph{double-precision Montgomery arithmetic}. The idea behind (single-precision) Montgomery arithmetic modulo $p$ is that instead of the ring $\Z/p\Z$, we can choose a parameter $r$ coprime to $p$, and use the isomorphic ring \begin{equation*} \M(p,r) = (\{0,1,\dots,p-1\},0,e,{\ominus},{\oplus},{\otimes}) \end{equation*} with the usual negation $\ominus x = (-x) \bmod p$, addition $x \oplus y = (x + y) \bmod p$ and additive identity $0$, but where multiplication is defined by $x \otimes y = xyr^{-1} \bmod p$ and, consequently, the multiplicative identity is defined by $e = r \bmod p$. The advantage of this is that, by choosing $r$ wisely, it is possible to arrange that the product $x \otimes y$ can always be computed without resorting to division by $p$. When the modulus $p$ is odd, one such choice is $r = 2^n > p$, which corresponds to the original idea of P. L. Montgomery~\cite{M85}. For double-precision Montgomery arithmetic, we use the same choice for $r$, but the modulus is now $p^2$. We also use a \emph{double-precision representation} for the elements of $\M(p^2,2^n)$, which consists in representing $x \in \{0,1,\dots,p^2-1\}$ as an ordered pair $(x_0,x_1)$ where $x = x_0 + px_1$ and $0 \leq x_0,x_1 \leq p-1$. Addition and subtraction of such pairs is straightforward. Multiplication is not so obvious, but it can still be done using without resorting to division. The following result shows how to multiply two double-precision elements of $\M(p^2,2^n)$ using only multiplication, addition, and subtraction of nonnegative $n$-bit integers. \begin{lemma} Given two odd numbers $p,q$ such that $0 < p,q < 2^n$ and $pq \equiv 1 \pmod{2^n}$ as parameters. The product of two double-precision elements of $\M(p^2,2^n)$ can be computed using at most 7 multiplications and 8 additions/subtractions of nonnegative $n$-bit integers. The square of a double-precision element of $\M(p^2,2^n)$ can be computed using at most 6 multiplications and 8 additions/subtractions of nonnegative $n$-bit integers. \end{lemma} \begin{proof} Let $x_0 + x_1p$ and $y_0 + y_1p$ be double-precision elements of $\M(p^2,2^n)$. Thus $0 \leq x_0,x_1,y_0,y_1 \leq p-1$. First compute: \begin{align*} t_0 + t_1 2^n &\coloneqq x_0 y_0, & \text{where } 0 &\leq t_0 \leq 2^n-1; \\ u_0 + u_1 2^n &\coloneqq qt_0, & \text{where } 0 &\leq u_0 \leq 2^n-1; \\ v_0 + v_1 2^n &\coloneqq pu_0, & \text{where } 0 &\leq v_0 \leq 2^n-1. \end{align*} This requires 3 multiplications. Note that $t_0 = v_0$ since $pq \equiv 1 \pmod{2^n}$, and so \begin{equation}\label{eq:montlo} x_0y_0 - pu_0 = 2^n(t_1 - v_1). \end{equation} Furthermore, note that $0 \leq t_1, v_1 \leq p-1$. Next compute: \begin{align*} t'_0 + t'_1 2^n &\coloneqq x_0y_1 + x_1y_0 + u_0, & \text{where } 0 &\leq t_0' \leq 2^n-1; \\ u'_0 + u'_1 2^n &\coloneqq qt'_0, & \text{where } 0 &\leq u_0' \leq 2^n-1; \\ v'_0 + v'_1 2^n &\coloneqq pu'_0, & \text{where } 0 &\leq v_0' \leq 2^n-1. \end{align*} This requires 4 multiplications and 3 additions (with carry). Again, $t'_0 = v'_0$ since $pq \equiv 1 \pmod{2^n}$, and so \begin{equation}\label{eq:monthi} (x_0 y_1 + x_1 y_0 + u_0) - pu'_0 = 2^n(t'_1 - v'_1). \end{equation} Furthermore, note that $0 \leq t'_1 \leq 2p-1$ and $0 \leq v'_1 \leq p - 1$. Combining \eqref{eq:montlo} and \eqref{eq:monthi}, we find that \begin{align*} (x_0 + x_1 p)(y_0 + y_1 p) &\equiv x_0 y_0 + (x_0 y_1 + x_1 y_0)p &\pmod{p^2} \\ &\equiv 2^n(t_1 - v_1) + 2^n (t'_1 - v'_1) p \end{align*} Since $-p < t_1 - v_1 < p$ and $-p < t'_1 - v'_1 < 2p$, with at most 5 more additions/subtractions, we can find $z_0, z_1$ such that \begin{equation*} (x_0 + x_1 p)(y_0 + y_1 p) \equiv 2^n(z_0 + z_1 p) \pmod{p} \end{equation*} and $0 \leq z_0, z_1 \leq p-1$. In total, this process requires at most 7 multiplications and 8 additions (possibly with carry). For squaring, we have $x_0 y_1 = x_1 y_0$, so we can save 1 multiplication by computing this product only once. \end{proof} \noindent For comparison, single-precision Montgomery multiplication in $\M(p^2,2^{2n})$ requires 3 multiplications and up to 2 additions/subtractions of nonnegative $2n$-bit integers. For small $n$, multiplication of $2n$-bit integers takes at least 3 times as long as multiplication of $n$-bit integers, and addition of $2n$-bit integers takes 2 times as long as addition of $n$-bit integers. Since multiplication is usually much slower than addition, double-precision Montgomery squaring results in approximately 30\% improvement over single-precision Montgomery squaring. To test whether $p$ is a Wieferich prime, we need to check whether $2^{p-1} \equiv 1 \pmod{p^2}$ or, equivalently, whether $2^{(p-1)/2} \equiv \pm1 \pmod{p^2}$, as suggested by Crandall, et al.~\cite{CDP97}. Our implementation used a standard binary powering ladder to accomplish this. \begin{theorem} Given a $n$-bit prime number $p$ ($n \geq 4$), we can test whether $2^{(p-1)/2} \equiv \pm1 \pmod{p^2}$ using at most $6n+2\lg(n)-10$ multiplications and $12n+\lg(n)-13$ additions/subtractions of nonnegative $n$-bit integers. \end{theorem} \begin{proof} We do the computations using double-precision Montgomery arithmetic in $\M(p^2,2^n)$. To get started, we need to compute the auxiliary parameter $q$ such that $pq \equiv 1 \pmod{2^n}$. This can be done in many ways; our implementation used a Newton iteration that requires $2\lg(n)-4$ multiplications and $\lg(n)-2$ subtractions (of $\leq n/2$-bit integers). To get started with the binary powering ladder, we need to compute the double-precision representation of the multiplicative unit $2^n$. Since $p$ has $n$-bits, this can be done with only 1 subtraction. Finally, the (left-to-right) binary powering ladder with the $(n-1)$-bit exponent $(p-1)/2$ requires $n-1$ squarings and at most $n-1$ doublings. Each squaring uses 6 multiplications and 8 additions/subtractions and each doubling uses up to 4 additions/subtractions. In total, this gives $6n-6$ multiplications and up to $12n-12$ additions/subtractions. \end{proof} \subsection{Faster Sieving} Previous searches of Wieferich primes employed a segmented sieve of Eratosthenes to completely sieve an interval for primes before testing those primes for being Wieferich. Using the testing methods described above, we found that testing a single number for being Wieferich was quite fast (about $1\mu\mathrm{s}$), and that therefore much of our computing time would be spent sieving. For traditional sieves, most of the sieving time is spent sieving small primes. Indeed, sieving an interval of length $\ell$ for the prime $p$ requires about $\ell/p$ memory write operations. So, for example, it takes about as much time to sieve an interval for the six primes $2,3,5,7,11,13$ as it takes to sieve the same interval for the primes $17,19,23,\dots,82139$. With this in mind, we began to look for better ways to sieve out multiples of very small primes. Since memory operations are usually much more costly than elementary arithmetic operations, we looked for sieves that require little or no memory. \subsubsection{The Magic Sieve} Sieving for a few small primes $p_1,\dots,p_k$ amounts to enumerating the elements of the unit group modulo $M = p_1 \cdots p_k$. Our first idea was to make better use of the structure of the unit group $(\Z/M\Z)^*$. After a few experiments, we found a special number that we ended up calling the \emph{Magic Modulus}: \begin{equation*} M = 2 \cdot 3 \cdot 5 \cdot 7 \cdot 11 \cdot 17 \cdot 23 \cdot 29 \cdot 47 \cdot 53 \cdot 59 \cdot 83 = 319514496269430. \end{equation*} This number was chosen because of the simple structure of the unit group mod $M$. Indeed, the group $(\Z/M\Z)^*$ has a large cyclic factor of order $6569843280$, generated by $31$, and a second factor of order $8192 = 2^{13}$. The integers in the interval $[kM+1, (k+1)M]$ that are coprime to $M$ can be listed in a rectangular array \begin{align*} x_{ij} &= kM + (a_i \cdot 31^j \bmod M) &(0 \leq i < 8192, 0 \leq j < 6569843280) \end{align*} where the numbers $a_i$ are chosen representatives of the cosets of the cyclic group $\langle 31 \rangle$ in $(\Z/M\Z)^*$. Given a coset representative $a_i$, listing the integers $x_{ij}$ ($0 \leq j < 6569843280$) only involves multiplying by $31$, which can be accomplished by five doubling and one subtraction operation modulo $M$ and essentially no memory access. While the Magic Sieve is somewhat less efficient than the Spin Sieve (described below), it has the advantage that it is easy to implement and requires essentially no memory storage. In fact, our implementation used only $104$ bytes of data to be stored in memory. The performance of the Magic Sieve was adequate for our purposes --- we only had to compute Fermat quotients for $16.8\%$ of the integers in an interval. \subsubsection{The Spin Sieve} For optimal results, a sieve should use the first few primes $2,3,\dots,p_k$. A weakness of the Magic Sieve is that it is sometimes preferable to omit a few small primes so that the unit group has a large cyclic factor. The Spin Sieve does away with the reliance on the structure of the unit group, while, like the Magic Sieve, requiring very little memory storage. The Spin Sieve was inspired by Pritchard's Wheel Sieve~\cite{P81}. A similar idea was independently discovered by Sorenson~\cite{S06}, from whom we borrowed some implementation ideas. The idea behind the Spin Sieve is to find a simple bijection between the set of $k$-tuples \begin{equation*} T_k = \{1\}\times\{1,2\}\times\{1,2,3,4\}\times\cdots\times\{1,\dots,p_k-1\} \end{equation*} and the set \begin{equation*} A_k = \{ x : 0 < x < M_k, (x,M_k) = 1 \} \end{equation*} where $M_k = 2 \cdot 3 \cdot 5 \cdots p_k$. This bijection $s_k:T_k \to A_k$ is defined recursively by $s_1(1) = 1$, and \begin{equation*} s_k(t_1,\dots,t_k) = s_{k-1}(t_1,\dots,t_{k-1}) + M_{k-1}((r_k(t_1,\dots,t_{k-1}) + t_k) \bmod p_k) \end{equation*} where $r_k:T_{k-1}\to\{0,\dots,p_k-1\}$ satisfies \begin{equation*} s_{k-1}(t_1,\dots,t_{k-1}) + M_{k-1} r_k(t_1,\dots,t_{k-1}) \equiv 0 \pmod{p_k}. \end{equation*} While the definition of $s_k$ is somewhat unwieldy, it is rather simple to compute the values of $s_k$ in succession with the lexicographic ordering of $T_k$. This ``spinning'' operation is the origin of the name of the sieve. There are many ways to use the values of $s_k$. In our implementation, for each value of $s_{10}$, we further sieve (in the traditional way) the arithmetic progression $s_{10}(t_1,\dots,t_{10}) \pmod{M_{10}}$ for the primes $31,\dots,65521$ and compute the Fermat quotients of the remaining values. On average, less than $15.8\%$ of numbers survive the Spin Sieve, then about $32\%$ of the remaining numbers survive the second sieve, so we only compute Fermat quotients for approximately $5\%$ of the numbers in a given interval. \end{section} \begin{section}{Other applications} \subsection{Base-$a$ Wieferich Primes} The definition of Wieferich prime uses base 2 for historical reasons, but mathematically there is no particular reason why we can't consider other bases as well. In an analogy to the base-2 case, define the \emph{Fermat coefficient base-$a$} of $p$ to be \begin{equation*} q_a(p) = \frac{a^{p-1}-1}{p}. \end{equation*} Primes for which $q_a(p)$ vanishes modulo $p$ are sometimes called \emph{base-a} Wieferich primes. Using the methods described above, we searched for solutions to this equation of $a=3, 5,$ and $7$. For each of these bases, we looked for solutions up to about $9.7 \times 10^{14}$ (or, more precisely, solutions not greater than $(2\cdot 3 \cdot 5 \cdot 7 \cdot 11 \cdot 13 \cdot 17\cdot 19\cdot 23\cdot 29) \cdot 1.5 \times 10^{5}$). The previous record for searches with $a=3$ and $a=5$ was held by Keller and Richstein~\cite{KR04}, who searched to $10^{11}$. The previous record for $a=7$ belonged to Montgomery~\cite{M93}, who searched to $2^{32}$. We found no new base-$a$ Wieferich primes for any of these bases, and therefore there are still only two known solutions for each base $2, 3,$ and $7$, and six solutions for base $5$. We do not know whether there is any significance to larger number of base-5 solutions; this may be simply a statistical aberration. We also found 203 primes $p$ with base-3 Fermat quotients less than 100, 179 with base-5, and 212 for base-7. These values, together with near-misses of larger Fermat quotient, are also available on the project web page, while primes with small relative Fermat coefficients base-$a$ are reported in Tables \ref{tbl:smallbase3}--\ref{tbl:smallbase7}. \subsection{Fibonacci-Wieferich primes} Another class of primes initially defined because of Fermat's Last Theorem are the \emph{ Fibonacci-Wieferich} primes, sometimes called \emph{Wall-Sun-Sun} primes. Let $F_u$ denotes the $u$th Fibonacci number, and $\left(\frac{p}{5}\right)$ denotes the Legendre symbol; that is, $$ \left(\frac{p}{5}\right) = \begin{cases} 1, & \text{ if } p \equiv \pm 1 \pmod 5;\\ -1, & \text{ if } p \equiv \pm 2 \pmod 5;\\ 0, & \text{ if } p \equiv 0 \pmod 5.\\ \end{cases} $$ Then although for any prime $p$, \begin{equation*} F_{p-\left(\frac{p}{5}\right)} \equiv 0 \pmod {p}, \end{equation*} there are no known solutions to \begin{equation}\label{fibthm} F_{p-\left(\frac{p}{5}\right)} \equiv 0 \pmod {p^2}. \end{equation} Any solution to \eqref{fibthm} is a Fibonacci-Wieferich prime. Previous searches for Fibonacci-Wieferich primes these primes have extended as far as to $2.0 \times 10^{14}$~\cite{MR07} by McIntosh and Roettger. By modifying the methods described above, we were also able to use our code in the search for these primes. We were able to search to about $9.7 \times 10^{14}$ (or, more precisely, solutions not greater than $(2\cdot 3 \cdot 5 \cdot 7 \cdot 11 \cdot 13 \cdot 17\cdot 19\cdot 23\cdot 29) \cdot 1.5 \times 10^{5}$), in which space no Fibonacci-Wieferich primes were found. We did, however, find several new ``near misses'' -- those primes for which $F_{p-\left(\frac{p}{5}\right)}$ is small $\pmod {p^2}$. These are reported in Table \ref{tbl:smallfib}. A more extensive list of near-misses for Fibonacci-Wieferich primes, together with the near misses for Wieferich primes of each base we studied, will be available on the project web page. After we completed this work, we discovered that McIntosh and Roettger had extended the search described their search~\cite{MR07} to $10^{15}$. Our values match theirs precisely. Because their new values have not yet been published, we include them here. \subsection{Computational considerations with Fibonacci-Wieferich primes} Searching for Fibonacci-Wieferich primes involves calculating $F_{p-\left(\frac{p}{5}\right)} \pmod {p^2}$, for (fairly) large $p$. Naturally space and memory considerations keep us from computing $F_{p-\left(\frac{p}{5}\right)}$ directly. Instead we do these calculations by recalling a well-known identity, namely that $$ \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}. $$ By doing all these computations modulo $n^2$ for some $n$, and the same type of binary ladder that we used above, we can compute our values quite quickly. In fact, we can save even more time by noting that all of our matrix calculations involve either squaring, or multiplying by $\begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}$. Although naively multiplying two $2 \times 2$ matrices requires eight multiplications and four additions, squaring a symmetric matrix requires only three squares $a^2, b^2, c^2$, one product $(a+c)\cdot b$, and 3 additions: $$ \begin{pmatrix} a & b\\ b & c \end{pmatrix}^2 = \begin{pmatrix} a^2 + b^2 & (a+c)\cdot b \\ (a+c)\cdot b & b^2 + c^2 \end{pmatrix}. $$ After implementing this algorithm, a test for a Fibonacci-Wieferich prime runs about four times slower than testing a Wieferich prime. In practice, since some of the multiplications can be done in parallel, the test takes a bit less than four times the time for Wieferich primes. \end{section} \begin{section}{The Computation} The largest part of our computation was the search for Wieferich primes to $6.7 \times 10^{15}$. This was performed on DISCOVERY cluster at Dartmouth College, a cluster of (at the time) about 500 AMD Opteron nodes with 64-bit processors. We ran our code on 24 processors for a period of about 200 days. For short periods of low cluster load, we utilized more processors, once reaching a total of 96. At other times, our computation was tabled for higher priority tasks. Altogether, the search used approximately 12000 CPU days, a value that compares well with that of the previous record search, which used (based on the information provided by Knauer and Richstein) used roughly 50000 CPU days -- although on slower computers. The searches for base-3, 5, and 7 Wieferich primes, along with the search for Fibonacci-Wieferich primes, were performed on the Condor cluster in the Dartmouth College Mathematics Department during low-load times over a period of many months. We used between one and eighteen 64-bit Linux machines of various architectures. Computation times were carefully recorded for these runs. In the following table, the times are recorded: For (base-2) Wieferich primes, the total time taken for the computation $[N, 6.7 \times 10^{15}]$, and for all other primes, the time taken for computing the range $[N, 150 000 N]$, where $N=2\cdot 3 \cdot 5 \cdot 7 \cdot 11 \cdot 13 \cdot 17\cdot 19\cdot 23\cdot 29 = 6469693230$. (Computations of the range $[0,N]$ were very short, and were implemented with unsophisticated methods.) \begin{table}[h] \begin{center} \caption{Time used in CPU days for calculations (Range values are approximate; see discussion above for precise values)} \label{tbl:ugly} \begin{tabular}{llr} \textrm{Search} & \textrm{Range} & \textrm{CPU Days}\\ \hline Wieferich primes & $[6.5 \times 10^9, 6.7 \times 10^{15}]$ & 12907.97\\ Base-3 Wieferich & $[6.5 \times 10^9, 9.7 \times 10^{14}]$ & 714.54\\ Base-5 Wieferich & $[6.5 \times 10^9, 9.7 \times 10^{14}]$ & 812.14\\ Base-7 Wieferich & $[6.5 \times 10^9, 9.7 \times 10^{14}]$ & 916.67\\ Fibonacci-Wieferich & $[6.5 \times 10^9, 9.7 \times 10^{14}]$ & 1978.49 \end{tabular} \end{center} \end{table} \end{section} \begin{section}{Results} \subsection{Wieferich and near-Wieferich primes} As we stated earlier, no new (base-2) Wieferich primes were found. It has become standard practice to report ``near-Wieferich'' primes; that is, those $p$ for which $2^{(p-1)/2} \equiv \pm 1 + Ap \pmod{ p^2}$, where $|A| \le 100$. However, as the magnitude of the primes under consideration grows, the density of the near-Wieferich primes diminishes, and there are fewer to report. We propose, therefore, a new definition: a \emph{near-Wieferich prime} is one for which the value of $A/p$ is small; say less than $10^{-13}$. Table~\ref{tbl:small} gives all such primes not greater than $6.7 \times 10^{15}$. For base-3 Wieferich primes, such a definition excludes some 200 previously unreported primes with $|A| < 100$, and simliar numbers are excluded for base-5 and base-7 primes. Rather than give a table of all these primes here, their values will be given on the project webpage. We might well ask whether our results were to be expected: that is, should we expect to find any Wieferich primes in the region $[1.25 \times 10^{15}, 6.7 \times 10^{15}]$? Certainly $A=A(p)$ can take on any of $p$ values (mod $p$). Assuming that $A$ takes these values randomly, the ``probability'' that $A$ takes any particular value (say, 0) is $1/p$. From this, a heuristic is given in \cite{CDP97} that the expected number of Wieferich primes in the interval $[x,y]$ is \begin{equation}\label{eq:conj} \sum_{x \le p \le y} \frac{1}{p} \approx \log \left( \frac{\log y}{\log x} \right) = \log\log y - \log\log x. \end{equation} From this we would conclude that the expected number of Wieferich primes in our interval is $.0472$, and therefore the lack of such primes there is not surprising. Because our program recorded all $p$ with ``small'' $A$ (that is, all those for which $|A| < 2^{24}$), we compiled a large data set which can be used to give more rigorous (experimental) confirmation of this conjecture. Indeed, our program recorded more than 2.1 million primes $p$ for which $A < 2^{24}$. Using this data, we checked the following conjecture, which follows from the same heuristic as does equation \eqref{eq:conj}: \begin{conjecture}[Crandall, Dilcher, and Pomerance]\label{conj:CDP} % Source? The number of prime $p \in [a,b]$ for which $A \in [K,L]$ is asymptotically \begin{equation*} (L-K) \cdot (\log\log b - \log\log a). \end{equation*} \end{conjecture} The table below presents a small snapshot of our experimental results confirming this conjecture. Complete data will be available on our webpage. In Table~\ref{tbl:stats}, the values in row $i$, column $k$ reflect the number of $p \in [i\times 10^{14}, (i+1)\times 10^{14}]$ with $A \in [4(k-1)\times 10^6, 4k \times 10^6)$. \begin{table}[h] \begin{center} \caption{} \label{tbl:stats} \begin{tabular}{|c|rrrrrrrr|r|} \hline $i$ & $k=1$ & $k=2$ & $k=3$ & $k=4$ & $k=5$ & $k=6$ & $k=7$ & $k=8$ & Expected \\ \hline 15 & 7451 & 7473 & 7435 & 7315 & 7511 & 7356 & 7361 & 7423 & 7380.8\\ 16 & 6779 & 6897 & 6999 & 6862 & 6858 & 6942 & 6879 & 6941 & 6920.8\\ 17 & 6449 & 6448 & 6545 & 6480 & 6391 & 6497 & 6420 & 6622 & 6514.2\\ 18 & 6000 & 6135 & 6099 & 6028 & 6071 & 6080 & 6256 & 6146 & 6152.1\\ 19 & 6053 & 5887 & 5839 & 5866 & 5854 & 5752 & 5911 & 5831 & 5827.7\\ \hline \end{tabular} \end{center} \end{table} Overall, the conjecture holds up very well. Indeed, in the strip given by $k=2$, the relative error between the conjectured and experimental values is never greater than $5.5 \%$. Furthermore, we can plot conjectured and actual numbers of near-Wieferich primes for different values of $i$. The result is a graph that looks remarkably like a straight line. Let $x_i$ be the expected number primes $p$ in the $i$th interval $[i\times 10^{14}, (i+1)\times 10^{14})$ for which $A < 4 \times 10^6$, and let $y_i$ be the actual number of such $p$ in the same interval (the values in column 1). We expect from Conjecture \ref{conj:CDP} that for any $i$, $y_i \approx x_i$. In fact, linear regression on the two data sets returns a best fit equation of \begin{equation*} y = 0.999958129 x + 9.7 \end{equation*} (with $R^2=.9992$), giving strong experimental agreement with the conjectured value. Similar tests using different parts of our data show no meaningful disagreement between the values of $A$ and what is expected heuristically. A data set of all near misses with Fermat quotients less than $2^{24}$ (comprising roughly 2 million primes) is available on the project web page. \end{section} \begin{section}{Acknowledgments} We would like to acknowledge the scientific support received from the Research Computing group at Dartmouth College and the generous amount of computational resources we were given on the DISCOVERY cluster. We thank Carl Pomerance, Giulio Genovese, and Chris Hall for valuable discussions which improved this work. \end{section} \begin{table} \caption{Primes $p < 6.7\times10^{15}$ for which $2^{p-1} \equiv 1 + Ap \pmod {p^2}$, with $|A/p| < 10^{-13}$. (Values of $|A/p|$ are given in multiples of $10^{-14}$.) } \label{tbl:small} \begin{equation*} \begin{array}{rrlcrll} \multicolumn{1}{c}{p} & \multicolumn{1}{c}{A} & \multicolumn{1}{c}{|A/p|} & \hspace{-.5 in} & \multicolumn{1}{c}{p} & \multicolumn{1}{c}{A} & \multicolumn{1}{c}{|A/p|}\\ 1093 & 0 & 0 \\ 3511 & 0 & 0\\ 765760560131939 & -76 & 9.925 & & 3723113065138349 & -36 & 0.967 \\ 993048728162299 & +81 & 8.157 & & 3925342714781797 & -139 & 3.541\\ 1302848719581529 & +76 & 5.833 & & 3948546628939699 & -186 & 4.711\\ 1515362530042687 & +149 & 9.833 & & 4032459967159163 & -172 & 4.265 \\ 1680898792774051 & -96 & 5.711 & & 4143792274787999 & +216 & 5.213\\ 1865546314599557 & +75 & 4.020 & & 4150209531584437 & +48 & 1.157\\ 1885825033325021 & +158 & 8.378 & & 5109286219780877 & -79 & 1.546\\ 2276306935816523 & +6 & 0.264 & & 5131427559624857 & -36 & 0.702\\ 2576594157291871 & -123 & 4.774 & & 5294488110626977 & -31 & 0.586\\ 2718566561783551 & +203 & 7.467 & & 5367369195612269 & +318 & 5.925\\ 2849352392161111 & +255 & 8.949 & & 5464249230405811 & +426 & 7.796\\ 3167939147662997 & -17 & 0.537 & & 5539428831517831 & +230 & 4.152\\ 3383577137448533 & +331 & 9.783 & & 5592905052127597 & +353 & 6.312\\ 3411159925463651 & +176 & 5.160 & & 5625021395769599 & -413 & 7.342\\ 3544715971857451 & +127 & 3.583 & & 5683778474515027 & +332 & 5.841\\ 3660747680296367 & -211 & 5.764 & & 5755502459289463 & +476 & 8.270\\ 3690728733648797 & -334 & 9.050 & & 6227907715670981 & -379 & 6.086\\ 3692386431182551 & +277 & 7.502 & & 6517506365514181 & +58 & 0.890\\ 6521780305210439 & -595 & 9.123\\ \end{array} \end{equation*} \end{table} \begin{table} \caption{Primes $p < 9.7\times10^{14}$ for which $3^{p-1} \equiv 1 + Ap \pmod {p^2}$, with $|A/p| < 10^{-13}$. (Values of $|A/p|$ are given in multiples of $10^{-14}$.) } \label{tbl:smallbase3} \begin{equation*} \begin{array}{rrl} \multicolumn{1}{c}{p} & \multicolumn{1}{c}{A} & \multicolumn{1}{c}{|A/p|}\\ 1006003 & 0 & 0.0\\ 39433103646379 & -1 & 2.536 \\ 61629351935149 & +2 & 3.245 \\ 191293826264479 & +19 & 9.932 \\ 229887238986217 & +11 & 4.785 \\ 431096201990017 & +12 & 2.784 \\ 481589141680567 & +13 & 2.699 \\ 566967768385507 & +12 & 2.117 \\ 631564776981199 & -59 & 9.342 \\ 638096726480497 & +40 & 6.269 \\ \end{array} \end{equation*} \end{table} \begin{table} \caption{Primes $p < 9.7\times10^{14}$ for which $5^{p-1} \equiv 1 + Ap \pmod {p^2}$, with $|A/p| < 10^{-13}$. (Values of $|A/p|$ are given in multiples of $10^{-14}$.) } \label{tbl:smallbase5} \begin{equation*} \begin{array}{rrl} \multicolumn{1}{c}{p} & \multicolumn{1}{c}{A} & \multicolumn{1}{c}{|A/p|}\\ 53471161 & 0 & 0.0\\ 1645333507 & 0 & 0.0\\ 6692367337 & 0 & 0.0\\ 188748146801 & 0 & 0.0\\ 319072335276077 &+2 & 0.627 \\ 419207873154803 & +26 & 6.202 \\ 817486743201059 & -59 & 7.21724\\ \end{array} \end{equation*} \end{table} \begin{table} \caption{Primes $p < 9.7\times10^{14}$ for which $7^{p-1} \equiv 1 + Ap \pmod {p^2}$, with $|A/p| < 10^{-13}$. (Values of $|A/p|$ are given in multiples of $10^{-14}$.) } \label{tbl:smallbase7} \begin{equation*} \begin{array}{rrl} \multicolumn{1}{c}{p} & \multicolumn{1}{c}{A} & \multicolumn{1}{c}{|A/p|}\\ 87121568306639 & +8 & 9.1826\\ \end{array} \end{equation*} \end{table} \begin{table} \caption{New primes for which $F_{p-\left(\frac{p}{5}\right)} \equiv Ap \pmod {p^2}$ satisfies $|A| < 100 \pmod {p^2}$ (for examples with $p < 2 \times 10^{14}$, see \cite{MR07}). } \label{tbl:smallfib} \begin{equation*} \begin{array}{rr} \multicolumn{1}{c}{p} & \multicolumn{1}{c}A\\ 267927950960309 & -9\\ 276225896955847 & 6\\ 299920665662731 & -49\\ 321208072276457 & -98\\ 331961404795379 & -98\\ 399729951985657 & -38\\ 481154641312217 & 31\\ 548865911671993 & -92\\ 549413206041731 & 62\\ 585297174492313 & 11\\ 635696842671829 & 52\\ 732698387434649 & -75\\ \end{array} \end{equation*} \end{table} \begin{thebibliography}{10} \bibitem{B22} N.~Beeger, On a new case of the congruence $2^{p-1} \equiv 1 \pmod{p^2}$, {\it Messenger of Mathematics} \textbf{51} (1922), 149--150. \bibitem{B39} N.~G. W.~H. Beeger, On the congruence {$2\sp {p-1}\equiv1\pmod {p\sp 2$}} and {F}ermat's last theorem, {\it Nieuw Arch. Wiskunde} \textbf{20} (1939), 51--54. \bibitem{BTW69} J.~Brillhart, J.~Tonascia, and P.~Weinberger, On the {F}ermat quotient, {\it Computers in Number Theory}, Proc. Sci. Res. Council Atlas Sympos., No.\ 2, Academic Press, 1971, pp.~213--222. \bibitem{BM2001} R.~Brown and R.~McIntosh, \newline \href{http://www.loaria.fr/~zimmerma/records/Wieferich.status}{\tt http://www.loaria.fr/\char'176zimmerma/records/Wieferich.status}, 2001. \bibitem{CDP97} Richard Crandall, Karl Dilcher, and Carl Pomerance, A search for {W}ieferich and {W}ilson primes, {\it Math. Comp.} \textbf{66} (1997), 433--449. \bibitem{C2002} J.~Crump, reported at \href{http://www.spacefire.com/numbertheory/Wieferich.htm}{\tt http://www.spacefire.com/numbertheory/Wieferich.htm}; website was not functioning at time of publication of this article, 2001. \bibitem{F58} C.~E. Fr{\"o}berg, Some computations with {W}ilson and {F}ermat remainders, {\it Math. Tables Aids Comput.} \textbf{12} (1958), 281. \bibitem{F68c} C.~E. Fr{\"o}berg, On some number-theoretical problems treated with computers, in {\it Computers in Mathematical Research}, North-Holland, 1968, pp.~84--88. \bibitem{KR04} Wilfrid Keller and J\"{o}rg Richstein, Solutions of the congruence $ a^{p-1} \equiv 1 \pmod {p^r}$, {\it Math. Comp.} \textbf{74} (2005), 927--936. \bibitem{KR2005} Joshua Knauer and J{\"o}rg Richstein, The continuing search for {W}ieferich primes, {\it Math. Comp.} \textbf{74} (2005), 1559--1563. \bibitem{K60} Sidney Kravitz, The congruence {$2\sp{p-1}\equiv 1 \ ({\rm mod}p\sp{2})$} for {$p<100,000$}, {\it Math. Comp.} \textbf{14} (1960), 378. \bibitem{L81} D.~H. Lehmer, On {F}ermat's quotient, base two, {\it Math. Comp.} \textbf{36} (1981), 289--290. \bibitem{L38} Emma Lehmer, On congruences involving {B}ernoulli numbers and the quotients of {F}ermat and {W}ilson, {\it Ann. of Math.} (2) \textbf{39} (1938), 350--360. \bibitem{MR07} R. J. McIntosh and E. L. Roettger, A search for Fibonacci-Wieferich and Wolstenholme primes, {\it Math. Comp.} \textbf{76} (2007), 2087--2094. \bibitem{M13} W.~Meissner, \"{U}ber die Teilbarkeit von $2^{p-2}$ durch das Quadrat der Primzahl $p=1093$, {\it Sitzungsberichte der Akademie der Wissenschaften, Berlin} {\bf 35} (1913), 663--667. \bibitem{M85} Peter~L. Montgomery, Modular multiplication without trial division, {\it Math. Comp.} \textbf{44} (1985), 519--521. \bibitem{M93} Peter L. Montgomery, New solutions of $a^{p-1} \equiv 1 \pmod {p^2}$, {\it Math. Comp.} {\bf 61} (1993), 361--363. \bibitem{P64} Erna~H. Pearson, On the congruences {$(p-1)!\equiv -1$} and {$2\sp{p-1}\equiv 1\,({\rm mod}\,p\sp{2})$}, {\it Math. Comp.} \textbf{17} (1964), 194--195. \bibitem{P81} Paul Pritchard, A sublinear additive sieve for finding prime numbers, {\it Comm. ACM} \textbf{24} (1981), 18--23. \bibitem{R64} Hans Riesel, Note on the congruence {$a\sp{p-1}\equiv 1$} {$({\rm mod}$} {$p\sp{2})$}, {\it Math. Comp.} \textbf{18} (1964), 149--150. \bibitem{S06} Jonathan~P. Sorenson, The pseudosquares prime sieve, in {\it Algorithmic Number Theory}, Lecture Notes in Comput. Sci., vol.\ 4076, Springer, 2006, pp.~193--207. \bibitem{SS92} Zhi~Hong Sun and Zhi~Wei Sun, Fibonacci numbers and {F}ermat's last theorem, {\it Acta Arith.} \textbf{60} (1992), 371--388. \bibitem{W09} A.~Wieferich, Zum letzten {F}ermat'schen {T}heorem, {\it J. Reine Angew. Math.} \textbf{136} (1909), 293--302. \end{thebibliography} \bigskip \hrule \bigskip \noindent 2010 {\it Mathematics Subject Classification}: Primary 11A41; Secondary 11Y16, 11Y11. \bigskip \hrule \bigskip \noindent \emph{Keywords:} Wieferich prime, Fibonacci-Wieferich prime, Wall-Sun-Sun prime, wheel sieve, magic sieve, Montgomery arithmetic. \noindent (Concerned with sequences \seqnum{A001220}, \seqnum{A014127}, \seqnum{A123692}, and \seqnum{A123693}.) \bigskip \hrule \bigskip \vspace*{+.1in} \noindent Received June 30 2011; revised version received September 12 2011. Published in {\it Journal of Integer Sequences}, October 16 2011. \bigskip \hrule \bigskip \noindent Return to \htmladdnormallink{Journal of Integer Sequences home page}{http://www.cs.uwaterloo.ca/journals/JIS/}. \vskip .1in \end{document}