\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}