Robin Larrieu

The Truncated Fourier Transform for mixed radices

The standard version of the Fast Fourier Transform (FFT) is applied to problems of size n = 2 k . For this reason, FFT-based evaluation/interpolation schemes often reduce a problem of size l to a problem of size n, where n is the smallest power of two with l n. However, this method presents "jumps" in the complexity at powers of two; and on the other hand, n -l values are computed that are actually unnecessary for the interpolation. To mitigate this problem, a truncated variant of the FFT was designed to avoid the computation of these unnecessary values. In the initial formulation [14], it is assumed that n is a power of two, but some use cases (for example in finite fields) may require more general values of n. This paper presents a generalization of the Truncated Fourier Transform (TFT) for arbitrary orders. This allows to benefit from the advantages of the TFT in the general case.

Introduction

Many basic arithmetic operations on polynomials can be performed efficiently using evaluation/interpolation techniques. A typical example is the multiplication of polynomials: let A, B be polynomials in K[X] for a field K, such that deg(AB) < n. The product AB can be computed by choosing n different values (x 1 , . . . , x n ) ∈ K n , and evaluating the A(x i ) and B(x i ) for all i. Then, a term-by-term multiplication leads to the values (AB)(x i ) for all i, and an interpolation on these values allows to retrieve the polynomial AB.

The Discrete Fourier Transform (DFT) is a way to perform these evaluations and interpolations on specifically chosen values: let ω be a primitive n-th root of unity and P ∈ K[X] of degree less than n; then the DFT of P is the n-tuple (P (1), P (ω), . . . , P (ω n-1 )). Symmetrically, the inverse DFT computes the coefficients of P from the values (P (ω i )) 0 i<n . When n is highly composite, the DFT and its inverse can be computed efficiently using the Fast Fourier Transform (FFT) algorithm. This method was known to Gauss around 1805 [START_REF] Gauss | Nachlass: Theoria interpolationis methodo nova tractata[END_REF], but it received little attention until after its rediscovery by Cooley and Tukey [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF]. For this reason, and because they were the first to use this method as a systematic computation algorithm, the modern formulation of the FFT is usually attributed to Cooley and Tukey.

Jump phenomenon

The standard version of FFT-based multiplication algorithms relies on a primitive n-th root of unity ω ∈ K, where n is a power of two. This requirement causes the following drawback: when a polynomial of degree less than d is considered (or when d evaluation points are needed) with d slightly larger than 2 k , one must perform a FFT of order 2 k+1 . This causes a significant overhead since up to twice as many values are computed as what is actually needed.

This jump phenomenon can be mitigated by allowing a more precise choice of n d. For example, instead of requiring n = 2 k , one can allow more general products of small primes such as n = 2 k 3 l 5 m . FFTs of such sizes reduce to DFTs of sizes 2, 3 and 5, for which efficient codelets are implemented e.g. in the FFTW3 library [START_REF] Frigo | The design and implementation of FFTW3[END_REF]. Alternatively, optimized radix-2 methods may be preferred because of their simplicity (fewer base cases to handle). For example, the FFT pruning [START_REF] Markel | FFT pruning[END_REF] aims to reduce the overhead for a zero-padded sequence. Another example is Crandall and Fagin's Discrete Weighted Transform [START_REF] Crandall | Discrete weighted transforms and large-integer arithmetic[END_REF], which reduces a problem of size d < 2 k to two problems of size 2 l and 2 m with d < 2 l + 2 m < 2 k .

Another elegant solution to this jump phenomenon is to use the Truncated Fourier Transform [START_REF] Van Der Hoeven | The truncated Fourier transform and applications[END_REF]. The TFT behaves as a usual FFT of order 2 k , but it performs a multipoint evaluation with exactly the desired length, while avoiding the computation of all intermediate values that are not needed for obtaining the output. Moreover, the interpolation can be performed with the same complexity using the inverse TFT. Improvements to this algorithm were made to reduce memory usage [START_REF] Harvey | An in-place truncated Fourier transform and applications to polynomial multiplication[END_REF], and improve cache friendliness [START_REF] Harvey | A cache-friendly truncated FFT[END_REF]. Mateer [START_REF] Mateer | Fast Fourier Transform Algorithms with Applications[END_REF]Chapter 6] also proposed a different formulation of the TFT inspired by Crandall and Fagin's reduction. He mentions that this alternative formulation can be used with a few adaptations when n = 3 k , or another prime power.

Goal of this paper

The methods discussed above rely on a choice of n with a very specific form. This requires the base field K to contain primitive n-th roots of unity for all such n. This is true for K = C, but in general, the choice of roots of unity is restricted. It is always possible to add virtual roots of unity (with certain restrictions on their order if the field has nonzero characteristic) as in the Schönhage-Strassen algorithm [START_REF] Schönhage | Fast multiplication of large numbers[END_REF], but this extension causes computational overhead.

Assume that the choice of roots of unity is restricted by both our base field K and practical considerations. Let S ⊂ N denote the set of orders n for the roots of unity that can be used in FFTs. For example, the use cases mentioned in previous section assume K = C, and only practical considerations are taken into account. This leads to sets of the form S = {2 k |k ∈ N}, or S = {2 k 3 l 5 m |k, l, m ∈ N}. In a finite field like F 2 l , there are roots of unity only for specific orders, so we would have S = {n|n divides 2 l -1}. A remarkable example is F 2 60 because there are efficient ways of computing in this field, and many roots of unity (with large, highly-composite order) are known [START_REF] Harvey | Fast polynomial multiplication over F 2 60[END_REF].

As discussed previously, there is a jump phenomenon at elements of S (more or less important depending on their distribution). This paper aims to reduce this jump phenomenon through a generalization of the Truncated Fourier Transform to an arbitrary n.

At first, we give a brief reminder of the Fast Fourier Transform in the general case, and some useful notations are introduced. In section 3 and 4, we present algorithms to compute the Truncated Fourier Transform and its inverse for any n. Finally, the complexity of the TFT is discussed and compared with the ordinary FFT.

2 The Cooley-Tukey FFT Let A = (A i ) 0 i<n be a vector in K n , and let ω be a primitive n-th root of unity. The Discrete Fourier Transform (DFT) of A with respect to ω is the vector  = ( Âi ) 0 i<n where Âi

= n-1 j=0 A j • ω ij . (2.1) Equivalently, if P A is the polynomial A 0 + A 1 X + • • • + A n-1 X n-1 , then Âi = P A (ω i ).
The Fast Fourier Transform (FFT) describes an algorithm to compute the DFT efficiently.

Fundamentals of the FFT algorithm

The Cooley-Tukey algorithm [START_REF] Cooley | An algorithm for the machine calculation of complex Fourier series[END_REF] relies on the following remark: assuming that n = n 1 n 2 is composite, the following holds for all k 1 < n 1 and

k 2 < n 2 Âk 1 +n 1 k 2 = n 2 -1 i=0 n 1 -1 j=0 A i+n 2 j • ω (i+n 2 j)(k 1 +n 1 k 2 ) ,
that is:

Âk 1 +n 1 k 2 = n 2 -1 i=0 ω ik 1 • n 1 -1 j=0 A i+n 2 j • (ω n 2 ) jk 1 • (ω n 1 ) ik 2 . (2.2)
Using formula (2.2), a Discrete Fourier Transform of length n = n 1 n 2 can be decomposed into n 2 DFTs of length n 1 (inner DFTs) followed by n 1 DFTs of length n 2 (outer DFTs). These smaller DFTs can be computed recursively using the same method as long as their size not prime. In the special case where n = 2 k , this leads to the well-known complexity bound O(n lg n).

Remark 2.1. Each inner DFT is applied to a subvector of the input (with offset i < n 1 and stride n 2 ). Instead, it can be preferable to reindex the working array before the inner DFTs, and once again before the outer DFTs. The main purpose of this reindexing is to perform the recursive calls on contiguous memory blocks. This approach may lead to more efficient implementation because of cache effects: sparing frequent data exchanges with the RAM during the inner DFT compensates for the cost of the reindexing.

Generalized bitwise mirror

In an in-place implementation of the Cooley-Tukey FFT, it is convenient to order the output differently (see for example [7, Section 2.1]). The purpose of this different order is to ensure that the result of the full FFT is simply the concatenation of the outputs from the outer DFTs. As a consequence, it spares a matrix transposition at the end of the algorithm, and therefore also at each recursion step. The v-mirror noted [i] v defined in this section is a notation for this new indexation: FFT(A) i = Â[i]v Let the vector v = (p 0 , . . . , p d-1 ) be a decomposition of n (n = d-1 j=0 p j ); with not necessarily prime p j 's. Then, any index i < n can be uniquely written under the form:

i = i 0 • p 1 • • • p d-1 + i 1 • p 2 • • • p d-1 + • • • + i d-2 • p d-1 + i d-1
, where for all j, i j < p j .

With these notations, the v-mirror of i is defined as

[i] v = i 0 + i 1 • p 0 + i 2 • p 0 p 1 + • • • + i d-1 • p 0 • • • p d-2 .
When all p j 's are 2, this definition coincides with the bitwise mirror introduced for the radix-2 TFT from [START_REF] Van Der Hoeven | The truncated Fourier transform and applications[END_REF]. Let h d, w 1 = (p 0 , . . . , p h-1 ), w 2 = (p h , . . . , p d-1 ), n 1 = h-1 j=0 p j and n 2 = d-1 j=h p j (so that n = n 1 n 2 ). It is easy to show the following basic properties of the v-mirror:

[[i] v ] v = i with v = (p d-1 , . . . , p 0 ), (2.3) [i] v = i if d = 1 i.e. v = (n), (2.4) [I + n 2 J] v = [J] w 1 + n 1 [I] w 2 for I < n 2 , J < n 1 .
(2.5)

Remark 2.2. As claimed at the beginning of this subsection, if FFT(A) i = Â[i]v , then the output of the top-level FFT is indeed the concatenation of the outputs from the outer FFTs. This comes from properties (2.5) and (2.2) combined, and it will be more formally discussed in the next subsection, once additional notations have been introduced.

Remark 2.3. As shown in [START_REF] Van Der Hoeven | The truncated Fourier transform and applications[END_REF] for n = 2 k , it is crucial for the functioning of the TFT algorithm to order the output according to the v-mirror. Actually, using the natural ordering would require a transposition after the outer FFTs (and also at each recursion depth). In the case of the TFT, such a transposition would mix the known and unknown values, while the TFT relies on a contiguous block of known values.

Further notations for the steps of the FFT algorithm

Now consider an execution of the FFT algorithm as in section 2.1. The recursive calls of the algorithm define a decomposition

n = p 0 • • • p d-1 of n, with n 1 = p 0 • • • p h-1 and n 2 = p h • • • p d-1 .
We reuse the notations from the previous section: v = (p 0 , . . . , p d-1 ), w 1 = (p 0 , . . . , p h-1 ) and w 2 = (p h , . . . , p d-1 ). As in [START_REF] Harvey | Fast polynomial multiplication over F 2 60[END_REF], we assume the output verifies FFT(A) i = Â[i]v (that is, the Discrete Fourier Transform that is returned is reordered according to the v-mirror).

Remark 2.4. Between the inner and outer DFTs, intermediates values are multiplied by a power of ω, more precisely ω i•k 1 as in formula (2.2). This factor is usually called the twiddle factor. With an indexation of the output according to the v-mirror, the twiddle factor must be ω

i[k 1 ]w 1
For clarity, we will introduce different input and output vectors. These vectors are actually just names to represent specific parts of the working vector at different steps of the algorithm, but no duplication of data should occur. Let A and B (size n) be the input and output vectors of the algorithm. Let α i and β i (i < n 2 , each vector has size n 1 ) be the input and output vectors of the i-th inner DFT (recursive call of the FFT algorithm). Let γ j and δ j (j < n 1 , each vector has size n 2 ) be the input and output vectors of the j-th outer DFT. These notations are illustrated in Figure 1. By definition, we have the following properties for all i, j:

(α i ) j = A i+n 2 j , (2.6) (β i ) j = FFT(α i ) j = ( α i ) [j]w 1 , (2.7)
(γ j ) i = ω i[j]w 1 • (β i ) j , (2.8) (δ j ) i = FFT(γ j ) i = ( γ j ) [i]w 2 ,
(2.9)

(δ j ) i = B i+n 2 j = Â[i+n 2 j]v . ( 2
. Let K 1 , k 1 < n 1 and K 2 , k 2 < N 2 such that k 1 = [K 1 ] w 1 and k 2 = [K 2 ] w 2 .
Then, using the property (2.5) of the v-mirror in relation (2.10), we get:

(α i ) j = A i+n 2 j , (β i ) K 1 = ( α i ) k 1 , (γ K 1 ) i = ω ik 1 • (β i ) K 1 (δ K 1 ) K 2 = ( γ K 1 ) k 2 and (δ K 1 ) K 2 = Âk 1 +n 1 k 2 .
This proves the correctness of the algorithm using formula (2.2).

Specification of the FFT algorithm and its inverse

This section formalizes FFT and inverse FFT as blackboxes to be used in TFT algorithms.

In particular, assumptions on the input and the output are precised.

Because of Remark 2.3, we assume the output is ordered according to the v-mirror. Hence, the FFT algorithm can be described as follows:

Algorithm 2.1. FFT • INPUT: an integer n ∈ N, a vector A = (a i ) i<n , a vector v = (p 0 , . . . , p d-1 )
such that n = d-1 i=0 p i , and a primitive n-th root of unity ω

• OUTPUT: the vector T = ( Â[i]v ) i<n
Symmetrically, the inverse transformation is computed by the IFFT (inverse FFT) algorithm: Algorithm 2.2. IFFT (inverse FFT)

• INPUT: an integer n ∈ N, a vector B = (b i ) i<n , a vector v = (p 0 , . . . , p d-1 )
such that n = d-1 i=0 p i , and a primitive n-th root of unity ω

• OUTPUT: the vector A = (a i ) i<n such that ∀i < n, b i = Â[i]v
Remark 2.6. The inverse FFT can be computed by reversing the FFT algorithm. Alternatively, using FFT as a blackbox, the following formula (and reordering of the input/output) can be used:

DFT ω -1 (DFT ω (A)) k = n-1 i=0 n-1 j=0 a i • ω (i-k)j = na k then, (DFT ω ) -1 = 1 n DFT ω -1 .
(2.11)

The Truncated Fourier Transform for arbitrary orders

This section generalizes the Truncated Fourier Transform (TFT) [START_REF] Van Der Hoeven | The truncated Fourier transform and applications[END_REF] for an arbitrary order n = p 0 p 1 • • • p d-1 . Given a vector A of length n, the TFT computes l n well chosen values of the Discrete Fourier Transform of A; that is the vector T = (T 0 , . . . , T l-1 ) = ( Âi 0 , . . . , Âi l-1 ). Note that these are not necessarily the first l values of Â. The algorithm presented in this section aims to perform less computation than by simply computing the DFT of A and discarding the unused values.

Atomic transforms

At first we consider the following base case: given (a 0 , . . . , a p-1 ) for p prime or reasonably small, we want to compute directly the TFT ( Â0 , . . . , Âl-1 ). To do so, one can naively apply Horner's rule for each value as in formula (2.1), which is especially efficient for small p (principle of a specialized codelet). For larger p, it becomes more interesting to compute the full DFT, then discard unused values. A full DFT of a such size can be computed using efficient transformations such as Rader's algorithm [START_REF] Rader | Discrete Fourier transforms when the number of data samples is prime[END_REF] and Bluestein's transform [START_REF] Bluestein | A linear filtering approach to the computation of discrete fourier transform[END_REF].

We assume that these considerations translate into the following algorithm:

Algorithm 3.1. atomicTFT

• INPUT: integers n ∈ N and l < n, a vector A = (a 0 , . . . , a n-1 ) and and a primitive n-th root of unity ω

• OUTPUT: the vector T = ( Â0 , . . . , Âl-1 )

Remark 3.1. In the base case, we do not use the indexation according to the v-mirror because of property (2.3).

General idea

Assume only l n values of the output are actually needed. The plain FFT algorithm can be modified to avoid computation of irrelevant intermediate values. As stated in Remark 2.3, the output of the DFT must be ordered according to the v-mirror.

If we want to return the tuple ( Â[i]v ) 0 i<l , then according to relation (2.10), the vectors δ j need to be computed only for j < m = l/n 2 . This means by definition (2.9) that only the γ j with j < m are needed. From formula (2.8), we conclude that for every i < n 2 , only the first m values of β i need to be computed.

Moreover, if q = (l quo n 2 ) < m, only the r = (l rem n 2 ) first values of δ q are needed (where quo and rem represent the quotient and remainder operations in the euclidean division). Figure 2 gives a visual representation of which values are actually needed.

Presentation of the algorithm

The previous discussion suggest that a TFT of order n = n 1 n 2 can be decomposed into n 2 TFTs of order n 1 followed by m TFTs of order n 2 (as for the usual FFT). If the top-level TFT has length l, then the inner TFTs have length m = l/n 2 . Most of the outer TFTs ) such that n = d-1 i=0 p i , and a primitive n-th root of unity ω

• OUTPUT: the vector T = ( Â[i]v ) i<l if d = 1 then return atomicTFT(n, l, A, ω) Algorithm 3.1 else if l = n then return FFT(n, A, v, ω)
Ordinary FFT (Algorithm 2.1) else choose h ∈ {1, . . . , d -1}; define n 1 , n 2 , w 1 , w 2 as in section 2.3 m ← l/n 2 ; q ← l quo n 2 ; r ← l rem n 2 for 0 i < n 2 do n 2 TFTs of size n 1 (α i ) j ← A i+n 2 j for all j < n 1

β i ← TFT(n 1 , m, α i , w 1 , ω n 2 ) (γ j ) i ← ω i[j]w 1 • (β i ) j for all j < m end for T ← ∅ for 0 j < m do n 1 TFTs of size n 2 if j < q then δ j ← FFT(n 2 , γ j , w 2 , ω n 1 ) Ordinary FFT (Algorithm 2.1) else δ j ← TFT(n 2 , r, γ j , w 2 , ω n 1 ) end if T ← (T, δ j ) end for return T end if Theorem 3.1.
Assuming correct implementations of Algorithms 2.1 and 3.1, the Algorithm 3.2 is correct.

Proof. Case l = n corresponds to a full FFT, then Algorithm 2.1 returns the expected result. If d = 1 (typically when n is prime), the Discrete Fourier Transform is computed directly. Then, unnecessary outputs are discarded. The results are ordered as expected because of property (2.4). In the other cases, the algorithm is called recursively, and its correctness results by induction from the discussion in section 3.2. Remark 3.2. As for the in-place Cooley-Tukey FFT presented in [START_REF] Harvey | Fast polynomial multiplication over F 2 60[END_REF], the result depends on the vector v, but not on the choice of h.

The inverse TFT for arbitrary orders

The formula (2.11) for the usual FFT cannot be used in the case of the TFT because not all values of the transform are known. As for the standard TFT (n = 2 k ) [START_REF] Van Der Hoeven | The truncated Fourier transform and applications[END_REF], we revert the algorithm computing the TFT instead.

The inversion of a TFT is to be understood as follows: assume that the values (B i ) 0 i<l of the output, and (A i ) l i<n of the input are known. Then, the goal is to retrieve the missing values (A i ) 0 i<l of the input. Typically, the values A i (for i l) are known to be 0 because of a simple analysis regarding the degree, but the coefficients of highest degree of a polynomial can also be deduced from a limit analysis.

At first, we provide a method to solve the base case of size p by direct computation (here p is not necessarily prime, but it should be reasonably small, so the problem can be solved without further decomposition). Then, we present a recursive algorithm that reduces the TFT inversion to such a base case.

Atomic inverse transforms

In this section, we consider the following skew butterfly problem: given Â0 , . . . , Âl-1 and a l , . . . , a p-1 , how can we compute the missing values Âl , . . . , Âp-1 and a 0 , . . . , a l-1 ?

The Discrete Fourier Transform is given by the following matrix-vector product:

     Â0 Â1 . . . Âp-1      =      1 1 • • • 1 1 ω • • • ω p-1 . . . . . . . . . . . . 1 ω p-1 • • • ω (p-1)(p-1)           a 0 a 1 . . . a p-1      , or, in a more compact form, Â = V ω • A
For any m p, we define the submatrices

V ω,m = (ω ij ) 0 i,j<m , Ṽω,m = (ω (i+m)(j+m ) 0 i,j<p-m , and
W ω,m = (ω i(j+m) ) i<m ; j<p-m .
Note that V ω,m has size m × m, Ṽω,m has size (p -m) × (p -m) and W ω,m has size m × (p -m). In other words:

V ω = V ω,m W ω,m W ω,m Ṽω,m (4.1)
The considered skew butterfly problem is equivalent to the resolution of the following matrix equation with parameters A 2 , B 1 and unknowns A 1 and B 2 :

B 1 B 2 = V ω,l W ω,l W ω,l Ṽω,l • A 1 A 2 (4.2)
V ω,l is a general Vandermonde matrix of determinant 0 i<j<l (ω i -ω j ) = 0, hence it is invertible. Therefore,

A 1 = (V ω,l ) -1 (B 1 -W ω,l • A 2 ) . ( 4

.3)

Once A 1 is known, it is easy to compute B 2 as Remark 4.2. Often, it is not necessary to compute B 2 entirely, but only specific values.

B 2 = W ω,l • A 1 + Ṽω,l • A 2 . ( 4
In this case, equation (4.4) reduces to a much smaller computation. For example, when composing atomic transforms for the inverse TFT, returning the first value of B 2 is actually sufficient. For our usage, we assume that the results from this section translate into the following algorithm:

= (A 1 , A 2 ) = (a i ) i<n

Recursive algorithm

In a similar way as in the case n = 2 k [14], intermediate results are not always computed in an order corresponding to the recursion depth. In a usual FFT, all outer FFTs are inverted, then the inner FFTs are inverted. On the contrary for the TFT, some of the inner TFTs must be inverted, and these inversions provide additional values that allow the outer TFT to be inverted. This behavior is illustrated in Figure 3). Therefore, the algorithm will return some of the missing output values (b i ) l i<n in addition to the desired input values (a i ) 0 i<l . It turns out that the outer TFT needs only one value from each inner TFT before it can be inverted, so returning b l is actually sufficient for the functioning of the recursive algorithm.

We assume that we have at our disposal an algorithm (Algorithm 2.2) to reverse a full FFT (case l = n). We also assume having another algorithm to reverse atomic TFTs (Algorithm 4.1). The design from previous paragraph translates into the following algorithm: Algorithm 4.2. ITFT (inverse TFT)

• INPUT: integers n ∈ N and l < n, vectors A 2 = (a i ) l i<n and B 1 = (b i ) i<l , a vector v = (p 1 , . . . , p d ) such that n = d i=1 p i , and a primitive n-th root of unity ω

• OUTPUT: the vector A 1 = (a i ) i<l and the value b l such that:

∀i l, b i = Â[i]v where A = (A 1 , A 2 ) = (a i ) i<n if l = 0 then return ∅, n-1 i=0 a i else if d = 1 then return atomicITFT(n, l, A 2 , B 1 , ω) Algorithm 4.1 else choose h ∈ {1, . . . , d -1}; define n 1 , n 2 , w 1 , w 2 as in section 2.3 q ← l quo n 2 , r ← l rem n 2 for j < q do "Step 1" (δ j ) i ← (B 1 ) i+n 2 j for all i < n 1 γ j ← IFFT(n 2 , δ j , w 2 , ω n 1 )
Algorithm 2.2 (β i,1 ) j ← ω -i[j]w 1 • (γ j ) i for all i < n 1 end for for r i < n 2 do "Step 2" (α i,2 ) j ← a i+n 2 j for all j q i + n 2 j l (α i,1 ), (

β i ) q ← ITFT(n 1 , q, α i,2 , β i,1 , w 1 , ω n 2 ) q < n 1 since l < n (γ q ) i ← ω i[q]w 1 • (β i ) q (A 1 ) i+n 2 j ← (α i,1 ) j for all j < q end for ∆ 1 ← ((B 1 ) i+n 2 q ) i<r "Step 3" Γ 2 ← ((γ q ) i ) r i<n 2 Γ 1 , (δ q ) r ← ITFT(n 2 , r, ∆ 1 , Γ 2 , w 2 , ω n 1 ) for i < r do "Step 4" (β i,1 ) q ← ω -i[q]w 1 • (Γ 1 ) i if q + 1 < n 1 then (α i,2 ) j ← a i+n 2 j for all j > q i + n 2 j l (α i,1 ), dummy ← ITFT(n 1 , q + 1, α i,2 , β i,1 , w 1 , ω n 2 ) else (α i,1 ) ← IFFT(n 1 , β i,1 , w 1 , ω n 2 ) Algorithm 2.2 end if (A 1 ) i+n 2 j ← (α i,1
) j for all j q end for return A 1 , (δ q ) r end if Theorem 4.1. Assuming correct implementations of Algorithms 2.2 and 4.1, the Algorithm 4.2 is correct.

Proof. We proceed by induction over d. The case l = 0 is clearly correct. For d = 1, the result is computed directly using Algorithm 4.1, which is supposed to be correct. As for Algorithm 3.2, the result is ordered as expected because of property (2.4) of the v-mirror.

In Step 1, the γ j are computed for all j < q using a full reverse FFT. This means γ j and δ j verify equation (2.9) for j < q. Then, the first part of every vector β i is computed according to equation (2.8).

At this point, the vectors α i and β i are partially known. More precisely, we know the values with j > q of α i and the values j < q of β i . Moreover, (α i ) q is known for i r (by definition, l = n 2 q + r). In Step 2, a recursive call computes the missing part of α i for these i, as well as the value (β i ) q . This means α i and ((β i ) j ) j q verify equation (2.7) for all i r (the recursive call is correct by the induction hypothesis).

At the end of Step 2, the second part of (γ q ) i (for i r) is computed according to equation (2.8). It is noted Γ 2 in the algorithm. The first part ∆ 1 of (δ q ) i (for i < r) is also known from the input. In Step 3, the missing part Γ 1 of γ q (that is the (γ q ) i for i < r) is computed as well as (δ q ) r through a recursive call. Then, γ q and ((δ q ) i ) i r verify equation (2.9).

Finally in Step 4, the (β i ) q are computed for i < r. The (α i ) j being given on input for j > q, and the (β i ) j being known from Step 1, the missing (α i ) j (j q) can be computed using a recursive call. Since this call is correct by the induction hypothesis, α i and ((β i ) j ) j q verify equation (2.7) for i < r (hence for all i because of Step 2.

All in all, the vectors (some truncated) α i , ((β i ) j ) j q , γ j (j q), δ j (j < q) and ((δ q ) i ) i r verify equations (2.6) to (2.10). This is sufficient to prove correctness as seen in section 3 (where l is replaced by l = l + 1)

Practical remarks

Remark 4.3. Algorithm 4.2 can be used to compute the unique polynomial P of degree less than n such that l evaluation points are given by the vector B 1 and the coefficients of degree at least l are given by A 2 : ∀i < l, P (ω [i]v ) = b i and ∀i l, P i = a i . In particular, it can be used to interpolate a polynomial of degree less than l by setting A 2 = (0, . . . , 0). Remark 4.4. In section 4.1, the order p may be composite (which can happen if an element of the vector v from Algorithm 4.2 is composite), but the resolution of equation (4.3) is not efficient if p is large. Algorithm 4.2 shows that the problem can be reduced to smaller sizes as long as its order is composite. However, it seems difficult to solve a skew butterfly problem if its size is a large prime.

For example, in F 2 60 , we have primitive roots of unity of order 2 60

-1 = 3 2 • 5 2 • 7 • 11 • 13 • 31 • 41 • 61 • 151 • 331 • 1321.
It is feasible to perform the inversion using linear algebra for size up to 13 efficiently, but direct computation may become too costly for p = 331 or 1321 for example. Remark 4.5. An inverse TFT of length 0 is actually very simple: it reduces to the computation of b 0 . For this reason, if (l rem n 2 ) = 0, then the recursive call in Algorithm 4.2 at Step 3 becomes trivial. This partially solves the problem mentioned in the previous remark: assume the prime factors of n are sorted in the vector v (in increasing order). If for example an inversion by direct computation is possible for p 0 , . . . , p k-1 but not for p k , . . . , p d-1 (because these primes are too large), then a TFT of length l can still be reverted

if ñ = p k × • • • × p d-1 divides l.
Remark 4.6. It is important that the recursive calls in Step 2 return the additional output value (β i ) q . However, it is not necessary that Algorithm 4.2 always returns the additional output value b l . For example, this value is simply discarded in the recursive calls from Step 4. Another typical case where this value is not needed is for the interpolation a polynomial of degree less than l from the l values ( Â[i]v ) i<l . It is possible to adapt Algorithm 4.2 to avoid this unnecessary computation when the value b l is not needed.

A remarkable duality for atomic inverse transforms

Direct resolution of the skew butterfly problem from section 4.1 requires the inversion of the matrix V ω,l of size l × l, which becomes expensive if l is large. In this section, we present a dual problem that can be solved through the inversion of matrix V ω -1 ,p-l , that has size (p -l) × (p -l). This duality ensures that the skew butterfly problem can always be solved through the inversion of a matrix V φ,m (and a few matrix-vector products), where m p/2 and φ is either ω or ω -1 .

Property (2.11) gives the inverse matrix

(V ω ) -1 = (1/p) • V ω -1 .
Then, with a decomposition of V ω -1 as in formula (4.1), equation (4.2) can be rewritten as follows:

A 1 A 2 = 1 p • V ω -1 ,l W ω -1 ,l W ω -1 ,l Ṽω -1 ,l B 1 B 2 .
We want to solve the equation above with parameters A 2 , B 1 and unknowns A 1 , B 2 . We have

B 2 = ( Ṽω -1 ,l ) -1 • (pA 2 -W ω -1 ,l B 1 ) . (4.5)
Once B 2 is known, it is easy to compute A 1 since

A 1 = 1 p (V ω -1 ,l B 1 + W ω -1 ,l B 2 ) . (4.6)
Except for the factor 1/p, these formulas are symmetric to the formulas (4.3) and (4.4). This shows the duality between these two problems.

The initial claim was that the problem could be solved through the inversion of V ω -1 ,p-l , but equation (4.5) involves the matrix Ṽω -1 ,l . Actually, it appears that Ṽω -1 ,l deduces from V ω -1 ,p-l by multiplication with diagonal matrices. To reduce notations, let φ = ω -1 . We introduce the diagonal matrices D φ,l = Diag (φ l(l+i) ) i<p-l and Dφ,l = Diag (φ li ) i<p-l of size (p -l) × (p -l). We have the following property:

( Ṽφ,l ) i,j = φ (l+i)(l+j) = φ l(l+j) • φ ij • φ li for all i, j , hence, Ṽφ,l = D φ,l • V φ,p-l • Dφ,l . (4.7)
Remark 4.7. The matrix V ω -1 deduces from V ω by a simple row (or column) permutation. For this reason, the dual problem presented here is equivalent (in terms of field operations) to the direct problem from section 4.1 up to scalar-vector multiplications and multiplication by diagonal matrices.

Remark 4.8. Assume that we have an efficient method to compute the function X → V φ • X. For large p, this can be done using Rader's [START_REF] Rader | Discrete Fourier transforms when the number of data samples is prime[END_REF] or Bluestein's [START_REF] Bluestein | A linear filtering approach to the computation of discrete fourier transform[END_REF] algorithm. Then, the required matrix-vector products (for example W ω -1 ,l • B 1 ) can be computed efficiently using the relation

dummy W ω -1 ,l • B 1 = V ω -1 • B 1 0 .

Complexity analysis

This section aims to evaluate the field operation count for a TFT of size n and length l n, and to compare it with the cost for a full FFT of size n. In the following, we note these costs T(l, n) and F(n) respectively. Asymptotic bounds involve the field operation count for common arithmetic operations on polynomials of degree n: M(n) for the multiplication and C(n) for the cyclic convolution (multiplication modulo X n -1).

Let n = p 0 • • • p d-1 be the size of the Discrete Fourier Transform (FFT or TFT) that is considered. If we develop completely the recursive calls of the FFT algorithm as in section 2.1, the execution decomposes into d successive transformations of a vector of length n (see Figure 4). At each row, the working vector is transformed using n/p i independent DFTs of size p i , which are computed directly.

Complexity of a full FFT

In a full FFT, all n/p i atomic DFTs are computed at each row. If we note f(p) = F(p)/p the normalized operation count (per intermediate value), then each of the atomic DFTs has an operation count of p i f(p i ). Summing for all i, we get the following result: Theorem 5.1 (Complexity of the usual Cooley-Tukey FFT). We have

F(n) = n • (f(p 0 ) + f(p 1 ) + • • • + f(p d-1 )) + O(nd) .
(5.1)

Remark 5.1. The term O(nd) is the cost for the operations between rows of atomic DFTs, that is the multiplications by twiddle factors. The multiplicative constant is actually small; in fact, with a precomputed table of twiddle factors, this represents n(d -1) multiplications, and even less if we take into account that some of the twiddle factors are equal to 1.

Remark 5.2. For small p, it is most efficient to compute the atomic DFTs using specialized codelets that perform naive matrix-vector products. This yields F(p) ∼ p 2 , or f(p) ∼ p. For larger p, methods like Rader's or Bluestein's algorithms are more efficient. In this case, F(p) ∼ p log p, or f(p) ∼ log p (with a larger constant factor than for the naive method).

Complexity of atomic TFTs

As discussed in section 3.1, there are two simple methods to compute an atomic TFT:

• for very small l, one can naively compute the l first values of the Fourier transform using Horner's rule (p additions and p multiplications for each evaluation point). This method has a cost of T(l, p) = 2lp.

• for l near p, it is interesting to compute the full atomic DFT and discard the last p -l values, at a cost of T(l, p) = F(p).

Efficient methods for intermediate values of l are more complex. As an example, we sketch briefly here a method that is adapted for p = O(l). This method is based on the following result from Bostan and Schost's work [START_REF] Bostan | Polynomial evaluation and interpolation on special sets of points[END_REF]:

Lemma 5.2. Let (1, ω, . . . , ω l-1 ) be a geometric progression such that the points ω i are pairwise distinct. Then, a polynomial of degree less than l can be evaluated on these points in M(l) + O(l) field operations.

In our case, the points ω i are indeed pairwise distinct because l p where p is the order of the root ω. Then, to evaluate the polynomial P of degree p, we write P (X) = P 0 (X) + X l P 1 (X) + • • • + X ql P q (X) with deg P i < l for all i, and q = p quo l .

Then, the multipoint evaluation of P can be done by evaluating (p quo l) polynomials of degree less than l, followed by operations on vectors of size l (addition and multiplication by diagonal matrices) This method to compute atomic TFTs has an operation cost of

T(l, p) = p l M(l) + O(p) .
As in the previous subsection, we normalize the operation count per intermediate value: t(l, p) = T(l, p)/l. We introduce the overhead for the atomic TFT as k(l, p) = t(l, p)/f(p). Since k(l, p) decreases with l, it is meaningful to also introduce

K(p) = sup l p l • k(l, p) .
By definition, we have: In particular, l k(l, p) p, which gives K(p) p.

t(l, p) = k(l, p)f(p) , T(l, p) = l • t(l, p) K(p)f(p) . ( 5
Remark 5.4. Actually, the case l = p corresponds to a full atomic DFT, so T(p, p) = F(p) = p f(p). Combined with the upper bound from the previous remark, we get K(p) = p. However, we keep the notation K(p) to respect the symmetry with the inverse TFT.

Complexity of atomic inverse TFTs

Similarly, we introduce the corresponding costs for the inverse TFT: let T * (l, p) be the cost for the inverse TFT, and t * (l, p) = T * (l, p)/l is the normalized cost. The corresponding overhead is k * (l, p) = t * (l, p)/f(p) and K * (p) = sup l p l • k * (l, p). We know from sections 4.1 and 4.4 that the resolution of the skew butterfly problem requires a few matrix-vector products and the inversion of a matrix V φ,m , where m p/2 and φ is either ω or ω -1 . Without loss of generality, we assume that φ = ω and m = l p/2 and we consider the direct method from section 4.1. (If l > p/2, then we reduce to the dual problem, which causes only O(p) additional operations because of Remark 4.7.) Actually, it is not necessary to compute (V ω,l ) -1 ; it is sufficient to compute the function Y → (V ω,l ) -1 Y , which is a polynomial interpolation on the points 1, ω, . . . , ω l-1 . The following result from Bostan and Schost [2, Section 5] gives an upper bound for the cost of this operation: Lemma 5.3. Let (1, ω, . . . , ω l-1 ) be a geometric progression such that the points ω i are pairwise distinct. Then, the interpolation of a polynomial of degree less than l on these points can be performed in 2M(l) + O(l) operations.

Since l p/2, a multiplication of size l can be seen as a cyclic convolution of size p, that is M(l) C(p). Moreover, Bluestein's transform is an efficient method to compute the DFT for a large p, which gives F(p) = C(p) + O(p) asymptotically. We are now able to bound the cost of atomic inverse TFTs: Proof. At first, we compute Y = B 1 -W ω,l A 2 . The transformation X → V ω X is a DFT that can be performed in C(p) + O(p) operations (using for example Bluestein's transform). Then, using Remark 4.5, the matrix-vector product W ω,l A 2 can be computed using C(p) + O(p). Then, we compute A 1 = (V ω,l ) -1 Y , which can be done using 2M(l) + O(l) base field operations from Lemma 5.3; and we have M(l) C(p) because l p/2.

Finally, the relation B = V ω A allows to retrieve B 2 in C(p) + O(p) operations.

Complexity of the TFT

When computing row i + 1 from row i in the TFT, some of the atomic DFTs are not performed, most of the remaining are full DFTs and only a few are TFTs. More precisely, we can isolate the i → i + 1 transform as follows: using the notations from section 3, we split the TFT at h = i + 1, then each of the n 2 inner TFTs are split at h = i. This shows that the i → i + 1 transform consists in π i l /p i full DFTs and π i atomic TFTs, where π i = p i+1 • • • p d-1 (= n 2 ) and l = l/π i . The atomic TFTs have length r i = (l rem p i ).

As a consequence, the complexity of the complete TFT is given by T(l, n) = d-1 i=0 π i l/π i p i p i f(p i ) + π i r i t(r i , p i ) + O(π i l/π i ) .

Figure 1 : 2 . 5 .

 125 Figure 1: Notations and steps of the FFT

  l

Figure 2 :

 2 Figure 2: Example of Truncated Fourier Transform (TFT)

. 4 ) 4 . 1 .

 441 RemarkIn the case l = 1, p = 2, we get the equations a 1 = b 1 -a 2 , then b 2 = b 1 -2a 2 , which corresponds to the inversion of the usual butterfly of size 2 as encountered when n = 2 k[START_REF] Van Der Hoeven | The truncated Fourier transform and applications[END_REF].

Algorithm 4 . 1 .

 41 atomicITFT (atomic inverse TFT) • INPUT: integers n ∈ N and l < n, vectors A 2 = (a i ) l i<n and B 1 = (b i ) i<l , and a primitive n-th root of unity ω • OUTPUT: the vector A 1 = (a i ) i<l and the value b l such that ∀i l, b i = Âi where A

Figure 3 :

 3 Figure 3: Example of inverse TFT. Black (resp. white) dots represent known (resp. unknown) values, ITFT* do not return additional output values

0 1 2 Figure 4 :

 24 Figure 4: Successive transforms in the FFT/TFT algorithm after complete development of the recursive calls

Lemma 5 . 4 .

 54 We have T * (l, p) 4C(p) + O(p) for large p .

  Then,K * (p) 4p(1 + o(1)) .

Since l /p i p i + r i = l (euclidean division), this rewrites

We have l/π i (l/π i ) + 1. Moreover, the term O(π i l/π i ) = O(l) + O(π i ) corresponds to the multiplication by twiddle factors at each step, as in Theorem 5.1. We then have the (otherwise abusive) simplification

This yields

By definition (5.3), we have r i k(r i , p i ) K(p i ), and it is clear that K(p i ) k(1, p i ) 1. Then, handling the cases r i = 0 and 1 r i p i -1 separately shows easily that r i (k(r i , p i ) -1) + 1 K(p i ). Finally, using equation (5.1), we get the bound: Theorem 5.5 (Complexity of the Truncated Fourier Transform (TFT)). The Truncated Fourier Transform can be performed using

field operations.

For the inverse TFT, all atomic inverse TFT (at a given row of atomic transform) do not necessarily have the same length. However the above reasoning still applies, and we get the bound: Theorem 5.6 (Complexity of the Inverse Truncated Fourier Transform (ITFT)). The Inverse Truncated Fourier Transform can be performed using

field operations.

Remark 5.5. The overhead for atomic inverse transforms K * (p) is larger than for atomic direct transforms K(p).

Remark 5.6. The bound (5.3) rewrites

Since K(p)f(p) increases with p, this shows that it is best to sort the prime factors of n in increasing order (p 0 p 1 • • • p d-1 ) to minimize the overhead. Moreover, assuming n is highly composite, p 0 • • • p i grows much faster than K(p i )f(p i ), so that the predominant term in the linear factor is K(p 0 )f(p 0 )/p 0 .