text:
\begin{document}

\maketitle

\begin{abstract}
In this paper, we analyze several variants of a simple method for
generating prime numbers with fewer random bits. To generate a prime $p$
less than $x$, the basic idea is to fix a constant $q\propto x^{1-\ve}$,
pick a uniformly random $a<q$ coprime to $q$, and choose $p$ of the form
$a+t\cdot q$, where only $t$ is updated if the primality test fails. We
prove that variants of this approach provide prime generation algorithms
requiring few random bits and whose output distribution is close to
uniform, under less and less expensive assumptions: first a relatively
strong conjecture by H.~Montgomery, made precise by Friedlander and
Granville; then the Extended Riemann Hypothesis; and finally fully
unconditionally using the Barban--Davenport--Halberstam theorem.

\medskip
We argue that this approach has a number of desirable properties compared
to previous algorithms. In particular:
\begin{itemize}
\item
it uses much fewer random bits than both the ``trivial algorithm''
(testing random numbers less than $x$ for primality) and Maurer's almost
uniform prime generation algorithm;

\item
the distance of its output distribution to uniform can be made
arbitrarily small, unlike algorithms like PRIMEINC (studied by Brandt and
Damg\aa{}rd), which we show exhibit significant biases;

\item
all quality measures (number of primality tests, output entropy,
randomness, etc.) can be obtained under very standard conjectures or even
unconditionally, whereas most previous nontrivial algorithms can only be
proved based on stronger, less standard assumptions like the
Hardy--Littlewood prime tuple conjecture.
\end{itemize}

\par\bigskip\noindent
\textbf{Keywords:} {Number Theory,
Cryptography,
Prime Number Generation.}
\end{abstract}

\allowdisplaybreaks

\section{Introduction}
\label{s:intro}

There are several ways in which we could assess the quality of a random
prime generation algorithm, such as its speed (time complexity), its
accuracy (the probability that it outputs numbers that are in fact
composite), its statistical properties (the regularity of the output
distribution), and the number of bits of randomness it consumes to
produce a prime number (as good randomness is crucial to key generation
and not easy to come by \cite{rfc4086}).

In a number of works in the literature, cryptographers have proposed
faster prime
generation algorithms \cite{DBLP:conf/asiacrypt/BrandtDL91,DBLP:conf/crypto/BrandtD92,DBLP:conf/ches/JoyePV00,DBLP:conf/ches/JoyeP06}
or algorithms providing a proof that the generated numbers are indeed prime
numbers \cite{DBLP:conf/eurocrypt/Maurer89,DBLP:journals/joc/Maurer95,DBLP:conf/crypto/Mihailescu94}.

A number of these works also prove lower bounds on the entropy of the
distribution of prime numbers they generate, usually based on very strong
conjectures on the regularity of prime numbers, such as the prime
$r$-tuple conjecture of Hardy--Littlewood~\cite{HL}. However, such bounds
on the entropy do not ensure that the resulting distribution is
statistically close to the uniform distribution: for example, they do not
preclude the existence of efficient distinguishers from the uniform
distribution, which can indeed be shown to exist in most cases.

But some cryptographic protocols (including most schemes based on the
Strong RSA assumption, such as Cramer--Shoup
signatures~\cite{DBLP:journals/tissec/CramerS00}) specifically require
uniformly distributed prime numbers for the security proofs to go
through.

Moreover, some cryptographers, like
Maurer~\cite{DBLP:conf/eurocrypt/Maurer89}, have argued that even for
more common uses of prime number generation, like RSA key generation, one
should preferably generate primes that are almost uniform, so as to avoid
biases in the RSA moduli $N$ themselves, even if it is not immediately
clear how such biases can help an adversary trying to factor $N$. This
view is counterbalanced by results of Mih\u{a}ilescu~\cite{MR1944733}
stating in particular that, provided the biases are not too large (a
condition that is satisfied by the algorithms with large output entropy
mentioned above, if the conjectures used to establish those entropy
bounds hold), then, asymptotically, they can give at most a
polynomial advantage to an adversary trying to factor $N$. This makes the
problem of uniformity in prime number generation somewhat comparable to
the problem of tightness in security reductions.

To the authors' knowledge, the only known prime generation algorithms for
which the statistical distance to the uniform distribution can be bounded
are the one proposed by
Maurer~\cite{DBLP:conf/eurocrypt/Maurer89,DBLP:journals/joc/Maurer95} on
the one hand, and the trivial algorithm (viz. pick a random odd integer
in the desired interval, return it if it is prime, and try again
otherwise) on the other hand. The output distribution of the trivial
algorithm is exactly uniform (or at least statistically close, once one
accounts for the compositeness probability of the underlying randomized
primality checking algorithm), and the same can be said for at least some
variants of Maurer's algorithm, but both of those algorithms have the
drawback of consuming a very large amount of random bits.

By contrast, the PRIMEINC algorithm studied by Brandt and
Damg{\aa}rd~\cite{DBLP:conf/crypto/BrandtD92} (basically, pick a random
number and increase it until a prime is found) only consumes roughly as
many random bits as the size of the output primes, but we can show that
its output distribution, even if it can be shown to have high entropy if
the prime $r$-tuple conjecture holds, is also provably quite far from
uniform, as we demonstrate in \S\ref{sec:why}. It is likely that
most algorithms that proceed deterministically beyond an initial random
choice, including those of Joye, Paillier and
Vaudenay~\cite{DBLP:conf/ches/JoyePV00,DBLP:conf/ches/JoyeP06},
exhibit similar distributional biases.

The goal of this paper is to achieve in some sense the best of both
worlds: construct a prime generation algorithm that consumes much fewer
random bits than the trivial algorithm while being efficient and having
an output distribution that is provably close to the uniform one.

We present such an algorithm in \S\ref{sec:method}: to generate a prime
$p$, the basic idea is to fix a constant $q\sim x^{1-\ve}$, pick a
uniformly random $a<q$ coprime to $q$, and choose $p$ of the form
$a+t\cdot q$, where only $t$ is updated if the primality test fails. We
prove that variants of this approach provide prime generation algorithms
requiring few random bits and whose output distribution is close to
uniform, under less and less expensive assumptions: first a relatively
strong conjecture by H.\,L.~Montgomery, made precise by Friedlander and
Granville; then the Extended Riemann Hypothesis; and finally fully
unconditionally using the Barban--Davenport--Halberstam theorem.

\section{Preliminaries}
\label{sec:prelims}

\subsection{Regularity measures of finite probability distributions}

In this subsection, we give some definitions on distances between random
variables and the uniform distribution on a finite set. We also provide some
relations which will be useful to bound the entropy of our prime
generation algorithms. These results can be found in~\cite{Shoup}.

\begin{defn}[Entropy and Statistical Distance]
Let $X$ and $Y$ be two random variables on a finite set $S$. The
\emph{statistical distance} between them is defined as the $\ell_1$
norm:
\footnote{An alternate definition frequently found in the literature
differs from this one by a constant factor $1/2$. That constant factor is
irrelevant for our purposes.}
\[ \Delta_1(X;Y)=\sum_{s\in S} \Big| \Pr[X=s] - \Pr[Y=s]\Big|. \]
We simply denote by $\Delta_1(X)$ the statistical distance between $X$
and the uniform distribution on $S$:
\[ \Delta_1(X)=\sum_{s\in S} \Big| \Pr[X=s] - \frac1{|S|}\Big|, \]
and say that $X$ is \emph{statistically close to uniform} when $\Delta_1(X)$ is
negligible.
\footnote{For this to be well-defined, we of course need a family of
random variables on increasingly large sets $S$. Usual abuses of language
apply.}

The \emph{squared Euclidean imbalance} of $X$ is the square of the
$\ell_2$ norm between $X$ and the uniform distribution on the same set:
\[\Delta^2_2(X)=\sum_{s \in S} \Big| \Pr[X=s]-1/|S|\Big|^2.\]
We also define the \emph{collision probability} of $X$ as:
\[\beta(X)=\sum_{s\in S} \Pr[X=s]^2,\]
and the \emph{collision entropy} (also known as the R\'enyi entropy)
of $X$ is then $H_2(X) = -\log_2 \beta(X)$. Finally, the
\emph{min-entropy} of $X$ is $H_\infty(X) = -\log_2\gamma(X)$, where
$\gamma(X)=\max_{s\in S}(\Pr[X=s])$.
\end{defn}

\begin{lemma}
Suppose $X$ is a random variable of a finite set $S$.
The quantities defined above satisfy the following relations:
\begin{align}
\gamma(X)^2 &\leq \beta(X) = 1/|S| + \Delta^2_2(X) \leq \gamma(X) \leq
1/|S|+\Delta_1(X), \label{eq:beta}\\
\Delta_1(X) &\leq \Delta_2(X)\sqrt{|S|}. \label{eq:d1d2}
\end{align}
\end{lemma}

\subsection{Prime numbers in arithmetic progressions}

All algorithms proposed in this paper are based on the key idea that, for
any given integer $q>1$, prime numbers are essentially equidistributed
among invertible classes modulo $q$. The first formalization of that idea
is de la Vall{\'e}e Poussin's \emph{prime number theorem for arithmetic
progressions}~\cite{ValleePoussin96}, which states that for any fixed
$q>1$ and any $a$ coprime to $q$, the number $\pi(x;q,a)$ of prime
numbers $p\leq x$ such that $p\equiv a\pmod q$ satisfies:
\begin{equation}
\label{eq:valleepoussin}
\pi(x;q,a) \mathop{\sim}_{x\to+\infty} \frac{\pi(x)}{\varphi(q)}.
\end{equation}

De la Vall{\'e}e Poussin established that estimate for constant $q$, but it
is believed to hold uniformly in a very large range for $q$. In fact,
H.\,L.~Montgomery conjectured~\cite{Montgomery71,Montgomery76} that for any
$\ve > 0$:
\footnote{As is usual in analytic number theory and related subjects,
we use the notations $f(u) \ll g(u)$ and $f(u) = O\big(g(u)\big)$
interchangeably. A subscripted variable on $\ll$ or $O$ means that the implied
constant depends only on that variable.}
\[ \left| \pi(x;q,a) - \frac{\pi(x)}{\varphi(q)} \right| \ll_\ve
   (x/q)^{1/2+\ve} \qquad \big( q<x,\ (a,q)=1 \big),
\]
which would imply that \eqref{eq:valleepoussin} holds uniformly for $q
\ll x/\log^{2+\ve} x$. However, Friedlander and Granville showed~\cite{FG89}
that conjecture to be overly optimistic, and proposed the following
corrected estimate.
\begin{conjecture}[Friedlander--Granville--Montgomery]
\label{conj:fgm}
For $q<x$, $(a,q) = 1$ and all $\ve > 0$, we have:
\[ \left| \pi(x;q,a) - \frac{\pi(x)}{\varphi(q)} \right| \ll_\ve
   (x/q)^{1/2}\cdot x^\ve.
\]
In particular, the estimate \eqref{eq:valleepoussin} holds uniformly for
$q \ll x^{1-3\ve}$.
\end{conjecture}

That conjecture is much more precise than what can be proved using
current techniques, however. The best unconditional result of the same
form is the Siegel--Walfisz theorem~\cite{Walfisz36}, which only implies
that \eqref{eq:valleepoussin} holds in the much smaller range $q \ll
(\log x)^A$ (for any $A>0$).

Stronger estimates can be established assuming the Extended Riemann
Hypothesis ({i.e.} the Riemann Hypothesis for $L$-functions of Dirichlet
characters), which
gives~\cite[p.~125]{Davenport}:
\[ \left| \pi(x;q,a) - \frac{\pi(x)}{\varphi(q)} \right| \ll
   x^{1/2}\log x \qquad \big( q<x,\ (a,q)=1 \big).
\]
This implies \eqref{eq:valleepoussin} in the range $q\ll
x^{1/2}/\log^{2+\ve} x$, which is again much smaller than the one from
Conjecture~\ref{conj:fgm}. The range can be extended using averaging,
however. The previous result under ERH is actually deduced from estimates
on the character sums $\pi(x,\chi) = \sum_{p\leq x} \chi(p)$ for
nontrivial Dirichlet characters $\chi \bmod q$, and more careful
character sum arguments allowed Tur{\'a}n to obtain the following
theorem.
\begin{theorem}[Tur{\'a}n~\cite{Turan37}]
The Extended Riemann Hypothesis implies that for all $q<x$:
\[ \sum_{a\in(\Z/q\Z)^*} \left| \pi(x;q,a) - \frac{\pi(x)}{\varphi(q)}
\right|^2 \ll x(\log x)^2
\]
where the implied constant is absolute.
\end{theorem}
That estimate is nontrivial in the large range $q\ll x/\log^{4+\ve}$,
and implies that \eqref{eq:valleepoussin} holds for all $q$ in that range
and \emph{almost all} $a\in(\Z/q\Z)^*$.

Averaging over the modulus as well, it is possible to obtain fully
unconditional estimates valid in a similarly wide range: this is a result
due to Barban~\cite{Barban66} and Davenport and Halberstam~\cite{DH66}.
We will use the following formulation due to
Gallagher~\cite{Gallagher67}, as stated in \cite[Ch.~29]{Davenport}.
\begin{theorem}[Barban--Davenport--Halberstam]
For any fixed $A > 0$ and any $Q$ such that $x(\log x)^{-A} < Q < x$, we
have:
\[ \sum_{q\leq Q} \sum_{a\in(\Z/q\Z)^*}
   \left| \pi(x;q,a) - \frac{\pi(x)}{\varphi(q)} \right|^2 \ll_A
   \frac{xQ}{\log x}.
\]
\end{theorem}

Finally, we will also need a few classical facts regarding Euler's totient
function (for example, \cite[Th. 328 \& 330]{hardy}).
\begin{lemma}
The following asymptotic estimates hold:
\begin{align}
\varphi(q) &\gg \frac{q}{\log\log q}, \\
\Phi(x)    &:=  \sum_{q\leq x} \varphi(q) = \frac{3x^2}{\pi^2} + O(x\log x).
\end{align}
\end{lemma}

\section{Close-to-uniform prime number generation with fewer random bits}
\label{sec:method}

\subsection{Basic algorithm}

A simple method to construct obviously uniformly distributed prime
numbers up to $x$ is to pick random numbers in $\{1,\dots,\lfloor
x\rfloor\}$ and retry until a prime is found. However, this method
consumes $\log_2 x$ bits of randomness per iteration (not counting the
amount of randomness consumed by primality testing), and hence an
expected amount of $(\log x)^2/\log 2$ bits of randomness to produce a
prime, which is quite large.

As mentioned in the introduction, we propose the following algorithm to
generate almost uniform primes while consuming fewer random bits: first
fix an integer
\begin{equation}
\label{eq:rangeq}
q\propto x^{1-\ve}
\end{equation}
and pick a random $a\in(\Z/q\Z)^*$. Then, search for prime numbers
$\leq x$ of the form $p = a + t\cdot q$.
This method, described as Algorithm~\ref{alg:basic}, only consumes
$\log_2 t = \ve\log_2 x$ bits of randomness per iteration, and the
probability of success at each iteration is
$\sim \frac{\pi(x;q,a)}{x/q}$. Assuming that Conjecture~\ref{conj:fgm} is
true, which ensure that~\eqref{eq:valleepoussin} holds in the
range~\eqref{eq:rangeq}, this probability is about $q/\big(\varphi(q)\log
x\big)$, and the algorithm should thus consume roughly:
\begin{equation}
\label{eq:basicrandbits}
\ve\cdot\frac{\varphi(q)}{q}\cdot\frac{(\log x)^2}{\log 2}
\end{equation}
bits of
randomness on average: much less than the trivial algorihm. Moreover, we
can also show, under the same assumption, that the output distribution
is statistically close to uniform and has close to maximal entropy.

\begin{algorithm}[t]
\caption{Our basic algorithm.}
\label{alg:basic}
\begin{algorithmic}[1]
  \State Fix $q \propto x^{1-\ve}$
  \State $a \randgets (\Z/q\Z)^*$
    \label{algstep:basic-a}
    \Comment{considered as an element of $\{1,\dots,q-1\}$}
  \Repeatx{forever}
    \State $t \randgets \{0,\dots,\lfloor \frac{x-a}q \rfloor\}$
    \State $p \gets a + t\cdot q$
    \Ifx{$p$ is prime}
      \label{algstep:test}
      \textbf{return} $p$
  \EndRepeatx
\end{algorithmic}
\end{algorithm}

We establish those results in \S\ref{sec:fgm}, and show in
\S\ref{sec:basic-erh} that Tur{\'a}n's theorem can be used to obtain
nearly the same results under the Extended Riemann Hypothesis. ERH is not
sufficient to prove that Algorithm~\ref{alg:basic} terminates almost
surely, or to bound the expectation of the number of random bits it
consumes, due to the possibly large contribution of negligibly few values
of $a$. We can avoid these problems by modifying the algorithm slightly,
as discussed in \S\ref{sec:as-erh}. Finally, in \S\ref{sec:bdh}, we show
that unconditional results of the same type can be obtained using the
Barban--Davenport--Halberstam theorem, for another slightly different
variant of the algorithm.

Before turning to these analyses, let us make a couple of remarks on
Algorithm~\ref{alg:basic}. First, note that one is free to choose $q$ in
any convenient way in the range~\eqref{eq:rangeq}. For example, one could
choose $q$ as the largest power of $2$ less than $x^{1-\ve}$, so as to
make Step~\ref{algstep:basic-a} very easy. It is preferable, however, to
choose $q$ as a (small multiple of a) primorial, to minimize the ratio
$\varphi(q)/q$, making it as small as $\propto 1/\log\log q \sim
1/\log\log x$; this makes the expected number of iterations and the
expected amount~\eqref{eq:basicrandbits} of consumed randomness
substantially smaller. In that case, Step~\ref{algstep:basic-a} becomes
slightly more complicated, but this is of no consequence.

Indeed, our second observation is that Step~\ref{algstep:basic-a} is
always negligible in terms of running time and consumed randomness
compared to the primality testing loop that follows. Indeed, even the
trivial implementation (namely, pick a random $a\in\{0,\dots,q-1\}$ and
try again if $\gcd(a,q)\neq 1$) requires $q/\varphi(q) \ll \log\log q$
iterations on average. It is thus obviously much faster than the
primality testing loop, and consumes $\ll \log x\log\log x$ bits of
randomness, which is negligible compared to~\eqref{eq:basicrandbits}.
Furthermore, an actual implementation would take advantage of the known
factorization of $q$ and use a unit generation algorithm such as the one
proposed by Joye and Paillier~\cite{DBLP:conf/ches/JoyeP06}, which we can
show requires only $O(1)$ iterations on average.

Finally, while we will not discuss the details of the primality test of
Step~\ref{algstep:test}, and shall pretend that it returns exact results
(as the AKS algorithm~\cite{AKS} would, for example), we note that it is
fine (and in practice preferable) to use a probabilistic compositeness
test such as Miller--Rabin~\cite{Rabin} instead, provided that the number
of rounds is set sufficiently large as to make the error probability
negligible. Indeed, the output distribution of our algorithm then stays
statistically close to uniform, and the number of iterations is never
larger.

\subsection{Analysis under the Friedlander--Granville--Montgomery
conjecture}
\label{sec:fgm}

As mentioned above, it is straightforward to deduce from the
Friedlander--Granville--Montgomery conjecture that
Algorithm~\ref{alg:basic} terminates almost surely, and to bound its
expected number of iterations and amount of consumed randomness.

\begin{thm}
\label{th:termfgm}
Assume that Conjecture~\ref{conj:fgm} holds. Then
Algorithm~\ref{alg:basic} terminates almost surely, requires
$(1+o(1))\varphi(q)/q\cdot\log x$ iterations of the main loop on average,
and consumes
$ \big(\ve + o(1)\big)\cdot\frac{\varphi(q)}{q}\cdot\frac{(\log x)^2}{\log 2}
$
bits of randomness on average.
\end{thm}
\begin{proof}
Indeed, fix $q\propto x^{1-\ve}$.
Conjecture~\ref{conj:fgm} implies, uniformly over $a\in(\Z/q\Z)^*$:
\[ \left| \pi(x;q,a) - \frac{\pi(x)}{\varphi(q)} \right| \ll
   (x/q)^{1/2}\cdot x^{\ve/4} \propto x^{3\ve/4},
\]
which is negligible compared to $\pi(x)/\varphi(q) \gg x^\ve/\log x$. As
a result, we get $\pi(x;q,a) = (1+o(1))\pi(x)/\varphi(q) =
(1+o(1))/\varphi(q)\cdot x/\log x$ uniformly over
$a$, and the success probability of the main loop becomes:
\[ \frac{\pi(x;q,a)}{1 + \lfloor\frac{x-a}q\rfloor} =
\frac{q}{\varphi(q)}\cdot\frac{1+o(1)}{\log x} \]
which implies the stated results immediately.\qed
\end{proof}

Now let $X$ be the output distribution of Algorithm~\ref{alg:basic},
{i.e.} the distribution on the set of prime numbers $\leq x$ such that
Algorithm~\ref{alg:basic} outputs a prime $p$ with probability exactly
$\Pr[X = p]$. Clearly, we have, for all $(a,q)=1$ and all $t$ such that
$a+t\cdot q\leq x$ is prime:
\[ \Pr[X = a + t\cdot q] = \frac1{\varphi(q)}\cdot\frac1{\pi(x;q,a)}. \]
As a result, the squared Euclidean imbalance of $X$ is:
\begin{align*}
   \Delta_2^2(X)
&= \sum_{a\in(\Z/q\Z)^*} \sum_{a+tq\leq x\text{ prime}}
      \Big|\Pr[X=a+tq] - \frac1{\pi(x)}\Big|^2 +
      \sum_{p|q}\frac1{\pi(x)^2} \\
&= \sum_{a\in(\Z/q\Z)^*} \pi(x;q,a)
      \Big|\frac1{\varphi(q)}\cdot\frac1{\pi(x;q,a)} - \frac1{\pi(x)}\Big|^2 +
      \sum_{p|q}\frac1{\pi(x)^2} \\
&= \frac1{\pi(x)^2} \sum_{a\in(\Z/q\Z)^*} \frac1{\pi(x;q,a)}
      \Big|\pi(x;q,a) - \frac{\pi(x)}{\varphi(q)}\Big|^2 +
      \sum_{p|q}\frac1{\pi(x)^2} \\
&\ll \frac1{\pi(x)^2} \sum_{a\in(\Z/q\Z)^*} \frac{\log x}{x^\ve}\cdot x^{3\ve/2}
 \ll \frac{\log^3 x}{x^2} \cdot \varphi(q)x^{\ve/2}
 \ll \frac{\log^3 x}{x^{1+\ve/2}} \ll \frac1{x^{1+\ve/3}}.
\end{align*}
We can then deduce the following.

\begin{thm}
\label{th:statfgm}
Assume that Conjecture~\ref{conj:fgm} holds. Then the output distribution
of Algorithm~\ref{alg:basic} is statistically close to uniform, and its
collision entropy is only negligibly smaller than that of the uniform
distribution.
\end{thm}
\begin{proof}
Indeed, by~\eqref{eq:d1d2}, the statistical distance to the uniform
distribution satisfies:
\[ \Delta_1(X) \leq \Delta_2(X)\sqrt{\pi(x)}
     \ll \frac{1}{x^{1/2+\ve/6}} \sqrt{\frac{x}{\log x}} \ll x^{-\ve/6},
\]
which is negligible. Moreover, the collision probability is:
\[ \beta(X) = \frac1{\pi(x)} + \Delta_2^2(X)
     = \frac1{\pi(x)}\left(1 + O\Big(\frac{\pi(x)}{x^{1+\ve/3}}\Big)\right)
     = \frac1{\pi(x)}\left(1 + o\big(x^{-\ve/3}\big)\right).
 \]
Hence:
\[ H_2(X) = \log_2\big(\pi(x)\big) - \log_2\big(1 + o(x^{-\ve/3})\big)
          = (H_2)_{\max} - o(x^{-\ve/3})
\]
as required.\qed
\end{proof}

\subsection{Analysis under the Extended Riemann Hypothesis}
\label{sec:basic-erh}

Assume the Extended Riemann Hypothesis, and denote by $\alpha$ the
fraction of all possible choices of $a\in(\Z/q\Z)^*$ such that the error
term $E(x;q,a) := \big| \pi(x;q,a) - \pi(x)/\varphi(q) \big|$ satisfies
$E(x;q,a) > x^{3\ve/4}$. Then, Tur{\'a}n's theorem asserts that:
\[ \sum_{a\in(\Z/q\Z)^*} E(x;q,a)^2 \ll x(\log x)^2, \]
and the left-hand side is greater or equal to
$\alpha \varphi(q)\cdot x^{3\ve/2}$ by
definition of $\alpha$. As a result, we get:
\[ \alpha \ll \frac{x^{1-3\ve/2}(\log x)^2}{\varphi(q)}
          \ll \frac{(\log x)^2\log\log x}{x^{\ve/2}}
\]
and hence $\alpha$ is negligible. Therefore, for all except at most a
negligible fraction of choices of $a\in(\Z/q\Z)^*$, we obtain that
$E(x;q,a) \leq x^{3\ve/4}$, and since $\pi(x)/\varphi(q) \gg x^\ve/\log
x$, this implies $\pi(x;q,a) = (1+o(1))\pi(x)/\varphi(q)$ as before. As a
result, under ERH, we obtain an analogue of Theorem~\ref{th:termfgm}
valid with overwhelming probability on the choice of $a$.
\begin{thm}
\label{th:termbasic-erh}
Assume ERH holds. Then
Algorithm~\ref{alg:basic} terminates with overwhelming probability.
Moreover, except for a negligible fraction of choices of the
class $a \bmod q$, it requires
$(1+o(1))\varphi(q)/q\cdot\log x$ iterations of the main loop on average,
and consumes
$ \big(\ve + o(1)\big)\cdot\frac{\varphi(q)}{q}\cdot\frac{(\log x)^2}{\log 2}
$
bits of randomness on average.
\end{thm}
Moreover, using Tur{\'a}n's theorem and the Cauchy--Schwarz inequality,
we can also establish under ERH alone the following analogue of
Theorem~\ref{th:statfgm}, regarding the output distribution of the
algorithm.
\begin{thm}
\label{th:stat1erh}
Assume ERH holds. Then the output distribution
of Algorithm~\ref{alg:basic} is statistically close to uniform, and its
collision entropy is no more than $O(\log\log x)$ bits smaller than that of the
uniform distribution.
\end{thm}
\begin{proof}
Algorithm~\ref{alg:basic} almost surely produces an output for a given
choice of $a$ if and only if $\pi(x;q,a)\neq0$, and this is no longer
certain under ERH. Therefore, the probability that the algorithm outputs
a prime $p = a + tq\leq x$ becomes:
\[ \Pr[X = a + tq] = \frac1{\varphi^*_x(q)}\cdot\frac1{\pi(x;q,a)}, \]
where $\varphi^*_x(q) = \#\{ a\in(\Z/q\Z)^*\ |\ \pi(x;q,a)\neq0 \}$. By
the previous discussion on the distribution of the values $\pi(x;q,a)$,
we know that $\varphi^*_x(q) = \varphi(q)\cdot\big(1-O(\alpha)\big)$.
As a result, a similar computation as in \S\ref{sec:fgm} gives:
\[ \Delta_2^2(X) = \frac1{\pi(x)^2}
     \sum_{\substack{a\in(\Z/q\Z)^* \\ \pi(x;q,a)\neq0}}
     \frac1{\pi(x;q,a)}\Big|\pi(x;q,a) - \frac{\pi(x)}{\varphi^*_x(q)}\Big|^2
   + \frac{\omega(q)}{\pi(x)^2},
\]
where $\omega(q)$ denotes as usual the number of prime factors of $q$.
Then, using the coarse lower bound $\pi(x;q,a)\geq1$ when
$\pi(x;q,a)\neq0$, we get:
\[ \Delta_2^2(X) \leq \frac{S^2 + \omega(q)}{\pi(x)^2} \]
where:
\begin{align*}
S  &= \sqrt{\sum_{a\in(\Z/q\Z)^*}
        \Big|\pi(x;q,a) - \frac{\pi(x)}{\varphi^*_x(q)}\Big|^2} \\
&\leq \sqrt{\sum_{a\in(\Z/q\Z)^*}
        \Big|\pi(x;q,a) - \frac{\pi(x)}{\varphi(q)}\Big|^2} +
      \sqrt{\sum_{a\in(\Z/q\Z)^*} \pi(x)^2
        \Big|\frac{\varphi(q) - \varphi^*_x(q)}{\varphi(q)\cdot
        \varphi^*_x(q)}\Big|^2} \\
&=    O(x^{1/2}\log x) + \sqrt{\varphi(q)\pi(x)^2
        \Big|\frac{\alpha}{1-\alpha}\cdot\frac1{\varphi(q)}\Big|^2} \\
&\ll  x^{1/2}\log x + \frac{\pi(x)}{\varphi(q)^{1/2}}\cdot\alpha
 \ll  x^{1/2}\log x + \frac{x\log x(\log\log x)^2}{x^{1/2-\ve/2}x^{\ve/2}}
 \ll  x^{1/2}\log x(\log\log x)^2
\end{align*}
by Tur{\'a}n's theorem again. Hence:
\[ \Delta_2^2(X) \ll \frac{x\log^{2+\ve} x}{\pi(x)^2} \ll
    \frac{\log^{3+\ve}
    x}{\pi(x)}.
\]
This is enough to obtain a relatively good bound on the collision entropy:
\[ H_2(X) = \log_2(\pi(x)) - \log_2(\log^{3+\ve} x) = (H_2)_{\max} - O(\log\log x) \]
but isn't sufficient for bounding the statistical distance. However, a
direct computation using Tur{\'a}n's theorem and the Cauchy--Schwarz
inequality is enough:
\begin{align*}
   \Delta_1(X)
&= \sum_{\substack{a\in(\Z/q\Z)^* \\ \pi(x;q,a)\neq0}} \pi(x;q,a)
      \Big|\frac1{\varphi^*_x(q)}\cdot\frac1{\pi(x;q,a)} - \frac1{\pi(x)}\Big| +
   \sum_{p|q} \frac1{\pi(x)} \\
&= \frac1{\pi(x)} \sum_{\substack{a\in(\Z/q\Z)^* \\ \pi(x;q,a)\neq0}}
      \Big|\pi(x;q,a) - \frac{\pi(x)}{\varphi^*_x(q)}\Big| +
   \frac{\omega(q)}{\pi(x)} \\
&\ll \frac1{\pi(x)}\cdot S\cdot\sqrt{\varphi^*_x(q)}
 \ll \frac{\log x}x\cdot x^{1/2}\log^2 x\cdot\sqrt{q}
 \ll \frac{\log^3 x}{x^{1/2}} \cdot x^{1/2-\ve/2}
 \ll \frac{\log^3 x}{x^{\ve/2}},
\end{align*}
which proves that the distribution is indeed statistically close to
uniform.\qed
\end{proof}

\subsection{Achieving almost sure termination under ERH}
\label{sec:as-erh}

Theorem~\ref{th:termbasic-erh} above is somewhat unsatisfactory, as we
have to ignore a negligible but possibly nonzero fraction of all values
$a\bmod q$ to obtain a bound on the average number of iterations and on
the randomness consumed by Algorithm~\ref{alg:basic} under ERH. But this
is unavoidable for that algorithm: as mentioned above, it is not known
whether ERH implies that for $q\propto x^{1-\ve}$, all $a\in(\Z/q\Z)^*$
satisfy $\pi(x;q,a)\neq 0$. And if an $a$ exists such that
$\pi(x;q,a)=0$, the choice of that $a$ in Step~\ref{algstep:basic-a} of
Algorithm~\ref{alg:basic}, however unlikely, is a case of
non-termination: as a result, the existence of such an $a$ prevents any
nontrivial bound on average running time or average randomness.

We propose to circumvent that problem by falling back to the trivial
algorithm (pick a random $p<x$, check whether it is prime and try again
if not) in case too many iterations of the main loop have been carried
out. This variant is presented as Algorithm~\ref{alg:erh}.

\begin{algorithm}[b]
\caption{A variant which terminates almost surely under ERH.}
\label{alg:erh}
\begin{algorithmic}[1]
  \State Fix $q \propto x^{1-\ve}$
  \State $a \randgets (\Z/q\Z)^*$
    \label{algstep:erh-a}
    \Comment{considered as an element of $\{1,\dots,q-1\}$}
  \Repeatx{$T=\log^2 x$ times}
    \State $t \randgets \{0,\dots,\lfloor \frac{x-a}q \rfloor\}$
    \State $p \gets a + t\cdot q$
    \Ifx{$p$ is prime}

      \textbf{return} $p$
  \EndRepeatx
  \Repeatx{forever}
    \State $p \randgets \{1,\dots,\lfloor x \rfloor\}$
    \Ifx{$p$ is prime}

      \textbf{return} $p$
  \EndRepeatx
\end{algorithmic}
\end{algorithm}

Clearly, since Algorithm~\ref{alg:erh} is the same as
Algorithm~\ref{alg:basic} except for the possible fallback to the trivial
algorithm, which has a perfectly uniform output distribution, the output
distribution of the variant is at least as close to uniform as the
original algorithm. In other words, the analogue of
Theorem~\ref{th:stat1erh} holds, with the same proof.
\begin{thm}
\label{th:stat2erh}
Assume ERH holds. Then the output distribution
of Algorithm~\ref{alg:erh} is statistically close to uniform, and its
collision entropy is no more than $O(\log\log x)$ bits smaller than that of the
uniform distribution.
\end{thm}
Moreover, as claimed above, we can obtain the following stronger analogue
of Theorem~\ref{th:termbasic-erh}.
\begin{thm}
\label{th:term2erh}
Assume ERH holds. Then
Algorithm~\ref{alg:erh} terminates almost surely, requires
$(1+o(1))\varphi(q)/q\cdot\log x$ iterations of the main loop on average,
and consumes
$ \big(\ve + o(1)\big)\cdot\frac{\varphi(q)}{q}\cdot\frac{(\log x)^2}{\log 2}
$
bits of randomness on average.
\end{thm}
\begin{proof}
Algorithm~\ref{alg:erh} terminates almost surely because the trivial
algorithm does. One can estimate its average number of iterations as
follows. Denote by $\varpi(t)$ the probability that
Algorithm~\ref{alg:erh} terminates after exactly $t$ iterations, and
$\varpi_a(t)$ the probability of the same event conditionally to $a$
being chosen in Step~\ref{algstep:erh-a}. We have:
\begin{align*}
  \varpi_a(t) &= \begin{cases}
     \displaystyle
     \bigg(1 - \frac{\pi(x;q,a)}{1+\lfloor\frac{x-a}q\rfloor}\bigg)^{t-1}
     \cdot\frac{\pi(x;q,a)}{1+\lfloor\frac{x-a}q\rfloor}
     & \text{for $t \leq T$;} \\
     \displaystyle
     \bigg(1 - \frac{\pi(x;q,a)}{1+\lfloor\frac{x-a}q\rfloor}\bigg)^{T}
     \bigg(1 - \frac{\pi(x)}{\lfloor x\rfloor}\bigg)^{t-T-1}
     \cdot\frac{\pi(x)}{\lfloor x\rfloor}
     & \text{otherwise.}
   \end{cases} \\
  \varpi(t) &= \frac1{\varphi(q)}\sum_{a\in(\Z/q\Z)^*} \varpi_a(t).
\end{align*}
Moreover, the expected number $N$ of iterations in
Algorithm~\ref{alg:erh} is given by $N = \sum_{t\geq 1} t\varpi(t)$.
We can denote by $N_a = \sum_{t\geq 1} t\varpi_a(t)$ the contribution of
a certain choice $a\in(\Z/q\Z)^*$.

Now, recall from \S\ref{sec:basic-erh} that $\pi(x;q,a)$ is within a
distance $\ll x^{3\ve/4}\log x$ of $\pi(x)/\varphi(q)$, except
for a fraction $\alpha\ll (\log x)^3/x^{\ve/2}$ of all possible
choices of $a$. If we
denote by $A$ the set of ``bad'' choices of $a$, we can write, for all
$a\in A$:
\[
  \varpi_a(t) \leq \begin{cases}
     \displaystyle
     1
     & \text{for $t \leq T$;} \\
     \displaystyle
     \Big(1 - \frac{\pi(x)}{\lfloor x\rfloor}\Big)^{t-T-1}
     \cdot\frac{\pi(x)}{\lfloor x\rfloor}
     & \text{otherwise.}
   \end{cases} \\
\]
Hence, if we let $\xi := \pi(x)/\lfloor x\rfloor$, we get:
\begin{align*}
 N_a &\leq \sum_{t=1}^T t + \sum_{t=T+1}^{+\infty} t(1-\xi)^{t-T-1} \xi
  = T(T+1)/2 + \sum_{k=1}^{+\infty} (T+k)(1-\xi)^{k-1}\xi \\
 N_a &\leq T(T+1)/2 + T\frac{\xi}{\xi}
+ \frac{\xi}{\xi^2} \leq T(T+3)/2 + \frac1\xi \ll \log^4 x.
\end{align*}
On the other hand, for $a\not\in A$, we have $\xi_a :=
\frac{\pi(x;q,a)}{1+\lfloor\frac{x-a}q\rfloor} =
\frac{q}{\varphi(q)}\cdot\frac{1+o(1)}{\log x}$. Therefore:
\begin{align*}
  N_a &= \sum_{t=1}^T t(1-\xi_a)^{t-1}\xi_a + (1-\xi_a)^T \sum_{t=T+1}^{+\infty}
t(1-\xi)^{t-T-1} \xi \\
  N_a &= \frac1{\xi_a} - \sum_{t=T+1}^{+\infty} t(1-\xi_a)^{t-1}\xi_a +
         (1-\xi_a)^T \sum_{t=T+1}^{+\infty} t(1-\xi)^{t-T-1} \xi \\
  \Big| N_a - \frac1{\xi_a} \Big| &\leq (1-\xi_a)^T
	 \sum_{k=1}^{+\infty} \Big[ (T+k)(1-\xi)^{k-1} \xi
	 + (T+k)(1-\xi_a)^{k-1} \xi_a \Big] \\
  \Big| N_a - \frac1{\xi_a} \Big| &\leq
	 \exp(-T\xi_a)\cdot(2T+1/\xi+1/\xi_a) \\
  \Big| N_a - \frac1{\xi_a} \Big| &\leq
	 \exp\Big(-(1+o(1))\frac{q}{\varphi(q)}\log x\Big)
	 \cdot(2T+1/\xi+1/\xi_a) \ll \frac1{x^{1-\ve}}.
\end{align*}
As a result, we obtain:
\begin{align*}
  N &= \frac1{\varphi(q)}\sum_{a\in(\Z/q\Z)^*} N_a
     = \Big(1-O(\log^3 x/x^{\ve/2})\Big)\cdot\Big(\frac1{\xi_a} +
       O(1/x^{1-\ve})\Big) + O(\log^3 x/x^{\ve/2})\cdot O(\log^4 x) \\
    &= \frac1{\xi_a} + O\Big(\frac{\log^7 x}{x^{\ve/2}}\Big)
     = (1+o(1))\frac{\varphi(q)}q\cdot \log x
 \end{align*}
as required. As for the expected number $R$ of random bits consumed by
the algorithm, it is given (ignoring the negligible amount necessary to
pick $a$) by:
\[ R = \frac{\log x}{\log 2}\bigg( \sum_{t=1}^T \ve t\cdot\varpi(t) +
\sum_{t=T+1}^{+\infty} (\ve T + t - T)\cdot\varpi(t) \bigg) \]
and the stated estimate is obtained by an exactly analogous computation.
\qed
\end{proof}

\subsection{An unconditional algorithm}
\label{sec:bdh}

Finally, we propose yet another variant of our algorithm for which both
almost sure termination and uniformity bounds can be established
unconditionally.
The idea is to no longer use a fixed modulus $q$, but to pick it
uniformly at random instead in the range $\{1,\dots,Q\}$ where $Q\propto
x(\log x)^{-A}$; uniformity bounds can then be deduced from the
Barban--Davenport--Halberstam theorem. Unfortunately, since $Q$ is only
polynomially smaller than $x$, we can no longer prove that the output
distribution is statistically close to uniform: the statistical distance
is polynomially small instead, with an arbitrarily large exponent
depending only on the constant $A$. On the other hand, termination is
obtained as before by falling back to the trivial algorithm after a
while, and since $q$ is often very close to $x$, we get an even better
bound on the number of consumed random bits.

\begin{algorithm}[t]
\caption{An unconditional variant.}
\label{alg:uncond}
\begin{algorithmic}[1]
  \State Fix $Q \propto x(\log x)^{-A}$ even
  \State $q \randgets \{Q/2+1,\dots,Q\}$
    \label{algstep:uncond_aq}
  \State $a \randgets \{0,\dots,q-1\}$
  \State \textbf{if} {$\gcd(a,q)\neq1$}
    \textbf{then goto} step \ref{algstep:uncond_aq}
    \label{algstep:uncond_aqcheck}

  \Repeatx{$T=\log^2 x$ times}
    \State $t \randgets \{0,\dots,\lfloor \frac{x-a}q \rfloor\}$
    \State $p \gets a + t\cdot q$
    \Ifx{$p$ is prime}

      \textbf{return} $p$
  \EndRepeatx
  \Repeatx{forever}
    \State $p \randgets \{1,\dots,\lfloor x \rfloor\}$
    \Ifx{$p$ is prime}

      \textbf{return} $p$
  \EndRepeatx
\end{algorithmic}
\end{algorithm}

Our proposed unconditional algorithm is described as
Algorithm~\ref{alg:uncond}. It picks the pair $(q,a)$ uniformly at random
among pairs of integers such that $q\in\{Q/2+1,\dots, Q\}$ and $a$ is a
standard representative of the classes in $(\Z/q\Z)^*$. There are:
\[ F(Q) := \sum_{Q/2<q\leq Q}\varphi(q) = \Phi(Q)-\Phi(Q/2) = \frac9{4\pi^2}Q^2 + O(Q\log Q) \]
possible such pairs, and we claim that for all except a polynomially
small fraction of them, $\pi(x;q,a)$ is close to $\pi(x)/\varphi(q)$.
Indeed, denote by $\alpha$ the fraction of all pairs $(q,a)$ such
that:
\[ E(x;q,a) := \Big|\pi(x;q,a) - \frac{\pi(x)}{\varphi(q)}\Big| >
   (\log x)^{3A/4}.
\]
Since $\pi(x)/\varphi(q) \gg x/Q \propto (\log x)^A$, we get
$\pi(x;q,a) = (1+o(1))\pi(x)/\varphi(q)$ for all pairs $(q,a)$ except a
fraction of at most $\alpha$. Moreover, we have the following trivial
lower bound:
\[ \sum_{Q/2 < q\leq Q} \sum_{a\in(\Z/q\Z)^*} E(x;q,a)^2 \geq
   \big[\alpha F(Q)\big]\cdot(\log x)^{3A/2}.
\]
On the other hand, the Barban--Davenport--Halberstam theorem ensures that
the sum on the left-hand side is $\ll xQ/\log Q$. As a result, we get:
\[
    \alpha
\ll \frac{(\log x)^{-3A/2}}{F(Q)} \cdot \frac{xQ}{\log Q}
\ll \frac{x(\log x)^{-3A/2}}{Q\log Q}
\ll \frac{x(\log x)^{-3A/2}}{x(\log x)^{-A+1}}
\ll \frac1{(\log x)^{A/2}}.
\]
This allows us to deduce the analogue of Theorem~\ref{th:term2erh} for
Algorithm~\ref{alg:uncond}.
\setcounter{section}{3}
\setcounter{subsection}{5}
\begin{thm}
\label{th:termuncond}
Algorithm~\ref{alg:uncond} with $A>6$ terminates almost surely, requires
$(1+o(1))\varphi(q)/q\cdot\log x$ iterations of the main loop on average,
and consumes:
\[ \big(A + o(1)\big)\cdot\frac{\varphi(q)}{q}\cdot\frac{\log x\log\log x}{\log 2}
\]
bits of randomness on average.
\end{thm}
\begin{proof}
Algorithm~\ref{alg:uncond} terminates almost surely because the trivial
algorithm does. One can estimate its average number of iterations as
follows. Denote by $\varpi(t)$ the probability that
Algorithm~\ref{alg:uncond} terminates after exactly $t$ iterations, and
$\varpi_{q,a}(t)$ the probability of the same event conditionally to the
pair $(q,a)$ being chosen in
Steps~\ref{algstep:uncond_aq}--\ref{algstep:uncond_aqcheck}. We have:
\begin{align*}
  \varpi_{q,a}(t) &= \begin{cases}
     \displaystyle
     \bigg(1 - \frac{\pi(x;q,a)}{1+\lfloor\frac{x-a}q\rfloor}\bigg)^{t-1}
     \cdot\frac{\pi(x;q,a)}{1+\lfloor\frac{x-a}q\rfloor}
     & \text{for $t \leq T$;} \\
     \displaystyle
     \bigg(1 - \frac{\pi(x;q,a)}{1+\lfloor\frac{x-a}q\rfloor}\bigg)^{T}
     \bigg(1 - \frac{\pi(x)}{\lfloor x\rfloor}\bigg)^{t-T-1}
     \cdot\frac{\pi(x)}{\lfloor x\rfloor}
     & \text{otherwise.}
   \end{cases} \\
  \varpi(t) &= \frac1{F(Q)}\sum_{Q/2< q\leq Q}\sum_{a\in(\Z/q\Z)^*}
\varpi_{q,a}(t).
\end{align*}
Moreover, the expected number $N$ of iterations in
Algorithm~\ref{alg:erh} is given by $N = \sum_{t\geq 1} t\varpi(t)$.
We can denote by $N_{q,a} = \sum_{t\geq 1} t\varpi_{q,a}(t)$ the contribution of
a certain choice $(q,a)$.

As we have just seen, $\pi(x;q,a) = (1+o(1))\pi(x)/\varphi(q)$, except
perhaps for a fraction $\alpha\ll (\log x)^{-A/2}$ of all choices of
$(q,a)$. If we denote by $A$ the set of ``bad'' choices of $(q,a)$, we
can write, as before, that for all $(q,a)\in A$:
\[
  \varpi_{q,a}(t) \leq \begin{cases}
     \displaystyle
     1
     & \text{for $t \leq T$;} \\
     \displaystyle
     \Big(1 - \frac{\pi(x)}{\lfloor x\rfloor}\Big)^{t-T-1}
     \cdot\frac{\pi(x)}{\lfloor x\rfloor}
     & \text{otherwise.}
   \end{cases} \\
\]
Hence, if we let $\xi = \pi(x)/\lfloor x\rfloor$, we again obtain:
\[
 N_{q,a} \leq \sum_{t=1}^T t + \sum_{t=T+1}^{+\infty} t(1-\xi)^{t-T-1} \xi
 \leq T(T+3)/2 + \frac1\xi \ll \log^4 x.
\]
On the other hand, for $(q,a)\not\in A$, we have $\xi_{q,a} :=
\frac{\pi(x;q,a)}{1+\lfloor\frac{x-a}q\rfloor} =
\frac{q}{\varphi(q)}\cdot\frac{1+o(1)}{\log x}$. Therefore, as before:
\[
  \Big| N_{q,a} - \frac1{\xi_{q,a}} \Big| \leq
	 \exp\Big(-(1+o(1))\frac{q}{\varphi(q)}\log x\Big)
	 \cdot(2T+1/\xi+1/\xi_{q,a}) \ll \frac1{x^{1-\ve}}.
\]
As a result, we get:
\begin{align*}
  N &= \frac1{F(Q)}\sum_{Q/2<q\leq Q}\sum_{a\in(\Z/q\Z)^*} N_{q,a}\\
    &= \Big(1-O\big((\log x)^{-A/2}\big)\Big)\cdot\Big(\frac1{\xi_{q,a}} +
       O(1/x^{1-\ve})\Big) + O\big((\log x)^{-A/2}\big)\cdot O(\log^2 x) \\
    &= \frac1{\xi_{q,a}} + O\big((\log x)^{4-A/2}\big)
     = (1+o(1))\frac{\varphi(q)}q\cdot \log x
\end{align*}
as required, since $4-A/2 < 1$.
As for the expected number $R$ of random bits consumed by
the algorithm, it is now given by:
\[ R = \frac{A\log\log x}{\log 2} \sum_{t=1}^T t\cdot\varpi(t) +
 \sum_{t=T+1}^{+\infty} \Big(T\cdot\frac{A\log\log x}{\log 2} + (t -
T)\cdot\frac{\log x}{\log 2}\Big)\cdot\varpi(t)
\]
where we have again ignored the random bits necessary to pick the pair $(q,a)$,
since we need only $\frac{Q(3Q+2)}{8F(Q)} \sim \pi^2/6$ iterations
of the loop from Step~\ref{algstep:uncond_aq} to
Step~\ref{algstep:uncond_aqcheck} to select it, and hence $O(\log x)$
random bits. The stated estimate is obtained by
essentially the same computation as for $N$.\qed
\end{proof}

We now turn to estimates on the uniformity of the output distribution of
the algorithm. For that purpose, we consider instead the output
distribution $X$ of Algorithm~\ref{alg:uncondnofall}, the variant of
Algorithm~\ref{alg:uncond} in which no bound is set to the number of
iterations of the main loop
(Steps~\ref{algstep:uncond-loopstart}--\ref{algstep:uncond-loopend}),
{i.e.} with no fallback to the trivial algorithm. Clearly, since the
trivial algorithm has a perfectly uniform output distribution, the
output distribution of Algorithm~\ref{alg:uncond} is at least as close to
uniform as $X$.

\begin{algorithm}[t]
\caption{An variant of Algorithm~\ref{alg:uncond} with no fallback.}
\label{alg:uncondnofall}
\begin{algorithmic}[1]
  \State Fix $Q \propto x(\log x)^{-A}$ even
  \State $q \randgets \{Q/2+1,\dots,Q\}$
  \State $a \randgets \{0,\dots,q-1\}$
  \State \textbf{if} {$\gcd(a,q)\neq1$}
    \textbf{then goto} step \ref{algstep:uncond_aq}
  \Repeatx{forever}
    \label{algstep:uncond-loopstart}
    \State $t \randgets \{0,\dots,\lfloor \frac{x-a}q \rfloor\}$
    \State $p \gets a + t\cdot q$
    \Ifx{$p$ is prime}

      \textbf{return} $p$
  \EndRepeatx
    \label{algstep:uncond-loopend}
\end{algorithmic}
\end{algorithm}

\begin{thm}
\label{th:statuncond}
The output distribution
of Algorithm~\ref{alg:uncondnofall} (and hence also
Algorithm~\ref{alg:uncond}) has a statistical distance $\Delta_1\ll
(\log x)^{(1-A)/2}$ to the uniform distribution. In particular, this
distance is polynomially small as long as $A>1$.
\end{thm}
\begin{proof}
Algorithm~\ref{alg:uncondnofall} produces an output (almost surely)
for exactly those choices of $(q,a)$ such that $\pi(x;q,a)\neq 0$. Let us
denote by $F^*_x(Q)$ the number of such choices. Clearly, with
notation from above, we have:
\[ F^*_x(Q) = \sum_{Q/2<q\leq Q} \varphi^*_x(q)\qquad\textrm{and}\qquad
   1-\alpha \leq \frac{F^*_x(Q)}{F(Q)} \leq 1.
\]
Then, the probability that Algorithm~\ref{alg:uncondnofall} outputs a
given prime $p\leq x$ can be written as:
\[ \Pr[X=p] = \frac1{F^*_x(Q)} \sum_{\substack{Q/2<q\leq Q\\ p\nmid q}}
     \frac1{\pi(x;q,p\bmod q)}.
\]
Therefore, we have:
\begin{align*}
   \Delta_1(X)
&= \sum_{p\leq x} \Big| \Pr[X=p] - \frac1{\pi(x)} \Big| \\
&= \sum_{p\leq x} \Bigg| \frac1{F^*_x(Q)} \sum_{\substack{Q/2<q\leq Q\\
   p\nmid q}} \frac1{\pi(x;q,p\bmod q)} - \frac1{\pi(x)} \Bigg| \\
&= \frac1{F^*_x(Q)} \sum_{p\leq x} \Bigg| \sum_{\substack{Q/2<q\leq Q\\
   p\nmid q}} \frac1{\pi(x;q,p\bmod q)} -
   \sum_{Q/2<q\leq Q}\frac{\varphi^*_x(q)}{\pi(x)} \Bigg| \\
&\leq \frac1{F^*_x(Q)} \sum_{p\leq x} \sum_{\substack{Q/2<q\leq Q\\
   p\nmid q}} \Big| \frac1{\pi(x;q,p\bmod q)} -
   \frac{\varphi^*_x(q)}{\pi(x)} \Big| + \frac1{F^*_x(Q)} \sum_{p\leq x}
   \sum_{\substack{Q/2<q\leq Q\\ p|q}} \frac{\varphi^*_x(q)}{\pi(x)} \\
&\leq \frac1{F^*_x(Q)} \sum_{Q/2<q\leq Q} \sum_{\substack{p\leq x\\
   p\nmid q}} \Big| \frac1{\pi(x;q,p\bmod q)} -
   \frac{\varphi^*_x(q)}{\pi(x)} \Big| + \frac1{F^*_x(Q)} \sum_{Q/2<q\leq
   Q} \frac{\omega(q)\varphi^*_x(q)}{\pi(x)} \\
&\leq \frac1{F^*_x(Q)\pi(x)} \sum_{\substack{Q/2<q\leq Q\\
   a\in(\Z/q\Z)^*
   }} \Big| \pi(x) -
   \varphi^*_x(q)\pi(x;q,a) \Big| + \frac{O(\log Q)}{\pi(x)}.
\end{align*}
We can then bound the sum over $(q,a)$ of $\big| \pi(x) -
\varphi^*_x(q)\pi(x;q,a) \big|$ as $D+D^*$, where:
\[ D = \sum_{\substack{Q/2<q\leq Q\\ a\in(\Z/q\Z)^* }}
         \Big| \pi(x) - \varphi(q)\pi(x;q,a) \Big|
   \quad\text{and}\quad
   D^* = \sum_{\substack{Q/2<q\leq Q\\ a\in(\Z/q\Z)^* }}
         \big| \varphi(q) - \varphi^*_x(q) \big|\cdot\pi(x;q,a).
\]
Now, on the one hand:
\begin{align*}
   D^*
&= \sum_{Q/2<q\leq Q} \big| \varphi(q) - \varphi^*_x(q) \big|
   \sum_{a\in(\Z/q\Z)^*} \pi(x;q,a) \\
&\leq \big( F(Q) - F^*_x(Q) \big) \cdot \pi(x)
 \leq \alpha F(Q)\pi(x)
 \ll  x^3(\log x)^{-5A/2-1},
\end{align*}
and on the other hand, applying the Cauchy--Schwarz inequality and the
Barban--Davenport--Halberstam theorem:
\begin{align*}
   D
&= \sum_{\substack{Q/2<q\leq Q\\ a\in(\Z/q\Z)^* }}
     \varphi(q)\cdot\Big| \pi(x;q,a) - \frac{\pi(x)}{\varphi(q)} \Big|
 \leq Q \sum_{\substack{Q/2<q\leq Q\\ a\in(\Z/q\Z)^* }}
          \Big| \pi(x;q,a) - \frac{\pi(x)}{\varphi(q)} \Big| \\
&\leq Q \cdot \sqrt{F(Q)} \cdot \sqrt{\frac{xQ}{\log Q}}
 \ll  x^3(\log x)^{-5A/2-1/2}.
\end{align*}
As a result, we obtain:
\[ \Delta_1(X) \ll \frac1{F^*_x(Q)\pi(x)}\cdot(D+D^*)
	\ll \frac{(\log x)^{2A+1}}{x^3}\cdot\frac{x^3}{(\log
	x)^{5A/2+1/2}} \ll \frac1{(\log x)^{(A-1)/2}}.
\]\qed

\end{proof}

\section{Comparison with other prime number generation algorithms}

In this section, we compare other prime number generation algorithms to
our method. In \S\ref{sec:why} we show that the output distribution of
Brandt and Damg\aa{}rd's PRIMEINC algorithm exhibits significant biases
(a fact that is intuitively clear, but which we make quite precise), and
in \S\ref{sec:maurer}, we discuss the advantages and drawbacks of our
method compared to that of Maurer, as the only previous work emphasizing
a close to uniform output distribution.

\subsection{PRIMEINC}
\label{sec:why}

Previous works on the generation of prime numbers, such as
\cite{DBLP:conf/crypto/BrandtD92,DBLP:conf/ches/JoyeP06},
provide a proof (based on rather strong assumptions) that the
output distribution of their algorithm has an entropy not much smaller
than the entropy of the uniform distribution. This is a reasonable
measure of the inability of an adversary to guess which particular prime
was output by the algorithm, but it doesn't rule out the possibility of
gaining some information about the generated primes. In particular, it
doesn't rule out the existence of an efficient distinguisher between the
output distribution and the uniform one.

Consider for example the PRIMEINC algorithm studied by Brandt and
Damg\aa{}rd in~\cite{DBLP:conf/crypto/BrandtD92}. In essence, it
consists in picking a random integer $y < x$ and returning the smallest
prime greater or equal to $y$. There are some slight technical
differences between this description and the actual PRIMEINC\footnote{To
wit, Brandt and Damg\aa{}rd restrict their attention to odd numbers
$x/2<y<x$, and set an upper bound to the number of iterations in the
algorithm}, but they have essentially no bearing on the following
discussion, so we can safely ignore them.

It is natural to suspect that the distribution of the output of this
algorithm is quite different from the uniform distribution: for example,
generating the second prime of a twin prime pair, i.e. a prime $p$ such
that $p-2$ is also prime, is abnormally unlikely. More precisely, if one
believes the twin prime conjecture, the proportion of primes of that form
among all primes up to $x$ should be $\sim 2c_2/\log x$, where
$c_2\approx 0.66$ is the twin prime constant. On the other hand, PRIMEINC
outputs such a $p$ if and only if it initially picks $y$ as $p$ or $p-1$.
Therefore, we expect the frequency of such primes $p$ in the output of
PRIMEINC to be much smaller, about $4c_2/(\log x)^2$. This provides an
efficient distinguisher between the output of PRIMEINC and the uniform
distribution.

More generally, it is easy to see that the method used
in~\cite{DBLP:conf/crypto/BrandtD92} to obtain the lower bound on the
entropy of the output distribution of PRIMEINC can also provide a
relatively large constant lower bound on the statistical distance to the
uniform distribution. Indeed, the method relies on the following result,
obtained by a technique first proposed by Gallagher~\cite{Gallagher}.

\def\mycite{\cite[Lemma 5]{DBLP:conf/crypto/BrandtD92}}
\begin{lemma}[\mycite]
\label{lem:gallagher}
Assume the prime $r$-tuple conjecture, and let $F_h(x)$ denote the number
of primes $p\leq x$ such that the largest prime $q$ less than $p$ satisfies
$p-q\leq h$. Then for any constant $\lambda$,
\[ F_{\lambda\log x}(x) = \frac{x}{\log
x}\big(1-e^{-\lambda}\big)(1+o(1)) \]
as $x\to+\infty$.
\end{lemma}

Now, the probability that the PRIMEINC algorithm outputs a fixed $p$ can
clearly be written as $d(p)/x$, where $d(p)$ is the distance between $p$
and the prime that immediately precedes it. Let $\Delta_1'$ be the
statistical distance between the output distribution of PRIMEINC and the
uniform distribution. We have:
\[ \Delta_1' = \sum_{p\leq x} \left|\frac{d(p)}{x} - \frac{1}{\pi(x)}\right|
     >\!\!\! \sum_{\substack{p\leq x \\ d(p)>2\log x}}
             \left|\frac{2\log x}{x} - \frac{\log x}{x}\right|
\]
for $x\geq 17$ in view of the already mentioned classical bound
$\pi(x)>x/\log x$ for $x\geq 17$. By Lemma~\ref{lem:gallagher}, this
gives:
\[ \Delta_1' > \frac{\log x}{x} F_{2\log x}(x) =
\big(1-e^{-2}\big)(1+o(1)) > 0.86 + o(1)
\]
as $x\to+\infty$, and in particular, $\Delta_1'$ admits a relatively large
constant lower bound (at least if one believes the prime $r$-tuple
conjecture on which the entropy bound is based).

Since the entropy lower bounds for other prime generation algorithms
like~\cite{DBLP:conf/ches/JoyeP06} are based on the same techniques, it
is very likely that their output distributions can similarly be shown to
be quite far from uniform.

\subsection{Maurer's algorithm}
\label{sec:maurer}

In~\cite{DBLP:conf/eurocrypt/Maurer89,DBLP:journals/joc/Maurer95}, Maurer
proposed a prime generation algorithm based on a similar principle as
Bach's technique to produce random numbers with known
factorization~\cite{DBLP:journals/siamcomp/Bach88}. This algorithm has
the major advantage of producing a primality certificate as part of its
output; as a result, it generates only provable primes, contrary to our
method.

Additionally, in the standard version of the algorithm, the distribution
of generated primes is heuristically close to uniform. More precisely,
this distribution would be exactly uniform if the following assertions
held:
\begin{enumerate}
\item \label{item:twoxplusone}
  The distribution of the relative sizes of the prime factors of
  an integer $x$ of given length \emph{conditional to} $2x+1$ being prime
  is the same as the distribution of the relative sizes of the prime
  factors of a uniformly random integer of the same length.

\item \label{item:dickson}
  That distribution is, moreover, the same as the asymptotic one (i.e. as
  the length goes to $+\infty$), as computed using Dickson's $\rho$
  function.
\end{enumerate}
Assertion~\ref{item:twoxplusone} is a rather nonstandard statement about
prime numbers, but Maurer points to some heuristic arguments that renders
it quite plausible \cite{Maurer92e}, at least asymptotically.
Assertion~\ref{item:dickson} is of course not exactly true, and it seems
difficult to quantify how much this affects the output distribution, but
probably not to a considerable extent.

That standard version of the algorithm, however, has a small probability
of not terminating, as discussed in
\cite[\S3.3]{DBLP:journals/joc/Maurer95}. To circumvent that problem,
Maurer suggests a modification to the algorithm which introduces biases
in the output distribution. If it is used, the algorithm will in
particular never generate primes $p$ such that $(p-1)/2$ has a large
prime factor. This provides an efficient distinguisher from the uniform
distribution (exactly analogous to the ``twin prime'' distinguisher of
the previous section), since e.g. safe primes have non negligible
density. There are other ways to avoid cases of non-termination that do
not affect the output distribution to such an extent (such as some form
of backtracking), but Maurer advises against them because of the large
performance penalty they may incur.

Furthermore, the paper \cite{DBLP:journals/joc/Maurer95} also describes
faster variants of the algorithm that are less concerned with the quality
of the output distribution, and they do indeed have outputs that are very
far from uniform, typically reaching only around 10\% of all prime
numbers.

More importantly, the main drawback of Maurer's algorithm compared to our
method is likely the amount of randomness it consumes: it is within a
small constant factor of the amount consumed by the trivial algorithm. In
fact, if we neglect everything but the last loop of the topmost step in
the recursion, and count as zero the cases when the recursion has more
than two steps, we see that the average amount of randomness used by the
algorithm is bounded below by $c$ times the corresponding amount for the
trivial algorithm, where (with the notations of
\cite[Appendix~1]{DBLP:journals/joc/Maurer95}):
\[ c \geq \int_{1/2}^1 \big(1-x\big)\, dF_1(x)
        = \int_{1/2}^1 \big(1-x\big)\, \frac{dx}{x}
        = \log 2 - \frac12
     \geq 0.19
\]
and we only counted a very small fraction of the actual randomness used!

As a result, in context where random bits are scarce and primality
certificates are not required, it seems that our method offers a rather
better trade-off than Maurer's algorithm: the uniformity of the
distribution is easier to estimate in our case, and we use much less
randomness. Moreover, Maurer finds
\cite[\S4.3]{DBLP:journals/joc/Maurer95} that an efficient implementation
of his algorithm is about 40\% slower than the trivial algorithm, whereas
Algorithm~\ref{alg:basic} is easily 5 times as fast in practice.

On the other hand, if proofs of primality are considered important, it is
clearly better to rely on Maurer's algorithm than to use our method and
then primality proving.

\bibliographystyle{abbrv}
\bibliography{biblio}

\iffalse
\subsection{Description}
\label{sec:methoddesc}

\subsubsection{Basic Method.}

A simple method to construct obviously uniformly distributed prime
numbers up to $2^n$ is to pick random numbers in $\{0,\dots,2^n-1\}$ and
retry until a prime is found. However, this method consumes $n$ bits of
randomness per iteration (not counting the amount of randomness consumed
by primality testing), and hence an expected amount of $n^2\log 2$
bits of randomness to produce a prime, which is quite large.

A simple way to reduce the amount of randomness consumed is to pick a
random odd number and to simply change the top $\ell$ bits until a prime
is found. This method, described in Algorithm~\ref{alg:basic}, consumes
only $\ell$ bits of randomness per iteration, and an expected
$(n-\ell)+(\ell n\log 2)/2$ bits overall. Despite the simplicity of this
method, we prove in the next section that its resulting distribution
achieves a statistical distance of less than about $2^{-\ell/2}$ to the
uniform distribution.

\begin{algorithm}
\caption{Our basic method}
\label{alg:basic}
\begin{algorithmic}[1]
  \State $b \randgets \{0,\ldots,2^{n-\ell-1}-1\}$
  \Repeat
  	\State $a \randgets \{0,\dots,2^{\ell}-1\}$
	\State $p \gets 2^{n-\ell}\cdot a + 2b + 1$
  \Until{$p$ is prime}
  \State \textbf{return} $p$
\end{algorithmic}
\end{algorithm}

Note that this method can be combined with efficient trial division to
speed up primality testing, since the reductions of the low-order bits
$(2b+1)$ modulo small primes can be precomputed. Thus, we only need to
carry out trial division on $a$ at each iteration of the \textbf{repeat}
loop.

\subsubsection{More efficient version.}

The previous method starts by fixing the reduction modulo $2^{n-\ell}$ of
the prime it searches for. But we can improve the odds of finding a prime
number by fixing the reduction $b$ modulo a number $m$ which is a product
of many small primes instead. The number $b$ should then be picked
uniformly at random among integers less than $m$ and coprime to $m$,
which can be achieved using the unit generation algorithm of Joye and
Paillier~\cite{DBLP:conf/ches/JoyeP06}.

Including the Joye-Paillier algorithm, the whole prime generation process
can be written as described in Algorithm~\ref{alg:efficient}. In this
setting, $m$ is usually picked as the product of all small primes up to a
certain bound $\beta$ chosen such that the bit size of $m$ is $n-\ell$
for some appropriate $\ell$.

\begin{algorithm}
\caption{More efficient method}
\label{alg:efficient}
\begin{algorithmic}[1]
  \State $b \randgets \{1,\ldots,m-1\}$
    \label{algstep:jpstart}
  \State $u \gets (1-b^{\lambda(m)}) \bmod m$
    \label{algstep:loopjp}
  \If{$u\neq 0$}
    \State $r\randgets \{1,\ldots,m-1\}$
    \State $b\gets b+ru \bmod m$
    \State \textbf{goto} step~\ref{algstep:loopjp}
  \EndIf
    \label{algstep:jpend}
  \Repeat
    \label{algstep:primstart}
    \State $a \randgets \{0,\dots,\lfloor 2^{n}/m \rfloor-1 \}$
    \State $p \gets am + b$
    \label{algstep:primsmallmult}
  \Until{$p$ is prime}
    \label{algstep:primend}
  \State \textbf{return} $p$
\end{algorithmic}
\end{algorithm}

In Algorithm~\ref{alg:efficient},
Steps~\ref{algstep:jpstart}--\ref{algstep:jpend} are the Joye-Paillier
unit generation technique: at the end of this loop, $b$ is distributed
exactly uniformly at random in $(\Z/m\Z)^*$. Then,
Steps~\ref{algstep:primstart}--\ref{algstep:primend} are exactly
analogous to Algorithm~\ref{alg:basic}, but are executed a much smaller number
of times on average.

More precisely, the probability that a number of the form $am+b<2^n$ is
prime is about $\frac{m}{\varphi(m)n\log 2}$ according to the prime
number theorem for arithmetic progressions. Therefore, the expected
number of primality testing loops in this case is only:
\[ n\log 2\cdot\frac{\varphi(m)}{m} =
   n\log 2\prod_{p\leq\beta}\left(1-\frac1p\right) \sim
   \frac{n\log 2\cdot e^{-\gamma}}{\log\beta}
\]
by Mertens' formula (where $\gamma$ is the Euler-Mascheroni constant).
Now, we have $\log m = \sum_{p\leq\beta}\log p \sim\beta$,
so that $\beta\approx(n-\ell)\log 2$. Thus, the number of primality
testing loops is down to about:
\[ \frac{n\log 2\cdot e^{-\gamma}}{\log(n-\ell)+\log\log 2} \]
which, for typical parameter sizes ($512\leq n\leq 3072$ and
$32\leq\ell\leq 256$, say), is around $5$ to $7$ times fewer than with the
previous approach. Accordingly, this part of the algorithm consumes $5$ to
$7$ times less randomness, namely about:
\[ \frac{\ell n\log 2\cdot e^{-\gamma}}{\log(n-\ell)+\log\log 2}\text{
bits overall.} \]
This technique can also lead to significantly faster prime generation
compared to an approach which uses trial division to reduce calls to the
more costly primality testing algorithm, depending on how efficient
modular reductions are in the environment of interest---for example, Joye
and Paillier report speed-ups of one order of magnitude on smart cards
using their algorithm~\cite[Fig. 2]{DBLP:conf/ches/JoyeP06}, which
is slightly more computationally demanding
than ours since it requires a large modular multiplication at each
iteration, instead of the small integer multiplication from
Step~\ref{algstep:primsmallmult}.

On the other hand, the unit generation part of the algorithm does consume
more randomness than the single-step generation of $b$ in
Algorithm~\ref{alg:basic}. This has limited influence, however, since the
number of iterations in this part is quite small. Indeed, it is easy to
see that $b$ is relatively prime to some $p|m$ as soon as one of the
random elements of $\Z/m\Z$ returned by the random number generator is
itself coprime to $p$ (more precisely: the assignment $b\gets b+ur$
preserves coprimality with $p$, since $(b,p)=1$ implies $u\equiv 0\pmod
p$; but when $b\equiv 0 \pmod p$, then $u$ is invertible mod $p$, and
thus $b$ becomes non zero mod $p$ if $r$ is). As a result, the
probability that this part of the algorithm takes $\leq k$ iterations
(including the initial random assignment of $b$) is exactly $z_k =
\prod_{p\leq\beta} (1-p^{-k})$, which converges rapidly to $1$. The
expected number of iterations is thus $E_\beta = \sum_{k\geq 1} k(z_k-z_{k-1})$
which we can evaluate numerically for any given $\beta$. An upper bound
is given by the limit as $\beta\to+\infty$, namely:
\begin{equation}
\label{eq:einfty}
 E_{\infty} = \sum_{k\geq 2} k\left(\frac1{\zeta(k)}-\frac1{\zeta(k-1)}\right)
   \leq 2.71
\end{equation}
This means in particular that the amount of randomness required by unit
generation, $(n-\ell)\cdot E_\beta$, is reasonably small (and typically small
enough that this method is better than the previous one in terms of
random bits).

An added benefit of this approach is that the bound on the statistical
distance to the uniform distribution is expressed in terms of
$\varphi(m)$ rather than $m$ itself: this is smaller when $m$ is a
product of many small primes than when $m$ is a power of two as before.

\subsection{Main result}

Both of the previous methods generate primes of the form $am+b$ where
$b$ is first picked uniformly at random among integers less than $m$ and
coprime to $m$: in the first method, $m$ is a power of two, whereas in
the second method, it is a product of small primes, but the principle
remains the same.

We now prove that the distribution of primes generated in this fashion is
close to the uniform distribution among all primes in the allowable
interval $\{1,\dots,x-1\}$ with $x = m\cdot\lfloor 2^{n}/m \rfloor
\approx 2^n$.

\begin{theorem}
\label{th:mainthm}
Assume the Extended Riemann Hypothesis. Suppose further that $x\geq 17$
and $x/m \geq e$ (the base of natural logarithms). Then
the statistical distance $\Delta$ between the uniform distribution over
primes less than $x$ and the distribution resulting from the previous
method is bounded as:
\[ \Delta < 3\log^2 x \cdot \sqrt{\frac{\varphi(m)}{x}} \]
\end{theorem}

The rest of this section is devoted to the proof of this theorem.
As usual, note $\pi(x)$ the number of primes up to $x$, and $\pi(x,m,b)$
the number of primes $p$ up to $x$ such that $p\equiv b\pmod m$.
Furthermore, for any Dirichlet character $\chi$ modulo $m$, let:
\[ \pi(x,\chi) = \sum_{p\leq x} \chi(p) \]
(where the sum is over prime numbers $p$, as will be the case for all
sums over $p$ henceforth). Then clearly, for any $b$ coprime to $m$, we
have:
\[ \pi(x,m,b) = \frac{1}{\varphi(m)}
                \sum_{\chi \bmod m} \overline{\chi(b)}\cdot\pi(x,\chi)
\]
and in this sum, the contribution of the trivial character $\chi_0$ is:
\[ \frac{1}{\varphi(m)} \pi(x,\chi_0) = \frac{\pi(x)-s(m)}{\varphi(m)} \]
where $s(m)$ is the sum of all distinct prime divisors of $m$. Hence:
\begin{equation}
\label{eq:sumnontrivchar}
\frac{1}{\varphi(m)} \sum_{\chi\neq\chi_0} \overline{\chi(b)}\cdot\pi(x,\chi)
= \pi(x,m,b) - \frac{\pi(x)-s(m)}{\varphi(m)}
\end{equation}
Therefore, we obtain the following result.

\begin{lemma}
\label{lem:statdistascharsum}
The statistical distance $\Delta$ satisfies:
\[ \Delta \leq \frac{1}{\pi(x)} \Bigg( 2s(m) + \frac{1}{\varphi(m)}
               \sum_{b\in(\Z/m\Z)^*} \bigg| \sum_{\chi\neq\chi_0}
               \overline{\chi(b)}\pi(x,\chi) \bigg| \Bigg)
\]
\end{lemma}
\begin{proof}
Indeed, the probability that a given prime $p\leq x$ is chosen by our
algorithm is $0$ if $p$ divides $m$, and $1/(\varphi(m)\pi(x,m,b))$ if
$p\equiv b\pmod m$ for $b$ coprime to $m$ ($b$ is first chosen uniformly
among the $\varphi(m)$ classes in $(\Z/m\Z)^*$, and $p$ is then chosen
uniformly among the $\pi(x,m,b)$ possible primes of the form $am+b$).
Thus:
\begin{align*}
  \Delta &= \frac{s(m)}{\pi(x)} + \sum_{b\in(\Z/m\Z)^*} \pi(x,m,b)\cdot
            \left|\frac1{\varphi(m)\pi(x,m,b)} - \frac1{\pi(x)}\right| \\
         &= \frac1{\pi(x)} \Bigg(s(m) + \sum_{b\in(\Z/m\Z)^*}
            \left|\frac{\pi(x)}{\varphi(m)} - \pi(x,m,b)\right| \Bigg) \\
         &= \frac1{\pi(x)} \Bigg(s(m) + \sum_{b\in(\Z/m\Z)^*}
            \bigg|\frac{s(m)}{\varphi(m)} - \frac1{\varphi(m)}
            \sum_{\chi\neq\chi_0} \overline{\chi(b)}\pi(x,\chi) \bigg|
            \Bigg)
         & \text{by \eqref{eq:sumnontrivchar}} \\
         &\leq \frac{1}{\pi(x)} \Bigg( 2s(m) + \frac{1}{\varphi(m)}
               \sum_{b\in(\Z/m\Z)^*} \bigg| \sum_{\chi\neq\chi_0}
               \overline{\chi(b)}\pi(x,\chi) \bigg| \Bigg)
\end{align*}
as required. \qed
\end{proof}

We can now estimate the sum over $b\in(\Z/m\Z)^*$ using the Cauchy-Schwarz
inequality. Indeed, let:
\[
  S(x) = \frac{1}{\varphi(m)}\sum_{b\in(\Z/m\Z)^*} \bigg|
           \sum_{\chi\neq\chi_0} \overline{\chi(b)}\pi(x,\chi) \bigg|
\]
Applying Cauchy-Schwarz, we get:
\begin{align*}
  S(x) &\leq \frac{1}{\varphi(m)}\cdot\sqrt{\varphi(m)}
             \sqrt{\sum_{b\in(\Z/m\Z)^*} \bigg|
             \sum_{\chi\neq\chi_0} \overline{\chi(b)}\pi(x,\chi)
             \bigg|^2} \\
       &\leq \frac{1}{\sqrt{\varphi(m)}} \sqrt{\sum_{b\in(\Z/m\Z)^*}
             \sum_{\chi,\chi'\neq\chi_0} \overline{\chi(b)}\chi'(b)\cdot
             \pi(x,\chi)\overline{\pi(x,\chi')}}
\end{align*}
But the sum over $b\in(\Z/m\Z)^*$ of $\overline{\chi(b)}\chi'(b)$ is
$\varphi(m)$ if $\chi=\chi'$ and vanishes otherwise. Thus:
\begin{equation}
  \label{eq:sbound}
  S(x) \leq \sqrt{\sum_{\chi\neq\chi_0} \big|\pi(x,\chi)\big|^2}
\end{equation}
We can thus derive a bound on $S(x)$, and hence on $\Delta$, from an
estimate of $\pi(x,\chi)$ for non trivial Dirichlet characters $\chi$.
Such an estimate can be obtained by locating the zeros of the
$L$-function $L(s,\chi)$ in the critical strip (see e.g. \cite[Chapter
20]{davenport}). In particular, a well-known consequence of the Extended
Riemann Hypothesis, due to Titchmarsh, is that $\pi(x,\chi) =
O\big(\sqrt{x}\log (mx)\big)$, where the implied constant is absolute.
One can actually give an effective version of this estimate:

\begin{lemma}[Oesterl\'e, \cite{oesterle}]
\label{lem:oesterle}
Assuming the Extended Riemann Hypothesis, we have, for all $x\geq 2$ and
all non trivial Dirichlet characters $\chi$ modulo $m$:
\[ \big|\pi(x,\chi)\big| \leq \sqrt{x} (\log x + 2\log m) \]
\end{lemma}

Therefore, under ERH, estimate~\eqref{eq:sbound} becomes:
\[ S(x) \leq \sqrt{\varphi(m) x}(\log x + 2\log m)
\]
since $m\leq x$. Putting this together with the result of
Lemma~\ref{lem:statdistascharsum}, we get:
\begin{align*}
  \Delta &\leq \frac{2s(m) + \big(3\log x - 2\log(x/m)\big)
                     \sqrt{\varphi(m) x}}{\pi(x)} \\
         &\leq \frac{3\log x\sqrt{\varphi(m) x}}{\pi(x)} +
               \frac{2}{\pi(x)}\big(s(m) - \sqrt{e\cdot\varphi(m)m}\big) \\
         &\leq \frac{3\log x\sqrt{\varphi(m) x}}{\pi(x)}
\end{align*}
since it is easy to check that $s(m)\leq \sqrt{e\cdot\varphi(m)m}$.
Finally, a classical estimate \cite[Th. 8.8.1]{bachshallit} states that
$\pi(x)>x/\log x$ for $x\geq 17$. Therefore, we obtain the stated bound
on $\Delta$:
\[ \Delta < 3\log^2 x \cdot \sqrt{\frac{\varphi(m)}{x}} \]
\qed

\begin{remark}
\begin{enumerate}
\item \label{rmk:termination}
Perhaps surprisingly, while the Extended Riemann Hypothesis provides
strong enough bounds to obtain this non trivial estimate for the
statistical distance, it is not actually sufficient to prove that the
algorithm as written above \emph{terminates} almost surely\footnote{The
statistical distance estimated above is between the uniform distribution
and the ``a priori'' output distribution of the algorithm, determined by the
probability over all random coin flips of returning a particular number
(or $\bot$ if a certain sequence of coin flips corresponds to a case of
non-termination).}.  Indeed, the
best bound that can be established with ERH (e.g. using
Lemma~\ref{lem:oesterle}) for the smallest prime $p$ in the arithmetic
progression $am+b$ is of the form $p < C\cdot(m\log m)^2$ for some
constant $C$. In particular, ERH isn't sufficient to
prove that $\pi(x,m,b)\neq 0$ unless $x \gg m^2$, which is not usually
the case in our setting, where we would like to choose $m$ only a few bits
shorter than $x$.

However, a conjecture of Heath-Brown~\cite{HeathBrown78}, for which
Wagstaff~\cite{wagstaff} gave several heuristic and experimental
arguments, ensures that the algorithm terminates on any
input provided that $m$ is chosen such that
$x>(2+\varepsilon)\cdot\varphi(m)\log^2 m$. This is satisfied, in
particular, whenever the bound given in Theorem~\ref{th:mainthm} is non
trivial!

Since we would like to avoid relying on a conjecture stronger than ERH,
though, we suggest another workaround in \S\ref{sec:termination}: we
simply modify the algorithm to fall back on the trivial random search
if no prime is found in a certain number $T$ of iterations. The modified
algorithm clearly terminates almost surely, its output distribution
is as least as close to uniform as before, and it doesn't use many more
random bits if $T$ is large enough.

\item
The analysis above assumes that the primality testing algorithm
used in the algorithm checks for exact primality. While a deterministic
primality test such as AKS could in principle be used,
efficiency reasons make probabilistic algorithms like Miller-Rabin, which
have a small probability of returning \texttt{true} on input of a
composite number, much more suitable for actual implementations. However,
the number of Miller-Rabin rounds carried out can be chosen to ensure that this
probability of erroneously recognizing a composite as prime is
negligible. This means that the actual output distribution is in fact
statistically close to the ideal distribution for which the primality
testing algorithm is perfect. The ``real'' $\Delta$ is thus at most
negligibly larger than the ideal one considered in this section.

\item
For simplicity, we have only considered the generation of primes smaller
than a given bound which is a multiple of $m$, but it is not difficult to
extend our algorithms to the generation of primes in a fixed interval,
using the same approach as~\cite{DBLP:conf/ches/JoyeP06}.
\end{enumerate}
\end{remark}

\subsection{Achieving termination under ERH}
\label{sec:termination}

As mentioned in Remark~\ref{rmk:termination}, ERH is not sufficient to
establish that our algorithm terminates almost surely. A simple
modification can solve that problem however: fall back on the trivial
algorithm (picking random numbers and checking them for primality) if no
prime has been found before a certain number $T$ of iterations.

Since the trivial algorithm terminates almost surely, the same is true
for this modified algorithm. Moreover, the statistical distance
$\Delta_T$ of the output distribution of the modified algorithm to
uniform is surely no larger than before, since we fall back to an
algorithm whose output distribution is exactly uniform. More precisely,
we have
\begin{align*}
  \Delta_T &= \frac{s(m)}{\pi(x)}\cdot \sum_{b\in(\Z/m\Z)^*}
  \frac{\varpi_T(b)}{\varphi(m)} + \sum_{b\in(\Z/m\Z)^*} \pi(x,m,b)\cdot
            \left|\frac{\varpi_T(b)}{\varphi(m)\pi(x,m,b)} +
\frac{1-\varpi_T(b)}{\pi(x)} - \frac1{\pi(x)}\right| \\
           &= \sum_{b\in(\Z/m\Z)^*} \varpi_T(b) \left(
\frac{s(m)}{\pi(x)\varphi(m)} + \pi(x,m,b)\cdot
            \left|\frac1{\varphi(m)\pi(x,m,b)} - \frac1{\pi(x)}\right|
\right)
\end{align*}
where we denote by $\varpi_T(b)$ the probability that a prime number is
found in $T$ steps or less in the class of $b\bmod m$, i.e.
\[ \varpi_T(b) = 1 - \left(1-\frac{\pi(x,m,b)}{x/m}\right)^T. \]
Since $\varpi_T(b)\leq 1$, it follows that $\Delta_T \leq \Delta$.
\fi

\iffalse
\section{Performance comparison}
\label{sec:perf}

Let us simplify and summarize the performance estimates of
\S\ref{sec:methoddesc} by giving asymptotic formulas for the
efficiency of the various algorithms mentioned above\footnote{Except
Maurer's algorithm, which has a very different structure than the other
techniques, so that trying to fit it in those tables seems to make little
sense. See \S\ref{sec:maurer} for an extended comparison.} in terms of number
of primality tests, random bits used, and statistical distance to the
uniform distribution. These formulas are collected in
Table~\ref{tab:effasympt}.

\begin{table}[p]
\begin{center}
\begin{tabular}{r|c|c|c|c|c}
 & Trivial & PRIMEINC & Joye-Paillier &
Algorithm~\ref{alg:basic} & Algorithm~\ref{alg:efficient} \\
\hline
Expected primality tests  & $\frac{\log 2}2\cdot n$  & $\frac{\log
2}2\cdot n$ & $e^{-\gamma}\log 2\cdot \frac{n}{\log n}$ & $\frac{\log 2}2\cdot n$
& $e^{-\gamma}\log 2\cdot \frac{n}{\log n}$ \\
\hline
Expected random bits & $\frac{\log 2}2\cdot n^2$ & $n$ & $E_\infty\cdot
n$ & $\frac{\log 2}2\cdot n\ell$ & $e^{-\gamma}\log 2\cdot
\frac{n\ell}{\log n}$ \\
\hline
Statistical distance & $= 0$ & $\gtrsim 0.86$ & $>\text{constant?}$
& $< 1.02 n^2\cdot 2^{-\ell/2}$ & $< \frac{1.09n^2}{\log n}\cdot
2^{-\ell/2}$
\end{tabular}
\end{center}
\caption{Comparison of the asymptotic efficiencies of the previous
algorithms (assuming $1\ll \ell\ll n$). The constant $E_\infty\approx
2.7$ is defined in Equation~\eqref{eq:einfty}.}
\label{tab:effasympt}
\end{table}

Concretely speaking, we can precisely compute those efficiency parameters
for the generation of, say, $512$ or $1024$-bit primes. We report those
numbers in Table~\ref{tab:eff512} and \ref{tab:eff1024}, together with
timings for straightforward implementations of these algorithms using
PARI/GP~\cite{pari}.

Note that there is, of course, a trade-off between the number of random
bits used and the uniformity of the distribution: both increase with
$\ell$ (resp. $m$). A reasonable compromise is probably to choose $\ell$
just under the number of bits produced by each call to the random number
generator of the system, so that a single word of randomness is produced
for each iteration of the algorithm. In our case, we picked $\ell=64$ in
Algorithm~\ref{alg:basic}, and chose $m$ as the product of all primes up
to $337$ or $691$ (for $512$ and $1024$ bits) in
Algorithm~\ref{alg:efficient} ($\lfloor 2^n/m \rfloor$ is then $62$ bit
long in both cases).

\begin{table}[p]
\begin{center}
\begin{tabular}{r|c|c|c|c|c}
 & Trivial & PRIMEINC & Joye-Paillier &
Algorithm~\ref{alg:basic} & Algorithm~\ref{alg:efficient} \\
\hline
Expected primality tests  & 177 & 177 & 35 & 177 & 35 \\
\hline
Expected random bits & 90,852 & 512 & 1098 & 11,805 & 3320 \\
\hline
Statistical distance & $= 0$ & $\gtrsim 0.86$ & $>\text{constant?}$
& $< 2^{-12.9}$ & $<2^{-13.3}$ \\
\hline
Average CPU time & 0.19~s & 0.20~s & 0.04~s & 0.19~s & 0.04~s
\end{tabular}
\end{center}
\caption{Practical efficiency of the various algorithms for the
generation of $512$-bit primes.}
\label{tab:eff512}
\end{table}

\begin{table}[p]
\begin{center}
\begin{tabular}{r|c|c|c|c|c}
 & Trivial & PRIMEINC & Joye-Paillier &
Algorithm~\ref{alg:basic} & Algorithm~\ref{alg:efficient} \\
\hline
Expected primality tests  & 355 & 355 & 61 & 355 & 61 \\
\hline
Expected random bits & 363,409 & 1024 & 2405 & 22,965 & 6205 \\
\hline
Statistical distance & $= 0$ & $\gtrsim 0.86$ & $>\text{constant?}$
& $< 2^{-10.9}$ & $<2^{-11.4}$ \\
\hline
Average CPU time & 2.03~s & 2.02~s & 0.44~s & 1.92~s & 0.41~s
\end{tabular}
\end{center}
\caption{Practical efficiency of the various algorithms for the
generation of $1024$-bit primes.}
\label{tab:eff1024}
\end{table}

As we can see, our algorithms are as efficient as comparable ones in
terms of speed, and achieve provable bounds on the statistical distance
at a much smaller cost than the only other such example, the trivial
uniform generator. In fact, Algorithm~\ref{alg:efficient} is as fast as
the Joye-Paillier algorithm and has provably good statistical properties
at a mere threefold increase in random bits requirement (as opposed to
thousandfold for the trivial uniform generator).
\fi

\end{document}


TITLE: Prove $ \int_0^t 2X_s \ dX_s = X_t^2-X_0^2-\langle X, X\rangle_t $ WITHOUT Ito's formula
QUESTION [1 upvotes]: Suppose $X$ is a continuous local martingale (on some complete filtered probability space). Then I want to prove that
$$
\int_0^t 2X_s \ dX_s = X_t^2-X_0^2-\langle X, X\rangle_t
$$
where $\langle X , X \rangle$ denotes the quadratic variation process of $X$.
One way to go about this is using Ito's formula but I want to see if I can give a direct proof of this just using properties of the stochastic integral.

What I have thought about so far is that by definition of the quadratic variation process the above yields that
$$
\int_0^t 2X_s \ dX_s - X_0^2
$$
is a continuous local martingale. Maybe this could be used along with the fact that the quadratic variation process is unique to show the relation?
It almost also looks like that $(a+b)^2=a^2+b^2+2ab$ could play a role. Maybe this would lead to too many terms but some of them could be indistinguashable from 0 (by showing they were both a continuous local martingale and a finite variation process).

The textbook I have is "Brownian Motion, Martingales, and Stochastic Calculus" by Jean-Francois Le Gall so any result prior to the section on Ito's formula is available.

REPLY [3 votes]: It sounds like you were almost there. Just use
\begin{align}\tag{1}
\int_0^tX_s\,dX_s&=\lim_{\max\limits_i|t_i-t_{i-1}|\to 0}\sum_{i=1}^n X_{t_{i-1}}(X_{t_i}-X_{t_{i-1}})
\end{align}
and
\begin{align}
\langle X,X\rangle_t&=\lim_{\max\limits_i|t_i-t_{i-1}|\to 0}\sum_{i=1}^n(X_{t_i}- X_{t_{i-1}})^2\\\tag{2}
&=\lim_{\max\limits_i|t_i-t_{i-1}|\to 0}\sum_{i=1}^nX_{t_i}(X_{t_i}- X_{t_{i-1}})-\int_0^tX_s\,dX_s\,.
\end{align}
The sum in the first term of (2) is
\begin{align}
\sum_{i=1}^nX_{t_i}(X_{t_i}- X_{t_{i-1}})&=\sum_{i=1}^nX_{t_i}^2- X_{t_{i-1}}^2+X_{t_{i-1}}^2-X_{t_i}X_{t_{i-1}}\\&=X_t^2-X_0^2-\sum_{i=1}^nX_{t_{i-1}}(X_{t_i}- X_{t_{i-1}})\,.
\end{align}
By (1) this is $X_t^2-X_0^2-\int_0^tX_s\,dX_s$ in the limit.
Altogether we have
$$
\langle X,X\rangle_t=X_t^2-X_0^2-2\int_0^tX_s\,dX_s\,.
$$


TITLE: $q$-series identity
QUESTION [4 upvotes]: I have to prove the following identity: $$\sum_{n\geq 0} (-1)^n(2n+1)q^{\frac{n(n+1)}{2}} = (q;q)_\infty^3$$ where $(a;q)_\infty = \prod_{i\geq 0}1-aq^i$ is the $q$-Pochhammer symbol. In my notes the identity is stated as a corollary of Jacobi triple product, whose statement is $$\sum_{n\in\mathbb Z} z^nq^{\frac{n(n+1)}{2}} = (-zq;q)_\infty(-1/z;q)_\infty(q;q)_\infty$$
What I have done so far is to change the range of $n$ to $\mathbb Z$, noting that
\begin{align*}\sum_{n\geq 0} (-1)^n(2n+1)q^{\frac{n(n+1)}{2}} &= \sum_{n\in\mathbb Z} (-1)^n n q^{\frac{n(n+1)}{2}} \\
&= \sum_{n\in\mathbb Z} (-1)^n \frac{2n+n^2-n^2}{2} q^{\frac{n(n+1)}{2}} \\
&= \sum_{n\in\mathbb Z} (-1)^n \frac{n^2+n}{2} q^{\frac{n(n+1)}{2}} + \sum_{n\in\mathbb Z} (-1)^n \frac{n-n^2}{2} q^{\frac{n(n+1)}{2}} \\
&= q\left(\sum_{n\in\mathbb Z} (-1)^n q^{\frac{n(n+1)}{2}}\right)' + \sum_{n\in\mathbb Z} (-1)^{n+1} \frac{n(n-1)}{2} q^{\frac{n(n+1)}{2}} \\
&= 0 + \sum_{n\in\mathbb Z} (-1)^{n+1} \frac{n(n-1)}{2} q^{\frac{n(n+1)}{2}} \end{align*}
where the first sum is zero, thanks to Jacobi triple product with $z=-1$.
However, I can't tell if I'm getting closer to my goal. Any hint would be appreciated at this point!
As a side question, I think I'm not getting the Jacobi triple product properly. In the LHS, there are negative powers of $z$. Does this mean that the equality is not happening in $\mathbb C[[z,q]]$ but only for values of $z,q$ where the sum converges? If so, what is the region of convergence of the series?

REPLY [7 votes]: The theorem is due to Jacobi.
An english-language demonstration can be found in
chapter 19.9, at theorem 357 (with typos) by

G. H. Hardy and E. M. Wright: An introduction to the theory of numbers.
Oxford University press, sixth edition 2008, ISBN 978-0-19-921985-8.
Note: First edition 1938.
The fourth edition seems to be preferable, the sixth has shameful OCR artefacts.

The proof goes like this.
Jacobi's triple product identity, as quoted in the question:
$$\sum_{n\in\mathbb{Z}} z^n q^{n(n+1)/2} =
(-zq;q)_\infty (-\tfrac{1}{z};q)_\infty (q;q)_\infty$$
To get somewhere near $(q;q)_\infty^3$, we want to set $z=-1$.
The problem is then the second product:
Its first factor is $(1+\tfrac{1}{z})$ which becomes zero for $z=-1$.
So let us pull this first factor out:
$$\begin{align}
 (-\tfrac{1}{z};q)_\infty &=
 (1+\tfrac{1}{z}) (-\tfrac{1}{z}q;q)_\infty
\\\therefore\quad
 (-zq;q)_\infty (-\tfrac{1}{z}q;q)_\infty (q;q)_\infty
 &= \frac{\sum_{n\in\mathbb{Z}} z^n q^{n(n+1)/2}}{1+\tfrac{1}{z}}
\\ &= \sum_{n=0}^\infty \frac{z^n+z^{-n-1}}{1+\tfrac{1}{z}} q^{n(n+1)/2}
\\ &= \sum_{n=0}^\infty z^{-n}\frac{z^{2n+1}+1}{z+1} q^{n(n+1)/2}
\\ &= \sum_{n=0}^\infty z^{-n}\,(z^{2n}-z^{2n-1}\pm\cdots-z+1)\,
 q^{n(n+1)/2}
\end{align}$$
Thus, the singularity at $z=-1$ can be removed.
Setting $z=-1$ now gives
$$(q;q)_\infty^3 = \sum_{n=0}^\infty (-1)^n (2n+1)\,q^{n(n+1)/2}$$
which is what we wanted to prove.


TITLE: Can we symmetrize connections on tangent bundles?
QUESTION [5 upvotes]: Let $M$ be a smooth manifold. Denote by $\mathcal{A}$ the space of all affine connections on $TM$, and by $\mathcal{S}$ the affine subspace of all the symmetric connections.
Is there a projection $P:\mathcal{A}\to \mathcal{S}$? (i.e $P^2=P$).
If there is more than one, is there some natural choice?
Note that $S$ is convex, so perhaps introducing some metric on $\mathcal{A}$ will be useful?

REPLY [2 votes]: The naturality statement in the end of the answer by @AnthonyCarapetis can indeed be made precise: There is a natural action of local diffemorphisms on affine connections (basically characterized by $(f^*\nabla)_{f^*\xi}(f^*\eta)=f^*(\nabla_\xi\eta)$). In this sense, forming the torsion is natural, i.e. if $\nabla$ has torsion $T$, then the torsion of $f^*\nabla$ is $f^*T$ (pullback of tensor fields). In particular, if $\nabla$ is symmetric, then so is $f^*\nabla$. On the other hand, naturality of forming the torsion also implies that the projection $P$ in Anthony's answer is compatible with the action of local diffeomorphism, i.e. $P(f^*\nabla)=f^*(P(\nabla))$ for any connection $\nabla$.
The projection $P$ is not quite uniquely determined by this naturality propery (plus maybe some assumptions like locality), since you can play algebraic games with the torsion tensor, compare with the beginning of section 25 of the book MR1202431 by Kolar, Michor and Slovak, which is available online here. It becomes unique if you require that the symmetrized connection has the same geodesics.


TITLE: Is there a limit for how "good" a numerical method can be?
QUESTION [0 upvotes]: Multiplying two matrices $A \cdot B$ of size $n \times n$ in the trivial way requires $n^3$ computations. However, more efficient algorithms such as the Strassen algorithm have a lower complexity of approximately $\mathcal{O}(n^{2.8})$. There's even the Coppersmith-Winograd algorithm of theoretical complexity $\mathcal{O}(n^{2.37})$.
It seems unlikely that there would be for example a linear function that computes the matrix product, but is there a theoretical limit as to how effective any algorithm to multiply matrices must be, would for example $\mathcal{O}(n^2)$ be possible?
I'm also interested in extending the question beyound matrix multiplication. Are there such "theoretical lower bounds on complexity" for accomplishing a certain task with numerical methods (like multiplying two matrices) in general?

REPLY [0 votes]: There is a survey paper by Virginia Vassilevska Williams who showed the currently best bound of $\omega \leq 2.373$.
In the introduction she says:

The best known lower bounds for the problem are of the form $c \cdot n^2$ for explicit constants $c \leq 3$. Perhaps because of this, a majority of experts believe that there should exist an essentially quadratic time algorithm for matrix multiplication. Some prominent experts, however,
  such as Strassen, have also expressed the opinion that matrix multiplication requires strictly more
  than quadratic time. The answer has evaded computer science researchers for more than 40 years!

Thus it seems that we don't know more than $2 \leq \omega \leq 2.373$.


TITLE: How would you integrate $\sin(1/x)$?
QUESTION [1 upvotes]: I have found out that if i use u-substitution here's what i get :
$$\int ( \sin(u) ) \, dx$$
\begin{align*}
u &= 1/x = x^{-1}, \\
du / dx &= -x^{-2}, \\
dx / du &= -x^2,
\end{align*}
$$\int ( -\sin (u)x^2) \, du$$
$$= \text{?}$$
as you can see this cannot be equal to $\cos u$ because there is the term in $x$ as well. i know that we should generally divide by the derivative of the function but here the derivative is not a constant. my question is : what to do when the derivative is not a constant?

REPLY [0 votes]: Here are a few basic ways to integrate. I will use sine and cosine properties like their expansions. The first step uses integration by parts and a series expansion switching it and the integral by integrating term by term:
$$\mathrm{\int sin\left(\frac1x\right)dx=x\, sin\left(\frac1x\right)+\int\frac{cos\left(\frac1x\right)}{x}dx= x\,sin\left(\frac1x\right)+\int\frac1x\sum_{n=0}^\infty \frac{(-1)^n\left(\frac1x\right)^{2n}}{(2n)!}= x\,sin\left(\frac1x\right)+\sum_{n=0}^\infty \frac{(-1)^n}{(2n)!}\int x^{-2n-1}dx= C+x\,sin\left(\frac1x\right)+\sum_{n=0}^\infty \frac{(-1)^n\left(\frac1x\right)^{2n}}{2n(2n)!}}$$
This matches the main definition given for the Cosine Integral Ci(y). This means we have the special function closed form of:
$$\mathrm{\int sin\left(\frac1x\right)dx=C+ x\,sin\left(\frac1x\right)-Ci\left(\frac1x\right)}$$
There is also a way to use a sine series right away without integration by parts:
$$\mathrm{\int sin\left(\frac1x\right)dx=C+ x\,sin\left(\frac1x\right)-Ci\left(\frac1x\right)=\int\sum_{n=0}^\infty \frac{(-1)^n\left(\frac1x\right)^{2n+1}}{(2n+1)!}dx=C-\sum_{n=0}^\infty \frac{(-1)^n}{2n(2n+1)!x^{2n}}}$$
Here is a beautiful complex plot of the function. Please correct me and give me feedback!


TITLE: If $n$ is an odd integer, then $4\mid n^2 − 1$.
QUESTION [0 upvotes]: How do I prove using direct proof:  if $n$ is an odd integer, then $4$ divides $n^2 − 1$.

REPLY [2 votes]: With $n$ odd, we have
$n = 2k + 1. \; k \in \Bbb Z; \tag 1$
then
$n^2 = 4k^2 + 4k + 1, \tag 2$
whence
$n^2 - 1 = 4k^2 + 4k = 4(k^2 + k) \Longrightarrow 4 \mid n^2 -1. \tag 3$
One may also write
$n + 1 = 2k + 2, \tag 4$
$n - 1 = 2k, \tag 5$
$n^2 - 1 = (n + 1)(n - 1) = (2k + 2)(2k) = 2(k + 1)(2k)$
$= 4k(k + 1) \Longrightarrow 4 \mid n^2 - 1. \tag 6$
You can pick whichever proof you like the most!