zoukankan      html  css  js  c++  java
  • FFTNTT总结

    学了好久,终于基本弄明白了

    推荐两个博客:
    戳我
    戳我
    再推荐几本书:
    《ACM/ICPC算法基础训练教程》
    《组合数学》(清华大学出版社)
    《高中数学选修》

    预备知识

    复数方面

    找数学老师去

    [i^{2}=-1,i为虚数的单位 ]

    坐标系上纵轴就是虚数轴,复数就是这上面的点
    三种表示法:
    $$一般:a + bi,a为实部,b为虚部$$
    $$指数:e^{i heta}坐标系上的模长$$
    $$三角:模长
    (cos heta + i sin heta)$$
    运算:
    加减法:实部虚部分别相加
    乘法:$$(a + bi) * (c + di) = ac + adi + bci + bdi^{2}
    = ac-bd+(ad+bc)i$$

    欧拉公式

    [e^{ix} = cosx + isinx(就是指数表示和三角表示) ]

    [特别的e^{ipi} = -1 ]

    多项式

    [系数表示法:A(x) = Sigma _{k=0}^{n - 1} a_kx^k ]

    [点值表示法:对于所有的x_k,求出它们对应的A(x),设为y_k ]

    [则可以用{(x_0, y_0), (x_1, y_1), ......, (x_n-1, y_n-1)} 表示这个多项式 并且是唯一确定的]

    单位复数根

    [n次单位复数根omega^{n} = 1,n次单位复数根刚好有n个对应e^{frac{2kpi i}{n}},其中k=0到n - 1 ]

    三个性质:
    消去引理:
    $$n, d, k为正整数,则omega^{dk}{dn}=omega^{k}{n}$$
    $$证明:套e^{frac{2kpi i}{n}} 即可$$
    折半引理:
    $$n为大于零的偶数,则(omega^{k+frac{n}{2}}{n})^{2}=omega^{2k+n}{n}=omega^{2k}{n}omega^{n}{n}=(omega^{k}{n})^{2}$$
    求和引理:
    大于1的整数n,和不被n整除的非负整数k,有
    $$Sigma^{n-1}
    {j=0}(omega^{k}_{n})^{j}=0$$
    证明可以用等比数列求和公式得到(很简单的,手推一遍就好)

    Rader排序

    其实就是二进制数位翻转

    正题

    DFT

    对于k=0~n-1,定义:

    [y_k=A(omega^{k}_{n}) = Sigma^{n-1}_{j=0} a_j(omega^{k}_{n})^j ]

    [得到的y称为a的离散傅里叶变换,记作y=DFT_n(a) (这里的y,a指的是所有的y_k, a_k,即向量y,a) ]

    逆DFT

    [就是DFT的逆变换,求出向量a,记为DFT^{-1} ]

    假设得到了向量y

    [对于y_k = Sigma^{n-1}_{i=0}a_i(omega^{k})^i ]

    [有a_k = frac{1}{n}Sigma^{n-1}_{i=0}y_i(omega^{-k})^i ]

    [证明:a_k=frac{1}{n}Sigma^{n-1}_{i=0}y_i(omega^{-k})^i=frac{1}{n}Sigma^{n-1}_{i=0}( Sigma^{n-1}_{j=0}a_j(omega^{k})^j)(omega^{-k})^i=frac{1}{n} Sigma^{n-1}_{i=0}a_i(Sigma^{n-1}_{j=0}(omega^{j-k})^i) ]

    [可以用等比数列求和出上面的就是a_k(当j e k是括号里的为0,当j=k时为1) ]

    FFT

    上面已经把DFT和逆DFT搞定了,两个几乎是一样的

    所以求多项式的积(卷积)可以用DFT转换成点值表示,就可以O(n),一一相乘,得到积的多项式的点值表示,最后用逆DFT得到系数表示

    复杂度瓶颈在于怎样快速求解DFT(逆DFT和DFT方法一样)

    FFT就是一个O(nlogn)求解DFT的方法

    首先把A(x)分成奇数项和偶数项记作

    [A^{[0]}(x) = a_0 + a_2x + a_4x^2 + ... + a_{n-2}x^{frac{n}{2} - 1} ]

    [A^{[1]}(x) = a_1 + a_3x + a_5x^2 + ... + a_{n-1}x^{frac{n}{2} - 1} ]

    [显然A(x) = A^{[0]}(x^2) + xA^{[1]}(x^2) ]

    那么

    [A(omega^k_n)=A^{[0]}((omega^k_n)^2) + omega^k_n A^{[1]}((omega^k_n)^2)=A^{[0]}(omega^k_{frac{n}{2}}) + omega^k_n A^{[1]}(omega^k_{frac{n}{2}}) ]

    [因为omega^{frac{n}{2}}_{n}=omega_{2}=e^{kpi i}= cos kpi + i sin kpi = -1 ]

    [所以A(omega^{k+frac{n}{2}}_n)=A^{[0]}(omega^k_{frac{n}{2}}) - omega^k_n A^{[1]}(omega^k_{frac{n}{2}}) ]

    这称为蝴蝶操作

    于是对每个y值的求解可以通过分组求出,若递归变成处理子任务,这样复杂度就成了O(nlogn)

    这样不停地分组,最后就相当于Rader排序了一番,所以也可以变成非递归的

    注意每次都要把多项式补成2的幂,便于FFT

    递归写可能好理解一些,但不好写

    还有一些东西什么的,其实记一记就好了其实自己说不清

    系统的复数complex代码

    # include <bits/stdc++.h>
    # define RG register
    # define IL inline
    # define Fill(a, b) memset(a, b, sizeof(a))
    using namespace std;
    typedef long long ll;
    const int _(3e6 + 10);
    const double Pi = acos(-1);
    
    IL ll Read(){
    	char c = '%'; ll x = 0, z = 1;
    	for(; c > '9' || c < '0'; c = getchar()) if(c == '-') z = -1;
    	for(; c >= '0' && c <= '9'; c = getchar()) x = x * 10 + c - '0';
    	return x * z;
    }
    
    int n, m, r[_], l;
    complex <double> a[_], b[_];
    
    IL void FFT(complex <double> *P, int opt){
    	for(RG int i = 0; i < n; ++i) if(i < r[i]) swap(P[i], P[r[i]]); //Rader排序
    	for(RG int i = 1; i < n; i <<= 1){
    		complex <double> W(cos(Pi / i), opt * sin(Pi / i));  //旋转因子
    		for(RG int p = i << 1, j = 0; j < n; j += p){
    			complex <double> w(1, 0);
    			for(RG int k = 0; k < i; ++k, w *= W){
    				complex <double> X = P[j + k], Y = w * P[j + k + i];
    				P[j + k] = X + Y; P[j + k + i] = X - Y;    //蝴蝶操作
    			}
    		}
    	}
    }
    
    int main(RG int argc, RG char *argv[]){
    	n = Read(); m = Read();
    	for(RG int i = 0; i <= n; ++i) a[i] = Read();
    	for(RG int i = 0; i <= m; ++i) b[i] = Read();
    	m += n;
    	for(n = 1; n <= m; n <<= 1) ++l;//补成2的幂
    	for(RG int i = 0; i < n; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));//Rader排序预处理
    	FFT(a, 1); FFT(b, 1); //DFT
    	for(RG int i = 0; i < n; ++i) a[i] = a[i] * b[i]; //点值直接相乘
    	FFT(a, -1);  //逆DFT
    	for(RG int i = 0; i <= m; ++i) printf("%d ", (int)(a[i].real() / n + 0.5));
    	return 0;
    }
    
    

    或者可以自己定义complex,用复数运算

    struct Complex{
    	double real, image;
    	IL Complex(){  real = image = 0;  }
    	IL Complex(RG double a, RG double b){  real = a; image = b;  }
    	IL Complex operator +(RG Complex B){  return Complex(real + B.real, image + B.image);  }
    	IL Complex operator -(RG Complex B){  return Complex(real - B.real, image - B.image);  }
    	IL Complex operator *(RG Complex B){  return Complex(real * B.real - image * B.image, real * B.image + image * B.real);  }
    }
    

    NTT(快速数论变换)

    前置技能原根
    (g)(p)(质数)的原根
    (e^{frac{2pi i}{n}}equivomega_nequiv g^{frac{p-1}{n}}(mod p))
    带进去就好了
    Reverse的那个不会证明
    (UOJ)的模板

    # include <bits/stdc++.h>
    # define RG register
    # define IL inline
    # define Fill(a, b) memset(a, b, sizeof(a))
    using namespace std;
    typedef long long ll;
    const int Zsy(998244353);
    const int Phi(998244352);
    const int G(3);
    const int _(4e5 + 5);
    
    IL ll Input(){
        RG ll x = 0, z = 1; RG char c = getchar();
        for(; c < '0' || c > '9'; c = getchar()) z = c == '-' ? -1 : 1;
        for(; c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) + (c ^ 48);
        return x * z;
    }
    
    int n, m, N, l, r[_], A[_], B[_];
    
    IL int Pow(RG ll x, RG ll y){
    	RG ll ret = 1;
    	for(; y; y >>= 1, x = x * x % Zsy)
    		if(y & 1) ret = ret * x % Zsy;
    	return ret;
    }
    
    IL void NTT(RG int *P, RG int opt){
    	for(RG int i = 0; i < N; ++i) if(r[i] < i) swap(P[r[i]], P[i]);
    	for(RG int i = 1; i < N; i <<= 1){
    		RG int W = Pow(G, Phi / (i << 1));
    		if(opt == -1) W = Pow(W, Zsy - 2);
    		for(RG int j = 0, p = i << 1; j < N; j += p){
    			RG int w = 1;
    			for(RG int k = 0; k < i; ++k, w = 1LL * w * W % Zsy){
    				RG int X = P[k + j], Y = 1LL * w * P[k + j + i] % Zsy;
    				P[k + j] = (X + Y) % Zsy, P[k + j + i] = (X - Y + Zsy) % Zsy;
    			}
    		}
    	}
    }
    
    int main(RG int argc, RG char* argv[]){
        n = Input(), m = Input();
    	for(RG int i = 0; i <= n; ++i) A[i] = Input();
    	for(RG int i = 0; i <= m; ++i) B[i] = Input();
    	for(n += m, N = 1; N <= n; N <<= 1) ++l;
    	for(RG int i = 0; i < N; ++i) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
    	NTT(A, 1); NTT(B, 1);
    	for(RG int i = 0; i < N; ++i) A[i] = 1LL * A[i] * B[i] % Zsy;
    	NTT(A, -1);
    	RG int inv = Pow(N, Zsy - 2);
    	for(RG int i = 0; i <= n; ++i) printf("%lld ", 1LL * A[i] * inv % Zsy);
        return 0;
    }
    
    
  • 相关阅读:
    hlgoj 1766 Cubing
    Reverse Linked List
    String to Integer
    Bitwise AND of Numbers Range
    Best Time to Buy and Sell Stock III
    First Missing Positive
    Permutation Sequence
    Next Permutation
    Gray Code
    Number of Islands
  • 原文地址:https://www.cnblogs.com/cjoieryl/p/8206721.html
Copyright © 2011-2022 走看看