FFT前置知识及FFT\FFT分治

首先依次介绍FFT(快速傅里叶变换)前置知识:
复数及单位根
复变函数欧拉定理
单位根三大引理
多项式的系数与点值表示法
多项式的求值与插值(顺介绍拉格朗日插值法)
向量卷积

复数及单位根

我们知道复数可以写成$z=a+bi$,其中$a$为实部,$bi$为虚部.而且一个复数可以在一个复平面上表示出来,而我们可以将其写成极坐标的形式$z=r(cos\theta +isin\theta )$而根据欧拉定理:
$$r{e}^{\theta i}==cos\theta +isin\theta$$
则有:$z=r{e}^{\theta i}$,这里给出乘法的例子:
$$z_1{z}_2==(r_1{e}^{{\theta }_{1}i})(r_2{e}^{{\theta }_{2}i})==r_1{r}_2{e}^{({\theta }_1+{\theta }_2)i}$$
下面介绍单位根:
如果存在复数${\omega }^n==1$,则其根被称作单位根,显然单位根的模为1,其中我们定义单位根:
$${\omega }_n^k==e^{i\frac{2k\pi }{n}}$$
其中定义主n次单位根(当k==1时):
$${\omega }_n==e^{i\frac{2\pi }{n}}$$
而在程序中我们一般写成三角函数形式:
$$cos\frac{2k\pi }{n}+isin\frac{2k\pi }{n}$$
这里给出复数运算模板代码:

struct Complex
{
    double r, i;
    Complex() {}//无参函数不要写东西!易超时
    Complex(double r, double i)
    :r(r), i(i) {}
}f[N], g[N];
Complex operator + (const Complex & a, const Complex & b)
{return Complex(a.r + b.r, a.i + b.i);}
Complex operator - (const Complex & a, const Complex & b)
{return Complex(a.r - b.r, a.i - b.i);}
Complex operator * (const Complex & a, const Complex & b)

下面给出单位根的三大引理
$${\omega }_{dn}^{dk}=={\omega }_n^k\text{\hspace{1em}(1)}$$
$${\omega }_n^{k+\frac{n}{2}}==({\omega }_n^k{)}^2=={\omega }_{\frac{n}{2}}^k\text{\hspace{1em}(2)}$$
$$\sum _{i=0}^{n-1}(w_n^k{)}^i==0\text{\hspace{1em}(3)}$$
多项式的两种表示方法:
第一种非常常见的系数表示法,及$f(x)={\sum }_{i=0}^n{a}_i{x}^i$,另外我们给出另一种表示方法,点值表示法,及用n+1种不同的取值代入$f(x)$,组成的集合:
A ( x ) = { ( x 0 , f ( x 0 ) ) , ( x 1 , f ( x 1 ) ) , . . . . ( x n , f ( x n ) ) } A(x)=\{(x_0,f(x_0)),(x_1,f(x_1)),....(x_n,f(x_n))\} A(x)={(x0​,f(x0​)),(x1​,f(x1​)),....(xn​,f(xn​))}
我们分析两种表示方法:
系数表达式:
求值:我们可以使用秦九昭算法$O(n)$得出,
求和:两个表达式的各项系数相加,是$O(n)$的
求积:一个表达式相乘的各项都需要和另一个的各项相乘,是$O(n^2)$的
点值表达式:
求值:我们使用拉格朗日插值法可以在$O(n^2)$的出系数表达式并求值(插值的定义一会给出)
这里给出拉格朗日插值法的公式(n项表达式):
$$A(x)=\sum _{k=0}^{n-1}{y}_k\frac{{\prod }_{j!=k}(x-x_j)}{{\prod }_{j!=k}(x_k-x_j)}$$
并有模板:
P4781 【模板】拉格朗日插值
下面是模板代码:

const int N = 2e3+5;
const ll mod = 998244353;
ll a[N], b[N];
ll inv(ll a, ll b)
{
    ll res = 1;
    while(b)
    {
        if (b&1) res = res * a % mod;
        b >>= 1;
        a = a * a % mod;
    }
    return res % mod;
}
int main()
{
    int n; ll k;
    cin >> n >> k;
    for (int i = 1; i <= n; i++)
        scanf("%d%d", &a[i], &b[i]);
    ll sum = 0;
    for (int i = 1; i <= n; i++)
    {
        ll ans = b[i];
        ll fa = 1, fb = 1;
        for (int j = 1; j <= n; j++)
        {
            if (j == i) continue;
            fa = (fa * (k - a[j] + mod)%mod)%mod;
            fb = (fb * (a[i] - a[j] + mod)%mod)%mod;
        }
        ans = (ans * fa % mod * inv(fb, mod-2)%mod)%mod;
        sum = (sum + ans)%mod;
    }
    printf("%lld\n", sum);
}

求和:两个表达式的各项点值的$y$值相加,是$O(n)$的
求积:两个表达式的各项点值的$y$值相乘,是$O(n)$的
我们不难发现,如果我们相求两个多项式的积$f,g$(系数向量分别是$a,b$),在点值表达式的情况下将十分高效.这时新表达式$h$(系数向量是c)的第i项系数为:
$$c_i=\sum _{k=0}^i{a}_k{b}_{j-k}$$
我们则称$c_i$为$a_i,b_i$的卷积,记作$a\otimes b$,两个多项式相乘,就是对应的两个系数向量的卷积,
我们将一个点值表达式转化为系数表达式称为插值.
由此可见:如何提高多项式求值与插值的效率.决定着两个表达式相乘的效率!

话讲完,正式下面介绍FFT

首先看DFT(离散傅里叶变换)
我们将单位根带入到系数表达式(方便起见下式中$\omega$用$w$替代)中:
$$y_i=A(w_n^i)=\sum _{j=0}^{n-1}{w}_n^{ij}{a}_j$$
记作$y=DFT(a)$
每次求值为$O(n)$所以离散傅里叶变换最终时间复杂度为$O(n^2)$
这样显然是不行的,毕竟直接相乘的时间复杂度都是$O(n^2)$
然后就有了FFT:
我们对于一个n项系数表达式(n为2的幂数):
$$A(x)\to a=[a_0,a_1.....a_{n-1}]$$
可以分成:
$$A^{[0]}(x)\to a^{[0]}=[a_0,a_2,a_4......a_{n-2}]$$
$$A^{[1]}(x)\to a^{[1]}=[a_1,a_3,a_5......a_{n-1}]$$
这时我们发现一个及其有趣的事实:
$$A(x)=A^{[0]}(x^2)+x{A}^{[1]}(x^2)$$
而根据(1),(2)两个引理:我们最终有:
$$\{\begin{array}{rl}A(w_n^k)& =A^{[0]}(w_{\frac{n}{2}}^k)+w_n^k{A}^{[1]}(w_{\frac{n}{2}}^k)\\ A(w_n^{k+\frac{n}{2}})& =A^{[0]}(w_{\frac{n}{2}}^k)-w_n^k{A}^{[1]}(w_{\frac{n}{2}}^k)\end{array}$$
于是问题得以递归解决:
$T(n)=O(\frac{n}{2})+2T(\frac{n}{2})$
时间复杂度为$O(nlogn)$

逆FFT

我们的多项式不可以FFT过去弄不过来了!比较简单,
因为:
$$y_k=\sum _{j=0}^{n-1}{a}_j{w}_n^{kj}$$
所以:
$$a_j=\frac{1}{n}\sum _{j=0}^{n-1}{y}_k{w}_n^{-kj}$$
证明略,总之就是证明原向量矩阵可逆
所以我们稍加改动就可以逆FFT,就是将单位根的指数取相反数,最后结果再乘1/n
所以两个多项式$f,g$的乘积,使其对应系数向量$a,b$的卷积!可以有:
$$a\otimes b=DF{T}_{2n}^{-1}(DF{T}_{2n}(a)DF{T}_{2n}(b))$$

不过我们很少书写递归代码,这里使用迭代写法;
这里还牵扯到两个点:位逆序置换 蝴蝶操作,不说了…直接上代码吧
P3803 【模板】多项式乘法（FFT）

const int N = 3e6+5;
const double pi = acos(-1);
struct Complex
{
    double r, i;
    Complex() {}
    Complex(double r, double i)
    :r(r), i(i) {}
}f[N], g[N];
Complex operator + (const Complex & a, const Complex & b)
{return Complex(a.r + b.r, a.i + b.i);}
Complex operator - (const Complex & a, const Complex & b)
{return Complex(a.r - b.r, a.i - b.i);}
Complex operator * (const Complex & a, const Complex & b)
{return Complex(a.r * b.r - a.i * b.i, a.i * b.r + a.r * b.i);}
int rev[N];
int len, lim = 1;
void init(int n)
{
    while(lim <= n) lim <<= 1, len++;
    for (int i = 0; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len-1));
}
void FFT(Complex *a, int op)//op为1是FFT运算,-1是FFT逆运算
{
    for (int i = 0; i < lim; i++)
    if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int dep = 1; dep <= log2(lim); dep++)
    {
        int m = 1 << dep;
        Complex wn = Complex(cos(2.0 * pi / m), op * sin(2.0 * pi / m));
        for (int k = 0; k < lim; k += m)
        {
            Complex w = Complex(1, 0);
            for (int j = 0; j < m / 2; j++)
            {
                Complex t = w * a[k + j + m / 2];
                Complex u = a[k + j];
                a[k + j] = u + t;
                a[k + j + m / 2] = u - t;
                w = w * wn;
            }
        }
    }
    if (op == -1) for (int i = 0; i < lim; i++)
        a[i].r /= lim;
}
int main()
{
    int n, m;
    cin >> n >> m;
    init(n + m);
    for (int i = 0; i <= n; i++) scanf("%lf", &f[i].r);
    for (int i = 0; i <= m; i++) scanf("%lf", &g[i].r);
    FFT(g, 1); FFT(f, 1);
    for (int i = 0; i <= lim; i++) f[i] = f[i] * g[i];
    FFT(f, -1);
    for (int i = 0; i <= n + m; i++)
        printf("%d ", (int)(f[i].r + 0.5));
    return 0;
}

FFT分治
【模板】分治 FFT
根据cdq分治思想,
我们在计算多项式的一个系数的时候，要统计之前的所有影响。我们区间折半递归解决，先解决左区间的数值，然后统计其对右区间的影响，累加到右区间。
洛谷模板需要取模，故模板代码使用NTT，值得注意的是，洛谷模板的模数998244353，符合NTT所要求的素数类型，且存在原根为3。
下面是模板代码：

const int N = 1e6+5;
const int G = 3;
const ll mod = 998244353;
ll f[N], g[N];
ll a[N], b[N];
int rev[N], n;
int len, lim = 1;
void init(int n)
{
    lim = 1; len = 0;
    while(lim <= n) lim <<= 1, len++;
    for (int i = 0; i < lim; i++) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (len-1));
}
ll _pow(ll a, ll b) {
    //....
}
inline void calcrev(int logn)
{
    rev[0]=0;
    for(int i = 1; i < (1 << logn); i++)
        rev[i]=(rev[i>>1]>>1)|((i&1)<<(logn-1));
}
void NTT(ll *a, int led, int op)//op为1是FFT运算,-1是FFT逆运算
{
    int lim = 1 << led;
    //......

}
void cdq_fft(int l, int r, int logn)
{
    
    if (logn <= 0) return ;
    if (l >= n) return ;
    int mid = (l + r) >> 1;
    ll t = _pow(r - l, mod-2);
    cdq_fft(l, mid, logn - 1);
    calcrev(logn);
    memset(a + (r - l) / 2, 0, sizeof(ll) * (r - l) / 2);
    memcpy(a, f+l, sizeof(ll) * (r - l) / 2);
    memcpy(b, g, sizeof(ll) * (r - l) );
    NTT(a, logn, 1); NTT(b, logn, 1);
    for (int i = 0; i  < r - l; i++) a[i] = a[i] * b[i] % mod;
    NTT(a, logn, -1);
    for (int i = (r - l) / 2; i < r - l; i++)
        f[l + i] = (f[l+i] + a[i]) % mod;
    cdq_fft(mid, r, logn-1);
}
int main()
{
    cin >> n;
    init(n);
    for (int i = 1; i < n; i++) scanf("%lld", &g[i]);
    f[0] = 1;
    cdq_fft(0, lim, len);
    for (int i = 0; i < n; i++)
        printf("%lld ", (f[i] + mod) % mod);
    puts("");
    return 0;
}