11institutetext:  University of Washington, Seattle, WA, USA
22institutetext: Microsoft Research, Redmond, WA, USA

# Finding twin smooth integers by solving Pell equations

Jan Buzek This work was started as an undergraduate research project at the
Washington eXperimental Mathematics Lab (WXML) of the Department of
Mathematics at the University of Washington, Seattle.11 Junaid Hasan 11 Jason
Liu 11 Michael Naehrig 22 Anthony Vigil 11

###### Abstract

Any pair of consecutive $B$-smooth integers for a given smoothness bound $B$
corresponds to a solution $(x,y)$ of the equation $x^{2}-2\Delta y^{2}=1$ for
a certain square-free, $B$-smooth integer $\Delta$ and a $B$-smooth integer
$y$. This paper describes algorithms to find such twin $B$-smooth integers
that lie in a given interval by using the structure of solutions of the above
Pell equation.

The problem of finding such twin smooth integers is motivated by the quest for
suitable parameters to efficiently instantiate recent isogeny-based
cryptosystems. While the Pell equation structure of twin $B$-smooth integers
has previously been used to describe and compute the full set of such pairs
for very small values of $B$, increasing $B$ to allow for cryptographically
sized solutions makes this approach utterly infeasible.

We start by revisiting the Pell solution structure of the set of twin smooth
integers. Instead of using it to enumerate all twin smooth pairs, we focus on
identifying only those that lie in a given interval. This restriction allows
us to describe algorithms that navigate the vast set of Pell solutions in a
more targeted way. Experiments run with these algorithms have provided
examples of twin $B$-smooth pairs that are larger and have smaller smoothness
bound $B$ than previously reported pairs. Unfortunately, those examples do not
yet provide better parameters for cryptography, but we hope that our methods
can be generalized or used as subroutines in future work to achieve that goal.

##  1 Introduction

Twin smooth integers are pairs of consecutive smooth integers and have been
studied for various reasons at least since the 18th century. They have been
mentioned by Gauss in the context of computing logarithms of integers, and
they can be used to find particular solutions to certain diophantine equations
[7].

######  Definition 1 (Smooth and twin smooth integers)

Let $B>1$ and $m\in\mathbb{Z}$. The integer $m$ is called _$B$ -smooth_ if no
prime factor of $m$ is larger than $B$. The pair $(m,m+1)$ is called a pair of
_twin $B$-smooth integers_ if $m\cdot(m+1)$ is $B$-smooth.

Størmer [12] proved that for a fixed smoothness bound $B>1$, there is a finite
number of twin $B$-smooth integers and that all such pairs can be described
via solutions of Pell equations of the form $x^{2}-Dy^{2}=1$. Lehmer [7]
improved upon Størmer’s work in utilizing the structure of the set of
solutions to the Pell equation and strongly reduced the number of Pell
equations that must be solved in order to find all twin $B$-smooth pairs.
Lehmer computed the full set for $B=41$. Unfortunately, the number of Pell
equations that must be solved is still exponential in $B$, and therefore
finding all $B$-smooth pairs is computationally prohibitive for $B$ larger
than 200.

In this paper, we focus on the problem of finding some B-smooth pairs that lie
in a certain interval of possibly very large integers, rather than finding all
pairs; for example, given a smoothness bound $B$ as well as a bitlength $b$,
finding twin $B$-smooth integers that lie in the interval $[2^{b-1},2^{b})$,
i.e. finding such pairs $(m,m+1)$, for which $m$ can be written with $b$ bits.
A slight variant of the problem only assumes the bitlength $b$ is fixed and
then asks to find twin $B$-smooth integers with $b$ bits and the lowest
possible smoothness bound $B$.

####  1.0.1 Cryptographic applications.

The latter variant of the problem is motivated by applications in
cryptography. It has occurred recently in the context of selecting parameters
for instantiating public key cryptography that uses isogenies between elliptic
curves. For example, the digital signature scheme SQISign [5] requires a large
prime number $p$ such that $p-1$ and $p+1$ have large smooth factors. For
security reasons, the prime $p$ (and therefore the integer $m$) needs to have
a specific minimal size, usually given as a minimal number of bits, to ensure
that the schemes cannot be easily broken.

If one aims to find $p$ such that both $p-1$ and $p+1$ are fully smooth, the
problem is to find twin smooth integers $m=(p-1)/2$ and $m+1=(p+1)/2$ with the
additional condition that they sum to a prime number $p=m+(m+1)$. While this
was an interesting scenario for the (now broken) key exchange scheme B-SIDH
[4], it is not strictly necessary. Allowing a small non-smooth cofactor in
$p-1$ and $p+1$ often leads to a lower smoothness bound and a more efficient
implementation of the cryptosystem. Still, algorithms for finding (fully)
$B$-smooth twins, typically with much shorter bitsize, can be used to
construct cryptographically sized primes.

This motivates to look for twin smooths $(m,m+1)$ in a prescribed interval, or
more specifically with a given bitlength $b$, i.e.
$\lfloor\log_{2}(m)\rfloor+1=b$ Furthermore, the smoothness bound $B$ for
which a suitable twin $B$-smooth pair $(m,m+1)$ of the given size needs to be
found, has a direct impact on the efficiency of cryptographic primitives. The
smaller $B$ is, the more efficient the protocol becomes.

####  1.0.2 Algorithms for finding twin smooth integers.

As mentioned above, Lehmer’s paper [7] suggests to enumerate all twin
$B$-smooth pairs for fixed $B$ by solving the Pell equation $x^{2}-2\Delta
y^{2}=1$ for an exponential number of values for $\Delta$. This becomes
infeasible for even moderately large smoothness bound $B$. Luca and Najman
[10, 9] succeeded in computing all twin 100-smooth pairs with a modified
version of Lehmer’s approach. This computation was very costly because of the
large number of Pell equations to be solved and the solutions’ potentially
huge sizes. Costello [4, §5.3] used this approach to search for B-SIDH
parameters and reports to have solved all $2^{30}$ Pell equations for $B=113$
(the 30-th prime) in a large computational effort, discovering the largest one
to be a 76-bit $113$-smooth pair. An interesting algorithm that approximates
this problem and appears to be much more efficient is presented in [3].

####  1.0.3 The Conrey-Holmstrom-McLaughlin approach.

Conrey, Holmstrom, and McLaughlin propose an algorithm in [3] that quickly
generates a large set of moderately-sized twin smooth integers for a given
smoothness bound $B$. Experiments for relatively small $B$ demonstrate that
their method succeeds in generating almost all $B$-smooth pairs, only missing
a few that are often picked up when running the algorithm for a slightly
larger $B$. The above example of finding 100-smooth twins was easily run, but
missed 37 pairs of the Luca-Najman results. A significantly more costly search
for 200-smooth twins recovered all but one of those.

The algorithm works by observing that two twin smooth integer pairs, $(m,m+1)$
and $(M,M+1)$ with $m<M$, generate a new pair of twin smooths $(\mu,\mu+1)$ if
$\mu$ is an integer and $\frac{\mu}{\mu+1}=\frac{m}{m+1}\frac{M+1}{M}$. This
can be used to iteratively expand any set of twin smooths. Begin with the set
$S$ of integers $m<B$, though $S$ can be any set of integers $m$ such that
$(m,m+1)$ forms a smooth pair. At each iteration, take all pairs $m<M$ in $S$
and collect all new $\mu$ of the above form in a set $S^{\prime}$; set
$S=S\cup S^{\prime}$; if $S^{\prime}=\phi$, the algorithm stops. The largest
member of $S^{\prime}$ is bounded by the square of the largest member of $S$,
so the best case scenario is doubling the number of bits at each iteration.

While much more efficient in constructing most twin smooth pairs, this
algorithm seems to quickly become infeasible when aiming for large integers
that require a larger smoothness bound. The main reason is that the
cardinality of $S$ balloons rapidly, which makes successive iterations
infeasible before reaching a large enough bitlength. Nonetheless, this
algorithm is a powerful tool for quickly generating large sets for
experimentation, and an interesting avenue for future research, as
demonstrated recently in [1], where it was used as a subroutine for
successfully constructing parameters for SQISign. The paper also reports that
the authors ran the full algorithm for $B=547$. The largest $547$-smooth twin
smooth pair found in these experiments has 122 bits.

Moreover, the following algorithms have been used in the context of finding
cryptographic parameters. They do not attempt to compute the complete finite
set of twin $B$-smooth integers, but instead have the goal of generating just
enough “good” examples of a specific bit size. Ultimately the goal is to find
a pair with prime sum, but the probability of this happening is often large
enough within a large enough set of pairs. Specifically, the Prime Number
Theorem gives an estimate that random numbers in the cryptographic range
(around 256 bits) that we are looking for have a chance on the order of one in
a hundred of being prime. In comparison, the chances of random consecutive
integers in that range being B-smooth for B at most $2^{19}$ are tiny. Thus
assuming there is no interaction between smoothness of $(m,m+1)$ and primality
of $2m+1$, an algorithm that efficiently finds smooth pairs and then checks if
their sum is prime is an efficient way to find cryptographic parameters,
whereas an algorithm that checks whether $(\frac{p-1}{2},\frac{p+1}{2})$ is a
smooth pair for large primes $p$ is not. Therefore, the main task is to find
the twin smooth pairs.

####  1.0.4 Sieving.

A simple algorithm to search a given interval is to use a variant of the sieve
of Eratosthenes to list all $B$-smooth integers $m$ in the interval and check
for each whether $m+1$ is also smooth. This is inefficient because for large
$m$, the probability that $m+1$ is also smooth is very low and the chance to
encounter such pairs is so low that searchable intervals are too small and
likely do not even contain a solution or the search space is too large to be
covered. Nevertheless, plain sieving has been successfully used as a
subroutine in other methods. For example, some search algorithms increase the
probability of hitting twin smooth integers by using polynomial
parameterizations.

####  1.0.5 Sieving with $x^{n}-1$.

There are various polynomial parameterizations leading to twin smooth integers
that can be exploited by sieving for patterns of smaller smooth integers. A
simple polynomial that was found to work well in [2], [4], and [1] is
$f(x)=x^{n}-1$ for small values of $n$. For example, if $n=2$, then
$(m,m+1)=(f(x),f(x)+1)$ is a $B$-smooth pair when $(x-1)$, $x$, and $(x+1)$
are all smooth by the factorization $x^{2}-1=(x-1)(x+1)$. One can then use the
sieve to find three consecutive smooth integers $x-1,x,x+1$ and obtain the
twin pair $(x^{2}-1,x^{2})$ of twice as many bits. Because each of these
factors are about the square root of the resulting smooth pair, it is easier
to find the smooth triple of small integers. This phenomenon extends similarly
to larger values of $n$ such as $4$ or $6$. We give further explanation of the
reasons that certain polynomials work well in §2 (see in particular the
polynomials $p_{n}(x)$).

####  1.0.6 Sieving with Solutions to the Prouhet-Tarry-Escott Problem.

The Prouhet-Tarry-Escott (PTE) method is a generalization of the above idea
that finding smooth pairs parameterized by polynomials can be easier because
the polynomial factors of $(m,m+1)$ required to be smooth are much smaller
than the $m$ and $m+1$ themselves. Solutions to the PTE problem correspond to
polynomials $f(x)$ such that both $f(x)$ and $f(x)+1$ factor into linear terms
over $\mathbb{Z}$. Given such polynomials, we can search for pairs
$(m=f(x),m+1=f(x)+1)$ by checking for smoothness of all of the linear factors
of $f$ and $f+1$ for the same value of $x$. This is again done using a sieve.
Because smaller integers are much more likely to be smooth, it is often more
efficient to search for many small integers that are all smooth than for a
single large pair directly. The PTE method has found many of the best-known
pairs of cryptographic size; for details see [2]. However, this method is
restricted in that only a finite number of solutions to the PTE problem are
known, and all have degree at most $12$. Furthermore, it only finds pairs of
this special form, and it is unclear what type of pairs might be missed.

####  1.0.7 Extended Euclidean Algorithm.

Another improvement, discussed for example in [4], uses the extended Euclidean
algorithm. It chooses two coprime, $B$-smooth integers $\alpha$, $\beta$
roughly the square root of the size of the desired pair and runs the Euclidean
algorithm. This produces an equality of the form $\alpha s+\beta t=1$, and we
take $|\alpha s|,|\beta t|$ as our candidate smooth pair. The benefit is that
the portions $\alpha$ and $\beta$ of these numbers are already known to be
smooth, increasing the probability that $|\alpha\beta st|$ will be smooth
overall. For more details and a detailed analysis of the probabilities for
various sizes of $B$ and $m$, see [2] and [4].

####  1.0.8 Contributions.

In this paper, we describe how Pell equations can be used to search for large
twin-smooth pairs. Pell equations have been believed to be mainly a
theoretical tool useful for ensuring that all pairs have been found; we show
that they can also provide an efficient way of searching. We give an algorithm
generalizing and explaining the polynomial sieving described above and another
based on heuristically choosing which Pell equations to solve. We give
experimental evidence for our algorithms through the pairs that they found in
our search of the interval $[2^{240},2^{256}]$, including pairs that are the
largest known for their smoothness bound.

####  1.0.9 Outline.

In Section 2, we explore the structure of the solutions to the Pell equation,
deducing some results on explicit polynomial parameterizations from the work
of Lehmer. Then we specialize to the case of searching for pairs in a
specified large interval and describe how solving Pell equations can be
improved tremendously with this assumption. In Section 3, we apply this
assumption and our results on the structure of solutions to choosing which
Pell equations and which of their solutions to consider.

##  2 Solutions of the Pell equation

For a fixed, non-square positive integer $D$, the Pell equation is a
Diophantine equation of the form

$x^{2}-Dy^{2}=1$ (1)

and its solutions are pairs of integers $(x,y)$ that satisfy (1). The Pell
equation was studied in the seventh century by Brahmagupta, and since at least
the seventeenth century by mathematicians including Euler, Fermat, and
Lagrange. We give a brief overview of the results that interest us, generally
omitting proofs as these are readily available in resources such as [7], [6]
and [8].

Equation (1) always has the trivial solutions $(\pm 1,0)$, but our main
interest lies in describing all non-trivial solutions $(x,y)$ and without loss
of generality, we focus on those with $x,y>0$. The following well-known result
describes the structure of all such solutions (see, e.g. [6, Chapter 1], [8,
§1] and [7, §2]).

###### Theorem 2.1

Let $D>0$ be an integer that is not a perfect square. Then, the Pell equation
$x^{2}-Dy^{2}=1$ has an infinite number of integer solutions $(x,y)$ with
$x,y>0$. The solution with the smallest value of $x+y\sqrt{D}$ among them is
called the _fundamental solution_ and is denoted by $(x_{1},y_{1})$. All
nontrivial solutions as above can be written as $(x_{n},y_{n})$ for some
integer $n\geq 1$, where

$x_{n}+y_{n}\sqrt{D}=(x_{1}+y_{1}\sqrt{D})^{n}.$ (2)

The solution $(x_{n},y_{n})$ is called the _$n$ -th solution_.

###  2.1 Expressing solutions in terms of the fundamental solution

Lehmer [7, §2] describes all solutions of (1) in terms of the Lucas function
$U_{n}$. Let $\alpha=x_{1}+y_{1}\sqrt{D}$ and $\beta=x_{1}-y_{1}\sqrt{D}$.
Then, $\alpha\beta=1$ and for all $n\geq 1$, we get
$\alpha^{n}+\beta^{n}=2x_{n}$, $\alpha^{n}-\beta^{n}=2y_{n}\sqrt{D}$. The
Lucas function $U_{n}$ is defined as

$U_{n}=(\alpha^{n}-\beta^{n})/(\alpha-\beta)=y_{n}/y_{1}.$ (3)

Hence $U_{1}=1$ and the above identities yield recursive formulas to compute
$U_{n}$ for all $n$, since $x_{2n}=2x_{n}^{2}-1$, $U_{2n}=2x_{n}U_{n}$,
$x_{n+1}=x_{n}x_{1}+Dy_{n}y_{1}=x_{n}x_{1}+U_{n}(x_{1}^{2}-1)$ and
$U_{n+1}=x_{n}+x_{1}U_{n}$. Therefore, $U_{n}$ assumes integer values and
$y_{1}\mid y_{n}$ for all $n\geq 1$. Also, writing out (2) gives the explicit
formulas

$x_{n}=\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\binom{n}{2i}(Dy_{1}^{2})^{i}x_{1}^{n-2i},\quad
U_{n}=\sum_{i=0}^{\lceil\frac{n}{2}\rceil-1}\binom{n}{2i+1}(Dy_{1}^{2})^{i}x_{1}^{n-2i-1}.$
(4)

Note that we can replace $Dy_{1}^{2}=x_{1}^{2}-1$ via the Pell equation to
obtain univariate polynomials in $x_{1}$. We write these as polynomials
$p_{n}(x)\in\mathbb{Z}[x]$ such that $x_{n}=p_{n}(x_{1})$, i.e. we have
$p_{n}=\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\binom{n}{2i}(x^{2}-1)^{i}x^{n-2i}$
and $\deg(p_{n})=n$. Expanding $(x^{2}-1)^{i}$ and reordering the terms gives

$p_{n}(x)=\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\sum_{j=0}^{i}\binom{n}{2i}\binom{i}{j}(-1)^{j}x^{n-2j}=\sum_{j=0}^{\lfloor\frac{n}{2}\rfloor}\left((-1)^{j}\sum_{i=j}^{\lfloor\frac{n}{2}\rfloor}\binom{n}{2i}\binom{i}{j}\right)x^{n-2j}.$
(5)

The first few $p_{n}$, $n\in\\{1,2,\dots,12\\}$ are given as follows:

$\displaystyle p_{1}(x)$ $\displaystyle=x,$ $\displaystyle p_{2}(x)$
$\displaystyle=2x^{2}-1,$ $\displaystyle p_{3}(x)$
$\displaystyle=x(4x^{2}-3),$ $\displaystyle p_{4}(x)$
$\displaystyle=8x^{4}-8x^{2}+1,$ $\displaystyle p_{5}(x)$
$\displaystyle=x(16x^{4}-20x^{2}+5),$ $\displaystyle p_{6}(x)$
$\displaystyle=(2x^{2}-1)(16x^{4}-16x^{2}+1),$ $\displaystyle p_{7}(x)$
$\displaystyle=x(64x^{6}-112x^{4}+56x^{2}-7),$ $\displaystyle p_{8}(x)$
$\displaystyle=128x^{8}-256x^{6}+160x^{4}-32x^{2}+1,$ $\displaystyle p_{9}(x)$
$\displaystyle=x(4x^{2}-3)(64x^{6}-96x^{4}+36x^{2}-3),$ $\displaystyle
p_{10}(x)$ $\displaystyle=(2x^{2}-1)(256x^{8}-512x^{6}+304x^{4}-48x^{2}+1),$
$\displaystyle p_{11}(x)$
$\displaystyle=x(1024x^{10}-2816x^{8}+2816x^{6}-1232x^{4}+220x^{2}-11),$
$\displaystyle p_{12}(x)$
$\displaystyle=(8x^{4}-8x^{2}+1)(256x^{8}-512x^{6}+320x^{4}-64x^{2}+1).$

Similarly, Equation (4) gives an explicit formula for polynomials
$u_{n}(x)\in\mathbb{Z}[x]$ of degree $n-1$ such that $U_{n}=u_{n}(x_{1})$,
namely
$u_{n}(x)=\sum_{i=0}^{\lceil\frac{n}{2}\rceil-1}\binom{n}{2i+1}(x^{2}-1)^{i}x^{n-2i-1}$.
As above, by expanding $(x^{2}-1)^{i}$ and reordering, $u_{n}(x)$ can be
rewritten as

$u_{n}(x)=\sum_{j=0}^{\lceil\frac{n}{2}\rceil-1}\left((-1)^{j}\sum_{i=j}^{\lceil\frac{n}{2}\rceil-1}\binom{n}{2i+1}\binom{i}{j}\right)x^{n-2j-1}.$
(6)

Here are the first $12$ polynomials $u_{n}(x)$:

$\displaystyle u_{1}(x)$ $\displaystyle=1,$ $\displaystyle u_{2}(x)$
$\displaystyle=2x,$ $\displaystyle u_{3}(x)$ $\displaystyle=(2x-1)(2x+1),$
$\displaystyle u_{4}(x)$ $\displaystyle=4x(2x^{2}-1),$ $\displaystyle
u_{5}(x)$ $\displaystyle=(4x^{2}-2x-1)(4x^{2}+2x-1),$ $\displaystyle u_{6}(x)$
$\displaystyle=2x(2x-1)(2x+1)(4x^{2}-3),$ $\displaystyle u_{7}(x)$
$\displaystyle=(8x^{3}-4x^{2}-4x+1)(8x^{3}+4x^{2}-4x-1),$ $\displaystyle
u_{8}(x)$ $\displaystyle=8x(2x^{2}-1)(8x^{4}-8x^{2}+1),$ $\displaystyle
u_{9}(x)$ $\displaystyle=(2x-1)(2x+1)(8x^{3}-6x-1)(8x^{3}-6x+1),$
$\displaystyle u_{10}(x)$
$\displaystyle=2x(4x^{2}-2x-1)(4x^{2}+2x-1)(16x^{4}-20x^{2}+5),$
$\displaystyle u_{11}(x)$
$\displaystyle=(32x^{5}-16x^{4}-32x^{3}+12x^{2}+6x-1)(32x^{5}+16x^{4}-32x^{3}-12x^{2}+6x+1),$
$\displaystyle u_{12}(x)$
$\displaystyle=4x(2x-1)(2x+1)(2x^{2}-1)(4x^{2}-3)(16x^{4}-16x^{2}+1).$

Using the above polynomials, one can easily write down the $n$-th solution to
the Pell equation (1) in terms of the fundamental solution $(x_{1},y_{1})$ as

$(x_{n},y_{n})=(p_{n}(x_{1}),y_{1}u_{n}(x_{1})).$ (7)

Because the $p_{n}$ and $u_{n}$ are directly derived from (2) and (1), we
obtain the following lemma stating some divisibility properties.

###### Lemma 1

Let $p_{n}(x),u_{n}(x)\in\mathbb{Z}[x]$ be the above polynomials such that
$x_{n}=p_{n}(x_{1})$, $U_{n}=u_{n}(x_{1})$ for $0<n\in\mathbb{Z}$,
$\deg(p_{n})=n$ and $\deg(u_{n})=n-1$.

  1. (a)

If $n$ is odd, then $x\mid p_{n}(x)$.

  2. (b)

If $n=km$ for $0<k,m\in\mathbb{Z}$, then $p_{n}(x)=p_{m}(p_{k}(x))$ and
therefore, if $m$ is odd, $p_{k}(x)\mid p_{n}(x)$.

  3. (c)

Furthermore, $u_{n}(x)=u_{k}(x)u_{m}(p_{k}(x))$. In particular, if $k\mid n$,
then $u_{k}\mid u_{n}$.

  4. (d)

The leading coefficients of $p_{n}(x)$ and $u_{n}(x)$ are both equal to
$2^{n-1}$.

  5. (e)

If $n=2^{k}$, $k\geq 1$, then $p_{n}(x)$ is irreducible in $\mathbb{Z}[x]$.

###### Proof

Statement (a) follows from the explicit formula for $x_{n}$ in (4). Note that
if $n$ is odd, the exponent $n-2i$ of $x_{1}$ is always at least $1$.

If $n$ factors as $n=km$, then the exponentiation in Equation (2) can be split
up into two parts,
$x_{n}+y_{n}\sqrt{D}=(x_{1}+y_{1}\sqrt{D})^{n}=((x_{1}+y_{1}\sqrt{D})^{k})^{m}=(x_{k}+y_{k}\sqrt{D})^{m}$.
Writing out the last exponentiation by $m$ in terms of $x_{k}$ and $y_{k}$
explicitly gives
$x_{n}+y_{n}\sqrt{D}=p_{m}(x_{k})+u_{m}(x_{k})y_{k}\sqrt{D}=p_{m}(x_{k})+u_{m}(x_{k})y_{1}u_{k}(x_{1})\sqrt{D}$.
Now writing $x_{k}$ in terms of $x_{1}$ gives
$p_{m}(p_{k}(x_{1}))+u_{m}(p_{k}(x_{1}))y_{1}u_{k}(x_{1})\sqrt{D}$. The exact
same calculations can be carried out as polynomials in $x$ instead of
expressions in $x_{1}$, which gives
$x_{n}(x)+y_{n}(x)\sqrt{D}=p_{m}(p_{k}(x))+y_{1}u_{k}(x)u_{m}(p_{k}(x))\sqrt{D}$,
which proves (b) and (c).

Expression (5) shows that the leading coefficient of $p_{n}(x)$ is
$\sum_{i=0}^{\lfloor\frac{n}{2}\rfloor}\binom{n}{2i}=\binom{n}{0}+\binom{n}{2}+\dots+\binom{n}{2\lfloor\frac{n}{2}\rfloor}=2^{n-1}$.
From (6), one obtains the leading coefficient of $u_{n}(x)$ as
$\sum_{i=0}^{\lceil\frac{n}{2}\rceil-1}\binom{n}{2i+1}=\binom{n}{1}+\binom{n}{2}+\dots+\binom{n}{2\lceil\frac{n}{2}\rceil-1}=2^{n-1}$,
and thus (d) holds.

To prove (e), we first take a look at the coefficients of $p_{2^{k}}(x)$ when
we write it as
$p_{2^{k}}(x)=\sum_{j=0}^{2^{k-1}}a_{j}x^{2^{k}-2j}=a_{0}x^{2^{k}}+a_{1}x^{2^{k}-2}+\dots+a_{2^{k-1}-1}x^{2}+a_{2^{k-1}}$
with

$a_{j}=(-1)^{j}\sum_{i=j}^{2^{k-1}}\binom{2^{k}}{2i}\binom{i}{j}$

as in (5). Note that $\lfloor n/2\rfloor=2^{k-1}$. The constant coefficient is
$a_{2^{k-1}}=1$ for $k>1$ and $a_{2^{k-1}}=a_{1}=-1$ for $k=1$. As we have
seen in (d), the leading coefficient is
$a_{0}=\sum_{i=0}^{2^{k-1}}\binom{2^{k}}{2i}=2^{2^{k}-1}$.

We now show by induction on $k$ that, for all coefficients $a_{j}$ with $1\leq
j\leq 2^{k-1}$, $2^{2^{k}-2j}$ divides $a_{j}$. For the polynomial
$p_{2}(x)=2x^{2}-1$ ($k=1$), this is true. For $k=2$, we have
$p_{4}(x)=8x^{4}-8x^{2}+1$ and $a_{1}=-8$ is divisible by $4=2^{2^{2}-2}$. For
the induction step, note that $p_{2^{k+1}}(x)=2(p_{2^{k}}(x))^{2}-1$ by the
recursive formula for $x_{2n}$ discussed right after Equation (3). If we write
$(p_{2^{k}}(x))^{2}=\sum_{j=0}^{2^{k}}b_{j}x^{2^{k+1}-2j}=b_{0}x^{2^{k+1}}+b_{1}x^{2^{k+1}-2}+\dots+b_{2^{k}-1}x^{2}+b_{2^{k}}$,
then

$b_{j}=\sum_{\begin{subarray}{c}0\leq l,m\leq 2^{k-1}\\\
l+m=j\end{subarray}}a_{l}a_{m}.$

If $lm\neq 0$, then $2^{2^{k}-2m}\mid a_{m}$ and $2^{2^{k}-2l}\mid a_{l}$, so
$2^{2^{k+1}-2(l+m)}\mid a_{l}a_{m}$, i.e. if $l+m=j$, then $2^{2^{k+1}-2j}\mid
a_{l}a_{m}$. If one of $l$ or $m$ is $0$, the term involves
$a_{0}=2^{2^{k}-1}$. In all $b_{j}$ with $j>0$, $a_{0}$ only occurs in
products with $a_{j}$, which appear twice. Therefore, in these sums, the two
terms can be combined to $2a_{0}a_{j}$ and it follows again that
$2^{2^{k+1}-2j}\mid 2a_{0}a_{j}$. Overall, this shows that $2^{2^{k+1}-2j}\mid
b_{j}$. Since
$p_{2^{k+1}}=2b_{0}x^{2^{k+1}}+2b_{1}x^{2^{k+1}-2}+\dots+2b_{2^{k}-1}x^{2}+(2b_{2^{k}}-1)$,
the divisibility conditions still hold, which concludes our induction proof.

Now consider the polynomial $q(x)=2p_{2^{k}}(x/2)$. The leading coefficient
cancels exactly, making $q(x)$ monic. The constant term in $p_{2^{k}}(x)$ is
equal to $1$, which means that the constant term in $q(x)$ is equal to $2$.
The property we just proved ensures that all other coefficients of
$p_{2^{k}}(x/2)$ are integers. Thus
$q(x)=x^{2^{k}}+c_{1}x^{2^{k}-2}+\dots+c_{2^{k-1}-1}x^{2}+2$ for even integers
$c_{1},\dots,c_{2^{k-1}-1}$. The Eisenstein irreducibility criterion with the
prime $p=2$ shows that $q(x)$ is irreducible in $\mathbb{Z}[x]$. Then,
$p_{2^{k}}(x/2)=q(x)/2$ is irreducible in $\mathbb{Q}[x]$, and so is
$p_{2^{k}}(x)$. Because $p_{2^{k}}(x)\in\mathbb{Z}[x]$ is primitive, it is
also irreducible in $\mathbb{Z}[x]$. ∎

This lemma shows that $p_{n}(x)$ is irreducible if and only if $n=2^{k}$ for
some $k\geq 1$.

###  2.2 Finding the fundamental solution via continued fractions

In light of the above, solving the Pell equation is reduced to finding its
fundamental solution $(x_{1},y_{1})$. This can be done via continued fraction
expansions as we explain next.

We focus on simple continued fractions (CF) for square roots, which are
expressions of the form

$\sqrt{D}=a_{0}+\cfrac{1}{a_{1}+\cfrac{1}{a_{2}+\cfrac{1}{a_{3}+\cfrac{1}{a_{4}+\dots}}}}$

for $D$ a squarefree integer and a sequence of integers $a_{0},a_{1},\dots$,
called the quotients, where $a_{i}>0$ for all $i$.

Every rational number can be expressed in two different ways via CF, and every
irrational number has a unique CF representation that is an infinite sequence
([13], 2.2); $\sqrt{D}$ falls in this second case. We form the sequence of
convergents $\frac{A_{i}}{B_{i}}$ to $\sqrt{D}$ by truncating the CF expansion
after a finite number of terms. Specifically, the $i$th convergent is given by

$a_{0}+\cfrac{1}{a_{1}+\cfrac{1}{\ddots+\cfrac{1}{a_{i}}}}.$

These convergents can also be computed by the recurrence relations

$\displaystyle A_{i}$ $\displaystyle=a_{i}A_{i-1}+A_{i-2},$ $\displaystyle
B_{i}$ $\displaystyle=a_{i}B_{i-1}+B_{i-2},$

beginning with the values $A_{0}=a_{0}$, $A_{1}=a_{1}a_{0}+1$, $B_{0}=1$, and
$B_{1}=a_{1}$ ([13], 2.4). An important observation is that the numerators and
denominators $A_{i}$ and $B_{i}$ of these convergents are monotonically
increasing; further because the $a_{i}$ are positive integers, these sequences
are bounded below by the Fibonacci sequence.

The quotients $a_{i}$ can be found almost directly from the definition by
recursively reciprocating and rounding down, however, the following recurrence
relations ([11, Theorem 10.19]) give a faster implementation.

$\displaystyle\alpha_{k}$
$\displaystyle=\left\lfloor\frac{P_{k}+\sqrt{D}}{Q_{k}}\right\rfloor,$
$\displaystyle P_{k+1}$ $\displaystyle=a_{k}Q_{k}-P_{k},$ $\displaystyle
Q_{k+1}$ $\displaystyle=\frac{d-P_{k+1}^{2}}{Q_{k}},$

where $Q_{0}=1,P_{0}=0$.

The connection between continued fractions and the Pell equation is given by
the following theorem ([6, Theorem 3.3]).

###### Theorem 2.2

Let $D>0$ be an integer that is not a perfect square. Then if $(x,y)$ is a
solution to the Pell equation $x^{2}-Dy^{2}=1$ then $\frac{x}{y}$ is a
convergent in the continued fraction expansion for $\sqrt{D}$.

This theorem suggests that to find the fundamental solution, we generate the
quotients $a_{i}$ with our recurrence relations, using these to generate
$A_{i}$ and $B_{i}$ for each $i$, and check whether $A_{i}$ and $B_{i}$
satisfy the Pell equation as $x$ and $y$. Because a fundamental solution
exists, this algorithm always terminates. We will only be interested in
solutions to the Pell equation that lie in a given range; this additional
restriction and our bounds on the monotonically increasing sequence $A_{i}$
will ensure that this algorithm is very fast. Specifically, if we give an
upper bound $\tilde{M}$ on the size of $x$ in the solution to the Pell
equation that we want to find, then we can cut off our algorithm after the
first convergent $A_{i}$ is generated that is larger than $\tilde{M}$. Because
the $A_{i}$ are bounded below by the Fibonnaci sequence, this means that our
algorithm will find any solution with $x\leq\tilde{M}$ in time proportional to
log($\tilde{M}$). In our search for cryptographically large twin smooth pairs
by solving Pell equations, $\tilde{M}$ is $2^{258}$, and a single core running
our C-code implementing this algorithm can find any solutions with
$x\leq\tilde{M}$ for tens of thousands of Pell equations a second.

##  3 Twin smooth integers as solutions of Pell equations

Lehmer states the correspondence between $B$-smooth twins and certain
solutions to Pell equations in Theorem 1 of [7]. The theorem is more general
because it allows a specific finite set of primes as the prime factors of the
twin smooth integers in question. But we are mainly interested in capturing
$B$-smoothness, so we state it here specialized to consider all primes not
exceeding $B$.

######  Theorem 3.1 (Lehmer [7], Thm. 1)

Let $B>2$ and $0<m\in\mathbb{Z}$. Then $(m,m+1)$ is a pair of twin $B$-smooth
integers if and only if there exists a $B$-smooth integer $y$ and a square-
free $B$-smooth integer $\Delta\neq 2$ such that $(x_{n},y_{n})$ with
$x_{n}=2m+1$, $y_{n}=y$ is the $n$-th solution to the Pell equation

$x^{2}-2\Delta y^{2}=1$ (8)

with $1\leq n\leq\max\\{3,(q+1)/2\\}$, where $q$ is the largest prime not
exceeding $B$.

###### Remark 1

If $\Delta$ in Equation (8) is odd then the equation is a Pell equation in the
sense of §2 with a square-free, positive integer $D=2\Delta$. But Theorem 3.1
allows $\Delta$ to be even, i.e. in which case $\Delta=2D$ for an odd, square-
free, positive $D$. In this case, a solution $(x_{n}^{\prime},y_{n}^{\prime})$
to Equation (8) corresponds to the solution
$(x_{n}=x_{n}^{\prime},y_{n}=2y_{n}^{\prime})$ of the Pell equation
$x^{2}-Dy^{2}=1$, i.e. here $D=\Delta/2$.

Furthermore, Equation (8) clearly implies that $x_{n}$ must be odd for any
solution $(x_{n},y_{n})$. So, if $x_{n}=2m+1$, then $\Delta
y_{n}^{2}=(x_{n}^{2}-1)/2=2m(m+1)$ is divisible by $4$. Since $\Delta$ is
square-free, $y_{n}$ must be even.

If $2=q_{1},q_{2},\dots,q_{t}\in\mathbb{Z}$ are all the prime numbers that are
not larger than $B$, listed consecutively in ascending order, then Lehmer’s
theorem shows that in order to find all twin $B$-smooth pairs, it suffices to
solve $2^{t}-1$ Pell equations of the form (8), because that is the number of
possible values $\Delta$ can take in Theorem 3.1. Indeed, following Lehmer’s
notation, let $Q$ denote the set of positive $B$-smooth integers, i.e.

$Q=\\{q_{1}^{\alpha_{1}}\cdot q_{2}^{\alpha_{2}}\cdot\ldots\cdot
q_{t}^{\alpha_{t}}\mid 0\leq\alpha_{i}\in\mathbb{Z},\ 1\leq i\leq t\\},$

and let $Q^{\prime}$ be the subset of square-free elements of $Q$ except $2$,
i.e.

$Q^{\prime}=\\{q_{1}^{\epsilon_{1}}\cdot q_{2}^{\epsilon_{2}}\cdot\ldots\cdot
q_{t}^{\epsilon_{t}}\mid\epsilon_{i}\in\\{0,1\\},\ 1\leq i\leq
t\\}\setminus\\{2\\}.$

Then $Q^{\prime}$ is the set of possible values of $\Delta$ that can occur in
Theorem 3.1 and has $2^{t}-1$ elements.

Theorem 3.1 gives an explicit correspondence between twin smooth pairs
$(m,m+1)$ and Pell equation coefficients $\Delta\in Q^{\prime}$ together with
solutions $(x,y)\in\mathbb{Z}_{>0}\times Q$. Define the set

$P=\\{(\Delta,x,y)\in Q^{\prime}\times\mathbb{Z}_{>0}\times Q\mid
x^{2}-2\Delta y^{2}=1\\}$

of triples $(\Delta,x,y)$ of such coefficient-solution triples and the set

$T=\\{m\in Q\mid m(m+1)\in Q\\}\subset Q$

of elements $m$ that describe all $B$-smooth twins. The correspondence is
given in the following lemma, which is simply a reformulation of the proof of
Theorem 3.1 as given in [7].

###### Lemma 2

Let the sets $P$ and $T$ be defined as above. The map $\delta:P\rightarrow T,\
(\Delta,x,y)\mapsto m=(x-1)/2$ is a bijection.

###### Proof

Let $(\Delta,x,y)\in P$. As noted in Remark 1, $x$ must be odd and $y$ must be
even. Hence both $m=(x-1)/2$ and $y/2$ are integers, as is $m+1=(x+1)/2$. This
means that $m\cdot(m+1)=(x^{2}-1)/4=2\Delta(y/2)^{2}\in Q$ because
$\Delta,y\in Q$, showing that $\delta(\Delta,x,y)$ lies in $T$.

Conversely, given $m\in T$ representing the twin $B$-smooth pair $(m,m+1)$, we
have $m(m+1)\in Q$ and we can uniquely write $m(m+1)=2q_{1}^{\alpha_{1}}\cdot
q_{2}^{\alpha_{2}}\cdot\ldots\cdot q_{t}^{\alpha_{t}}$, where
$\alpha_{i}=\epsilon_{i}+2\beta_{i}$ for $\epsilon_{i}\in\\{0,1\\}$ and
integers $\alpha_{i},\beta_{i}\geq 0$. Setting $x=2m+1$,
$y=2q_{1}^{\beta_{1}}\cdot q_{2}^{\beta_{2}}\cdot\ldots\cdot
q_{t}^{\beta_{t}}$ and $\Delta=q_{1}^{\epsilon_{1}}\cdot
q_{2}^{\epsilon_{2}}\cdot\ldots\cdot q_{t}^{\epsilon_{t}}$, we see that
$x^{2}-1=4m(m+1)=2\Delta y^{2}$. Therefore, given $m$, we obtain $x=2m+1$ and
unique $\Delta$ and $y$ by the square-free and square parts of $4m(m+1)$ via
$x^{2}-1=2\Delta y^{2}$, so we can define $\delta^{-1}$ by
$\delta^{-1}(m)=(\Delta,x,y)\in P$. We have $\delta\circ\delta^{-1}(m)=m$ for
all $m\in T$ and $\delta^{-1}\circ\delta(\Delta,x,y)=(\Delta,x,y)$ for all
$(\Delta,x,y)\in P$, which implies the lemma.∎

Any solution $(x,y)$ to Equation (8) is the $n$-th solution for some positive
integer $n$ as defined in Theorem 2.1. Lehmer [7, Thm. 1] showed that $n\leq
M=\max\\{3,(q_{t}+1)/2\\}$. Therefore, the set $P$ is finite. An algorithm to
compute all twin $B$-smooth pairs is thus to enumerate all coefficients
$\Delta\in Q^{\prime}$, compute the fundamental solution of the corresponding
Pell equation $x^{2}-2\Delta y^{2}=1$ and then to compute all $n$-th solutions
from that for $n\leq M$. To check whether any of those solutions yields a pair
of twin $B$-smooth integers, it is sufficient to check whether the
$y$-component of the solution is $B$-smooth. Clearly, this requires solving
$2^{t}-1$ Pell equations and quickly becomes infeasible for even moderate
values of $B$.

###  3.1 Twin smooth integers from $n$-th solutions for $n>1$

Next we consider for which $n$ the $n$-th solution to a Pell equation leads to
a twin smooth pair via the above correspondence. We start by showing that if
the $n$-th solution leads to a twin smooth pair then so does the fundamental
solution. In other words, an $n$-th solution for $n>1$ can only correspond to
a twin smooth pair if the corresponding fundamental solution corresponds to a
twin smooth pair.

###### Lemma 3

Let $\Delta\in Q^{\prime}$, $(x_{1},y_{1})$ be the fundamental solution and
$(x_{n},y_{n})$ the $n$-th solution to $x^{2}-2\Delta y^{2}=1$. Let
$m_{1}=(x_{1}-1)/2$ and $m_{n}=(x_{n}-1)/2$. If $m_{n}\in T$, then so is
$m_{1}$.

###### Proof

Since $m_{n}\in T$, we have $\delta^{-1}(m_{n})=(\Delta,x_{n},y_{n})\in P$ and
therefore, $y_{n}\in Q$. We have seen in §2 that $y_{1}\mid y_{n}$, so
$y_{1}\in Q$, which means that $(\Delta,x_{1},y_{1})\in P$. Then
$m_{1}=\delta(\Delta,x_{1},y_{1})\in T$. ∎

This lemma shows that when solving Pell equations to find twin smooth pairs,
an equation can be discarded if the fundamental solution does not yield a twin
smooth pair and there is no need to check any of the successive solutions
$(x_{n},y_{n})$, $n>1$ for that particular value of $\Delta$. From now on, let
$m_{n}=(x_{n}-1)/2$ for $n\geq 1$, where $(x_{n},y_{n})$ is the $n$-th
solution to the Pell equation (8).

Whenever the fundamental solution $(x_{1},y_{1})$ yields a twin smooth pair
$(m_{1},m_{1}+1)$, successive solutions might yield one as well. To decide
whether for a fixed $(\Delta,x_{1},y_{1})\in P$ we also have
$(\Delta,x_{n},y_{n})\in P$, it is sufficient to check whether $U_{n}\in Q$
because then $y_{n}=U_{n}y_{1}\in Q$. For that purpose, we can simply write
$U_{n}=v_{n}(m_{1})$ as a polynomial in $m_{1}$, where $v_{n}(x)=u_{n}(2x+1)$.
Then $U_{n}=u_{n}(x_{1})=u_{n}(2m_{1}+1)=v_{n}(m_{1})$. These expressions for
$n\leq 12$ are shown here:

$\displaystyle v_{1}(m)$ $\displaystyle=1,$ $\displaystyle v_{2}(m)$
$\displaystyle=2(2m+1),$ $\displaystyle v_{3}(m)$
$\displaystyle=(4m+1)(4m+3),$ $\displaystyle v_{4}(m)$
$\displaystyle=4(2m+1)(8m^{2}+8m+1),$ $\displaystyle v_{5}(m)$
$\displaystyle=(16m^{2}+12m+1)(16m^{2}+20m+5),$ $\displaystyle v_{6}(m)$
$\displaystyle=2(2m+1)(4m+1)(4m+3)(16m^{2}+16m+1),$ $\displaystyle v_{7}(m)$
$\displaystyle=(64m^{3}+80m^{2}+24m+1)(64m^{3}+112m^{2}+56m+7),$
$\displaystyle v_{8}(m)$
$\displaystyle=8(2m+1)(8m^{2}+8m+1)(128m^{4}+256m^{3}+160m^{2}+32m+1),$
$\displaystyle v_{9}(m)$
$\displaystyle=(4m+1)(4m+3)(64m^{3}+96m^{2}+36m+1)(64m^{3}+96m^{2}+36m+3),$
$\displaystyle v_{10}(m)$
$\displaystyle=2(2m+1)(16m^{2}+12m+1)(16m^{2}+20m+5)(256m^{4}+512m^{3}+304m^{2}+48m+1),$
$\displaystyle v_{11}(m)$
$\displaystyle=(1024m^{5}+2304m^{4}+1792m^{3}+560m^{2}+60m+1)$
$\displaystyle\quad\cdot(1024m^{5}+2816m^{4}+2816m^{3}+1232m^{2}+220m+11),$
$\displaystyle v_{12}(m)$
$\displaystyle=4(2m+1)(4m+1)(4m+3)(8m^{2}+8m+1)(16m^{2}+16m+1)$
$\displaystyle\quad\cdot(256m^{4}+512m^{3}+320m^{2}+64m+1).$

For the $n$-th solution to give a twin $B$-smooth pair, $U_{n}=v_{n}(m_{1})$
must be smooth, where $(m_{1},m_{1}+1)$ is the twin smooth pair given by the
fundamental solution. This means that each of the irreducible terms in the
polynomial representation of $v_{n}$ needs to evaluate to a $B$-smooth value
at $m_{1}$. In the same way as above, we may write polynomials for the
potential twin-smooth pair $(m_{n},m_{n}+1)$ corresponding to the $n$-th
solution, namely, we have that

$m_{n}=(p_{n}(2m_{1}+1)-1)/2.$ (9)

In the proof of Lemma 2 we have seen that
$m_{1}(m_{1}+1)=(x_{1}^{2}-1)/4=2\Delta(y_{1}/2)^{2}$ and
$m_{n}(m_{n}+1)=(x_{n}^{2}-1)/2=2\Delta(y_{n}/2)^{2}$. Since
$y_{n}=y_{1}U_{n}$, this means that

$m_{n}(m_{n}+1)=m_{1}(m_{1}+1)U_{n}^{2}=m_{1}(m_{1}+1)v_{n}(m_{1})^{2}.$ (10)

The following lemma gives a more specific description about how $m_{1}$ and
$m_{1}+1$ divide $m_{n}$ and $m_{n}+1$.

###### Lemma 4

Let $\Delta\in Q^{\prime}$, $(x_{1},y_{1})$ be the fundamental solution and
$(x_{n},y_{n})$ the $n$-th solution to $x^{2}-2\Delta y^{2}=1$. Let
$m_{1}=(x_{1}-1)/2$ and $m_{n}=(x_{n}-1)/2$. Then $m_{1}\mid m_{n}$. If $n$ is
even, then $(m_{1}+1)\mid m_{n}$, if it is odd, then $(m_{1}+1)\mid(m_{n}+1)$.

###### Proof

Note that $x_{1}^{2}-1=4m_{1}(m_{1}+1)$ and therefore, the explicit formula
for $p_{n}$ in §2 shows that $p_{n}(2m_{1}+1)\equiv 1\mod(2m_{1})$. It follows
that $2m_{1}\mid p_{n}(2m_{1}+1)-1$ and thus
$m_{1}\mid(p_{n}(2m_{1}+1)-1)/2=(x_{n}-1)/2=m_{n}$.

We also see that $p_{n}(2m_{1}+1)\equiv(-1)^{n}\mod(2(m_{1}+1))$ such that
$2(m_{1}+1)\mid p_{n}(2m_{1}+1)-1$ if $n$ is even and $2(m_{1}+1)\mid
p_{n}(2m_{1}+1)+1$ if $n$ is odd. ∎

Below are the expressions for $(m_{n},m_{n}+1)$ with $1\leq n\leq 12$.

$\displaystyle(m_{2},m_{2}+1)$ $\displaystyle=(4m(m+1),\ (2m+1)^{2}),$
$\displaystyle(m_{3},m_{3}+1)$ $\displaystyle=(m(4m+3)^{2},\
(m+1)(4m+1)^{2}),$ $\displaystyle(m_{4},m_{4}+1)$
$\displaystyle=(16m(m+1)(2m+1)^{2},\ (8m^{2}+8m+1)^{2}),$
$\displaystyle(m_{5},m_{5}+1)$ $\displaystyle=(m(16m^{2}+20m+5)^{2},\
(m+1)(16m^{2}+12m+1)^{2}),$ $\displaystyle(m_{6},m_{6}+1)$
$\displaystyle=(4m(m+1)(4m+1)^{2}(4m+3)^{2},\ (2m+1)^{2}(16m^{2}+16m+1)^{2}),$
$\displaystyle(m_{7},m_{7}+1)$ $\displaystyle=(m(64m^{3}+112m^{2}+56m+7)^{2},\
(m+1)(64m^{3}+80m^{2}+24m+1)^{2}),$ $\displaystyle(m_{8},m_{8}+1)$
$\displaystyle=(64m(m+1)(2m+1)^{2}(8m^{2}+8m+1)^{2},$
$\displaystyle\qquad(128m^{4}+256m^{3}+160m^{2}+32m+1)^{2}),$
$\displaystyle(m_{9},m_{9}+1)$
$\displaystyle=(m(4m+3)^{2}(64m^{3}+96m^{2}+36m+3)^{2},$
$\displaystyle\qquad(m+1)(4m+1)^{2}(64m^{3}+96m^{2}+36m+1)^{2}),$
$\displaystyle(m_{10},m_{10}+1)$
$\displaystyle=(4m(m+1)(16m^{2}+12m+1)^{2}(16m^{2}+20m+5)^{2},$
$\displaystyle\qquad(2m+1)^{2}(256m^{4}+512m^{3}+304m^{2}+48m+1)^{2}),$
$\displaystyle(m_{11},m_{11}+1)$
$\displaystyle=(m(1024m^{5}+2816m^{4}+2816m^{3}+1232m^{2}+220m+11)^{2},$
$\displaystyle\qquad(m+1)(1024m^{5}+2304m^{4}+1792m^{3}+560m^{2}+60m+1)^{2}),$
$\displaystyle(m_{12},m_{12}+1)$
$\displaystyle=(16m(m+1)(2m+1)^{2}(4m+1)^{2}(4m+3)^{2}(16m^{2}+16m+1)^{2},$
$\displaystyle\qquad(8m^{2}+8m+1)^{2}(256m^{4}+512m^{3}+320m^{2}+64m+1)^{2}).$

###### Remark 2

Lemma 1(b) shows that as soon as $n$ has an odd factor, the polynomial
$p_{n}(x)$ used to represent $x_{n}=p_{n}(x_{1})$ is reducible. This also
means that $2m_{n}+1=p_{n}(2m_{1}+1)$ as a polynomial in $m_{1}$ is reducible
and cannot assume a prime value for any $m_{1}$. Only $p_{2^{k}}(x)$ is
irreducible for any $k\geq 1$ according to Lemma 1(e). Therefore, when we are
interested in finding twin smooth pairs $(m_{n},m_{n}+1)$, where
$p=m_{n}+(m_{n}+1)$ is prime, only those solutions need to be checked where
$n=2^{k}$ is a power of $2$.

Since pairs $(m_{n},m_{n}+1)$ are polynomial expressions in the corresponding
$m_{1}$, one can easily obtain bounds on the size of $m_{n}$ in terms of the
size of $m_{1}$. For example, equations (9) and (10) allow to deduce a simple
lower bound for $m_{n}$ in terms of the leading coefficient of $p_{n}$ or
$v_{n}$. We obtain the following lemma.

###### Lemma 5

Let $\Delta\in Q^{\prime}$, $(x_{1},y_{1})$ be the fundamental solution and
$(x_{n},y_{n})$ the $n$-th solution to $x^{2}-2\Delta y^{2}=1$, and let
$m_{1}=(x_{1}-1)/2$, $m_{n}=(x_{n}-1)/2$, $U_{n}=y_{n}/y_{1}$ and
$v_{n}(x)\in\mathbb{Z}[x]$ the above polynomial such that
$U_{n}=v_{n}(m_{1})$. Then

$m_{1}v_{n}(m_{1})\leq m_{n}<(m_{1}+1)v_{n}(m_{1}).$

It follows that $m_{n}>4^{n-1}m_{1}^{n}$.

###### Proof

We use Equation (10) to see that
$m_{n}^{2}<m_{n}(m_{n}+1)=m_{1}(m_{1}+1)v_{n}(m_{1})^{2}<(m_{1}+1)^{2}v_{n}(m_{1})^{2}$,
which yields $m_{n}<(m_{1}+1)v_{n}(m_{1})$. Similarly,
$(m_{n}+1)^{2}>m_{n}(m_{n}+1)=m_{1}(m_{1}+1)v_{n}(m_{1})^{2}>m_{1}^{2}v_{n}(m_{1})^{2}$
and therefore $m_{n}\geq m_{1}v_{n}(m_{1})$. Because $v_{n}(x)=u_{n}(2x+1)$,
$\deg(u_{n})=n-1$ and the leading coefficient of $u_{n}$ is $2^{n-1}$
according to Lemma 1, the leading coefficient of $v_{n}$ is
$(2^{n-1})^{2}=4^{n-1}$. The explicit formulas for $u_{n}$ in Section 2 show
that all coefficients of $v_{n}(x)$ are positive and therefore,
$v_{n}(m_{1})>4^{n-1}m_{1}^{n-1}$, which concludes the proof.∎

The last inequality in the lemma provides an upper bound on $m_{1}$ if we are
given an upper bound on $m_{n}$. Assume, we are interested in twin smooth
integers given by $m_{n}$ corresponding to an $n$-th solution to the Pell
equation with $n>1$ with an upper bound $m_{n}<R$. Such pairs must then come
from a fundamental solution pair given by $m_{1}$ such that
$4^{n-1}m_{1}^{n}<R$, or, equivalently, $m_{1}<\frac{1}{4}(4R)^{1/n}$. For
example, to find twin smooth integers where $m_{n}$ has at most $b$ bits, we
have $R=2^{b}$, i.e. $m_{1}<R_{b}:=2^{(b-2n+2)/n}$. Table 1 shows concrete
values for the maximal possible bitsize of $m_{1}$ when $b=256$ and
$n\in\\{1,2,\dots,12\\}$.

$n$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $8$ | $9$ | $10$ | $11$ | $12$
---|---|---|---|---|---|---|---|---|---|---|---|---
$\lceil 258/n\rceil-2$ | $256$ | $127$ | $84$ | $63$ | $50$ | $41$ | $35$ | $31$ | $27$ | $24$ | $22$ | $20$
Table 1: Maximal bitsizes $\lceil\log(R_{256})\rceil=\lceil 258/n\rceil-2$ of
$m_{1}$ such that the corresponding $m_{n}$ has at most $256$ bits for
$n=1,2,\dots,12$.

Further we have the following absolute upper bound on $n$ for which we can
have twin smooth pairs corresponding to the $n$-th solution to a Pell equation
(see [7], Theorem 1).

###### Theorem 3.2

If $(m,m+1)$ is a pair of twin $B$-smooth integers corresponding to the $n$-th
solution to a Pell equation, then $n\leq\frac{B+1}{2}$.

This upper bound is in many ways much weaker than the results of Lemma 5
because Lemma 5 shows that for large $n$, the $n$-th solutions to Pell
equations under our bound all correspond to fundamental solutions that are
exponentially smaller. Together, this suggests an algorithm for finding all
smooth pairs in an interval corresponding to $n$-th solutions to Pell
equations for large $n$. The algorithm works as shown in Algorithm 1.

1:Fix a lower bound $s$ on $n$.

2:for $n\in\\{s,s+1,\dots,\frac{B+1}{2}$} do

3: Compute $T=\lceil M/n\rceil-2$ as in Table 1.

4: for $w\leq T$ do

5: if $(w,w+1)$ is a twin $B$-smooth pair then

6: if $v_{n}(w)$ is $B$-smooth then

7: Record $(m_{n},m_{n}+1)$ as a twin $B$-smooth pair.

8: end if

9: end if

10: end for

11:end for

Algorithm 1 Algorithm to search for twin $B$-smooth pairs of size at most $M$
bits.

The key benefit of this algorithm is that it finds smooth pairs of size
roughly $2^{M}$ while only searching through integers up to $2^{T}$, an
exponentially smaller range. In our running example of searching for pairs in
the range $[2^{240},2^{256}]$, we ran the above algorithm with a lower bound
of $s=6$. Then the most costly iteration above with $n=6$ only required
checking approximately $2^{41}$ possibilities for $w$. Notice further that as
$n$ gets larger, the possibilities for $w$ grow extremely restricted, and the
algorithm terminates quickly. The overall runtime is dominated by the
iteration with $n=s$.

For example, our run with $s=6$ and a maximum of $M=2^{275}$ found the
following smooth pairs. A $16759$-smooth, $245$-bit pair $(m,m+1)$ with

$\displaystyle m$
$\displaystyle=44746808406030847930450201970587971020922341276429366152081686798603000000$
$\displaystyle=2^{6}\cdot 3^{4}\cdot 5^{6}\cdot 7^{5}\cdot 11^{2}\cdot
13^{2}\cdot 41\cdot 43^{2}\cdot 53\cdot 97\cdot 241^{2}\cdot 337\cdot 509\cdot
673\cdot 4703^{2}$ $\displaystyle\quad\cdot 5981\cdot 9413^{2}\cdot
13669^{2}\cdot 16759^{2},$ $\displaystyle m+1$ $\displaystyle=31^{2}\cdot
157^{4}\cdot 181^{2}\cdot 251^{2}\cdot 349^{2}\cdot 359^{2}\cdot 457^{2}\cdot
1427^{2}\cdot 2617^{2}\cdot 9109^{2}\cdot 9649^{2}\cdot 10253^{2}$

and a $24551$-smooth, $260$-bit pair with

$\displaystyle m$
$\displaystyle=1248045507865502270977250845951694434798578493856490782548653674169732908101560$
$\displaystyle=2^{3}\cdot 3^{5}\cdot 5\cdot 7^{2}\cdot 17\cdot 19\cdot 29\cdot
31\cdot 43^{2}\cdot 53^{2}\cdot 149^{2}\cdot 211^{2}\cdot 227^{2}\cdot
233\cdot 827\cdot 919^{2}\cdot 2659\cdot 4723^{2}$ $\displaystyle\quad\cdot
6907^{2}\cdot 10831\cdot 16691^{2}\cdot 24551,$ $\displaystyle m+1$
$\displaystyle=11^{2}\cdot 13^{2}\cdot 23^{2}\cdot 71^{2}\cdot 107^{2}\cdot
263^{2}\cdot 587^{2}\cdot 1021^{2}\cdot 4297^{2}\cdot 6491^{2}\cdot
7309^{2}\cdot 8089^{2}\cdot 9049^{2}\cdot 19009^{2}.$

This second pair is larger than any previously known for its smoothness bound.

###  3.2 Lower order solutions

The above algorithm finds all twin-smooth pairs corresponding to $n$-th
solutions for $n$ large. In order to complete our approach, we now consider
algorithms for finding pairs corresponding to Pell equation solutions with $n$
small. In this case, the bounds from Lemma 5 are too weak, and the algorithm
from the previous section is intractable. Therefore, we search for such smooth
pairs directly by solving Pell equations via continued fractions to obtain the
fundamental solution, and then checking the successive solutions for small $n$
with our polynomials $v_{n}$.

The obvious difficulty with the approach is that there are exponentially many
Pell equations which could have corresponding smooth pairs, and we cannot
solve all of them. While we believe (see discussion) that this difficulty may
be intrinsic, and that finding all pairs may not be tractable, we can still
cleverly choose which Pell equations to solve in order to find some large
twin-smooth pairs.

Recall that each of the $2^{t}-1$ Pell equations $x^{2}-Dy^{2}=1$ with
potential corresponding smooth pairs has a unique coefficient $D=2\Delta$
where $\Delta$ is a smooth, square-free integer and $\Delta\neq 2$. Given any
particular $\Delta$, solving the corresponding Pell equation can be done
extremely quickly (see 2.2). Choosing which Pell equations to solve means
choosing the values of $\Delta$. Generally, the difficulty is in finding
$\Delta$s for which the size of the solution to the Pell equation is in the
desired range. Fundamental solutions to Pell equations can be extremely large
(exponential in $\Delta$), and such large solutions $(x,y)$ are very unlikely
to have smooth $y$, and are thus not helpful to our search. It is therefore
beneficial to choose coefficients $\Delta$ for which we expect the resulting
pairs to generally be small.

###  3.3 The smallest coefficient heuristic

A very simple heuristic for choosing $\Delta$ is to choose the smallest values
possible. This is beneficial for a simple combinatorial reason. Note that if
$(x,y)$ is a solution to the Pell equation

$x^{2}-1=2\Delta y^{2}$

then $y=\sqrt{\frac{x^{2}-1}{2\Delta}}$ is an integer, and as we have seen in
Section 3, it must be smooth if $(x,y)$ corresponds to a smooth pair. Now, if
$x$ lies in a certain range, such as $[2^{240},2^{256}]$ in our running
example, then a small $\Delta$ means that there is a larger range for possible
values of $y$. Thus we expect that having smaller $\Delta$ increases the
probability that some $y$ solving the Pell equation will meet our criteria.
Assuming that this is the only property of $\Delta$ affecting the chance of
having a solution in the desired range, the probability of having a solution
scales like $\frac{1}{\sqrt{\Delta}}$. Since $\Delta$ can range from $3$ to
the square of the maximum bound on the size of smooth pairs we are interested
in ($2^{512}$ in our example), the effect of this combinatorial heuristic is
very strong.

This reasoning agrees with experiments. In our search for pairs in the range
$[2^{240},2^{256}]$, we solved all Pell equations with $\Delta\leq 2^{44}$.
This computation ran in parallel on a 64-core computer for 3 days and found a
$19949$-smooth, $215$-bit pair $(m,m+1)$ corresponding to a fundamental
solution with

$\displaystyle m$
$\displaystyle=51963397732665557125190357543988479960188331933248699616266017360$
$\displaystyle=2^{4}\cdot 3^{5}\cdot 5\cdot 7\cdot 11^{2}\cdot 647\cdot
911\cdot 919^{2}\cdot 1103^{2}\cdot 2099^{2}\cdot 2203^{2}\cdot 2423^{2}\cdot
5279^{2}\cdot 8641^{2}\cdot 19949,$ $\displaystyle m+1$
$\displaystyle=19^{2}\cdot 23^{2}\cdot 47^{2}\cdot 277^{2}\cdot 359^{2}\cdot
419^{2}\cdot 541^{2}\cdot 887^{2}\cdot 1993^{2}\cdot 4549^{2}\cdot
4813^{2}\cdot 12721^{2}.$

While this pair is comparatively small to many others found, it is different
in that it corresponds directly to a fundamental solution. Certain
cryptographic applications require smooth pairs $(m,m+1)$ for which the sum
$m+(m+1)=2m+1$ is prime. This property is not achievable by pairs
corresponding to $n$-th solutions unless $n=2^{k}$ for some integer $k$ by
Lemma 1.

Checking the $n$-th solutions of the pairs found using the previously
discussed methods for $n\in\\{1,2,\dots,12\\}$ also found a pair

$\displaystyle m$
$\displaystyle=13751398930221343029252446890555634360426232035770955214339102967231667741706727080$
$\displaystyle=2^{3}\cdot 3^{5}\cdot 5\cdot 59\cdot 103^{2}\cdot 113^{2}\cdot
1697\cdot 2381\cdot 2383^{2}\cdot 2399^{2}\cdot 3623^{2}\cdot 4733\cdot
9151^{2}\cdot 10607^{2}$ $\displaystyle\quad\cdot 16267\cdot 18059^{2}\cdot
18289\cdot 23603,$ $\displaystyle m+1$ $\displaystyle=11^{4}\cdot 17^{2}\cdot
53^{2}\cdot 83^{2}\cdot 109^{2}\cdot 227^{2}\cdot 263^{2}\cdot 347^{2}\cdot
373^{2}\cdot 599^{2}\cdot 2341^{2}\cdot 7883^{2}\cdot 10223^{2}\cdot
10883^{2}\cdot 12511^{2},$

which has $273$ bits and is $23603$-smooth. It corresponds to the sixth
solution to a Pell equation. This pair, though larger than the target
cryptographic size of 256 bits, is strong evidence for the efficiency of this
algorithm. It is both larger and more smooth than the second pair presented in
3.1, and is far larger than any previously known for its smoothness bound.
This pair is also a nice example of how the methods in 3.1 can be combined
well with other algorithms to yield $B$-smooth pairs larger than either
algorithm could individually.

###  3.4 The small primes heuristic

While the above heuristic of choosing small values of $\Delta$ works well, a
limitation of it is that there are relatively few small $\Delta$, and the
probability of each Pell equation resulting in smooth pairs in the desired
range decreases strongly as $\Delta$ increases. Furthermore, once we relax the
restrictions on the maximum size of $\Delta$, exhausting all of its possible
values becomes computationally infeasible. Another heuristic that narrows the
search space is to choose values of $\Delta$ with smaller prime factors.

Among coefficients $\Delta$ of roughly the same size, those having more (and
hence on average smaller) prime factors correspond to Pell equations that are
more likely to result in smooth pairs in the desired range. We again give a
combinatorial argument. If $p\mid\Delta$, then from the Pell equation, we have

$x^{2}-1\equiv 0\pmod{p}\iff(x-1)(x+1)\equiv 0\pmod{p}\iff x\equiv\pm
1\pmod{p}$

because the integers modulo $p$ form a field. Suppose that $\Delta$ is
divisible by $n$ primes, say $\Delta=p_{1}\cdot p_{2}\cdot\ldots\cdot p_{n}$.
Any integer $x$ that solves the Pell equation must satisfy the above
condition, so we can associate to it a sequence $s\in\\{-1,+1\\}^{n}$, where
the $i$th element of the sequence is $x\text{ mod }p_{i}$.

Now suppose we fix such a sequence $s$, and consider all $x$ in our desired
range that are associated to $s$. By the Chinese remainder theorem, there is
exactly one such $x$ in every interval of length $\Delta$. This gives further
evidence for the above heuristic of choosing $\Delta$ small, because there are
more $x$s in our desired range with the potential of solving the Pell
equation.

Further, notice that $x$ in a solution to the Pell equation with coefficient
$\Delta$ can correspond to any of the length $n$ sequences. Since there are
$2^{n}$ such sequences and there are (up to $\pm 1$) the same number of $x$’s
in the range corresponding to each, the implicit search space in which we look
for solutions when solving Pell equations with coefficients $\Delta$ of a
specified size is much larger when $\Delta$ has more prime factors.

Small searches with this method over coefficients larger than $2^{42}$ yielded
a $201$-bit, $29881$-smooth pair

$\displaystyle m$
$\displaystyle=2317395102010090979961970844394697249269453177012017537196644$
$\displaystyle=2^{2}\cdot 13\cdot 53^{2}\cdot 113\cdot 139^{2}\cdot
269^{2}\cdot 347^{2}\cdot 509^{4}\cdot 569^{2}\cdot 743^{2}\cdot
12823^{2}\cdot 14149\cdot 29881,$ $\displaystyle m+1$
$\displaystyle=3^{9}\cdot 5\cdot 23^{2}\cdot 43^{2}\cdot 167^{2}\cdot
179^{2}\cdot 227\cdot 349^{2}\cdot 613^{2}\cdot 661^{2}\cdot 2927^{2}\cdot
4099^{2}\cdot 6421^{2}.$

It corresponds to $\Delta=13\cdot 113\cdot 14149\cdot 29881$. This pair is
smaller than those given previously, which is expected given the increase in
the sizes of $\Delta$ searched over. However, it is still close to
cryptographic size, indicating that more precise heuristics could extend our
methods to search effectively over large values of $\Delta$.

##  4 Discussion

We have seen how solutions to Pell equations give both a theoretical and
practical approach to finding twin smooth integers, including several new
algorithms that have produced new pairs. We conclude by discussing two
observations for exploration in further research.

Based on experiments, our algorithms are more likely to find smooth pairs
$(m,m+1)$ with $2m+1$ composite, which is undesirable for some applications.
Based on Lemma 1, the only pairs that can sum to a prime are those
corresponding to $n$-th solutions with $n=2^{k}$ for integer $k$. Generally
such pairs are easy to exhaust with the algorithms given in 3.1 for $k$ even
moderately large; for our example we found that no pairs in
$[2^{240},2^{256}]$ exist with $k\geq 3$. Therefore, it is interesting to
consider whether we can improve our search through fundamental solutions to
directly find large pairs more often.

One approach is to consider the structure of the continued fractions
expansions. Because it was not useful to our algorithms, we have omitted
discussion of many remarkable structures in the continued fractions of square
roots, such as the periodicity of the quotients. Because the size of
fundamental solutions is determined by the length of the period and size of
the quotients of $\sqrt{\Delta}$, one idea is to choose coefficients $\Delta$
by fixing the length of quotients in the corresponding continued fraction
expansion.

Unfortunately, while this approach guarantees that we have pairs $(x,y)$
satisfying the Pell equation $x^{2}-1=2\Delta y^{2}$ with $x$ in a given
range, however, it does not guarantee smoothness of $y$ or $D$. It is
interesting to consider whether a clever way of choosing the period lengths
and continued fractions quotients could result in a high probability of $y$
and $D$ being smooth, thus resulting in an efficient algorithm for finding
twin smooth pairs. This is a promising approach because of the strong and
well-studied structure of continued fractions. An overview of the theory of
continued fractions can be found in [6], chapters 3 and 6.

The tension described above between searching for solutions to Pell equations
(or equivalently smooth pairs) with a fixed smoothness bound and searching for
solutions (or smooth pairs) in a given size range is common across many
algorithms for finding twin smooth integers. This difficulty seems intrinsic,
as the problem of predicting the size of the solution to a Pell equation is
related to questions about computing class numbers of quadratic fields, which
remain highly mysterious.

Solving the Pell equation is equivalent to finding the fundamental unit in the
associated real quadratic field $\mathbb{Q}(\sqrt{D})$. Because of the
structure of units in its ring of integers, given the fundamental unit
$\epsilon_{D}=t+u\sqrt{D}$, then either $x=t,y=u$ or $x=t^{2}+u^{2},y=2tu$
solve the Pell equation $x^{2}-Dy^{2}=1$. The two cases arise from the fact
that the fundamental unit can have norm $\pm 1$. For simplicity we focus on
the first case; the analysis is analogous for the second. Since
$x\approx\sqrt{D}y$, $\epsilon_{D}\approx 2x$.

Dirichlet’s class number formula for quadratic fields relates the size of the
fundamental unit to that of the fields’ class number $h(D)$,

$h(D)=\frac{\mathcal{L}(\chi,1)\sqrt{D}}{2\text{log}(\epsilon_{D})}.$

Based on both experimental evidence and theoretical bounds,
$\mathcal{L}(\chi,1)$ grows extremely slowly with $D$, see for example the
following plot.

Therefore if we want to find coefficients $D$ with corresponding $x\approx M$,
we want

$h(D)\approx\frac{\sqrt{D}}{\text{log}(M)}.$

Since $M$ is fixed, choosing smooth, squarefree coefficients $D$ with Pell
equation solutions in the right range is a problem equivalent to choosing $D$
with class number roughly $\sqrt{D}$. This problem is likely very difficult,
since little is know about the sizes of class numbers of real quadratic
fields. For example, Gauss conjectured in 1801 that there are infinitely many
real quadratic fields with class number 1; this conjecture remains open. In
the absence of new breakthroughs in predicting the sizes of class numbers,
heuristics such as those presented in 3.3 and 3.4 appear to be the best way to
search for large twin smooth integers corresponding to fundamental solutions.

## References

  * [1] Giacomo Bruno, Maria Corte-Real Santos, Craig Costello, Jonathan Komada Eriksen, Michael Naehrig, Michael Meyer, and Bruno Sterner.  Cryptographic smooth neighbors.  Cryptology ePrint Archive, Paper 2022/1439, 2022.  https://eprint.iacr.org/2022/1439.
  * [2] M. Naehrig C. Costello, M. Meyer.  Sieving for twin smooth integers with solutions to the prouhet-tarry-escott problem.
  * [3] J. B. Conrey, M. A. Holmstrom, and T. L. McLaughlin.  Smooth neighbors.  Experimental Mathematics, 22(2):195–202, 2013.
  * [4] C. Costello.  B-SIDH: supersingular isogeny Diffie-Hellman using twisted torsion.  In S. Moriai and H. Wang, editors, ASIACRYPT 2020, volume 12492 of Lecture Notes in Computer Science, pages 440–463. Springer, 2020\.
  * [5] L. De Feo, D. Kohel, A. Leroux, C. Petit, and B. Wesolowski.  SQISign: Compact post-quantum signatures from quaternions and isogenies.  In Shiho Moriai and Huaxiong Wang, editors, ASIACRYPT 2020, volume 12491 of Lecture Notes in Computer Science, pages 64–93. Springer, 2020.
  * [6] M. J. Jacobson and H. C. Williams.  Solving the Pell Equation.  Springer, 2009.
  * [7] D. H. Lehmer.  On a problem of Störmer.  Illinois Journal of Mathematics, 8(1):57–79, 1964.
  * [8] H. W. Lenstra.  Solving the Pell equation.  Notices Amer. Math. Soc., 49(2), 2002.
  * [9] F. Luca and F. Najman.  On the largest prime factor of $x^{2}-1$.  Math. Comp., 80(273):429–435, 2011.
  * [10] F. Luca and F. Najman.  Errata to ”On the largest prime factor of $x^{2}-1$”.  Math. Comp., 83(285):337, 2014.
  * [11] K. H. Rosen.  Elementary Number Theory and Its Applications.  Addison Wesley Publishing Company, 1984.
  * [12] C. Størmer.  Quelques théorèmes sur l’équation de Pell $x^{2}-dy^{2}=\pm 1$ et leurs applications.  Christiania Videnskabens Selskabs Skrifter, Math. Nat. Kl, (2):48, 1897.
  * [13] J. Unger.  Solving Pell’s equation with continued fractions.  2009\.

##  5 Code

The code can be found at
https://github.com/lightningwolf111/twin_smooth_integers. Large searches were
done using the parallelized C implementations and use the GMP library.
Experiments with many other algorithms can be found in the repository as well,
and were implemented in python using the SAGE library.