学习笔记::fft

上次学fft还是5月份，昨天发现已经忘记怎么推导了，代码也看不懂了，就又学习了一发，大概是看menci的博客

0.fft可以进行多项式乘法，朴素的乘法跟手算一样是O(n^2)，fft可以通过分治做到nlogn

1.点值表示:首先我们平常看见的多项式都是系数表示，类似于a0+a1*x^1+a2*x^2......,然而通过这种形式我们是不可能降低复杂度的，怎么搞都是O(n^2),于是我们换成点值表示。点值表示是什么呢?设多项式A(x),那么我们带进去一个x可以得到一个y，y=A(x),然后我们取b个不同的x,也就得到n个不同的y，这两个n维向量就是点值表示(x0,x1,x2,...,xn-1),(y0,y1,y2,...,yn-1),其实也可以看成一个函数上取了n个不同的点。每个点值表示对应了唯一的多项式。

两个点值表达式可以相乘，而且可以O(n)相乘，就是对应项乘对应项，设两个多项式A(x),B(x),乘出来是C(x),假设点值表达式分别是(1,1),(2,2);(2,1),(2,2);那么乘出来就是(2,1),(2,4)。

所以现在的问题就是如何把一个系数表示的多项式变换为点值表达式，和从系数表示变成点值表示，分别叫多项式的求值和插值，考虑朴素的求值就是带进去一个一个算，复杂度是O(n^2)的，插值朴素也很慢，fft可以将这两个过程通过分治优化到O(nlogn),于是我们可以在O(n+nlogn)的时间复杂度解决多项式乘法

也就是说是这样的过程o（nlogn)求值(dft)->O(n)乘法->O(nlogn)插值(idft)完成多项式乘法

2.求值(dft)与插值(idft):为了分治，我们先倍增到2^k，这样很方便。考虑带进去的n个x具体是什么数，我们用单位复数根。单位复数根是一个虚数，就是在复平面上画一个单位圆，满足w^n=1的w就是单位负数根，也就是在单位圆上和x实轴夹角为2*pi/n的那个角。

这个东西有很多好的性质。现在考虑fft，设要求的多项式为A(x),然后我们按下奇偶分类，A(x)=a0+a1*x^1+a2*x^2+...+an-1*x^n-1,这里n都默认为2^k,那么A1(x)=a0+a2*x^2+...+an-2*x^n-2,A2(x)=a1*x+a3*x^3+...+an-1*x^n-1,那么A(x)=A1(x^2)+x*A2(x^2)

带进去单位复数根，A(w(n,k))=A1(w(n,k)^2)+w(n,k)* A2(w(n,k)^2)，根据单位负数根的性质，w(n,k)=w(n/2,k/2),那么A(w(n,k))=A1(w(2*n,k))+w(n, k) *A2(w(2*n,k)) = A1(w(n/2,k))+w(n,k)* A2(w(n/2,k))

A(w(n,k+n/2))=A1(w(n,k+n/2))+ w(n,k) *A2(w(n,k+n/2)) 因为W(n,n)=1,w(n,k+n/2)^2=w(n,2*k+n)=w(n,2*k)=w(n/2,k),而w(n,k+n/2)=-w(n,k),因为w(n,n/2)=-1，那么A(w(n,k))=A1(w(n/2,k))- w(n,k)*A2(w(n/2,k)),于是枚举k∈[0,n/2)就可以得到[0,n-1)的所有值，也就是说我们知道了A(w(n/2,0))->A(w(n/2,n/2-1))的所有值后就可以用得出A(w(n,0))->A(w(n,n-1)),那么每次的复杂度是T(n)=2*T(n/2)+O(n)=O(nlogn)，解决了求值的问题，至于插值没有搞清楚，就是把-1带入，然后每项除以n，具体没有搞清楚。

然而上面这样是递归的形式，比较慢，我们想可以直接迭代求。

首先我们发现最后每个数最终的位置是原来二进制位置的对称也就是如果01就变成10,也就是2号位和1号位交换，那么我们可以预处理出最终的位置。考虑合并的过程，我们先知道在进行了n=l的合并后，每个数组位置存的是当前这段[i,i+l)的求值后的值，也就是A1和A2的值，那么现在我们把长度乘2合并。

考虑合并的时候的两个式子

B(w(n,k)) = A1(w(n/2,k))+w(n,k)* A2(w(n/2,k))

B(w(n,k+n/2))=A1(w(n/2,k))- w(n,k)*A2(w(n/2,k))

我们要多开一个b保存，其实不用，我们记录x=A1(w(n/2,k)),y=w(n,k)* A2(w(n/2,k)),然后直接A1=x+y,A2=x-y就行了。

板子：记住循环要循环到n，否则最后取不到

#include<bits/stdc++.h>
using namespace std;
#define pi acos(-1)
const int N = 3e5 + 5;
int n1, n2, n, k;
complex<double> a[N], b[N];
void fft(complex<double> *a, int f)
{
    for(int i = 0; i < n; ++i) 
    {
        int t = 0;
        for(int j = 0; j < k; ++j) if(i >> j & 1) t |= 1 << (k - j - 1);
        if(i < t) swap(a[i], a[t]);
    }
    for(int l = 2; l <= n; l <<= 1) 
    {
        int m = l >> 1;
        complex<double> w(cos(pi / m), f * sin(pi / m));
        for(int i = 0; i < n; i += l) 
        {
            complex<double> t(1, 0);
            for(int k = 0; k < m; ++k, t *= w)
            {
                complex<double> x = a[i + k], y = t * a[i + m + k];
                a[i + k] = x + y;
                a[i + m + k] = x - y;
            }
        }
    }    
}
int main()
{
    scanf("%d%d", &n1, &n2);
    for(int i = 0; i <= n1; ++i) 
    {
        int x;
        scanf("%d", &x);
        a[i] = x;
    }
    for(int i = 0; i <= n2; ++i) 
    {
        int x;
        scanf("%d", &x);
        b[i] = x;
    }
    n = 1;
    while(n <= n1 + n2) n <<= 1, ++k;
    fft(a, 1);
    fft(b, 1);
    for(int i = 0; i < n; ++i) a[i] = a[i] * b[i];
    fft(a, -1);
    for(int i = 0; i <= n1 + n2; ++i) printf("%d ", (int)(a[i].real() / n + 0.5));
    return 0;
}

 ntt

取模怎么办呢

我们要用ntt

我们就不能用单位根了，这样不能取模，于是我们用原根，原根和单位跟性质很想。

单位根w(n,k)

1.w(n,n)=1

2.w(n,n/2)=-1->w(n.k+n/2)=-w(n,k)

3.w(n,k)=w(n/2,k/2)

我们要原根有这些性质

定义原根g满足g^0,g^1,...g^n-2%P各不相同

设一个质数为p=n*q+1,其中n=2^k,则g^(n*q)=1(mod p),这是由费马小定理 a^(p-1) = 1 (mod p)

那么变成单位根

单位根w(n,k)相当于把一个圆分成n分取其中k份的值，也就是一个单位是1/n个圆

那么这里我们定义g(n,k)表示一份是g^q，取k个，g(n,n)就是g^n*q

又因为费马小定理所以g(n,n)=1 (mod p)

并且g(n,n)=g(n,n/2)^2所以g(n,n/2)=1或-1.又因为%p各不相同，并且g^0%p=1,所以g(n,n/2)%p=-1,满足单位根性质