# Fast Fourier Transform
###### tags: `Summer Project` `FFT`

[TOC]

## Discrete Fourier Transform
The DFT sequence $X_0,X_2,\ldots,X_{N-1}$ of sequence $x_0,x_1,\ldots,x_{N-1}\in\mathbb{C}$ is defined by
$$
X_k=\sum_{n=0}^{N-1}x_ne^{-i2\pi k\frac{n}{N}},\quad k=0,1,\ldots,N-1
$$
Let $w=e^{-i\frac{2\pi}{N}}$, then
$$
X=\begin{pmatrix}
X_0\\X_1\\\vdots\\X_{N-1}
\end{pmatrix}=
\begin{pmatrix}
1&1&\dots&1\\
1&w&\dots&w^{N-1}\\
\vdots&\vdots&\ddots&\vdots\\
1&w^{N-1}&\dots&w^{(N-1)(N-1)}
\end{pmatrix}
\begin{pmatrix}
x_0\\x_1\\\vdots\\x_{N-1}
\end{pmatrix}
=Wx=\mathcal F_N(x)
$$
where $W$ is called DFT matrix, $\mathcal F_N$ is called $N$-points Fourier transform.

### Inverse Discrete Fourier Transform
$$
x_n=\frac{1}{N}\sum_{k=0}^{N-1}X_ke^{i2\pi k\frac{n}{N}},\quad n=0,1,\ldots,N-1
$$

Remark that $W$ is invertible and $W^{-1}=\frac{1}{N}W^\star$, thus $x=\frac{1}{N}W^\star X$.

### Two Dimension and More
Let $x=\{x_{ij}\mid 0\leq i<N,0\leq j<M\}$ then the DFT of $x$ is defined by
$$
X_{uv}=\sum_{n=0}^{N-1}\sum_{m=0}^{M-1}x_{nm}e^{-i2\pi\left(u\frac{n}{N}+v\frac{m}{M}\right)},\quad 0\leq u<N,0\leq v<M
$$
and IDFT
$$
x_{nm}=\frac{1}{NM}\sum_{u=0}^{N-1}\sum_{v=0}^{M-1}X_{uv}e^{i2\pi\left(n\frac{u}{N}+m\frac{v}{M}\right)},\quad 0\leq n<N,0\leq m<M
$$
For higher dimension, let $I=\{(n_1,n_2,\ldots,n_M)\mid0\leq n_i<N_i\}$, we define
$$
X_{k}=\sum_{n\in I}x_ne^{-i2\pi k\cdot\frac{n}{N}},\quad k\in I
$$

### Example
- Let $x_n=c$ for $n=0,1,\ldots,N-1$, then
$$
X_k=
\begin{cases}
cN&,k=0\\
0&,k\neq 0
\end{cases}
$$
- Let $y_n=cx_n+dz_n$ for $n=0,1,\ldots,N-1$, then
$$
\begin{split}
Y_k=&\sum_{n=0}^{N-1}y_nW_N^{kn}\\
=&\sum_{n=0}^{N-1}(cx_n+dz_n)W_N^{kn}\\
=&cX_k+dZ_k
\end{split}
$$
- Let $y_n=x_{n-c}$ for $n=0,1,\ldots,N-1$ and $c\in\mathbb N$, then
$$
\begin{split}
Y_k=&\sum_{n=0}^{N-1}y_nW_N^{kn}=\sum_{n=0}^{N-1}x_{n-c}W_N^{kn}=\sum_{n=-c}^{N-c-1}x_nW_N^{k(n+c)}\\
=&W_N^{kc}\sum_{n=0}^{N-1}x_nW_N^{kn}=W_N^{kc}X_k
\end{split}
$$

## Time-Frequency Perspective
Given a function $x_c:\mathbb R\to\mathbb C$, the *sampling* of $x_c$ is
$$
x_n=x_c(nT),\quad n=0,1,\ldots,L-1
$$
where
- $T$ is called the *sampling period*
- $F_s=\frac{1}{T}$ is called *sampling rate*
- $\Omega_s=\frac{2\pi}{T}$ is called *sampling frequency*
- $L$ is called the *length* of signal

### Example
Let $x_c(t)=\cos(2\pi ft)$, which is periodic with period $\frac{1}{f}$, then the DFT of $x_n$ is
$$
\begin{split}
X_k=&\sum_{n=0}^{L-1}x_ne^{-i2\pi k\frac{n}{L}}\\
=&\sum_{n=0}^{L-1}\cos(2\pi fnT)e^{-i2\pi fTn}\quad\mbox{when}\quad k=fTL=\frac{f}{F_s}L\\
=&\sum_{n=0}^{L-1}\frac{e^{i2\pi fnT}+e^{-i2\pi fnT}}{2}e^{-i2\pi fTn}\\
=&\frac{L}{2}+\frac{1}{2}\sum_{n=0}^{L-1}(e^{ir})^n\quad\mbox{where}\quad r=-4\pi fT\\
=&\frac{L}{2}+\frac{1}{2}\left(\frac{1-e^{irL}}{1-e^{ir}}\right)
\end{split}
$$
which is a peak on real line intuitively when $k\approx\frac{f}{F_s}L$. Suppose it is not a peak, that means
$$
-L=\frac{1-e^{irL}}{1-e^{ir}}\\
-L+L\cos(r)+iL\sin(r)=L-Le^{ir}=1-e^{irL}=1-\cos(rL)-i\sin(rL)
$$
thus
$$
L\cos(r)+\cos(rL)=L+1,\quad \cos(r)=\cos(rL)=1
$$
only when
$$
-4\pi fT=r=2\pi n,\quad n\in\mathbb Z\\2\frac{f}{F_s}=2fT=-n
$$
Hence, we need sampling rate $F_s> 2f$ to avoid this, which is called *reconstruction* condition.

Notice that when $k=(1-\frac{f}{F_s})L$,
$$
\begin{split}
X_k=&\sum_{n=0}^{L-1}\cos(2\pi fnT)e^{-i2\pi (1-fT)n}\\
=&\sum_{n=0}^{L-1}\frac{e^{i2\pi fnT}+e^{-i2\pi fnT}}{2}e^{i2\pi fTn}\\
=&\frac{L}{2}+\frac{1}{2}\sum_{n=0}^{L-1}e^{irn}\quad\mbox{where}\quad r=4\pi fT\\
=&\frac{L}{2}+\frac{1}{2}\left(\frac{1-e^{irL}}{1-e^{ir}}\right)
\end{split}
$$
which implies $X_k$ is symmetrical and it is called *double-sided amplitude spectrum*. If we change
$$
Y_k=\begin{cases}
X_k,&\quad k=0\vee k=L/2+1\\
2X_k,&\quad k=1,2,\ldots,\lfloor L/2\rfloor
\end{cases}
$$
then $Y_k$ is called *single-sided amplitude spectrum*.

Fs = 1000;
T = 1/Fs;  
L = 300;
t = (0:L-1)*T;
xc = @(x) 0.7*sin(2*pi*50*x) + sin(2*pi*120*x) + rand(size(x));
X = xc(t);
plot(1000*t,X)
xlabel('t (milliseconds)')
ylabel('X(t)')


Y = fft(X);
P2 = abs(Y/L);
f = Fs*(0:L-1)/L; % k=f/Fs*L, f=Fs*f/L
plot(f,P2) 
xlabel('f (Hz)')
ylabel('|P2(f)|')


P1 = P2(1:L/2+1); % L is even
P1(2:end-1) = 2*P1(2:end-1); % L is even
f = Fs*(0:(L/2))/L;
plot(f,P1) 
xlabel('f (Hz)')
ylabel('|P1(f)|')


### Zero Padding
Given sequence $x'_n$ with length $N$, padding zero with length $M$ is
$$
x_n=
\begin{cases}
x'_n,&0\leq n<N\\
0,&N\leq n<N+M
\end{cases}
$$
then the DFT is
$$
X_k=\sum_{n=0}^{N+M-1}x_nW_{N+M}^{kn}=\sum_{n=0}^{N-1}x_nW_{N+M}^{kn}
$$
Similar to above, the peak appears at $k=fT(N+M)$ and $k=(1-fT)(N+M)$. This will

- increase the resolution of frequency
- extend the length so that it performs FFT fast

| |Original | Padding |
|-|-------- | -------- |
|$x_n$ || |
| $\mid P_1(f)\mid$ |||

## Fast Fourier Transform
Assume $N=2^k,k\in\mathbb N\cup\{0\}$. Let $W_N=e^{-i\frac{2\pi}{N}}$
$$
\begin{split}
X_k=&\sum_{n=0}^{N-1}x_nW_N^{kn}=\sum_{n=0}^{N/2-1}x_{2n}W_N^{2nk}+\sum_{n=0}^{N/2-1}x_{2n+1}W_N^{(2n+1)k}\\
=&\sum_{n=0}^{N/2-1}x_{2n}W_N^{2nk}+W_N^k\sum_{n=0}^{N/2-1}x_{2n+1}W_N^{2nk}\\
=&\sum_{n=0}^{N/2-1}x_{2n}W_{N/2}^{nk}+W_N^k\sum_{n=0}^{N/2-1}x_{2n+1}W_{N/2}^{nk}
\end{split}
$$
which is calculated by the DFT of sequence $x_0,x_2,\ldots,x_{N-2}$ and $x_1,x_3,\ldots,x_N-1$ with length $N/2$. Thus we do *Divide and Conquer*.

Note that if $N/2\leq k<N$, it is not defined in sequence of length $N/2$, then with

- $W_{N/2}^{n(k+N/2)}=(e^{-i\frac{2\pi}{N/2}})^{n(k+N/2)}=(e^{-i\frac{2\pi}{N/2}})^{nk}e^{-i2\pi n}=(e^{-i\frac{2\pi}{N/2}})^{nk}=W_{N/2}^{nk}$
- $W_N^{k+N/2}=(e^{-i\frac{2\pi}{N}})^{k+N/2}=e^{-i\frac{2\pi}{N}k}e^{-i\pi}=-e^{-i\frac{2\pi}{N}k}=-W_N^k$

thus

$$
\begin{split}
X_{k+N/2}=&\sum_{n=0}^{N/2-1}x_{2n}W_{N/2}^{n(k+N/2)}+W_N^{k+N/2}\sum_{n=0}^{N/2-1}x_{2n+1}W_{N/2}^{n(k+N/2)}\\
=&\sum_{n=0}^{N/2-1}x_{2n}W_{N/2}^{nk}-W_N^{k}\sum_{n=0}^{N/2-1}x_{2n+1}W_{N/2}^{nk}
\end{split}
$$

this method is called *radix-2 DIT Cooley–Tukey Algorithm*

### Pseudo Code

funtion fft(N, x = x[0],x[1],...,x[N-1]): // N=2^k
    if N == 1:
        return x
    Xe = fft(N/2, (x[0],x[2],...,x[N-2])
    Xo = fft(N/2, (x[1],x[3],...,x[N-1])
    w = exp(-2i * pi / N)
    for k from 0 to N/2 - 1:
        X[k] = Xe[k] + w^k * Xo[k]
        X[k + N/2] = Xe[k] - w^k * Xo[k]
    return X

and to perform IFFT, it remains to change `w=exp(2i*pi/N)` and lastly to be divided by `N`.

### Time Complexity
$$
O(\mathcal F_1)=O(1)\\
\begin{split}
O(\mathcal F_N)=&O(\mathcal F_{\frac{N}{2}}+\mathcal F_{\frac{N}{2}}+N)=O(2\mathcal F_{\frac{N}{2}}+N)\\
=&O(2(2\mathcal F_{\frac{N}{4}}+\frac{N}{2})+N)=O(4\mathcal F_{\frac{N}{4}}+2N)\\
&\vdots\\
=&O(2^k\mathcal F_1+kN)=O(N+N\log_2N)\\
=&O(N\log N)
\end{split}
$$

### Matlab Implementation

function y = yafft(x)
    N = length(x);
    if N == 1
        y = x;
        return
    end
    
    xe = zeros(1, N/2);
    xo = zeros(1, N/2);
    for k = 1:N/2
        xe(k) = x(2*k-1);
        xo(k) = x(2*k);
    end
    
    ye = yafft(xe);
    yo = yafft(xo);
    w = exp(-2i*pi/N);
    
    y = zeros(1, N);
    for k = 1:N/2
        y(k) = ye(k) + w^(k-1) * yo(k);
        y(k+N/2) = ye(k) - w^(k-1) * yo(k);
    end
end

## General Cooley–Tukey Algorithm
Let $N=M\cdot L$ where $M\neq 1\neq L$, then
$$
\begin{split}
X_k=&\sum_{n=0}^{N-1}x_nW_N^{kn}
=\sum_{l=0}^{L-1}\sum_{m=0}^{M-1}x_{mL+l}W_N^{k(mL+l)}\\
=&\sum_{l=0}^{L-1}W_N^{kl}\sum_{m=0}^{M-1}x_{mL+l}W_N^{kmL}=\sum_{l=0}^{L-1}W_N^{kl}\sum_{m=0}^{M-1}x_{mL+l}W_M^{km}
\end{split}
$$
which splits the original problem into $L$ FFT problems of length $M$. Similarly, if $k\geq M$, we have

- $W_M^{(aM+k)m}=W_M^{aMm}W^{km}=e^{-i2\pi lm\frac{M}{M}}W_M^{km}=W_M^{km}$

thus
$$
\begin{split}
X_{aM+k}=&\sum_{l=0}^{L-1}W_N^{(aM+k)l}\sum_{m=0}^{M-1}x_{mL+l}W_M^{(aM+k)m}\\
=&\sum_{l=0}^{L-1}W_N^{(aM+k)l}\sum_{m=0}^{M-1}x_{mL+l}W_M^{km}
\end{split}
$$

### Pseudo Code

function fft(N, x = x[0],x[1],...,x[N-1]): \\ N=M*L
    for l from 0 to L-1:
        for m from 0 to M-1:
            z[l][m] = x[m*L+l]
    for l from 0 to L-1:
        Z[l] = fft(M, z[l][:])
    w = exp(-2i*pi/N)
    for a from 0 to L-1:
        for k from 0 to M-1:
            X[a*M+k] = 0
            for l from 0 to L-1:
                X[a*M+k] += w^((a*M+k)*l) * Z[l][k]
    return X

### Time Complexity
$$
O(\mathcal F_N)=O(LM+L\mathcal F_M+L^2M)=O(NL+L\mathcal F_{\frac{N}{L}})
$$

## Rader's Algorithm
Let $N>2$ be a prime, then $\mathbb Z_N^\times$ is a cyclic group, and there exists primitive root $g$. That is, $f_g:\mathbb Z_{N-1}\to\mathbb Z_{N}\setminus\{0\}$ by
$$
f_g(n)=g^n\pmod N
$$
is a bijection. Recall we have
$$
X_0=\sum_{n=0}^{N-1}x_n\\
X_k=x_0+\sum_{n=1}^{N-1}x_nW_N^{kn},\quad k=1,2,\ldots,N-1
$$
and we change the order
$$
X_k=x_0+\sum_{q=0}^{N-2}x_{g^q}W_N^{kg^q},\quad k=1,2,\ldots,N-1\\
X_{g^{-p}}=x_0+\sum_{q=0}^{N-2}x_{g^q}W_N^{g^{q-p}},\quad p=0,1,\ldots,N-2
$$
Now let $a_k=x_{g^{k}},b_k=W_N^{g^{-k}}$, then
$$
X_{g^{-p}}-x_0=\sum_{q=0}^{N-2}a_{q}b_{p-q},\quad p=0,1,\ldots,N-2
$$
is called the *cyclic convolution* $b\star a$ of $a$ and $b$.

### Cyclic Convolution
Calculating the cyclic convolution of two sequence $a_k$ and $b_k$ is another problem that can be solved by FFT.

:::info
**Property.**
Let $a=(a_0,\ldots,a_{N-1}),b=(b_0,\ldots,b_{N-1})$, then
$$
\mathcal F_N(a\star b)=\mathcal F_N(a)\cdot\mathcal F_N(b)
$$
**Proof.**
Let $c=(c_0,\ldots,c_{N-1})$ where $c_n=\sum_{k=0}^{N-1}b_ka_{n-k}$ and $A=\mathcal F_N(a),B=\mathcal F_N(b)$, then
$$
\begin{split}
c_n=&\sum_{k=0}^{N-1}b_ka_{n-k}\\
=&\sum_{k=0}^{N-1}b_k\left(\frac{1}{N}\sum_{m=0}^{N-1}A_{m}W_N^{-m(n-k)}\right)\\
=&\frac{1}{N}\sum_{k=0}^{N-1}b_k\left(\sum_{m=0}^{N-1}A_{m}W_N^{-mn}\right)W_N^{mk}\\
=&\frac{1}{N}\sum_{m=0}^{N-1}A_{m}\left(\sum_{k=0}^{N-1}b_kW_N^{mk}\right)W_N^{-mn}\\
=&\frac{1}{N}\sum_{m=0}^{N-1}A_{m}B_mW_N^{-mn}
\end{split}
$$
:::
Also, with proper way of extending the length, we can have the same convolution result.
:::info
**Property.**
Let $a,b\in\mathbb {C}^N$,, $a',b',c'\in\mathbb{C}^M$ where $M\geq 2N-1$ and

- $a'_k=\begin{cases}a_k,&\quad k=0,1,\ldots,N-1\\0,&\quad k=N,\ldots,M-1\end{cases}$
- $b'_k=\begin{cases}b_k,&\quad k=0,1,\ldots,N-1\\ 0,&\quad k=N,\ldots,M-N\\b_{k+N-M},&\quad k=M-N+1,\ldots,M-1\end{cases}$

namely,

$$
(a'\star b')^T=
\begin{pmatrix}
\begin{matrix}
a_0&0&\dots&0\\
a_1&a_0&\dots&0\\
\vdots&\vdots&\ddots&\vdots\\
a_{N-1}&a_{N-2}&\dots&a_0
\end{matrix}&\dots&
\begin{matrix}
a_{N-1}&a_{N-2}&\dots&a_1\\
0&a_{N-1}&\dots&a_2\\
\vdots&\vdots&\ddots&\vdots\\
0&0&\dots&a_{N-1}
\end{matrix}\\
\vdots&&\vdots
\end{pmatrix}
\begin{pmatrix}
b_0\\b_1\\\vdots\\b_{N-1}\\0\\\vdots\\0\\b_1\\\vdots\\b_{N-1}
\end{pmatrix}
$$

Then
$$
a\star b=a'\star b'\mid_{n=0,1,\ldots,N-1}
$$
namely, the first $N$ items of $a'\star b'$.

**Proof.**
Let $c=a\star b$, $c'=a'\star b'$, then for $0\leq k< N$,
$$
\begin{split}
c'_k=&\sum_{n=0}^{M-1}b'_na'_{k-n}\\
=&\sum_{n=0}^{k}b'_na'_{k-n}+\sum_{n=k+1}^{M-N+k}b'_na'_{k-n}+\sum_{n=M-N+k+1}^{M-1}b'_na'_{k-n}\\
=&\sum_{n=0}^{k}b_na_{k-n}+\sum_{n=M-N+k+1}^{M-1}b'_na'_{k-n}\\
=&\sum_{n=0}^{k}b_na_{k-n}+\sum_{n=k+1}^{N-1}b'_{M-N+n}a'_{M-N+k-n}\\
=&\sum_{n=0}^{k}b_na_{k-n}+\sum_{n=k+1}^{N-1}b_na_{k-n}=\sum_{n=0}^{N-1}b_na_{k-n}=c_k
\end{split}
$$
:::

function cyclic_convolution(N, a, b):
    M = minimum 2-power at least 2N-1
    for k from 0 to N-1:
        A[k] = a[k], B[k] = b[k]
    for k from N to M-N:
        A[k] = 0, B[k] = 0
    for k from M-N+1 to M-1:
        A[k] = 0, B[k] = b[k+N-M]
    A2 = fft(M, A)
    B2 = fft(M, B)
    return ifft(N, A2.*B2)[0:N-1]

and the time complexity of cyclic convolution is
$$
O(M+3M\log M)=O(N\log N)
$$

### Pseudo Code
The naive way to find the primitive root is by checking all $g$ and all powers, but it can be speed up due to the proterty.
:::info
**Property.**
Let $N$ be prime, $p_1,p_2,\ldots,p_k$ be prime factors of $N-1$ then $g$ is primitive root of $N$ if and only if
$$
g^\frac{N-1}{p_i}\neq 1\pmod N,\quad i=1,2,\ldots,k
$$

**Proof.**
Assume $g$ is not a primitive root, then $g^p=1\pmod N$ for some smallest $0<p<N-1$. Since $g^{N-1}=1\pmod N$, we have $p\mid N-1$. And thus $p_i\mid \frac{N-1}{p}$ for some $p_i$, and
$$
g^{\frac{N-1}{p_i}}=g^{p(N-1)/(pp_i)}=(g^p)^{(N-1)/pp_i}=1\pmod N
$$
:::

function primitive_root(N): // N is prime
    S = prime factors of N-1
    n = length(S)
    for g from 2 to n-1:
        for k in S:
            if g^((n-1)/k) mod N == 1:
                break
        else:
            return g

function fft(N, x = x[0],x[1],...,x[N-1]): // N is prime
    X[0] = sum(x)
    g = primitive_root(N)
    w = exp(-2i * pi / N)
    for k from 0 to N-2:
        a[k] = x[g^k]
        b[k] = w^(g^{-k})
    c = cyclic_convolution(N-1, a, b)
    for p from 0 to N-2:
        X[g^(-p)] = x[0] + c[p]
    return X

## Bluestein's algorithm
Note that
$$
nk=\frac{-(k-n)^2}{2}+\frac{k^2}{2}+\frac{n^2}{2}
$$
thus
$$
\begin{split}
X_k=&\sum_{n=0}^{N-1}x_nW_N^{kn}\\
=&\sum_{n=0}^{N-1}x_nW_N^{-(k-n)^2/2}W_N^{k^2/2}W_N^{n^2/2}\\
=&W_N^{k^2/2}\sum_{n=0}^{N-1}\left(x_nW_N^{n^2/2}\right)W_N^{-(k-n)^2/2}\\
=&W_N^{k^2/2}\sum_{n=0}^{N-1}a_nb_{k-n}\quad\mbox{where}\quad a_k=x_kW_N^{k^2/2},b_k=W_N^{-k^2/2}
\end{split}
$$

function fft(N, x): // for all N
    w = exp(-2i*pi/N)
    for k from 0 to N-1:
        a[k] = x[k] * w^(k^2/2)
        b[k] = w^(-k^2/2)
    c = cyclic_convolution(N, a, b)
    for k from 0 to N-1:
        X[k] = w^(k^2/2) * c[k]
    return X

and the time complexity is
$$
O(\mathcal F_N)=O(N+N\log N+N)=O(N\log N)
$$

## Application

### Multiplication of Polynomial
Let $p(x)=a_0+a_1x+\ldots+a_nx^n,q(x)=b_0+b_1x+\ldots+b_mx^m$, $a_n\neq0\neq b_m$. By interpolation, we find

$$
(W_{n+m+1}^k,p(W_{n+m+1}^k))\quad\mbox{and}\quad (W_{n+m+1}^k,q(W_{n+m}^k)),\quad k=0,1,\ldots,n+m
$$

then $r(x)$ has points
$$
(W_{n+m}^k,p(W_{n+m}^k)q(W_{n+m}^k)),\quad k=0,1,\dots,n+m
$$

function padding(a, N, v):
    return [a [v]*(N-length(a))]

function mul(a, b): 
    n = length(a), m = length(b)
    A = fft(n+m+1, padding(a, n + m + 1, 0))
    B = fft(n+m+1, padding(b, n + m + 1, 0))
    return ifft(n+m+1, A.*B)

### Solution of Circulant matrix
Let $A=(a_{ij})\in\mathbb{M}_{n\times n}(\mathbb{C})$, $A$ is called *circulant* if
$$
a_{i,j}=a_{i+1,j+1},\quad i,j\in\mathbb{Z}_n
$$
namely, in the form
$$
\begin{pmatrix}
a_0&a_{n-1}&\dots&a_2&a_{1}\\
a_1&a_0&\dots&a_3&a_2\\
\vdots&\vdots&\ddots&\vdots&\vdots\\
a_{n-2}&a_{n-3}&\dots&a_0&a_{n-1}\\
a_{n-1}&a_{n-2}&\dots&a_1&a_0
\end{pmatrix}
$$
To solve linear system
$$
Ax=b
$$
where $A$ is circulant, by the fact that it is equivalent to solve $a\star x=b$ where $a=(a_0,a_1,\ldots,a_{n-1})$ and $\star$ is the cyclic convolution. Thus
$$
\mathcal F_n(a\star x)=\mathcal F_n(a)\cdot\mathcal F_n(x)=\mathcal F_n(b)
$$
and hence
$$
\mathcal F_n(x)=\frac{\mathcal F_n(b)}{\mathcal F_n(a)},\quad x=\mathcal F_n^{-1}(\frac{\mathcal F_n(b)}{\mathcal F_n(a)})
$$

function solve(a, b):
    n = length(a)
    A = fft(n, a)
    B = fft(n, b)
    return ifft(n, A./B)

### aplusb
[Kattis](https://open.kattis.com/problems/aplusb)
Given $N$ integers $a_0,a_1,\ldots,a_{N-1}$ in the range $[-S,S]$, find the number of tuple $(i,j,k)$ such that

- $a_i+a_j=a_k$
- $i,j,k$ are distinct

with some limitations

- $1\leq N\leq 2\times10^5$
- $S=5\times10^4$

First, shift $b_i=a_i+S$ so that they are non-negtive. Let $M=2S$ and
$$
p(x)=c_0+c_1x+\ldots+c_{M}x^{M}
$$ where $c_i$ is the number of $i$ in $\{b_i\}$. Now define
$$
q(x)=p^2(x)=d_0+d_1x+\ldots+d_{2M}x^{2M}
$$
then observe that
$$
\begin{split}
d_{b_k+S}=&|\{(i,j)\mid b_i+b_j=b_k+S\}|\\
=&|\{(i,j)\mid a_i+a_j=a_k\}|
\end{split}
$$
and thus the solution is $\sum_{i=0}^{N-1}d_{b_i+S}$ excluding

1. $a_i+a_j=a_k$, but $a_i=0$ and $j=k$
2. $a_i+a_j=a_k$, but $a_j=0$ and $i=k$
3. $a_i+a_j=a_k$, but $i=j$

function solve(N, a):
    S = 5*10^4, M = 2*S, Z = 0
    for i from 0 to N-1:
        if a[i] == 0:
            Z++
        a[i] += S
        b[a[i]]++
    
    B = fft(2*M, padding(b, 2*M, 0))
    C = ifft(2*M, B.*B)
    
    for i from 0 to N-1:
        C[2*b[i]]--
    
    ans = 0
    for i from 0 to N-1:
        ans += C[b[i]+S]
    return ans-2*Z*(n-1)