多项式基础1 FFT NTT FWT

本着不在NOIP前碰多项式的原则，我现在开始学FFT。

多项式的表示

系数表示：$f(x)={\sum }_{i=0}^n{a}_i{x}^i$，则$a$是$f$的系数表示。

点值表示：一个长为$n+1$的数组，其中元素为$(d_i,f(d_i))$，$d$是元素互不相同的数组。

可以使用DFT和IDFT在$O(n\log n)$复杂度内将系数表示与点值表示互相转换。

$k$次单位根

$k$次单位根是其$k$次幂为$1$的一些复数，共有$k$个，分别为$e^{\frac{2\pi ij}{k}}(j\in [0,k-1])$。记${\omega }_k=e^{\frac{2\pi i}{k}}$，则$k$次单位根可分别表示为$w_k^i(i\in [0,k-1])$。

在复平面上，$k$次单位根是单位圆上均匀分布的$k$个点，其中第$0$个在$(1,0)$，接下来的依次按逆时针绕一圈。容易想到单位根有以下性质：

• ${\omega }_k^i={\omega }_{\frac{k}{2}}^{\frac{i}{2}}$
• ${\omega }_k^i=-{\omega }_k^{i+\frac{k}{2}}$
• ${\omega }_k^i=\cos \frac{i}{k}2\pi +i\cdot \sin \frac{i}{k}2\pi$

DFT

用于将系数表示转化为点值表示。

先将$f(x)$的阶向上扩充为$2^k-1$（多余部分的$a_i=0$）。设此时$f(x)$的阶为$n-1$，考虑将$n$次单位根作为$d$数组，即$d_i={\omega }_n^i$。

设$f_0(x)=a_0{x}^0+a_2{x}^1+\cdots +a_{n-1}{x}^{\frac{n}{2}}$，$f_1(x)=a_1{x}^0+a_3{x}^1+\cdots +a_{n-2}{x}^{\frac{n}{2}}$。显然$f(x)=f_0(x^2)+x{f}_1(x^2)$。

于是$f({\omega }_n^i)=f_0({\omega }_n^{2i})+{\omega }_n^i{f}_1({\omega }_n^{2i})=f_0({\omega }_{\frac{n}{2}}^i)+{\omega }_n^i{f}_1({\omega }_{\frac{n}{2}}^i)$。

同理可以得到$f({\omega }_n^{i+\frac{n}{2}})=f(-{\omega }_n^i)=f_0({\omega }_{\frac{n}{2}}^i)-{\omega }_n^i{f}_1({\omega }_{\frac{n}{2}}^i)$。

看到没？$f_0$和$f_1$中带入的值是一样的，因此我们可以分治递归解决问题。复杂度为$O(n\log n)$。

IDFT

用于将点值表示转化为系数表示。该过程等价于用${\omega }_n^{i+\frac{n}{2}}$替换${\omega }_n^i$做DFT，再将结果除以$n$。复杂度为$O(n\log n)$。

DFT和IDFT统称为FFT（Fa Fa Tower）。实际实现中，应用雷德（Rader）算法来优化常数，用迭代模拟递归的过程。即先枚举层，再枚举每个“计算块”。数组下标需要进行变换成如 { 0 , 4 , 2 , 6 , 1 , 5 , 3 , 7 } \{0,4,2,6,1,5,3,7\} {0,4,2,6,1,5,3,7}这样。

const double pi = 3.14159265358979323846264338327950288419716939937510582097494459230781640628620899;
struct complex {
    double r, v;
    complex() {}
    complex(double rr, double vv) {r = rr; v = vv;} 
};
inline complex operator + (complex a, complex b) {
    return complex(a.r+b.r, a.v+b.v);
}
inline complex operator - (complex a, complex b) {
    return complex(a.r-b.r, a.v-b.v);
}
inline complex operator - (complex x) {
    return complex(-x.r, -x.v);
}
inline complex operator * (complex a, complex b) {
    return complex(a.r*b.r-a.v*b.v, a.r*b.v+a.v*b.r);
}

const int maxn = 1<<17;
void fft(complex a[], int n, bool typ) {
    for(int i = 1, j = n>>1; i < n; i++) {
        if(i < j) swap(a[i], a[j]);
        for(int k = n>>1; (j^=k) < k; k >>= 1);
    }
    for(int i = 1; i < n; i <<= 1) for(int j = 0; j < n; j += i<<1) {
        complex w = complex(1, 0), wn = complex(cos(pi/(double)i), (typ?1:-1)*sin(pi/(double)i));
        for(int k = 0; k < i; k++) {
            complex u = a[j+k], v = w * a[j+i+k];
            a[j+k] = u + v;
            a[j+i+k] = u - v;
            w = w * wn;
        }
    }
    if(!typ) for(int i = 0; i < n; i++) a[i].r /= n;
}
complex tmp[maxn<<1];
void Mul(complex a[], int an, complex b[], int bn) {
    if(an <= 48 || bn <= 48) {
        for(int i = 0; i < an+bn-1; i++) tmp[i] = complex(0, 0);
        for(int i = 0; i < an; i++) for(int j = 0; j < bn; j++) tmp[i+j] = tmp[i+j] + a[i] * b[j];
        for(int i = 0; i < an+bn-1; i++) a[i] = tmp[i];
        return;
    }
    int n = 1;
    while(n < an+bn-1) n <<= 1;
    fft(a, n, 1); fft(b, n, 1);
    for(int i = 0; i < n; i++) a[i] = a[i] * b[i];
    fft(a, n, 0);
}

卷积

即多项式乘法。多项式$f$与$g$的卷积为$h$，设$f$的阶为$n$，$g$的阶为$m$，则$h$的阶为$n+m$。设$f$、$g$的系数数组分别为$a$、$b$，则$h$的系数数组$c$满足$c_i={\sum }_{j=0}^i{a}_j{b}_{i-j}$。也可理解为$a_i$和$b_j$对$c_{i+j}$有$+a_i{b}_j$的贡献。

可以看到朴素的卷积算法复杂度为$O(nm)$。但若先将$f$与$g$转化为点值表示后再相乘，则乘法的复杂度直接降为$O(n+m)$。因此卷积的复杂度等于FFT的复杂度，即$O(n\log n)$。

NTT

用于解决模意义下的系数表示与点值表示相转化。设模数为$P$。

由于${\omega }_n^n\equiv g^{P-1}\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}P)$，因此${\omega }_n^i\equiv g^{\frac{i\cdot (P-1)}{n}}$。将原来的单位根换为原根即可。

const int maxn = 1<<17;
int w[2][1<<19];
void nttinit() {
    for(int i = 0; i <= 18; i++) {
        w[1][1<<i] = w[0][1<<i] = 1;
        int wn = qpow(gmd, md-1>>i+1), invwn = inv(wn);
        for(int j = (1<<i)+1; j < 1<<i+1; j++) {
            w[1][j] = mul(w[1][j-1], wn);
            w[0][j] = mul(w[0][j-1], invwn);
        }
    }
}
void ntt(int a[], int n, bool typ) {
    for(int i = 1, j = n>>1; i < n; i++) {
        if(i < j) swap(a[i], a[j]);
        for(int k = n>>1; (j^=k) < k; k >>= 1);
    }
    for(int i = 1; i < n; i <<= 1) for(int j = 0; j < n; j += i<<1) for(int k = 0; k < i; k++) {
        int u = a[j+k], v = mul(w[typ][i+k], a[j+i+k]);
        a[j+k] = add(u, v);
        a[j+i+k] = sub(u, v);
    }
    if(!typ) {
        int invn = inv(n);
        for(int i = 0; i < n; i++) f[i] = mul(f[i], invn);
    }
}
int tmp[maxn<<1];
void Mul(int a[], int an, int b[], int bn) {
    if(an <= 48 || bn <= 48) {
        memset(tmp, 0, (an+bn-1)*sizeof(int));
        for(int i = 0; i < an; i++) for(int j = 0; j < bn; j++) Add(tmp[i+j], mul(a[i], b[j]));
        memcpy(a, tmp, (an+bn-1)*sizeof(int));
        return;
    }
    int n = 1;
    while(n < an+bn-1) n <<= 1;
    ntt(a, n, 1); ntt(b, n, 1);
    for(int i = 0; i < n; i++) a[i] = mul(a[i], b[i]);
    ntt(a, n, 0);
}

采用vector的版本：

typedef vector<int> pol;

void ntt(pol& f, int typ) {
    int n = f.size();
    for(int i = 1, j = n>>1; i < n; i++) {
        if(i < j) swap(f[i], f[j]);
        for(int k = n>>1; (j ^= k) < k; k >>= 1);
    }
    for(int i = 1; i < n; i <<= 1) for(int j = 0; j < n; j += i<<1) for(int k = 0; k < i; k++) {
        int u = f[j+k], v = mul(w[typ][i+k], f[j+i+k]);
        f[j+k] = add(u, v);
        f[j+i+k] = sub(u, v);
    }
    if(!typ) {
        int invn = inv(n);
        for(int i = 0; i < n; i++) f[i] = mul(f[i], invn);
    }
}

FWT

用于求位运算（集合）卷积。同样我们先对原来的两个数组进行处理，再$O(n)$相乘，再逆变换回去。

或FWT

考虑对数组的一种变换$FWT(A)=B$，满足$B[i]={\sum }_{i|j=i}A[j]$。这个可以通过分治在$O(n\log n)$时间内做到。

将数组的长度向上扩充到$2^k$，设其为$n$。假设当前我们在处理区间$[l,r]$，并且已经计算完了$[l,mid]$与$[mid+1,r]$自己区间内部的影响，考虑合并。我们发现$l+i$这个数及其二进制子集一定是$mid+i+1$这个数的二进制子集，因此对于所有$i\in [0,\frac{r-l+1}{2}]$将$l+i$的答案累加到$mid+i+1$上即可。即$FWT(L:R)=FWT(L):(FWT(L)+FWT(R))$，其中$:$表示连接。不难推出其逆变换为$UFWT(L:R)=UFWT(L):(UFWT(R)-UFWT(L))$。

与FWT

此时需要实现的变换是$B[i]={\sum }_{i\&j=i}A[j]$。

与或运算类似，容易推出$FWT(L:R)=(FWT(L)+FWT(R)):FWT(R)$，$UFWT(L:R)=(UFWT(L)-UFWT(R)):UFWT(R)$。

以上两变换又称FMT（快速莫比乌斯变换），又称高维前缀和。

异或FWT

这个比较奇怪，反正我想不到。此时需要实现的变换是$B[i]={\sum }_{|i\cap j|\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}2)}A[j]-{\sum }_{|i\cap j|\equiv 1(\mathrm{m}\mathrm{o}\mathrm{d}2)}A[j]$。

变换是$FWT(L:R)=(FWT(L)+FWT(R)):(FWT(L)-FWT(R))$，$UFWT(L:R)=\frac{UFWT(L)+UFWT(R)}{2}:\frac{UFWT(L)-UFWT(R)}{2}$（和差问题）。手推可验证正确。
*2018.11~12