{'fraction_non_alphanumeric': 0.05512105549510337, 'fraction_numerical': 0.017784276387377583, 'mean_word_length': 4.157663977551736, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Fast Partial Fourier Transform

Given a time series vector, how can we efficiently compute a specified part of Fourier coefficients? Fast Fourier transform (FFT) is a widely used algorithm that computes the discrete Fourier transform in many machine learning applications. Despite its pervasive use, all known FFT algorithms do not provide a fine-tuning option for the user to specify one's demand, that is, the output size (the number of Fourier coefficients to be computed) is algorithmically determined by the input size. This matters because not every application using FFT requires the whole spectrum of the frequency domain, resulting in an inefficiency due to extra computation. In this paper, we propose a fast Partial Fourier Transform (PFT), a careful modification of the Cooley-Tukey algorithm that enables one to specify an arbitrary consecutive range where the coefficients should be computed. We derive the asymptotic time complexity of PFT with respect to input and output sizes, as well as its numerical accuracy. Experimental results show that our algorithm outperforms the state-of-the-art FFT algorithms, with an order of magnitude of speedup for sufficiently small output sizes without sacrificing accuracy.


INTRODUCTION
How can we efficiently compute a specified part of Fourier coefficients for a given time series vector? Discrete Fourier transform (DFT) is a crucial task in several application areas, including anomaly detection (Hou & Zhang (2007); Rasheed et al. (2009);Ren et al. (2019)), data center monitoring (Mueen et al. (2010)), and image processing (Shi et al. (2017)). Notably, in many such applications, it is well known that the DFT results in strong "energy-compaction" or "sparsity" in the frequency domain. That is, the Fourier coefficients of data are mostly small or equal to zero, having a much smaller support compared to the input size. Moreover, the support can often be specified in practice (e.g., a few low-frequency coefficients around the origin). These observations arouse a great interest in an efficient algorithm capable of computing only a specified part of Fourier coefficients.
Fast Fourier transform (FFT) is an algorithm that rapidly computes the DFT of a vector, which reduces the arithmetic complexity from naïve O(N 2 ) to O(N log N ), where N is data size. Nevertheless, all known FFT algorithms do not provide the aforementioned fine-tuning option for the user, i.e., the output size (the number of Fourier coefficients to be computed) is algorithmically determined by the input size. Such a lack of flexibility is often followed by just discarding many unused coefficients, not to mention the inefficiency due to the extra computation.
In this paper, we propose a fast Partial Fourier Transform (PFT), an efficient algorithm for computing a part of Fourier coefficients. Specifically, we consider the following problem: given a complex-valued vector a of size N , a non-negative integer M , and an integer µ, estimate the Fourier coefficients of a for the interval [µ − M, µ + M ]. The resulting algorithm is of remarkably simple structure, composed of several "smaller" FFTs combined with linear pre-and post-processing steps; consequently, we achieve O(N + M log M ) complexity of PFT. When M N , this is a significant improvement compared to the conventional FFT which never benefits from the information of M , resulting in O(N log N ) complexity, consistently.
To the best of our knowledge, PFT is the first DFT algorithm that enables one to control the output interval, providing great versatility and efficiency on the computation. There have been studies for estimating the top-k (the k largest in magnitude) Fourier coefficients of a given vector (Hassanieh et al. (2012); Iwen et al. (2007)), yet none of them grants a "freedom" of specifying the output interval. Furthermore, PFT does not require the input size to be a power of 2, unlike many other variants of FFT. This is because the idea of PFT derives from a modification of the Cooley-Tukey algorithm (Cooley & Tukey, 1965), which also makes it straightforward to extend the idea to a higher dimensionality (indeed, we present 2-dimensional PFT in Appendix B).
Through experiments, we show that PFT outperforms the state-of-the-art FFT libraries, FFTW by Frigo & Johnson (2005) and Intel Math Kernel Library (MKL), with an order of magnitude of speedup for sufficiently small output sizes without sacrificing accuracy.

RELATED WORK
We describe related works on Fast Fourier transform (FFT), as well as its applications.
Fast Fourier Transform. Cooley and Tukey proposed by far the most commonly used FFT algorithm (Cooley & Tukey, 1965), in which a DFT is recursively broken down into several smaller DFTs provided that the original data is of composite size. Johnson & Frigo (2006) reduced the arithmetic complexity of FFT to the state-of-the-art ∼ 34 9 N log N , where N is data size. Meanwhile, a great interest in an algorithm that efficiently computes only a small number of Fourier coefficients has grown because of the frequently observed energy-compaction property. A few techniques have been proposed by Hassanieh et al. (2012) and Iwen et al. (2007), which estimate the top-k Fourier coefficients of a given vector. However, all known FFT algorithms lack the ability to efficiently compute only a specified part of Fourier coefficients.
Applications of FFT. Fast Fourier transform has been widely used for anomaly detection (Hou & Zhang (2007); Rasheed et al. (2009); Ren et al. (2019)). Hou & Zhang (2007) and Ren et al. (2019) detect anomaly points of a given data by extracting a compact representation with FFT. Rasheed et al. (2009) use FFT to detect local spatial outliers which have similar patterns within a region but different patterns from the outside. Several works (Pagh (2013); Pham & Pagh (2013); Malik & Becker (2018)) exploit FFT for efficient operations. Pagh (2013) leverages FFT to efficiently compute a polynomial kernel used with support vector machines (SVMs). Malik & Becker (2018) propose an efficient tucker decomposition method using FFT. In addition, FFT has been used for fast training of convolutional neural networks (Mathieu et al. (2014); Rippel et al. (2015)) and an efficient recommendation model on a heterogeneous graph (Jin et al. (2020)).

OVERVIEW
We propose PFT, an efficient algorithm for computing a specified part of Fourier coefficients. The main challenges and our approaches are as follows: 1. How can we extract essential information for a specified output? Considering the fact that only a specified part of Fourier coefficients should be computed, we need to find an algorithm requiring fewer operations than the direct use of conventional FFT. This is achievable by carefully modifying the Cooley-Tucky algorithm, finding twiddle factors (trigonometric constants) with small oscillations, and approximating those factors using polynomial functions (Section 3.2.1).
2. How can we decrease approximation costs? The approach given above involves an approximating process, which would be computationally demanding. We propose using a base exponential function, by which all data-independent constants can be precomputed, so that one can bypass the approximation problem during the run-time (Section 3.2.2 and Section 3.3).
3. How can we further reduce numerical computation? We carefully reorder operations and factorize terms in order to alleviate the complexity of PFT. Such techniques separate all data-independent factors from data-dependent factors, allowing further precomputation. The arithmetic cost of the resulting algorithm has an asymptotic upper bound where N and M are input and output size descriptors, respectively (Section 3.4 and Section 3.5.1).
We describe details of PFT from Section 3.2 to Section 3.4, The time complexity and the approximation bound of PFT is analyzed in Section 3.5.

APPROXIMATION OF TWIDDLE FACTORS
The key of our algorithm is to approximate a part of twiddle factors with relatively small oscillations by using polynomial functions, which reduces the computational complexity of DFT due to the mixture of many twiddle factors. Using polynomial approximation also allows one to carefully control the degree of polynomial (or the number of approximating terms), enabling fine-tuning the output range and the approximation bound of the estimation. Our first goal is to find a collection of twiddle factors with small oscillations. This can be achieved by slightly adjusting the summand of DFT and splitting the summation as in the Cooley-Tukey algorithm (Section 3.2.1). Next, using a proper base exponential function, we give an explicit form of polynomial approximation to the twiddle factors (Section 3.2.2).

TWIDDLE FACTORS WITH SMALL OSCILLATIONS
Recall that the DFT is defined as follows: where a is a complex-valued vector of size N , and [ν] denotes {0, 1, · · · , ν −1} for a positive integer ν. Assume that N is a composite, so there exist p, q > 1 such that N = pq. The Cooley-Tukey algorithm re-expresses (1): a qk+l e −2πiml/N · e −2πimk/p , yielding two collections of twiddle factors, namely {e −2πiml/N } l∈ [q] and {e −2πimk/p } k∈ [p] . Consider the problem of computingâ m for −M ≤ m ≤ M , where M ≤ N/2 is a nonnegative integer. In this case, note that the exponent of e −2πiml/N ranges from −2πiM (q − 1)/N to +2πiM (q − 1)/N and that the exponent of e −2πimk/p ranges from −2πiM (p − 1)/p to +2πiM (p − 1)/p. Here (q−1)/N (p−1)/p ∼ 1 p , meaning that the first collection contains twiddle factors with smaller oscillations compared to the second one. Typically, a function with smaller oscillation results in a better approximation via polynomials. In this sense, it is reasonable to approximate the first collection of twiddle factors in (2) with polynomial functions, thereby reducing the complexity of the computation due to the mixture of two collections of twiddle factors. Indeed, one can further reduce the complexity of approximation. We slightly adjust the summand in (1) and split it: a n e −2πim(n−q/2)/N · e −πim/p = k∈[p] l∈ [q] a qk+l e −2πim(l−q/2)/N · e −2πimk/p · e −πim/p .
In (3), we observe that the range of exponents of the first collection {e −2πim(l−q/2)/N } l∈[q] of twiddle factors is [−πiM/p, +πiM/p], a contraction by a factor of around 2 when compared with [−2πiM (q − 1)/N, +2πiM (q − 1)/N ], hence the twiddle factors with even smaller oscillations. There is an extra twiddle factor e −πim/p in (3). Note that, however, it depends on neither k nor l, so the amount of the additional computation is relatively small.

BASE EXPONENTIAL FUNCTION
The first collection of twiddle factors in (3) consists of q distinct exponential functions. One can apply approximation process for each function in the collection; however, this would be timeconsuming. A more plausible approach is to 1) choose a base exponential function e uix for some fixed u ∈ R, 2) approximate e uix by using a polynomial, and 3) exploit a property of exponential functions: the laws of exponents. Specifically, suppose that we obtained a polynomial P(x) that approximates e uix on |x| ≤ |ξ|, where u, ξ are non-zero real numbers. Consider another exponential function e vix where v = 0. Since e vix = e ui(vx/u) , the re-scaled polynomial function P(vx/u) approximates e vix on |x| ≤ |uξ/v|. This observation indicates that once we find an approximation P to e uix on |x| ≤ |ξ| for properly selected u and ξ, all elements belonging to {e −2πim(l−q/2)/N } l∈ [q] can be approximated by re-scaling P. Fixing a base exponential function also enables precomputing a polynomial that approximates it, so that one can avoid solving the approximation problem during the run-time. We further elaborate this idea in a rigorous manner after giving a few definitions (see Definitions 3.1 and 3.2) and present a theoretical approximation bound in Theorem 4.
Let · R be the uniform norm (or supremum norm) restricted to a set R ⊆ R, that is, f R = sup{|f (x)| : x ∈ R} and P α be the set of polynomials on R of degree at most α. Definition 3.1. Given a non-negative integer α and non-zero real numbers ξ, u, we define a polynomial P α,ξ,u as the best approximation to e uix out of the space P α under the restriction |x| ≤ |ξ|: and P α,ξ,u = 1 when ξ = 0 or u = 0. 2 Smirnov & Smirnov (1999) proved the unique existence of such a polynomial. Also, a few techniques called minimax approximation algorithms for computing P α,ξ,u are reviewed in Fraser (1965). Definition 3.2. Given a tolerance > 0 and a positive integer r ≥ 1, we define ξ( , r) to be the scope about the origin such that the exponential function e πix can be approximated by a polynomial of degree less than r with approximation bound : We express the corresponding polynomial as In Definition 3.2, we choose e πix as a base exponential function. The rationale behind is as follows. First, using a minimax approximation algorithm, we precompute ξ( , r) and {w ,r−1,j } j∈ [r] for several tolerance 's (e.g. 10 −1 , 10 −2 , · · · ) and positive integer r's (typically 1 ≤ r ≤ 25). When N, M, p and are given, we find the minimum r satisfying ξ( , r) ≥ M/p. Then, by the preceding argument, it follows that the re-scaled polynomial function P r−1,ξ( ,r),π (−2x(l − q/2)/N ) . Therefore, we obtain a polynomial approximation on |m| ≤ M for each twiddle factor in {e −2πim(l−q/2)/N } l∈[q] , namely {P r−1,ξ( ,r),π (−2m(l − q/2)/N )} l∈ [q] . Then, it follows from (3) that Thus, we obtain an estimation ofâ m for −M ≤ m ≤ M by approximating the first collection of twiddle factors in (3).

ARBITRARILY CENTERED TARGET RANGES
In the previous section, we have focused on the problem of calculatingâ m for m belonging to [−M, M ]. We now consider a more general case: let us use the term target range to indicate the range where the Fourier coefficients should be calculated, and Note that the previously given method works only when our target range is centered at µ = 0. A slight modification of the algorithm allows the target range to be arbitrarily centered. One possible approach is as follows: given a complex-valued vector x of size N , we define y as y n = x n · e −2πiµn/N . Then, the Fourier coefficients of x and y satisfy the following relationship: x n · e −2πiµn/N · e −2πimn/N = n∈ [N ] x n · e −2πi(m+µ)n/N =x m+µ .
Therefore, the problem of calculatingx m for m ∈ R µ,M is equivalent to calculatingŷ m for m ∈ R 0,M , to which our previous method can be applied. This technique, however, requires extra N multiplications due to the computation of y.
A better approach, where one can bypass the extra process during the run-time, is to exploit the following lemma (see Appendix A.1 for the proof). Lemma 1. Given a non-negative integer α, non-zero real numbers ξ, u, and any real number µ, the following equality holds: This observation implies that, in order to obtain a polynomial approximating e uix on |x − µ| ≤ |ξ|, we first find a polynomial P approximating e uix on |x| ≤ |ξ|, then translate P by −µ and multiply it with the scalar e uiµ . Applying this process to the previously obtained approximation polynomi- . We substitute these polynomials for the twiddle factors {e −2πim(l−q/2)/N } l∈ [q] in (3), which gives the following estimation ofâ m for m ∈ R µ,M : where k ∈ [p], l ∈ [q], and j ∈ [r].

EFFICIENT SUMMATIONS
We have found that three main summation steps (each being over j, k and l) take place when computing the partial Fourier coefficients. Note that in (5), the innermost summation j is moved to the outermost position, and the term −2(m − µ)(l − q/2)/N is factorized into two independent terms, (m − µ)/p and 1 − 2l/q. Interchanging the order of summations and factorizing the term result in a significant computational benefit; we elucidate what operator we should utilize for each summation and how we can save the arithmetic costs from it. As we will see, the innermost sum over l corresponds to a matrix multiplication, the second sum over k can be viewed as multiple DFTs, and the outermost sum over j is an inner product.
For the first sum, let A = (a kl ) = a qk+l and B = (b lj ) = e −2πiµ(l−q/2)/N w ,r−1,j (1 − 2l/q) j , so that (5) can be written as follows: Here, note that the matrix B is data-independent (not dependent on a), and thus can be precomputed. Indeed, we have already seen that {w ,r−1,j } j∈ [r] can be precomputed. The other factors e −2πiµ(l−q/2)/N and (1 − 2l/q) j composing the elements of B can also be precomputed if (N, M, µ, p, ) is known in advance. Thus, as long as the setting (N, M, µ, p, ) is unchanged, we can reuse the matrix B for any input data a once the configuration phase of PFT is completed (Algorithm 1). We shall denote the multiplication A × B as C = (c kj ): c kj · e −2πimk/p · ((m − µ)/p) j · e −πim/p .
For each j ∈ [r], the summationĉ j;m = k∈[p] c kj · e −2πimk/p is a DFT of size p. We perform FFT r times for this computation and denote the corresponding Fourier coefficient asĉ j;m , which yields the following estimation ofâ m : Note thatĉ j;m is a periodic function of period p with respect to m, so we use the coefficient at m modulo p when m < 0 or m ≥ p. Therefore, the m th Fourier coefficient of a can be calculated approximately by the inner product ofĉ j;m and ((m − µ)/p) j with respect to j, followed by a multiplication with the extra twiddle factor e −πim/p (we also precompute ((m − µ)/p) j and e −πim/p ). The full computation is outlined in Algorithm 2.
Algorithm 1: Configuration phase of PFT input : Input size N , output descriptors M and µ, divisor p, and tolerance output: Matrix B, divisor p, and numbers of rows and columns, q and r

THEORETICAL ANALYSIS
We give theoretical analysis regarding the time complexity of PFT as well as its approximation bound.

TIME COMPLEXITY
We analyze the time complexity of PFT. Theorem 3 shows that the time cost T (N, M ) of PFT, where N and M are input and output size descriptors, respectively, is bounded by O(N + M log M ), provided N has sufficiently many divisors. 1 When M N , this is a significant improvement compared to the conventional FFT which never benefits from the information of M , resulting in O(N log N ) time cost, consistently. Before presenting the theorem, we consider an explicit condition which guarantees the "sufficiently many divisors" property. A positive integer is called b-smooth if none of its prime factors is greater than b. For example, the 2-smooth integers are equivalent to the powers of 2. Lemma 2 says that if N is a smooth number, then given any 0 < M ≤ N , one can always find a divisor of N that is tightly bounded by Θ(M ). We leave the proofs of the lemma and theorem in Appendices A.2 and A.3. Lemma 2. Let b ≥ 2. If N is b-smooth and M ≤ N is a positive integer, then there exists a positive divisor p of N satisfying M/ √ b ≤ p < √ bM .

APPROXIMATION BOUND
We now give a theoretical approximation bound of the estimation via the polynomial P. We denote the estimated Fourier coefficient of a as E(â). The following theorem states that the approximation bound is data-dependent of the total weight a 1 of the original vector, where · 1 denotes the 1 norm, and the given tolerance (see Appendix A.4 for the proof). Theorem 4. Given a tolerance > 0, the following inequality holds: where R µ,M is the target range.

EXPERIMENTS
Through experiments, the following questions should be answered: • Q1. Run-time cost (Section 4.2). How quickly does PFT compute a part of Fourier coefficients compared to other competitors without sacrificing accuracy?
• Q2. Effect of hyper-parameter p (Section 4.3). How the different choices of divisor p of input size N affect the overall performance of PFT?
• Q3. Anomaly detection (Section 4.4). How well does PFT work for a practical application using FFT (anomaly detection)?

EXPERIMENTAL SETUP
Machine. All experiments are performed on a machine equipped with Intel Core i7-6700HQ @ 2.60GHz and 8GB of RAM.
Datasets. We use both synthetic and real-world datasets listed in Table 1.
Competitors. We compare PFT with two state-of-the-art FFT algorithms, FFTW and MKL. All of them including PFT are implemented in C++.
1. FFTW: FFTW 4 is one of the fastest public implementation for FFT, offering a hardwarespecific optimization. We use the optimized version of FFTW 3.3.5, and do not include the pre-processing for the optimization as the run-time cost.
2. MKL: Intel Math Kernel Library 5 (MKL) is a library of optimized math routines including FFT, and often shows a better run time result than the FFTW. All the experiments are conducted with an Intel processor for the best performance. Run time vs. Input size (Target range: R 0, 2 9 ) (a)  Figure 1: (a) Run time vs. input size for target range R 0,2 9 with {S n } 22 n=12 datasets, and (b) run time vs. output size for S 22 . We observe that PFT outperforms both FFTW and MKL if the output size is small enough compared to the input size. When the output is sufficiently smaller than the input, the performance gain is significant: an order of magnitude of speedup.

Measure.
In all experiments, we use single-precision floating-point format, and the parameters p and are chosen so that the relative 2 error is strictly less than 10 −6 , which ensures that the overall estimated coefficients have at least 6 significant figures. Explicitly, whereâ is the actual coefficient, E(â) is the estimated coefficient, and R is the target range.

RUN-TIME COST
Run time vs. input size. We fix the target range to R 0,2 9 and evaluate the run time of PFT vs. input sizes N : 2 12 , 2 13 , · · · , 2 22 . The results are averaged over 10 thousands runs. Figure 1(a) shows how the three competitive algorithms scale with varying input size, wherein PFT outperforms both FFTW and MKL provided that the output size is small enough (< 10%) compared to the input size. Consequently, PFT achieves up to 21× speedup compared to its competitors. Due to the overhead of the O(N ) pre-and O(M ) post-processing steps, we see that PFT runs slower than FFT when M is close to N , so the time complexity tends to O(N + N log N ).
Run time vs. output size. In this experiment, we fix the input size to N = 2 22 and evaluate the run time of PFT vs. target ranges R 0,2 9 , R 0,2 10 , · · · , R 0,2 18 . The result is illustrated as a run time vs. output size plot (recall that |R 0,M | 2M ) in Figure 1(b), where each point is an average over 10 thousands runs. Note that the run times of FFTW and MKL are consistent because they do not benefit from the information of the output size M . We find that when the output size is sufficiently smaller than the input size, PFT is up to 13.8× faster than the other competitors.
Real-world data. When it comes to real-world data, it is not generally the case that the size of an input vector is a power of 2. Notably, PFT still shows a promising performance regardless of the fact that the input size is not a power of 2 or not even a highly composite 6 number: a strong indication that our proposed technique is robust for many different applications in real-world.
• Urban Sound dataset contains various sound recording vectors of size N = 32000 = 2 8 × 5 3 . We evaluate the run time of PFT vs. output sizes: 100, 200, 400, 800, 1600, 3200, and 6400. Figure 2(a) illustrates the average run times of the three competitive algorithms. We see that PFT outperforms both FFTW and MKL if the output size is small enough compared to the input size.   • Air Condition dataset is composed of time series vectors of size N = 19735 = 5 × 3947.
Note that N has only two non-trivial divisors, namely 5 and 3947, forcing one to choose p = 3947 in any practical settings; if we choose p = 5, the value M/p often turns out to be too large, which results in a poor performance. We evaluate the run time of PFT vs. output sizes: 125, 250, 500, 1000, 2000, 4000, 8000, and 16000, as shown in Figure 2(b). It is noteworthy that PFT still outperforms its competitors even in such pathological examples, which implies the robustness of our algorithm for various real-world situations.

EFFECT OF HYPER-PARAMETER p
To investigate the effect of different choices of p, we fix N = 2 22 and vary the ratio M/p from 1/32 to 4 for different target ranges: R 0,2 9 , R 0,2 10 , · · · , R 0,2 18 . Table 2 shows the resulting run time for each setting, where the bold highlights the best choice of M/p for each M , and the missing entries are due to worse performance than the FFT. One crucial observation is as follows: with the increase of output size, the best choice of the ratio M/p also increases or, equivalently, the optimal value of p tends to remain stable. Intuitively, this is the consequence of "balancing" the three summation steps (Section 3.4): when M N , the most computationally expensive operation is the matrix multiplication with O(rN ) time complexity, and thus, M/p should be small so that the number r of approximating terms decreases, despite the sacrifices in the batch FFT step requiring O(rp log p) operations (Appendix A.3). As the M becomes larger, however, more concern is needed regarding the batch FFT and post-processing steps, so the parameter p should not change rapidly. This observation, even though we do not provide an explicit formulation, indicates the possibility that the optimal value of p can be algorithmically auto-selected given a setting (N, M, µ, ), which we leave as a future work.

ANOMALY DETECTION
We demonstrate an example of how PFT is applied to practical applications. Here is one simple but fundamental principle: replace the "perform FFT and discard unused coefficients" procedure with "just perform PFT". Considering the anomaly detection method proposed in Rasheed et al. (2009), where one first performs FFT and then inverse FFT with only a few low-frequency coefficients to obtain an estimated fitted curve, we can directly apply the principle to the method. To verify this experimentally, we use a time series vector from Air Condition dataset, and set the target range as R 0,125 ( 250 low-frequency coefficients). Note that, in this setting, PFT results in around 8× speedup compared to the conventional FFT (see Figure 2(b)). The top-20 anomalous points detected from the data are presented in Figure 3. In particular, we found that replacing FFT with PFT does not change the result of top-20 anomaly detection, with all its computational benefits.

CONCLUSIONS
In this paper, we propose PFT (fast Partial Fourier Transform), an efficient algorithm for computing a specified part of Fourier coefficients. PFT approximates some of twiddle factors with relatively small oscillations using polynomial functions, reducing the computational complexity of DFT due to the mixture of many twiddle factors. Experimental results show that our algorithm outperforms the state-of-the-art FFT algorithms, FFTW and MKL, with an order of magnitude of speedup for sufficiently small output sizes without sacrificing accuracy. Future works include optimizing the implementation of PFT; for example, the optimal divisor p of input size N might can be algorithmically auto-selected. We also believe that hardware-specific optimizations (similar to FFTW or MKL) would further increase the performance of PFT.

A PROOFS
A.1 PROOF OF LEMMA 1 Proof. Recall that the polynomial P α,ξ,u is defined by arg min P ∈Pα P (x)−e uix |x|≤|ξ| . If P (x) ∈ P α , it is clear that e uiµ · P (x − µ) ∈ P α , because translation and non-zero scalar multiplication on a polynomial do not change its degree. Thus, we may re-express the definition of P α,ξ,u as follows: where the third equality holds since |e uiµ | = 1, and hence the proof.
A.2 PROOF OF LEMMA 2 Proof. Suppose that none of N 's divisors belongs to [M/ √ b, √ bM ). Let 1 = p 1 < p 2 < · · · < p d = N be the enumeration of all positive divisors of N in increasing order. It is clear that p 1 < √ bM and M/ √ b < p d since b ≥ 2 and 1 ≤ M ≤ N . Then, there exists an i ∈ {1, 2, · · · , d − 1} so that p i < M/ √ b and p i+1 ≥ √ bM . Since N is b-smooth and p i < N , at least one of 2p i , 3p i , · · · , bp i must be a divisor of N . However, this is a contradiction because we have p i+1 /p i > ( , bp i can be a divisor of N , which completes the proof.

A.3 PROOF OF THEOREM 3
Proof. Following the convention in counting FFT operations, we assume that all data-independent elements such as configuration results B, p, q, r and twiddle factors are precomputed, and thus not included in the run-time cost. We begin with construction of the matrix A. For this, we merely interpret a as an array representation for A of size p × q = N . Also, recall that the matrix B can be precomputed as described in Section 3.4. For the two matrices A of size p × q and B of size q × r, standard matrix multiplication algorithm has running time of O(pqr) = O(r · N ). Next, the expression (6) contains r DFTs of size p. We use FFT multiple times for the computation, then it is easy to see that the time cost is given by O(r · p log p). Finally, there are 2M + 1 coefficients to be calculated in (7) Note that r is only dependent of and M/p by its definition. Therefore, when is fixed, T (N, M ) is dependent of the choice of p. By Lemma 2, we can always find a p = Θ(M ). In this case, M/p is bounded, and thus, so is r. Then, from (8), we obtain the following asymptotic upper bound with respect to N and M : hence the proof.