arXiv:1712.03693v1 [cs.SC] 11 Dec 2017

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION
USING CYCLOTOMIC COEFFICIENT RINGS
DAVID HARVEY AND JORIS VAN DER HOEVEN
Abstract. We present√an algorithm that computes the product of two n-bit
∗
integers in O(n log n (4 2)log n ) bit operations. Previously, the best known
∗
log
n
bound was O(n log n 6
). We also prove that for a fixed prime p, polynomi∗
als in Fp [X] of degree n may be multiplied in O(n log n 4log n ) bit operations;
∗
the previous best bound was O(n log n 8log n ).

1. Introduction
In this paper we present new complexity bounds for multiplying integers and
polynomials over finite fields. Our focus is on theoretical bounds rather than practical algorithms. We work in the deterministic multitape Turing model [19], in
which time complexity is defined by counting the number of steps, or equivalently,
the number of ‘bit operations’, executed by a Turing machine with a fixed, finite
number of tapes. The main results of the paper also hold in the Boolean circuit
model.
The following notation is used throughout. For x ∈ R, we denote by log∗ x the
iterated logarithm, that is, the least non-negative integer k such that log◦k x 6 1,
where log◦k x := log · · · log x (iterated k times). For a positive integer n, we define
lg n := max(1, ⌈log2 n⌉); in particular, expressions like lg lg lg n are defined and
take positive values for all n > 1. We denote the n-th cyclotomic polynomial by
φn (X) ∈ Z[X], and the Euler totient function by ϕ(n).
All absolute constants in this paper are in principle effectively computable. This
includes the implied constants in all uses of O(·) notation.
1.1. Integer multiplication. Let M(n) denote the number of bit operations required to multiply two n-bit integers. For over 35 years, the best known bound for
M(n) was that achieved by the Schönhage–Strassen algorithm [23], namely
M(n) = O(n lg n lg lg n).

(1.1)

In 2007, Fürer described an asymptotically faster algorithm that achieves
∗

M(n) = O(n lg n KZlog

n

)

(1.2)

for some unspecified constant KZ > 1 [11, 12]. His algorithm reduces a multiplication of size n to a large collection of multiplications of size exponentially smaller
∗
than n; these smaller multiplications are handled recursively. The KZlog n term may
be understood roughly as follows: the number of recursion levels is log∗ n + O(1),
and the constant KZ measures the amount of ‘data expansion’ that occurs at each
level, due to phenomena such as zero-padding.
Immediately following Fürer’s work, De, Kurur, Saha and Saptharishi described
a variant based on modular arithmetic [10], instead of the approximate complex

2

DAVID HARVEY AND JORIS VAN DER HOEVEN

arithmetic used by Fürer. Their algorithm also achieves (1.2), again for some
unspecified KZ > 1.
The first explicit value for KZ was given by Harvey, van der Hoeven and Lecerf,
who described an algorithm that achieves (1.2) with KZ = 8 [17]. Their algorithm
borrows some important ideas from Fürer’s work, but also differs in several respects.
In particular, their algorithm has no need for the ‘fast’ roots of unity that were the
cornerstone of Fürer’s approach (and of the variant of [10]). The main presentation
in [17] is based on approximate complex arithmetic, and the paper includes a sketch
of a variant based on modular arithmetic that also achieves KZ = 8.
In a recent preprint, the first author announced that in the complex arithmetic
case, the constant may be reduced to KZ = 6, by taking advantage of new methods
for truncated integer multiplication [14]. This improvement does not seem to apply
to the modular variants.
The first main result of this paper is the following further improvement.
Theorem 1.1. There is an integer multiplication algorithm that achieves
√
∗
M(n) = O(n lg n (4 2)log n ).
√
In other words, (1.2) holds with KZ = 4 2 ≈ 5.657. The proof is given in
Section 5. The new algorithm works with modular arithmetic throughout; curiously,
there is no obvious analogue based on approximate complex arithmetic.
There have been several proposals in the literature for algorithms that achieve
KZ = 4 under various unproved number-theoretic conjectures: see [17, §9], [15],
and [8]. Whether KZ = 4 can be reached unconditionally remains an important
open question.
1.2. Polynomial multiplication over finite fields. For a prime p, let Mp (n)
denote the number of bit operations required to multiply two polynomials in Fp [X]
of degree less than n. The optimal choice of algorithm for this problem depends
very much on the relative size of n and p.
If n is not too large compared to p, say lg n = O(lg p), then a reasonable choice
is Kronecker substitution: one lifts the polynomials to Z[X], packs the coefficients
of each polynomial into a large integer (i.e., evaluates at X = 2b for b := 2 lg p +
lg n), multiplies these large integers, unpacks the resulting coefficients to obtain the
product in Z[X], and finally reduces the output modulo p. This leads to the bound
log∗ (n lg p)

Mp (n) = O(M(n lg p)) = O(n lg p lg(n lg p)KZ

),

(1.3)

where KZ is any admissible constant in (1.2). To the authors’ knowledge, this is
the best known asymptotic bound for Mp (n) in the region lg n = O(lg p).
When n is large compared to p, the situation is starkly different. The Kronecker
substitution method leads to poor results, due to coefficient growth in the lifted
product: for example, when p is fixed, Kronecker substitution yields
∗

Mp (n) = O(M(n lg n)) = O(n(lg n)2 KZlog

n

).

For many years, the best known bound in this regime was that achieved by the
algebraic version of the Schönhage–Strassen algorithm [23, 22], namely
Mp (n) = O(n lg n lg lg n lg p + n lg n M(lg p)).

(1.4)

The first term arises from performing O(n lg n lg lg n) additions in Fp , and the
second term from O(n lg n) multiplications in Fp . (In fact, this sort of bound holds

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

3

for polynomial multiplication over quite general rings [7].) For fixed p, this is faster
than the Kronecker substitution method by a factor of almost lg n. The main
reason for its superiority is that it exploits the modulo p structure throughout the
algorithm, whereas the Kronecker substitution method forgets this structure in the
very first step.
After the appearance of Fürer’s algorithm, it was natural to ask whether a Fürertype bound could be proved for Mp (n), in the case that n is large compared to p.
This question was answered in the affirmative by Harvey, van der Hoeven and
Lecerf, who gave an algorithm that achieves
∗

Mp (n) = O(n lg p lg(n lg p) 8log

(n lg p)

),

uniformly for all n and p [18]. This is a very elegant bound; however, written in this
way, it obscures the fact that the constant 8 plays two quite different roles in the
complexity analysis. One source of the value 8 is the constant KZ = 8 arising from
the integer multiplication algorithm mentioned above, but there is also a separate
constant KF = 8 arising from the polynomial part of the algorithm. There is no
particular reason to expect that KZ = KF , and it is somewhat of a coincidence that
they have the same numerical value in [18].
To clarify the situation, we mention that one may derive a complexity bound
for the algorithm of [18] under the assumption that one has available an integer
multiplication algorithm achieving (1.2) for some KZ > 1, where possibly KZ 6= 8.
Namely, one finds that
max(0,log∗ n−log∗ p)

Mp (n) = O(n lg p lg(n lg p) KF

∗

KZlog

p

)

(1.5)

where KF = 8 (we omit the proof). The second main result of this paper, proved
in Section 4, is the following improvement in the value of KF .
Theorem 1.2. Let KZ > 1 be any √
constant for which (1.2) holds (for example, by
Theorem 1.1, one may take KZ = 4 2). Then there is a polynomial multiplication
algorithm that achieves
Mp (n) = O(n lg p lg(n lg p) 4max(0,log

∗

n−log∗ p)

∗

KZlog

p

),

(1.6)

uniformly for all n > 1 and all primes p.
In other words, (1.5) holds with KF = 4. In particular, for fixed p, one can
∗
multiply polynomials in Fp [X] of degree n in O(n lg n 4log n ) bit operations.
Theorem 1.2 may be generalised in various ways. We briefly mention a few possibilities along the lines of [18, §8] (no proofs will be given). First, we may obtain analogous bit complexity bounds for multiplication in Fpa [X] and (Z/pa Z)[X] for a > 1,
and in (Z/mZ)[X] for arbitrary m > 1 (see Theorems 8.1–8.3 in [18]). We may
also obtain complexity bounds for polynomial multiplication in various algebraic
complexity models. For example, we may construct a straight-line program that
multiplies two polynomials in A[X] of degree less than n, for any Fp -algebra A, us∗
∗
ing O(n lg n 4log n ) additions and scalar multiplications and O(n 2log n ) nonscalar
multiplications (compare with [18, Thm. 8.4]).
1.3. Overview of the new algorithms. To explain the new approach, let us first
recall the idea behind the polynomial multiplication algorithm of [18].
Consider a polynomial multiplication problem in Fp [X], where the degree n is
very large compared to p. By splitting the inputs into chunks, we convert this

4

DAVID HARVEY AND JORIS VAN DER HOEVEN

to a bivariate multiplication problem in Fp [Y, Z]/(f (Y ), Z m − 1), for a suitable
integer m and irreducible polynomial f ∈ Fp [Y ]. This bivariate product is handled
by means of DFTs (discrete Fourier transforms) of length m over Fp [Y ]/f . The key
innovation of [18] was to choose deg f so that pdeg f − 1 is divisible by many small
primes, so many, in fact, that their product is comparable to n, even though deg f
itself is exponentially smaller. This is possible thanks to a number-theoretic result
of Adleman, Pomerance and Rumely [1], building on earlier work of Prachar [21].
Taking m to be a product of many of these primes, we obtain m | pdeg f − 1, and
hence Fp [Y ]/f contains a root of unity of order m. As m is highly composite, each
DFT of length m may be converted to a collection of much smaller DFTs via the
Cooley–Tukey method. These in turn are converted into multiplication problems
using Bluestein’s algorithm. These multiplications, corresponding to exponentially
smaller values of n, are handled recursively.
The recursion continues until n becomes comparable to p. The number of recursion levels during this phase is log∗ n − log∗ p + O(1), and the constant KF = 8
represents the expansion factor at each recursion level. When n becomes comparable to p, the algorithm switches strategy to Kronecker substitution combined with
∗
ordinary integer multiplication. This phase contributes the KZlog p term.
It was pointed out in [16, §8] that the value of KF can be improved to KF = 4
if one is willing to accept certain unproved number-theoretic conjectures, including
Artin’s conjecture on primitive roots. More precisely, under these conjectures, one
may find an irreducible f of the form f (Y ) = Y α−1 + · · · + Y + 1, where α is prime,
so that Fp [Y, Z]/(f (Y ), Z m − 1) is a direct summand of Fp [Y, Z]/(Y α − 1, Z m − 1).
This last ring is isomorphic to Fp [X]/(X αm − 1), and one may use this isomorphism
to save a factor of two in zero-padding at each recursion level. These savings lead
directly to the improved value for KF .
To prove Theorem 1.2, we will pursue a variant of this idea. We will take f to be
a cyclotomic polynomial φα (Y ) for a judiciously chosen integer α (not necessarily
prime). Since φα | Y α − 1, we may use the above isomorphism to realise the same
economy in zero-padding as in the conjectural construction of [16, §8]. However,
unlike [16], we do not require that f be irreducible in Fp [Y ]. Thus Fp [Y ]/f is no
longer in general a field, but a direct sum of fields. The situation is reminiscent of
r
Fürer’s algorithm, in which the coefficient ring C[Y ]/(Y 2 + 1) is not a field, but a
direct sum of copies of C. The key technical contribution of this paper is to show
that we have enough control over the factorisation of φα in Fp [Y ] to ensure that
Fp [Y ]/φα contains suitable principal roots of unity. This approach avoids Artin’s
conjecture and other number-theoretic difficulties, and enables us to reach KF = 4
unconditionally. The construction of α is the subject of Section 3, and the main
polynomial multiplication algorithm is presented in Section 4.
Let us now outline how we go about proving Theorem 1.1 (the integer case).
The algorithm is heavily dependent on the polynomial multiplication algorithm
just sketched. We take the basic problem to be multiplication in Z/(2n − 1)Z, for
arbitrary positive n. We choose a collection of small primes p, each having around
2 lg lg n bits, and whose product P has (lg n)1+o(1) bits. By cutting the input integers into many small chunks, we convert to a multiplication in (Z/P Z)[X]/(X N −1)
for a suitable N ≈ 2n/ lg P . One technical headache is that n is not necessarily
divisible by N ; following [17], we deal with this by adapting an idea of Crandall

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

5

and Fagin [9]. Next, by the Chinese remainder theorem, we reduce to multiplying in Fp [X]/(X N − 1) for each p separately. This is reminisicent of Pollard’s
algorithm [20], but instead of using three primes, here the number of primes grows
with n. At this stage, the coefficient size lg p is doubly exponentially smaller than N .
We perform these multiplications in Fp [X]/(X N − 1) by applying two recursion levels of the polynomial multiplication algorithm of Theorem 1.2. This reduces the
problem to a collection of multiplication problems in Fp [X], each doubly exponentially smaller than the original problem. Using Kronecker substitution, these are
′
converted back to multiplications in Z/(2n − 1)Z, where n′ is doubly exponentially
smaller than n, and the algorithm is applied recursively.
In effect, each recursive call in the new integer multiplication algorithm corresponds to two recursion levels of the existing Fürer-type algorithms. The speedup
relative to [17] may be understood as follows. In the algorithm of [17], at each
recursion level we incur a factor of two in overhead due to the zero-padding that
occurs when we split the inputs into small chunks. In the new algorithm, the passage from Z/P Z to Fp manages the same exponential size reduction without any
zero-padding. This roughly corresponds to saving a factor of two √
at every second
∗
recursion level of the algorithm of [17], and explains the factor of ( 2)log n overall
speedup.
2. Preliminaries
2.1. Logarithmically slow functions. Let x0 ∈ R, and let Φ : (x0 , ∞) → R
be a smooth increasing function. We recall from [17, §5] that Φ is said to be
logarithmically slow if there exists an integer ℓ > 0 such that
(log◦ℓ ◦ Φ ◦ exp◦ℓ )(x) = log x + O(1)
5

as x → ∞. For example, the functions log(5x), 5 log x, (log x)5 , and 2(log log x) are
logarithmically slow, with ℓ = 0, 1, 2, 3 respectively.
We will always assume that x0 is chosen large enough to ensure that Φ(x) 6 x−1
for all x > x0 . According to [17, Lemma 2], this is possible for any logarithmically
slow function, and it implies that the iterator Φ∗ (x) := min{k > 0 : Φ◦k (x) 6 x0 }
is well-defined on R. It is shown in [17, Lemma 3] that this iterator satisfies
Φ∗ (x) = log∗ x + O(1)

(2.1)

as x → ∞. In other words, logarithmically slow functions are more or less indistinguishable from log x, as far as iterators are concerned.
As in [17] and [18], we will use logarithmically slow functions to measure size
reduction in multiplication algorithms. The typical situation is that we have a
function T (n) measuring the (normalised) cost of a certain multiplication algorithm
for inputs of size n; we reduce the problem to a collection of problems of size
ni < Φ◦κ (n) for some κ > 1, leading to a bound for T (n) in terms of the various
T (ni ). Applying the reduction recursively, we wish to convert these bounds into an
explicit asymptotic estimate for T (n). This is achieved via the following ‘master
theorem’.
Proposition 2.1. Let K > 1, B > 0, and let ℓ > 0 and κ > 1 be integers. Let
x0 > exp◦ℓ (1), and let Φ : (x0 , ∞) → R be a logarithmically slow function such that
Φ(x) 6 x − 1 for all x > x0 . Assume that x1 > x0 is chosen so that Φ◦κ (x) is

6

DAVID HARVEY AND JORIS VAN DER HOEVEN

defined for all x > x1 . Then there exists a positive constant C (depending on x0 ,
x1 , Φ, K, B, ℓ and κ) with the following property.
Let σ > x1 and L > 0. Let S ⊆ R, and let T : S → R> be any function
satisfying the following recurrence. First, T (y) 6 L for all y ∈ S, y 6 σ. Second,
for all y ∈ S, y >P
σ, there exist y1 , . . . , yd ∈ S with yi 6 Φ◦κ (y), and weights
γ1 , . . . , γd > 0 with i γi = 1, such that

X
d
B
T (y) 6 K 1 +
γi T (yi ) + L.
log◦ℓ y i=1
For all y ∈ S, y > σ, we then have

T (y) 6 CL(K 1/κ )log

∗

y−log∗ σ

.

Proof. The special case κ = 1, x1 = x0 is exactly [17, Prop. 8]. We indicate
briefly how the proof of [17, Prop. 8] must be modified to obtain this more general
statement.
The first two paragraphs of the proof of [17, Prop. 8] may be read verbatim. In
the third paragraph, the inductive statement is changed to
T (y) 6 E1 · · · Ej L(K ⌈j/κ⌉ + · · · + K + 1),

where j := Φ∗σ (y). The inductive step is modified slightly: for 0 < j < κ we use
the fact that yi 6 σ, and for j > κ the fact that Φ∗σ (yi ) 6 Φ∗σ (Φ◦κ (y)) = Φ∗σ (y) − κ.
With these changes, the proof given in [17] goes through without difficulty.

2.2. Discrete Fourier transforms. Let n > 1 and let R be a commutative ring
in which n is invertible. A P
principal n-th root of unity is an element ω ∈ R such
n−1 ij
that ω n = 1 and such that j=0
ω = 0 for i = 1, 2, . . . , n − 1. If m is a divisor
n/m
of n, then ω
is easily seen to be a principal m-th root of unity.
Fix a principal n-th root of unity ω. The discrete Fourier transform (DFT) of the
sequence (a0 , . . . , an−1 ) ∈ Rn with respect to ω is the sequence (â0 , . . . , ân−1 ) ∈ Rn
Pn−1
Pn−1
defined by âj := i=0 ω ij ai . Equivalently, âj = A(ω j ) where A = i=0 ai X i ∈
R[X]/(X n − 1).
The inverse DFT recovers (a0 , . . . , an−1 ) from (â0 , . . . , ân−1 ). Computationally
it corresponds to a DFT with respect to ω −1 , followed by a division by n, because
n−1
n−1 n−1
1 X −kj
1 X X (i−k)j
ω
âj =
ω
ai = ak ,
n j=0
n i=0 j=0

k = 0, . . . , n − 1.

DFTs may be used to implement cyclic convolutions. Suppose that we wish to
compute C := AB where A, B ∈ R[X]/(X n −1). We first perform DFTs to compute
A(ω i ) and B(ω i ) for i = 0, . . . , n − 1. We then compute C(ω i ) = A(ω i )B(ω i ) for
each i, and finally perform an inverse DFT to recover C ∈ R[X]/(X n − 1).
This strategy may be generalised to handle a multidimensional cyclic convolution, that is, to compute C := AB for
A, B ∈ R[X1 , . . . , Xd ]/(X1n1 − 1, . . . , Xdnd − 1).

For this, we require that each nk be invertible in R, and that R contain a principal
nk -th root of unity ωk for each k. Let n := n1 · · · nd . We first perform multidimensional DFTs to evaluate A and B at the n points {(ω1j1 , . . . , ωdjd ) : 0 6 jk < nk }.
We then multiply pointwise, and finally recover C via a multidimensional inverse
DFT.

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

7

Each multidimensional DFT may be reduced to a collection of one-dimensional
DFTs as follows. We first compute A(X1 , . . . , Xd−1 , ωdj ) ∈ R[X1 , . . . , Xd−1 ] for
each j = 0, . . . , nd − 1; this involves n/nd DFTs of length nd . We then recursively
jd−1
). Altogether,
evaluate each of these polynomials at the n/nd points (ω1j1 , . . . , ωd−1
this strategy involves computing n/nk DFTs of length nk for each k = 1, . . . , d.
Finally, we briefly recall Bluestein’s method [5] for reducing a (one-dimensional)
DFT to a convolution problem (see also [17, §2.5]). Let n > 1 be odd and let ω ∈ R
be a principal n-th root of unity. Set ξ := ω (n+1)/2 , so that ξ 2 = ω and ξ n = 1.
Then computing the DFT of a given sequence (a0 , . . . , an−1 ) ∈ Rn with respect
to ω reduces to computing the product of the polynomials
f (Z) :=

n−1
X
i=0

2

ξ i ai Z i ,

g(Z) :=

n−1
X

2

ξ −i Z i

i=0

n

in R[Z]/(Z − 1), plus O(n) auxiliary multiplications in R. Notice that g(Z) is
fixed and does not depend on the input sequence.
2.3. The Crandall–Fagin algorithm. Consider the problem of computing a
‘cyclic’ integer product of length n, that is, a product uv where u, v ∈ Z/(2n − 1)Z.
If N and P are positive integers such that N | n and lg P > 2n/N + lg N , then
we may reduce the given problem to multiplication in (Z/P Z)[X]/(X N − 1), by
cutting up the integers into chunks of n/N bits. In this section we briefly recall a
variant [17, §9.2] of an algorithm due to Crandall and Fagin [9] that achieves the
same reduction without the assumption that N | n.
Assume that N 6 n and lg P > 2⌈n/N ⌉ + lg N + 1, and that we have available
some θ ∈ Z/P Z with θN = 2. (This θ plays the same role as the real N -th
root of 2 in the original Crandall–Fagin algorithm.) Set ei := ⌈ni/N ⌉ and ci :=
N ei − ni. Observe that ei+1 − ei = ⌊n/N ⌋ or ⌈n/N ⌉ for each i. Decompose
PN −1
PN −1
the inputs as u = i=0 2ei ui and v = i=0 2ei vi where 0 6 ui , vi < 2ei+1 −ei
P −1 ci
(i.e., a decomposition with respect to a ‘variable base’). Set U (X) := N
i=0 θ ui
PN −1 ci
and V (X) := i=0 θ vi , regarded as polynomials in (Z/P Z)[X]/(X N − 1), and
compute the product W (X) := U (X)V (X). Then one finds (see [17, §9.2]) that
PN −1
the product uv may be recovered by the formula uv = i=0 2ei wi (mod 2n − 1),
PN −1
where the wi are integers in [0, P ) defined by W (X) = i=0 θci wi .
To summarise, the problem of computing uv reduces to computing a product
in (Z/P Z)[X]/(X N − 1), together with O(N ) auxiliary multiplications in Z/P Z,
and O(N (lg n)2 + N lg P ) bit operations to compute the ei and to handle the final
overlap-add phase (again, see [17, §9.2] for details).
2.4. Data layout. In this section we discuss several issues relating to the layout
of data on the Turing machine tapes.
Integers will always be stored in the standard binary representation. If n is a
positive integer, then elements of Z/nZ will always be stored as residues in the
range 0 6 x < n, occupying lg n bits of storage.
If R is a ring and f ∈ R[X] is a polynomial of degree n > 1, then an element
of R[X]/f (X) will always be represented as a sequence of n coefficients in the
standard monomial order. This convention is applied recursively, so for rings of the
type (R[Y ]/f (Y ))[X]/g(X), the coefficient of X 0 is stored first, as an element of
R[Y ]/f (Y ), followed by the coefficient of X 1 , and so on.

8

DAVID HARVEY AND JORIS VAN DER HOEVEN

A multidimensional array of size nd × · · · × n1 , whose entries occupy b bits
each, will be stored as a linear array of bn1 · · · nd bits. The entries are ordered
lexicographically in the order (0, . . . , 0, 0), (0, . . . , 0, 1), . . . , (nd − 1, . . . , n1 − 1). In
particular, an element of (· · · (R[X1 ]/f1 (X1 )) · · · )[Xd ]/fd (Xd ) is represented as an
nd × · · · × n1 array of elements of R. We will generally prefer the more compact
notation R[X1 , . . . , Xd ]/(f1 (X1 ), . . . , fd (Xd )).
There are many instances where an n×m array must be transposed so that its entries can be accessed efficiently either ‘by columns’ or ‘by rows’. Using the algorithm
of [6, Lemma 18], such a transposition may be achieved in O(bnm lg min(n, m)) bit
operations, where b is the bit size of each entry. (The idea of the algorithm is to split
the array in half along the short dimension, and transpose each half recursively.)
One important application is the following result, which estimates the data rearrangement cost associated to the the Agarwal–Cooley method [2] for converting
between one-dimensional and multidimensional convolution problems (this is closely
related to the Good–Thomas DFT algorithm [13, 26]).
Lemma 2.2. Let n, m > 2 be relatively prime, and let R be a ring whose elements
are represented using b bits. There exists an isomorphism
∼ R[Y, Z]/(Y n − 1, Z m − 1)
R[X]/(X nm − 1) =
that may be evaluated in either direction in O(bnm lg min(n, m)) bit operations.
Proof. Let c := m−1 mod n, and let
β : R[X]/(X nm − 1) → R[Y, Z]/(Y n − 1, Z m − 1)

denote the homomorphism that maps X to Y c Z, and acts as the identity
on R. SupPnm−1
pose that we wish to compute β(F ) for some input polynomial F = k=0 Fk X k ∈
R[X]/(X nm − 1). Interpreting the list (F0 , . . . , Fnm−1 ) as an n × m array, the
(i, j)-th entry corresponds to Fim+j . After transposing the array, which costs
O(bnm lg min(n, m)) bit operations, we have an m × n array, whose (j, i)-th entry is Fim+j . Now for each j, cyclically permute the j-th row by (jc mod n)
slots; altogether this uses only O(bnm) bit operations. The result is an m × n
array whose (j, i)-th entry is F(i−jc mod n)m+j , which is exactly the coefficient of
Y ((i−jc)m+j)c Z (i−jc)m+j = Y i Z j in β(F ). The inverse map β −1 may be computed
by reversing this procedure.

Corollary 2.3. Let n1 , . . . , nd > 2 be pairwise relatively prime, let n := n1 · · · nd ,
and let R be a ring whose elements are represented using b bits. There exists an
isomorphism
n
R[X]/(X n − 1) ∼
= R[X1 , . . . , Xd ]/(X n1 − 1, . . . , X d − 1)
1

d

that may be evaluated in either direction in O(bn lg n) bit operations.

Proof. Using Lemma 2.2, we may construct a sequence of isomorphisms
∼ R[X1 , W2 ]/(X n1 − 1, W n2 ···nd − 1)
R[X]/(X n1 ···nd − 1) =
2

1

∼
= R[X1 , X2 , W3 ]/(X1n1 − 1, X2n2 − 1, W3n3 ···nd − 1)
···
n
∼
= R[X1 , . . . , Xd ]/(X n1 − 1, . . . , X d − 1),
1

d

the i-th
P of which may be computed in O(bn lg ni ) bit operations. The overall cost

is O( i bn lg ni ) = O(bn lg n) bit operations.

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

9

3. Cyclotomic coefficient rings
The aim of this section is to construct certain coefficient rings that play a central
role in the multiplication algorithms described later. The basic idea is as follows.
Suppose that we want to multiply two polynomials in Fp [X], and that the degree
of the product is known to be at most n. If N is an integer with N > n, then by
appropriate zero-padding, we may embed the problem in Fp [X]/(X N − 1). Furthermore, if we have some factorisation N = αm, where α and m are relatively prime,
then there is an isomorphism
Fp [X]/(X N − 1) ∼
= Fp [Y, Z]/(Y α − 1, Z m − 1),
and the latter ring is closely related to

Fp [Y, Z]/(φα (Y ), Z m − 1) ∼
= (Fp [Y ]/φα )[Z]/(Z m − 1)

(recall that φα (Y ) denotes the α-th cyclotomic polynomial). In particular, computing the product in (Fp [Y ]/φα )[Z]/(Z m − 1) recovers ‘most’ of the information
about the product in Fp [X]/(X N − 1).
In this section we show how to choose N , α and m with the following properties:
(1) N is not much larger than n, so that not too much space is ‘wasted’ in the
initial zero-padding step;
(2) ϕ(α) (= deg φα ) is not much smaller than α, so that we do not lose much
information by working modulo φα (Y ) instead of modulo Y α − 1 (this
missing information must be recovered by other means);
(3) the coefficient ring Fp [Y ]/φα contains a principal m-th root of unity, so that
we can multiply in (Fp [Y ]/φα )[Z]/(Z m − 1) efficiently by means of DFTs
over Fp [Y ]/φα ;
(4) m is a product of many integers that are exponentially smaller than n, so
that the DFTs of length m may be reduced to many small DFTs; and
(5) α is itself exponentially smaller than n.
The last two items ensure that the small DFTs can be converted to multiplication
problems of degree exponentially smaller than n, to allow the recursion to proceed.
Definition 3.1. An admissible tuple is a sequence (q0 , q1 , . . . , qe ) of distinct primes
(e > 1) satisfying the following conditions. First,
2

(lg N )3 < qi < 2(lg lg N ) ,

i = 0, . . . , e,

(3.1)

where N := q0 · · · qe . Second, qi − 1 is squarefree for i = 1, . . . , e, and
2

λ(q0 , . . . , qe ) := LCM(q1 − 1, . . . , qe − 1) < 2(lg lg N ) .

(3.2)

(Note that q0 − 1 need not be squarefree, and q0 does not participate in (3.2).)
An admissible length is a positive integer N of the form N = q0 · · · qe where
(q0 , . . . , qe ) is an admissible tuple.
If N is an admissible length, we treat (q0 , . . . , qe ) and λ(N ) := λ(q0 , . . . , qe ) as
auxiliary data attached to N . For example, if an algorithm takes N as input, we
implicitly assume that this auxiliary data is also supplied as part of the input.
Example 3.2. For n = 10100000 , there is a nearby admissible length
N = 1000000000000000000156121 . . .(99971 digits omitted) . . . 26353
= q0 q1 · · · q6035

10

DAVID HARVEY AND JORIS VAN DER HOEVEN

where
q0 = 206658761261792645783,
q1 = 36658226833235899 = 1 + 2 ·3 ·11 ·17 ·23 ·29 ·37 ·53 ·59 ·67 ·71 ·89,

q2 = 36658244723486119 = 1 + 2 ·3 ·17 ·29 ·47 ·59 ·67 ·73 ·83 ·101 ·109,
q3 = 36658319675739343 = 1 + 2 ·3 ·7 ·17 ·29 ·31 ·41 ·47 ·53 ·61 ·89 ·103,

q4 = 36658428883190467 = 1 + 2 ·3 ·11 ·31 ·43 ·61 ·71 ·73 ·107 ·109 ·113,
···

and

q6035 = 37076481100386859 = 1 + 2 ·3 ·13 ·29 ·31 ·59 ·83 ·97 ·101 ·103 ·107

λ(N ) = 2 ·3 ·5 · · · 113 = 31610054640417607788145206291543662493274686990.
Definition 3.3. Let p be a prime. An admissible length N is called p-admissible
if N > p2 and p ∤ N (i.e., p is distinct from q0 , . . . , qe ).
The following result explains how to choose a p-admissible length close to any
prescribed target.
Proposition 3.4. There is an absolute constant z1 > 0 with the following prop2
erty. Given as input a prime p and an integer n > max(z1 , p2 ), in 2O((lg lg n) ) bit
operations we may compute a p-admissible length N in the interval


1
n.
(3.3)
n<N < 1+
lg n
The key ingredient in the proof is the following number-theoretic result of Adleman, Pomerance and Rumely.
Lemma 3.5 ([1, Prop. 10]). There is an absolute constant C1 > 0 with the following
property. For all x > 10, there exists a positive squarefree integer λ0 < x2 such that
X
1 > exp(C1 log x/ log log x).
q prime
q−1|λ0

2

2

Proof of Proposition 3.4. Let λmax := ⌈2 9 (lg lg n) ⌉, and for λ > 1 define f (λ) to be
the number of primes q in the interval (lg n)4 < q 6 λmax + 1 such that q − 1 | λ
and q 6= p. We claim that, provided n is large enough, there exists some squarefree
λ0 ∈ {1, . . . , λmax } such that f (λ0 ) > lg n. To see this, apply Lemma 3.5 with
2
1
x := 2 9 (lg lg n) ; for large n we then have
C1 log x/ log log x > 15(log2 x)1/2 = 5 lg lg n,
so Lemma 3.5 implies that there exists a positive squarefree integer λ0 < x2 6 λmax
for which
X
1 > exp(5 lg lg n) > (lg n)5
q prime
q−1|λ0

and hence
f (λ0 ) =

X

(lg n)4 <q6λmax +1
q prime, q 6= p
q−1|λ0

1 > (lg n)5 − (lg n)4 − 1 > lg n.

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

11

We may locate one such λ0 by means of the following algorithm (adapted from
the proof of [18, Lemma 4.5]). First use a sieve to enumerate the primes q in
the interval (lg n)4 < q 6 λmax + 1, and to determine which λ = 1, . . . , λmax are
squarefree, in (λmax )1+o(1) bit operations. Now initialise an array of integers cλ := 0
for λ = 1, . . . , λmax . For each q 6= p, scan through the array, incrementing those cλ
for which λ is squarefree and divisible by q − 1, and stop as soon as one of the cλ
reaches lg n. We need only allocate O(lg lg n) bits per array entry, so each pass
through the array costs O(λmax lg lg n) bit operations. The number of passes is
2
O(λmax ), so the total cost of finding a suitable λ0 is O(λ2max lg lg n) = 2O((lg lg n) )
bit operations. Within the same time bound, we may also easily recover a list of
primes q1 , q2 , . . . , qlg n for which qi − 1 | λ0 .
Next, compute the partial products q1 , q1 q2 , . . . , q1 q2 · · · qlg n , and determine the
2
1
smallest integer e > 1 for which q1 · · · qe > n/2 2 (lg lg n) . Such an e certainly exists,
as q1 · · · qlg n > 2lg n > n. Since each qi occupies O((lg lg n)2 ) bits, this can all be
done in (lg n)O(1) bit operations. Also, as
2

2

1

qe 6 λ0 + 1 6 2 9 (lg lg n) + 1 < 2 4 (lg lg n)

2

2

1

and q1 · · · qe−1 6 n/2 2 (lg lg n) , we find that
1

2

2 4 (lg lg n) <

2
n
1
< 2 2 (lg lg n)
q1 · · · qe

for large n.
Let q0 be the least prime that exceeds n/(q1 · · · qe ) and that is distinct from p.
According to [4], the interval [x − x0.525 , x] contains at least one prime for all
sufficiently large x; therefore
0.6

n
n
q0 <
+
q1 · · · qe
q1 · · · qe




2
1
n
1
n
< 1+
< 1 + (2 4 (lg lg n) )−0.4
q1 · · · qe
lg n q1 · · · qe
2

for n sufficiently large. We may find q0 in 2O((lg lg n) ) bit operations, by using trial
division to test successive integers for primality.
Set N := q0 q1 · · · qe . Then (3.3) holds, and certainly N > p2 and p ∤ N . Let us
check that (q0 , . . . , qe ) is admissible, provided n is large enough. For i = 1, . . . , e
we have
2

1

2

(lg N )3 < (lg n)4 < qi 6 λ0 + 1 < 2 4 (lg lg n) < 2(lg lg N ) ,
and also
3

(lg N ) < 2

2
1
4 (lg lg n)



2
2
1
1
< q0 < 1 +
2 2 (lg lg n) < 2(lg lg N ) ;
lg n
1

2

this establishes (3.1). Also, as q0 > 2 4 (lg lg n) > qi for i = 1, . . . , e, we see that q0
is distinct from q1 , . . . , qe . Finally, (3.2) holds because
2

2

2

LCM(q1 − 1, . . . , qe − 1) | λ0 6 2 9 (lg lg n) < 2(lg lg N ) .
2

This also shows that we may compute the auxiliary data λ(q0 , . . . , qe ) in 2O((lg lg n) )
bit operations.


12

DAVID HARVEY AND JORIS VAN DER HOEVEN

Remark 3.6. Example 3.2 was constructed by enumerating the smallest primes
q1 , q2 , . . . exceeding (lg n)3 for which qi − 1 | 2 · 3 · 5 · · · 113, halting just before
their product reached n, and then choosing q0 to make N as close to n as possible.
The proof of Proposition 3.4 goes a different way: rather than choosing λ first,
the proof constructs q1 , . . . , qe and λ simultaneously. In particular, one cannot
guarantee that λ will be a product of an initial segment of primes, as occurred
in the example. Indeed, the proof of [1, Prop. 10] (and of its predecessor [21])
yields very little information at all about the prime factorisation of λ. For further
discussion, see [1, Remark 6.2].
Definition 3.7. Let p be a prime and let N = q0 · · · qe be a p-admissible length.
A p-admissible divisor of N is a positive divisor α of N , with q0 | α, such that the
ring Fp [Y ]/φα (Y ) contains a principal (q1 · · · qe )-th root of unity, and such that
lg N < α < 2(lg lg N )
and


ϕ(α) > 1 −

1
lg N



4

α.

(3.4)
(3.5)

The next result shows how to construct a p-admissible divisor for any sufficiently
large p-admissible length N . The idea behind the construction is as follows. Let
ordn p denote the order of p in the multiplicative group of integers modulo n. For
any α > 1, not divisible by p, the ring Fp [Y ]/φα (Y ) is a direct sum of fields of
order pr , where r = ordα p [27, Lemma 14.50]. Our goal is to ensure that pr − 1
is divisible by q1 · · · qe , so that Fp [Y ]/φα (Y ) contains the desired principal root of
unity. One way to force qi to divide pr − 1 is simply to choose α divisible by qi ,
as this implies that ordqi p | r. The difficulty is that we cannot do this for all qi ,
because then α would become too large, violating (3.4). Fortunately, we can take
advantage of the fact that the qi − 1 share a small common multiple λ = λ(N ); this
enables us to take α to be a product of a small subset of the qi , in such a way that
still every one of q1 , . . . , qe divides pr − 1.
Proposition 3.8. There is an absolute constant z2 > 0 with the following property.
Given as input a prime p and a p-admissible length N > z2 , we may compute a
p-admissible divisor α of N , together with the cyclotomic polynomial φα ∈ Fp [Y ]
4
and a principal (q1 · · · qe )-th root of unity in Fp [Y ]/φα , in 2O((lg lg N ) ) p1+o(1) bit
operations.
Proof. We are given as input an admissible tuple (q0 , . . . , qe ) with N = q0 · · · qe ,
and the squarefree integer λ := λ(q0 , . . . , qe ). Let L be the set of primes dividing λ.
By (3.2) we have |L| 6 log2 λ < (lg lg N )2 , and we may compute L in λO(1) =
2
2O((lg lg N ) ) bit operations.
We start by computing a table of values of ordqi p for i = 1, . . . , e; note that
p 6= qi by hypothesis, so ordqi p is well-defined. We have qi − 1 | λ and hence
ordqi p | λ for each i. To compute ordqi p, we first compute p mod qi in O(lg qi lg p)
bit operations, and then repeatedly multiply by p modulo qi until reaching 1. Since
ordqi p 6 λ, and there are e = O(lg N ) primes qi , the total cost to compute the
table is
2
O((λ lg 2 qi + lg qi lg p) lg N ) = (2(lg lg N ) lg p)O(1)
bit operations.

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

13

Using the above table, we construct a certain vector σ = (σ1 , . . . , σe ) ∈ {0, 1}e
as follows. Initialise the vector as σ := (0, . . . , 0). For each ℓ ∈ L, search for the
smallest i = 1, . . . , e such that ℓ | ordqi p. If such an i is found, set σi := 1; if no i
is found, ignore this ℓ. The cost of computing σ is O(|L|e(lg λ)2 ) = (lg N )O(1) bit
operations.
Q
Set α := q0 i:σi =1 qi . To establish (3.4), note that the number of i for which
σi = 1 is at most |L|, so (3.1) implies that
2

2

2

4

lg N < q0 6 α < (2(lg lg N ) )|L|+1 6 (2(lg lg N ) )(lg lg N ) = 2(lg lg N ) .
For (3.5), first observe that
ϕ(α)
=
α

(lg lg N )2
 Y 
 

1
1
1
1−
> 1−
.
1−
q0 i:σ =1
qi
(lg N )3
i

Since − log(1 − ε) < 2ε for any ε ∈ (0, 21 ), we obtain
− log

1
2(lg lg N )2
ϕ(α)
<
<
α
(lg N )3
lg N

and hence ϕ(α)/α > exp(−1/ lg N ) > 1 − 1/ lg N for sufficiently large N .
Now compute the cyclotomic polynomial φα ∈ Fp [Y ] (i.e., the reduction modulo p of φα (Y ) ∈ Z[Y ]). This can be done in (α lg p)O(1) bit operations, using for
example [27, Algorithm 14.48]. We may then determine the factorisation of φα
into irreducibles in Fp [Y ], say φα = f1 · · · fk , in αO(1) p1/2+o(1) bit operations [24,
Thm. 1]. Since p ∤ α, the fj are distinct, and each fj has degree r := ordα p [27,
Lemma 14.50]. In other words, Fp [Y ]/φα is isomorphic to a direct sum of k copies
of Fpr .
We claim that qh | pr − 1 for all h = 1, . . . , e. For this, it suffices to prove that
ordqh p | r for each h. Since λ is squarefree, it suffices in turn to show that every
prime ℓ dividing ordqh p also divides r. But for every such ℓ, the procedure for
constructing σ must have succeeded in finding some i for which ℓ | ordqi p (since at
least one value of i works, namely i = h). Then σi = 1 for this i, so qi | α. This
implies that ordqi p | ordα p = r, and hence that ℓ | r.
We conclude that q1 · · · qe | pr − 1, so each Fp [Y ]/fj contains a primitive root
of unity of order q1 · · · qe . As the factorisation of q1 · · · qe is known, we may locate
one such primitive root in each Fp [Y ]/fj in αO(1) p1+o(1) bit operations [25] (see
also [17, Lemma 3.3]). Combining these primitive roots via the Chinese remainder
theorem, we obtain the desired principal (q1 · · · qe )-th root of unity in Fp [Y ]/φα in
another (α lg p)O(1) bit operations.

Remark 3.9. The p1+o(1) term in Proposition 3.8 arises from the best known deterministic complexity bounds for factoring polynomials and finding primitive roots.
If we permit randomised algorithms, then p1+o(1) may be replaced by (lg p)O(1) .
This has no effect on the main results of this paper.
Example 3.10. Continuing with Example 3.2, let us take p = 3. In the notation
of the proof of Proposition 3.8, we have L = {2, 3, 5, . . . , 113}. For each ℓ ∈ L, let

14

DAVID HARVEY AND JORIS VAN DER HOEVEN

us write q (ℓ) for the smallest qi for which ordqi 3 is divisible by ℓ. Then we have
q (2) = q1 ,

q (3) = q1 ,

q (5) = q5 ,

q (7) = q9 ,

q (11) = q1 ,

q (13) = q5 ,

q (17) = q1 ,

q (19) = q5 ,

q (23) = q1 ,

q (29) = q1 ,

q (31) = q3 ,

q (37) = q1 ,

q (41) = q3 ,

q (43) = q4 ,

q (47) = q2 ,

q (53) = q1 ,

q (59) = q1 ,

q (61) = q3 ,

q (67) = q1 ,

q (71) = q1 ,

q (73) = q2 ,

q (79) = q6 ,

q (83) = q2 ,

q (89) = q1 ,

q (97) = q5 ,

q (101) = q2 ,

q (103) = q3 ,

q (107) = q4 ,

q (109) = q2 ,

q (113) = q4 .

Therefore σi = 1 for i = 1, 2, 3, 4, 5, 6, 9, and we have
α = q0 q1 q2 q3 q4 q5 q6 q9
≈ 1.8385309928916569681 × 10136 ,

ϕ(α) ≈ 1.8385309928916566171 × 10136 ,

r = ordα 3 = 2 · 3 · 5 · · · 109 · 113 · 883 · 9041 · 327251 · 39551747.

The ring F3 [Y ]/φα is isomorphic to a direct sum of ϕ(α)/r copies of F3r . The
extraneous factors in r (namely 883, 9041, 327251 and 39551747) arise from the
auxiliary prime q0 . Let
m := N/α = q7 q8 q10 q11 · · · q6035 ≈ 5.439125 × 1099863 ;
then since m | q1 · · · qe | 3r − 1, each copy of F3r contains a primitive m-th root of
unity, so F3 [Y ]/φα contains a principal m-th root of unity. Thus it is possible to
multiply in the ring F3 [Y, Z]/(φα (Y ), Z m − 1) by using DFTs over F3 [Y ]/φα .
Remark 3.11. In Example 3.10, every ℓ ∈ L divides ordqi p for some i. It seems
likely that this always occurs (at least for large n), but we do not know how to prove
this. If it fails for some ℓ, then r may turn out not to be divisible by ℓ, but the
proof of Proposition 3.8 shows that we still have qh | pr − 1 for every h = 1, . . . , e.
4. Faster polynomial multiplication
The goal of this section is to prove Theorem 1.2. We assume that KZ > 1 is a
constant for which we have available an integer multiplication algorithm achieving
∗
M(n) = O(n lg n KZlog n ).
We will describe a recursive routine PolynomialMultiply, that takes as input
integers r, t > 1, a prime p, and polynomials U1 , . . . , Ut , V ∈ Fp [X]/(X r − 1), and
computes the products U1 V, . . . , Ut V . Its running time is denoted by Cpoly (t, r, p).
Note that the input polynomials U1 , . . . , Ut , V are expected to be supplied consecutively on the input tape (first U1 , then U2 , and so on), and the outputs U1 V, . . . , Ut V
should also be written consecutively to the output tape.
The role of the parameter t is to allow us to amortise the cost of transforming
the fixed operand V across t products. This optimisation (borrowed from [17] and
[18]) saves a constant factor in time at each recursion level of the main algorithm.
Altogether the algorithm will perform 2t + 1 transforms: t + 1 forward transforms
for U1 , . . . , Ut and V , followed by t inverse transforms to recover the products
U1 V, . . . , Ut V .

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

15

To simplify the analysis, it is convenient to introduce the normalisation
C⋆poly (r, p) := sup
t>1

Cpoly (t, r, p)
.
(2t + 1)r lg p lg(r lg p)

We certainly have Mp (n) < Cpoly (1, 2n, p) + O(n lg p), so to prove Theorem 1.2 it
is enough to show that
C⋆poly (r, p) = O(4max(0,log

∗

r−log∗ p)

∗

KZlog

p

).

(4.1)

The algorithms presented in this section perform many auxiliary multiplications
and divisions involving ‘small’ integers and polynomials. We assume that all auxiliary divisions are reduced to multiplication via Newton’s method [27, Ch. 9], so
that the cost of a division (by a monic divisor) is at most a constant multiple of the
cost of a multiplication of the same bit size. We also assume that, unless otherwise
specified, all auxiliary multiplications are handled using the integer and polynomial variants of the Schönhage–Strassen algorithm, whose complexities are given
by (1.1) and (1.4).
We first discuss a subroutine Transform that handles DFTs over rings of the
form Rp,α := Fp [Y ]/φα (Y ), where p is a prime and α > 1. It takes as input p and α,
positive integers t and n such that n is odd and relatively prime to α, a principal
n-th root of unity ω ∈ Rp,α , and t input sequences (as,0 , . . . , as,n−1 ) ∈ Rnp,α for
s = 1, . . . , t. Its output is the sequence of transforms (âs,0 , . . . , âs,n−1 ) ∈ Rnp,α with
respect to ω, for s = 1, . . . , t. Just like PolynomialMultiply, the input and
output sequences are stored consecutively on the tape.
Let T(t, n, α, p) denote the running time of Transform. The following result
shows how to reduce the DFT problem to an instance of PolynomialMultiply.
Proposition 4.1. We have
T(t, n, α, p) < Cpoly (t, nα, p) + O(tnα lg α lg lg α lg p lg lg p lg lg lg p).
Proof. Let R := Rp,α . We use Bluestein’s method to reduce each DFT to the
problem of computing a certain product fs (Z)g(Z) in R[Z]/(Z n − 1), plus O(n)
multiplications in R, where fs (Z) and g(Z) are defined as in Section 2.2. By (1.1)
and (1.4), each multiplication in R costs
O((α lg α lg lg α)(lg p lg lg p lg lg lg p))

(4.2)

bit operations. To handle the products fs (Z)g(Z), we first lift the polynomials from
Fp [Y, Z]/(φα (Y ), Z n −1) to Fp [Y, Z]/(Y α −1, Z n −1) (for example, by zero-padding
in Y up to degree α). We then compute their images under the isomorphism
Fp [Y, Z]/(Y α − 1, Z n − 1) ∼
= Fp [X]/(X nα − 1)
provided by Lemma 2.2; this costs altogether O(tnα lg α lg p) bit operations. We
call PolynomialMultiply to compute the products in Fp [X]/(X nα − 1), at a cost
of Cpoly (t, nα, p) bit operations. We evaluate the inverse of the above isomorphism
to bring the products back to Fp [Y, Z]/(Y α − 1, Z n − 1). Finally, we reduce modulo
φα (Y ) to obtain the desired products in R[Z]/(Z n − 1); the cost of each of these
divisions is given by (4.2).

We now return to multiplication in Fp [X]/(X r − 1). Our implementation of
PolynomialMultiply chooses one of two algorithms, depending on the size of r

16

DAVID HARVEY AND JORIS VAN DER HOEVEN

relative to p. For r 6 p2 it uses the straightforward Kronecker substitution method
described in Section 1. By (1.3) this yields the bound
log∗ (r lg p)

Cpoly (t, r, p) = O(t M(r lg p)) = O(tr lg p lg(r lg p)KZ

)

and hence
log∗ (p2 lg p)

C⋆poly (r, p) = O(KZ

) = O(KZlog

∗

p

),

r 6 p2 .

(4.3)

Therefore (4.1) holds in this case.
For r > p2 , most of the work will be delegated to a subroutine AdmissibleMultiply, which is defined as follows. It takes as input an integer t > 1, a prime p,
a p-admissible length N , and polynomials U1 , . . . , Ut , V ∈ Fp [X]/(X N − 1), and
computes the products U1 V, . . . , Ut V . In other words, it has the same interface as
as PolynomialMultiply, but it only works for p-admissible lengths. We denote
its running time by Cad (t, N, p). As above we also define the normalisation
C⋆ad (N, p) := sup
t>1

Cad (t, N, p)
.
(2t + 1)N lg p lg(N lg p)

The reduction from PolynomialMultiply to AdmissibleMultiply in the case
r > p2 is given by the following proposition.
Proposition 4.2. There is an absolute constant z3 > 0 with the following property. For any prime p and any integer r > max(z3 , p2 ), there exists a p-admissible
length N in the interval


1
2r
(4.4)
2r < N < 1 +
lg r
such that


O(1)
⋆
C⋆ad (N, p) + O(1).
(4.5)
Cpoly (r, p) < 2 +
lg r
Proof. Given as input U1 , . . . , Ut , V ∈ Fp [X]/(X r − 1), our goal is to compute the
products U1 V, . . . , Ut V . For sufficiently large r we may apply Proposition 3.4 with
n := 2r to find a p-admissible length N such that (4.4) holds. Since N > 2r, we
may simply zero-pad to reduce each problem to multiplication in Fp [X]/(X N − 1).
This yields
2

Cpoly (t, r, p) < Cad (t, N, p) + O(tr lg p) + 2O((lg lg r) ) ,
where the tr lg p term arises from the reduction modulo X r − 1, and the last term
from Proposition 3.4. Dividing by (2t + 1)r lg p lg(r lg p) and taking suprema over
t > 1, we find that
C⋆poly (r, p) <

N lg(N lg p) ⋆
C (N, p) + O(1).
r lg(r lg p) ad

Finally, since lg(N lg p) 6 lg(r lg p) + 2 we obtain



N lg(N lg p)
2
O(1)
1
1+
<2+
<2 1+
.
r lg(r lg p)
lg r
lg(r lg p)
lg r



The motivation for defining admissible lengths is the following result, which
shows how to implement AdmissibleMultiply in terms of a large collection of
exponentially smaller instances of PolynomialMultiply.

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

17

Proposition 4.3. There is an absolute constant z4 > 0 with the following property.
Let p be a prime and let N > z4 be a p-admissible length. Then there exist integers
r1 , . . . , rd in the interval
6

7

2(lg lg N ) < ri < 2(lg lg N ) ,
P
and weights γ1 , . . . , γd > 0 with i γi = 1, such that
C⋆ad (N, p) <

 d

O(1) X
γi C⋆poly (ri , p) + O(1).
2+
lg lg N i=1

(4.6)

(4.7)

Proof. We are given as input a prime p, a p-admissible length N = q0 · · · qe and
polynomials U1 , . . . , Ut , V ∈ Fp [X]/(X N − 1). Our goal is to compute the products
U1 V, . . . , Ut V . We will describe a series of reductions that converts this problem
to a collection of exponentially smaller multiplication problems, plus overhead of
O(tN lg N lg p) bit operations incurred during the reductions.
Step 1 — reduce to products over cyclotomic coefficient ring. Invoking Proposition 3.8, we compute a p-admissible divisor α of N , the cyclotomic polynomial
φα ∈ Fp [Y ], and a principal (q1 · · · qe )-th root of unity ω ∈ Fp [Y ]/φα . As p2 < N ,
4
this requires at most 2O((lg lg N ) ) p1+o(1) < N 1/2+o(1) bit operations.
Set ψα := (Y α − 1)/φα ∈ Fp [Y ]. Since Y α − 1 has no repeated factors in
Fp [Y ], we have (φα , ψα ) = 1. Using the Euclidean algorithm, compute polynomials
χ1 , χ2 ∈ Fp [Y ] of degree at most α such that χ1 φα + χ2 ψα = 1; this costs at most
(α lg p)O(1) < N o(1) bit operations.
Let m := N/α. As m and α are coprime, Lemma 2.2 provides an isomorphism
Fp [X]/(X N − 1) ∼
= Fp [Y, Z]/(Y α − 1, Z m − 1)
that may be evaluated in either direction in O(mα lg α lg p) bit operations. By (3.4)
this simplifies to O(N (lg lg N )4 lg p) = O(N lg N lg p) bit operations. Next, since
(φα , ψα ) = 1, there is an isomorphism
Fp [Y ]/(Y α − 1) ∼
= (Fp [Y ]/φα ) ⊕ (Fp [Y ]/ψα ).
Using the precomputed polynomials χ1 and χ2 , we may evaluate the above isomorphism in either direction in
O((α lg α lg lg α)(lg p lg lg p lg lg lg p)) = O(α(lg lg N )5 (lg lg lg N )2 lg p)
= O(α lg N lg p)
bit operations (here we have again used (3.4) and the fact that p2 < N ). This
isomorphism induces another isomorphism
Fp [Y ]/(Y α − 1, Z m − 1) ∼
= (Fp [Y ]/φα )[Z]/(Z m − 1) ⊕ (Fp [Y ]/ψα )[Z]/(Z m − 1)
by acting on the coefficient of each Z i separately; it may be evaluated in either
direction in O(mα lg N lg p) = O(N lg N lg p) bit operations. Chaining these isomorphisms together, we obtain an isomorphism
Fp [X]/(X N − 1) ∼
= (Fp [Y ]/φα )[Z]/(Z m − 1) ⊕ (Fp [Y ]/ψα )[Z]/(Z m − 1)
that may be evaluated in either direction in O(N lg N lg p) bit operations.

18

DAVID HARVEY AND JORIS VAN DER HOEVEN

We now use the following algorithm. First, at a cost of O(tN lg N lg p) bit
operations, apply the above isomorphism to U1 , . . . , Ut and V to obtain polynomials
U1′ , . . . , Ut′ , V ′ ∈ (Fp [Y ]/φα )[Z]/(Z m − 1),

Ũ1′ , . . . , Ũt′ , Ṽ ′ ∈ (Fp [Y ]/ψα )[Z]/(Z m − 1).

Second, compute the products Ũ1′ Ṽ ′ , . . . , Ũt′ Ṽ ′ : since deg ψα < α/ lg N by (3.5),
each of these products may be converted, via Kronecker substitution, to a product
of univariate polynomials in Fp [X] of degree O(mα/ lg N ) = O(N/ lg N ) (i.e., map
Y to X and Z to X 2 deg ψα ). The cost of these multiplications is
O(t((N/ lg N ) lg N lg lg N )(lg p lg lg p lg lg lg p)) = O(tN lg N lg p)
bit operations. Third, compute the products U1′ V ′ , . . . , Ut′ V ′ , using the method
explained in Step 2 below. Finally, at a cost of O(tN lg N lg p) bit operations, apply
the inverse isomorphism to the pairs (Us′ V ′ , Ũs′ Ṽ ′ ) to obtain the desired products
U1 V, . . . , Ut V .
Step 2 — convert to multidimensional convolutions. Let R := Fp [Y ]/φα . In this
step our goal is to compute the products U1′ V ′ , . . . , Ut′ V ′ , where U1′ , . . . , Ut′ , V ′ ∈
R[Z]/(Z m − 1). We do this by converting each problem to a multidimensional
convolution of size md × · · · × m1 , for a suitable decomposition m = m1 · · · md . For
the subsequent complexity analysis, it is important that the mi are chosen to be
somewhat larger than the coefficient size. To achieve this we proceed as follows.
Let m = ℓ1 · · · ℓu be the prime factorisation of m. The ℓj form a subset of
{q1 , . . . , qe }, so by (3.1) we have
(lg N )3 < ℓj < 2(lg lg N )
⌊ 25 (lg lg N )5 ⌋.

for each j. Let w :=
as (4.8) and (3.4) imply that
u>

2

(4.8)

We certainly have u > w for large enough N ,

log2 N − log2 α
log2 N − (lg lg N )4
log2 m
=
>
≫ (lg lg N )5 .
2
2
(lg lg N )
(lg lg N )
(lg lg N )2

Therefore we may take
m1 := ℓ1 · · · ℓw ,

m2 := ℓw+1 · · · ℓ2w ,
···

md−1 := ℓ(d−2)w+1 · · · ℓ(d−1)w ,

md := ℓ(d−1)w+1 · · · ℓdw ℓdw+1 · · · ℓu ,

where d := ⌊u/w⌋ > 1. Each mi is a product of exactly w primes, except possibly md , which is a product of at least w and at most 2w − 1 primes. For large N
we have
2
7
4
mi < (2(lg lg N ) )2w 6 2 5 (lg lg N )
(4.9)
and
mi > ((lg N )3 )w > (2lg lg N −1 )3w > 2(lg lg N )

6

(4.10)

for all i, and hence
d6

log2 m
lg N
6
.
(lg lg N )6
(lg lg N )6

(4.11)

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

19

Computing the decomposition m = m1 · · · md requires no more than (lg N )O(1) bit
operations.
As the mi are pairwise relatively prime, Corollary 2.3 furnishes an isomorphism
m
R[Z]/(Z m − 1) ∼
= R[Z1 , . . . , Zd ]/(Z m1 − 1, . . . , Z d − 1)
1

d

that may be computed in either direction in O((m lg m)(α lg p)) = O(N lg N lg p)
bit operations. Therefore we may use the following algorithm. First, at a cost of
O(tN lg N lg p) bit operations, compute the images
U1′′ , . . . , Ut′′ , V ′′ ∈ R[Z1 , . . . , Zd ]/(Z1m1 − 1, . . . , Zdmd − 1)

of U1′ , . . . , Ut′ , V ′ under the above isomorphism. Next, as explained in Step 3 below,
compute the products U1′′ V ′′ , . . . , Ut′′ V ′′ . Finally, apply the inverse isomorphism to
recover the products U1′ V ′ , . . . , Ut′ V ′ ; again this costs O(tN lg N lg p) bit operations.
Step 3 — reduce to DFTs over R. In this step our goal is to compute the products U1′′ V ′′ , . . . , Ut′′ V ′′ , where U1′′ , . . . , Ut′′ and V ′′ are as above. Let ωi := ω q1 ···qe /mi
for i = 1, . . . , d, where ω is the principal (q1 · · · qe )-th root of unity in R computed
in Step 1. According to the discussion in Section 2.2, the desired multidimensional
convolutions may be computed by performing t+1 multidimensional m-point DFTs
with respect to the evaluation points (ω1j1 , . . . , ωdjd ), followed by tm pointwise multiplications in R, and then t multidimensional m-point inverse DFTs and tm divisions
by m. The total cost of the pointwise multiplications and divisions is
O(tm(α lg α lg lg α)(lg p lg lg p lg lg lg p)) = O(tN lg N lg p)
bit operations.
Each of the 2t + 1 multidimensional DFTs may be converted to a collection
of one-dimensional DFTs of lengths m1 , . . . , md by the method explained in Section 2.2. Note that the inputs must be rearranged so that the data to transform
along each dimension may be accessed sequentially. Let 1 6 i 6 d, and consider the
transforms of length mi . Treating each input vector as a sequence of mi+1 · · · md
arrays of size mi × (m1 · · · mi−1 ), we must transpose each array into an array of size
(m1 · · · mi−1 ) × mi , perform m/mi DFTs of length mi , and then transpose back to
the original ordering. The total cost of all these transpositions is
P
O(tmα lg p i lg mi ) = O(tN lg p lg m) = O(tN lg N lg p)
bit operations.
The one-dimensional DFTs over R are handled by the Transform subroutine.
Combining the contributions from Steps 1, 2 and 3 shows that
d

m
X
, mi , α, p + O(tN lg N lg p).
T
Cad (t, N, p) < (2t + 1)
mi
i=1

This concludes the description of the algorithm; it remains to establish the overall
complexity claim. First, Proposition 4.1 yields
d
X
i=1

T

m

mi



, mi , α, p <

d
X

Cpoly

i=1

m

mi

, mi α, p



+ O(dmα lg α lg lg α lg p lg lg p lg lg lg p).
By (4.11), the last term lies in
O(dN (lg lg N )5 (lg lg lg N )2 lg p) = O(N lg N lg p).

20

DAVID HARVEY AND JORIS VAN DER HOEVEN

Setting ri := mi α for i = 1, . . . , d, we obtain
Cad (t, N, p) < (2t + 1)

d
X

Cpoly

i=1

N
ri


, ri , p + O(tN lg N lg p).

Notice that (4.6) follows immediately from (4.9), (4.10) and (3.4) (for large N ).
For the normalised quantities, we have

d
X
Cpoly N
ri , ri , p
⋆
Cad (N, p) <
+ O(1)
N lg p lg(N lg p)
i=1
<

d 
X
2N
i=1

ri

+1

 r lg(r lg p)
i
i
C⋆ (ri , p) + O(1).
N lg(N lg p) poly

Now observe that
log2 mi + log2 α + lg lg p + O(1)
log2 mi + O((lg lg N )4 )
lg(ri lg p)
<
<
.
lg(N lg p)
log2 N
log2 m
P
Put γi := log2 mi / log2 m, so that i γi = 1. Then (4.10) implies that




lg(ri lg p)
O(1)
O((lg lg N )4 )
γi .
γi < 1 +
< 1+
lg(N lg p)
log2 mi
lg lg N
Moreover, from (4.6) we certainly have
 2N
r
ri
O(1)
i
+1
=2+
<2+
.
ri
N
N
lg lg N
The desired bound (4.7) follows immediately.



Combining Proposition 4.2 and Proposition 4.3, we obtain the following recurrence inequality for C⋆poly (r, p). (This is identical to Theorem 7.1 of [18], but with
the constant 8 replaced by 4.)
Proposition 4.4. There are absolute constants z5 , C2 , C3 > 0 and a logarithmically
slow function Φ : (z5 , ∞) → R with the following property. For any prime p and
any integer r > max(z5 , p2P
), there exist positive integers r1 , . . . , rd < Φ(r), and
weights γ1 , . . . , γd > 0 with i γi = 1, such that

X
d
C2
⋆
Cpoly (r, p) < 4 +
γi C⋆poly (ri , p) + C3 .
(4.12)
lg lg r i=1
Proof. We first apply Proposition 4.2 to construct a p-admissible length N such
that (4.4) and (4.5) both hold; then we apply Proposition 4.3 to construct integers r1 , . . . , rd and weights γ1 , . . . , γd satisfying (4.6) and (4.7). Define Φ(x) :=
8
7
2(log log x) ; then certainly ri < 2(lg lg 3r) < Φ(r) for large r. The bound (4.12)
follows immediately by substituting (4.7) into (4.5).

Now we may prove our main result for multiplication in Fp [X]. The proof is very
similar to that of [18, Thm. 1.1].
∗

Proof of Theorem 1.2. We have already noted that C⋆poly (r, p) = O(KZlog p ) in the
region r 6 p2 (see (4.3)). To handle the case r > p2 , let z5 , C2 , C3 and Φ(x) be as

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

21

in Proposition 4.4. Increasing z5 if necessary, we may assume that z5 > exp(exp(1))
and that Φ(x) 6 x − 1 for all x > z5 . For each prime p, set σp := max(z5 , p2 ) and
Lp := max(C3 , max C⋆poly (r, p)) = O(KZlog
26r6σp

∗

p

).

Now apply Proposition 2.1 with K = 4, B = C2 /4, S = {1, 2, . . .}, ℓ = 2, κ = 1,
x0 = x1 = z5 , σ = σp , L = Lp , and T (r) = C⋆poly (r, p). The first part of the
recurrence for T (y) is satisfied due to the definition of Lp , and the second part due
∗
∗
to Proposition 4.4. We conclude that C⋆poly (r, p) = O(Lp 4log r−log σp ) for r > p2 .
Since log∗ σp = log∗ p + O(1) and Lp = O(KZlog
∗
∗
∗
C⋆poly (r, p) = O(4log r−log p KZlog p ) for r > p2 .

∗

p

), we obtain the desired bound


5. Faster integer multiplication
The goal of this section is to prove Theorem 1.1. We will describe a recursive
routine IntegerMultiply, that takes as input positive integers n and t, and integers u1 , . . . , ut , v ∈ Z/(2n − 1)Z, and computes the products u1 v, . . . , ut v. We
denote its running time by Cint (t, n). As in Section 4, it is convenient to define the
normalisation
Cint (t, n)
C⋆int (n) := sup
.
t>1 (2t + 1)n lg n
We certainly have M(n) < Cint (1, 2n) + O(n), so to prove Theorem 1.1 it is enough
to prove that
√
∗
(5.1)
C⋆int (n) = O((4 2)log n ).
We begin by revisiting the polynomial multiplication algorithm from Section 4.
Recall that to handle a multiplication problem in Fp [X]/(X r − 1) for r 6 p2 , we
used Kronecker substitution to convert it to an integer multiplication problem of
size O(r lg p) (see (4.3)). This approach is suboptimal because it ignores the cyclic
structure of Fp [X]/(X r − 1).
To exploit this structure, we introduce new routines RefinedPolynomialMultiply and RefinedAdmissibleMultiply. They have exactly the same interface as PolynomialMultiply and AdmissibleMultiply. Their running times
are denoted by C̃poly (t, r, p) and C̃ad (t, N, p), with corresponding normalisations
C̃⋆poly (r, p) and C̃⋆ad (N, p). The implementation of RefinedAdmissibleMultiply
is exactly the same as AdmissibleMultiply, except that calls to PolynomialMultiply are replaced by calls to RefinedPolynomialMultiply. Similarly, the
implementation of RefinedPolynomialMultiply for r > p2 is exactly the same
as PolynomialMultiply, except that calls to AdmissibleMultiply are replaced
by calls to RefinedAdmissibleMultiply. Therefore these routines satisfy analogues of Proposition 4.2 and Proposition 4.3, with C⋆poly and C⋆ad replaced by C̃⋆poly
and C̃⋆ad .
Where the new routines differ is in the implementation of RefinedPolynomialMultiply for the case r 6 p2 , which is described in the proof of the following
result. The idea is to exploit the cyclic structure by using IntegerMultiply to
handle the (cyclic) integer multiplication. This device saves a constant factor at
each recursion level of the main algorithm.
Proposition 5.1. For any prime p, and for any positive integer r satisfying
(lg lg p)2 < lg r < (lg p)1/2 ,

22

DAVID HARVEY AND JORIS VAN DER HOEVEN

there exists an integer n in the interval


1
2r lg p
2r lg p < n < 1 +
lg r

(5.2)

such that
C̃⋆poly (r, p) <



O(1)
C⋆int (n) + O(1).
2+
lg r

(Note that this bound does not hold over the whole range r 6 p2 ; to obtain the
constant 2, we need to restrict to a smaller range of r.)
Proof. We are given as input U1 , . . . , Ut , V ∈ Fp [X]/(X r − 1), and we wish to
compute the products U1 V, . . . , Ut V .
We use the following algorithm. Lift the inputs to polynomials U1′ , . . . , Ut′ , V ′ ∈
Z[X]/(X r − 1), whose coefficients lie in the interval 0 6 x < p. Evaluate these
polynomials at X = 2b , where b := 2 lg p + lg r; that is, pack the coefficients
together to obtain integers us := Us (2b ) and v := V (2b ) in Z/(2rb − 1)Z. Call
IntegerMultiply with n := rb to compute the cyclic integer products ws := us v.
Then we have ws = Ws′ (2b ) where Ws′ := Us′ V ′ ∈ Z[X]/(X r − 1). Observe that the
coefficients of Ws′ lie in the interval 0 6 x 6 r(p − 1)2 < rp2 − 1, and since 2b > rp2 ,
we may unpack ws to recover the coefficients of Ws′ unambiguously. Finally, by
reducing the coefficients of Ws′ modulo p, we arrive at the desired products Ws ∈
Fp [X]/(X r − 1).
As n = 2r lg p + r lg r, the bound (5.2) follows by taking into account the hypothesis that lg r < (lg p)1/2 . For the complexity we have
C̃poly (t, r, p) < Cint (t, n) + O(tr lg p lg lg p lg lg lg p),
where the last term covers the divisions by p at the end of the algorithm (and also
the linear-time packing and unpacking steps). Dividing by (2t + 1)r lg p lg(r lg p)
and taking suprema over t > 1, we obtain


n lg n
lg lg p lg lg lg p
⋆
⋆
C̃poly (r, p) <
.
C (n) + O
r lg p lg(r lg p) int
lg(r lg p)
The last term lies in O(1) thanks to the assumption lg r > (lg lg p)2 . Moreover,
(5.2) implies that lg n 6 lg(r lg p) + 2, so we find that



1
2
O(1)
n lg n
<2 1+
.

1+
<2+
r lg p lg(r lg p)
lg r
lg(r lg p)
lg r
Now we describe the implementation of IntegerMultiply. It chooses one of
two algorithms, depending on the size of n. For small n, it calls any convenient
basecase multiplication algorithm, such as the Schönhage–Strassen algorithm. For
large n, it uses the algorithm described in the proof of Proposition 5.2 below.
This algorithm reduces the problem to a collection of instances of RefinedAdmissibleMultiply, one for each prime p ∈ P(n), where P(n) is defined to be
the set consisting of the smallest lg n primes that exceed 12 (lg n)2 (we will see in
the proof below that these primes satisfy lg p = 2 lg lg n + O(1)). For example,
P(105 ) = {149, 151, . . . , 233} (the first 14 primes after 144.5).

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

23

Proposition 5.2. There is an absolute constant z7 > 0 with the following property.
For all n > z7 , there exists an admissible length N in the interval


n
3
n
<N < 1+
,
(5.3)
lg n lg lg n
lg lg n lg n lg lg n
such that N is p-admissible for all p ∈ P(n), and such that

 X
O(1)
1 ⋆
⋆
Cint (n) < 2 +
C̃ (N, p) + O(1).
lg lg n
lg n ad

(5.4)

p∈P(n)

Proof. We are given as input u1 , . . . , ut , v ∈ Z/(2n − 1)Z, and we wish to compute
the products u1 v, . . . , ut v.
Step 1 — choose parameters. In this preliminary step we compute a number of
parameters that depend only on n.
The prime number theorem (see for example [3, p. 9]) implies that the number
of primes between 21 (lg n)2 and (lg n)2 is asymptotically
(lg n)2
≫ lg n.
2 log((lg n)2 )
Therefore, for large n we certainly have
2
1
2 (lg n)

for all p ∈ P(n). Define P :=

Q

p∈P(n)

< p < (lg n)2

(5.5)

p; then

(2 log2 lg n − 1) lg n 6 log2 P 6 (2 log2 lg n) lg n,
so
2 lg n lg lg n − 3 lg n 6 lg P 6 2 lg n lg lg n.
Clearly we may compute P(n) and P within (lg n)O(1) bit operations.
Let


2n
n′ :=
;
lg P − lg n − 3

(5.6)

(5.7)

this makes sense for large n, as lg P ≫ lg n. Using Proposition 3.4 (with p = 2),
construct an admissible length N in the interval


1
n′ < N < 1 +
n′ .
lg n′
′ 2

The invocation of Proposition 3.4 costs 2O((lg lg n ) ) = O(n) bit operations. Let us
check that (5.3) holds for this choice of N . In one direction, by (5.6) we have
N > n′ >

2n
2n
n
>
>
.
lg P − lg n − 3
lg P
lg n lg lg n

(5.8)

24

DAVID HARVEY AND JORIS VAN DER HOEVEN

For the other direction we have



2n
1
+
1
N < 1+
lg n′
2 lg n lg lg n − 4 lg n − 3



1
4 lg n + 3
n
lg n lg lg n
= 1+
1
+
+
′
lg n
2 lg n lg lg n − 4 lg n − 3
n
lg n lg lg n



2 + o(1)
n
1
1+
< 1+
lg n′
lg lg n
lg n lg lg n


3
n
< 1+
lg lg n lg n lg lg n
for large n.
Finally, let us verify that N is p-admissible for all p ∈ P(n) (for large n). First,
by (5.8) and (5.5) we have N > (lg n)4 > p2 . Also, by (3.1) and (5.8), every prime
divisor of N = q0 · · · qe satisfies qi > (lg N )3 > (lg n)2 > p; in particular, p ∤ N .
Step 2 — convert to polynomial product modulo P . In this step we apply
the Crandall–Fagin algorithm from Section 2.3. For this, we require that lg P >
2⌈n/N ⌉ + lg N + 1; this follows from (5.7) as
lnm
2n
2n
lg P > ′ + lg n + 3 >
+ lg N + 1.
+ lg n + 3 > 2
n
N
N
We also require an element θ ∈ Z/P Z such that θN = 2. To construct θ, we first
compute ap := N −1 (mod p − 1) for each p ∈ P(n). This modular inverse exists
because N is a product of primes that are all greater than p, and hence relatively
prime to p − 1. Then we put θp := 2ap (mod p), so that (θp )N = 2 (mod p). Using
the Chinese remainder theorem, we compute θ ∈ Z/P Z such that θ = θp (mod p)
for all p ∈ P(n); then θN = 2 as desired. All of this can be effected in (lg n)O(1) bit
operations.
According to Section 2.3, the problem of computing u1 v, . . . , ut v reduces to computing U1 V, . . . , Ut V for certain polynomials U1 , . . . , Ut , V ∈ (Z/P Z)[X]/(X N − 1),
plus auxiliary operations amounting to
O(t(N lg P + N (lg n)2 + N lg P lg lg P lg lg lg P )) = O(tn lg n/ lg lg n)
bit operations.
Step 3 — reduce to products modulo small primes. In this step we convert each
multiplication problem in (Z/P Z)[X]/(X N − 1) into a collection of products in
Fp [X]/(X N − 1), for p ∈ P(n).
We start with the isomorphism Z/P Z ∼
= ⊕p∈P(n) Fp , which may be computed in
either direction in O(lg P (lg lg P )2 lg lg lg P ) bit operations using fast simultaneous
modular reduction and fast Chinese remaindering algorithms [27, §10.3]. It induces
an isomorphism
M
Fp [X]/(X N − 1),
(Z/P Z)[X]/(X N − 1) ∼
=
p∈P(n)

which may be computed in either direction, for all s = 1, . . . , t, in
O(tN lg P (lg lg P )2 lg lg lg P ) = O(tn(lg lg n)2 lg lg lg n)
bit operations. Note that the isomorphism Z/P Z ∼
= ⊕p Fp must be applied to
the coefficient of each X i independently, but the subroutine for multiplying in

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

25

Fp [X]/(X N − 1) needs sequential access to all of the residues for a single prime p.
The required data rearrangement corresponds to transposing a tN × |P(n)| array,
which costs only O(tN |P(n)| lg |P(n)| maxp lg p) = O(tn lg lg n) bit operations. Finally, for each p ∈ P(n), the products in Fp [X]/(X N − 1) may be computed by
calling RefinedAdmissibleMultiply, since N is a p-admissible length.
Combining the contributions from Steps 1, 2 and 3, we obtain
X
C̃ad (t, N, p) + O(tn lg n/ lg lg n).
Cint (t, n) <
p∈P(n)

Dividing by (2t + 1)n lg n and taking suprema over t > 1 yields
C⋆int (n) <

X

C̃⋆ad (N, p)

p∈P(n)

N lg p lg(N lg p)
+ O(1).
n lg n

By (5.5) we have lg p 6 2 lg lg n, so (5.3) implies that


O(1)
n
N lg p < 2 +
,
lg lg n lg n
and also lg(N lg p) 6 lg n, for large n. The bound (5.4) follows immediately.



We may now glue together the various pieces to obtain a doubly-exponential
recurrence for C⋆int (n).
Proposition 5.3. There are absolute constants z9 > z8 > 0 and C5 , C6 > 0,
and a logarithmically slow function Ψ : (z8 , ∞) → R such that Ψ(z9 ) > z8 , with
the following property. For any n > z9 , there
P exist positive integers n1 , . . . , nd <
Ψ(Ψ(n)), and weights γ1 , . . . , γd > 0 with i γi = 1, such that
C⋆int (n) <


32 +

C5
lg lg lg n

X
d

γi C⋆int (ni ) + C6 .

(5.9)

i=1

Proof. Step 1 — top-level call to IntegerMultiply. Applying Proposition 5.2,
we obtain an admissible N in the interval


n
n
3
<N < 1+
,
(5.10)
lg n lg lg n
lg lg n lg n lg lg n
such that N is p-admissible for all p ∈ P(n), and such that
 X

1 ⋆
O(1)
C̃ (N, p) + O(1).
C⋆int (n) < 2 +
lg lg n
lg n ad

(5.11)

p∈P(n)

In what follows, we frequently use the estimates lg N = lg n + O(lg lg n) and
lg lg N = lg lg n + O(1), which follow from (5.10). Also, from (5.5) we have lg p =
2 lg lg n + O(1) for all p ∈ P(n).
Step 2 — first call to RefinedAdmissibleMultiply. In this step we use (the
refined analogue of) Proposition 4.3 to estimate the C̃⋆ad (N, p) term in (5.11), for a
fixed p ∈ P(n). We obtain integers rp,1 , . . . , rp,dp such that
6

2(lg lg N ) < rp,i < 2(lg lg N )

7

(5.12)

26

DAVID HARVEY AND JORIS VAN DER HOEVEN

P
and weights γp,1 , . . . , γp,dp > 0 with i γp,i = 1, such that
 dp

O(1) X
γp,i C̃⋆poly (rp,i , p) + O(1).
C̃⋆ad (N, p) < 2 +
lg lg N i=1
Substituting into (5.11), and using lg lg N = lg lg n + O(1), yields
 X X

dp
γp,i ⋆
O(1)
C̃
(rp,i , p) + O(1).
C⋆int (n) < 4 +
lg lg n
lg
n poly
i=1

(5.13)

p∈P(n)

Step 3 — first call to RefinedPolynomialMultiply. In this step we use (the
refined analogue of) Proposition 4.2 to estimate the C̃⋆poly (rp,i , p) term in (5.13), for
a fixed p ∈ P(n) and i ∈ {1, . . . , dp }. The precondition rp,i > max(z3 , p2 ) holds for
large n, because by (5.12) we have
lg(p2 ) = 4 lg lg n + O(1) < (lg lg N )6 6 lg rp,i .
Thus there exists a p-admissible length Np,i in the interval


1
2rp,i
2rp,i < Np,i < 1 +
lg lg n
such that


O(1)
⋆
C̃poly (rp,i , p) < 2 +
C̃⋆ad (Np,i , p) + O(1).
lg lg n
Substituting into (5.13) yields
 X X

dp
γp,i ⋆
O(1)
⋆
C̃ (Np,i , p) + O(1).
Cint (n) < 8 +
lg lg n
lg
n ad
i=1

(5.14)

(5.15)

p∈P(n)

Step 4 — second call to RefinedAdmissibleMultiply. In this step we use
Proposition 4.3 again, to estimate the C̃⋆ad (Np,i , p) term in (5.15), for a fixed p ∈
P(n) and i ∈ {1, . . . , dp }. We obtain integers rp,i,1 , . . . , rp,i,dp,i such that
6

7

2(lg lg Np,i ) < rp,i,j < 2(lg lg Np,i )
P
and weights γp,i,1 , . . . , γp,i,dp,i > 0 with j γp,i,j = 1, such that
C̃⋆ad (Np,i , p)


< 2+

O(1)
lg lg Np,i

X
dp,i

(5.16)

γp,i,j C̃⋆poly (rp,i,j , p) + O(1).

j=1

We have lg lg Np,i > lg lg rp,i > lg lg lg n, so substituting into (5.15) yields
 X X

dp dp,i
X γp,i γp,i,j
O(1)
C̃⋆poly (rp,i,j , p) + O(1). (5.17)
C⋆int (n) < 16 +
lg lg lg n
lg
n
i=1 j=1
p∈P(n)

Step 5 — second call to RefinedPolynomialMultiply. In this step we use
Proposition 5.1 to estimate the C̃⋆poly (rp,i,j , p) term in (5.17), for a fixed p ∈ P(n),
i ∈ {1, . . . , dp } and j ∈ {1, . . . , dp,i }. The precondition
(lg lg p)2 < lg rp,i,j < (lg p)1/2

holds for large n, as (5.16), (5.14) and (5.12) imply that
(6 lg lg lg n + O(1))6 < lg rp,i,j < (7 lg lg lg n + O(1))7 ,

FASTER INTEGER AND POLYNOMIAL MULTIPLICATION

27

whereas (lg lg p)2 = (lg lg lg n + O(1))2 and (lg p)1/2 = (2 lg lg n + O(1))1/2 . We thus
obtain an integer np,i,j in the interval


1
2rp,i,j lg p
(5.18)
2rp,i,j lg p < np,i,j < 1 +
lg lg lg n
such that
C̃⋆poly (rp,i,j , p) <


2+

O(1)
lg lg lg n



C⋆int (np,i,j ) + O(1).

Substituting into (5.17) produces

 X X
dp dp,i
X γp,i γp,i,j
O(1)
C⋆int (n) < 32 +
C⋆int (np,i,j ) + O(1).
lg lg lg n
lg
n
i=1 j=1
p∈P(n)

The weights γp,i γp,i,j / lg n sum to 1, so after appropriate reindexing we obtain the
desired bound (5.9).
8
Finally, for the logarithmically slow function Ψ(x) := 2(log log x) , let us verify
that for large n, we have np,i,j < Ψ(Ψ(n)) for all p, i and j. First, since lg lg Np,i >
lg lg lg n, we have lg p < 3 lg lg n 6 3 · 2lg lg Np,i , and hence, by (5.18) and (5.16),
np,i,j < 3rp,i,j lg p < 9 · 2(lg lg Np,i )

7

+lg lg Np,i

< Ψ(Np,i ).

Then by (5.14) and (5.12) we have

7

7

Np,i < 3rp,i < 3 · 2(lg lg N ) 6 3 · 2(lg lg n) < Ψ(n).

Since Ψ(x) is increasing, we get the desired inequality np,i,j < Ψ(Ψ(n)).



Now we may prove the main theorem for integer multiplication.
Proof of Theorem
√ 1.1.∗ We have already noted that it suffices to establish that
C⋆int (n) = O((4 2)log n ) (see (5.1)). Let z8 , z9 , C5 , C6 and Ψ(x) be as in Proposition 5.1. Increasing z8 if necessary, we may assume that z8 > exp(exp(exp(1)))
and that Ψ(x) 6 x − 1 for all x > z8 . Applying Proposition 2.1 with K = 32,
B = C5 /32, S = {1, 2, . . .}, ℓ = 3, κ = 2, x0 = z8 , x1 = σ = z9 , L =
max(C6 , max16n6z9 C⋆int (n)), and T (n) = C⋆int (n) leads immediately to the desired
bound.

Acknowledgments
The authors thank Grégoire Lecerf for his comments on a draft of this paper.
The first author was supported by the Australian Research Council (DP150101689
and FT160100219).
References
1. L. M. Adleman, C. Pomerance, and R. S. Rumely, On distinguishing prime numbers from
composite numbers, Ann. of Math. (2) 117 (1983), no. 1, 173–206. MR 683806 (84e:10008)
2. R. Agarwal and J. Cooley, New algorithms for digital convolution, IEEE Transactions on
Acoustics, Speech, and Signal Processing 25 (1977), no. 5, 392–410.
3. T. M. Apostol, Introduction to analytic number theory, Springer-Verlag, New York-Heidelberg,
1976, Undergraduate Texts in Mathematics. MR 0434929
4. R. C. Baker, G. Harman, and J. Pintz, The difference between consecutive primes. II, Proc.
London Math. Soc. (3) 83 (2001), no. 3, 532–562. MR 1851081
5. L. I. Bluestein, A linear filtering approach to the computation of discrete Fourier transform,
IEEE Transactions on Audio and Electroacoustics 18 (1970), no. 4, 451–455.

28

DAVID HARVEY AND JORIS VAN DER HOEVEN

6. A. Bostan, P. Gaudry, and É. Schost, Linear recurrences with polynomial coefficients and
application to integer factorization and Cartier-Manin operator, SIAM J. Comput. 36 (2007),
no. 6, 1777–1806. MR 2299425 (2008a:11156)
7. D. G. Cantor and E. Kaltofen, On fast multiplication of polynomials over arbitrary algebras,
Acta Inform. 28 (1991), no. 7, 693–701. MR 1129288 (92i:68068)
8. S. Covanov and E. Thomé, Fast integer multiplication using generalized Fermat primes,
http://arxiv.org/abs/1502.02800, 2016.
9. R. Crandall and B. Fagin, Discrete weighted transforms and large-integer arithmetic, Math.
Comp. 62 (1994), no. 205, 305–324. MR 1185244
10. A. De, P. Kurur, C. Saha, and R. Saptharishi, Fast integer multiplication using modular
arithmetic, SIAM Journal on Computing 42 (2013), no. 2, 685–699.
11. M. Fürer, Faster integer multiplication, STOC’07—Proceedings of the 39th Annual ACM
Symposium on Theory of Computing, ACM, New York, 2007, pp. 57–66. MR 2402428
(2009e:68124)
12.
, Faster integer multiplication, SIAM J. Comput. 39 (2009), no. 3, 979–1005.
MR 2538847 (2011b:68296)
13. I. J. Good, The interaction algorithm and practical Fourier analysis, J. Roy. Statist. Soc. Ser.
B 20 (1958), 361–372. MR 0102888 (21 #1674)
14. D. Harvey, Faster truncated integer multiplication, https://arxiv.org/abs/1703.00640,
2017.
15. D. Harvey and J. van der Hoeven, Faster integer multiplication using plain vanilla FFT
primes, https://arxiv.org/abs/1611.07144, to appear in Mathematics of Computation,
2016.
16. D. Harvey, J. van der Hoeven, and G. Lecerf, Faster polynomial multiplication over finite
fields, technical report, http://arxiv.org/abs/1407.3361, 2014.
, Even faster integer multiplication, J. Complexity 36 (2016), 1–30. MR 3530637
17.
, Faster polynomial multiplication over finite fields, J. ACM 63 (2017), no. 6, 52:1–
18.
52:23.
19. C. H. Papadimitriou, Computational complexity, Addison-Wesley Publishing Company, Reading, MA, 1994. MR 1251285 (95f:68082)
20. J. M. Pollard, The fast Fourier transform in a finite field, Math. Comp. 25 (1971), 365–374.
MR 0301966
21. K. Prachar, Über die Anzahl der Teiler einer natürlichen Zahl, welche die Form p − 1 haben,
Monatsh. Math. 59 (1955), 91–97. MR 0068569
22. A. Schönhage, Schnelle Multiplikation von Polynomen über Körpern der Charakteristik 2,
Acta Informat. 7 (1976/77), no. 4, 395–398. MR 0436663
23. A. Schönhage and V. Strassen, Schnelle Multiplikation grosser Zahlen, Computing (Arch.
Elektron. Rechnen) 7 (1971), 281–292. MR 0292344 (45 #1431)
24. V. Shoup, On the deterministic complexity of factoring polynomials over finite fields, Inform.
Process. Lett. 33 (1990), no. 5, 261–267. MR 1049276 (91f:11088)
25.
, Searching for primitive roots in finite fields, Math. Comp. 58 (1992), no. 197, 369–
380. MR 1106981
26. L. H. Thomas, Using computers to solve problems in physics, Applications of digital computers
458 (1963), 42–57.
27. J. von zur Gathen and J. Gerhard, Modern computer algebra, third ed., Cambridge University
Press, Cambridge, 2013. MR 3087522
E-mail address: d.harvey@unsw.edu.au
School of Mathematics and Statistics, University of New South Wales, Sydney NSW
2052, Australia
E-mail address: vdhoeven@lix.polytechnique.fr
Laboratoire d’informatique, UMR 7161 CNRS, École polytechnique, 91128 Palaiseau
Cedex, France

