• 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;
    }
    
    
  • 相关阅读:
    二、云计算openstack共享组件--时间同步服务ntp
    一、云计算openstack介绍
    五、Kvm虚拟机迁移
    四、Kvm虚拟化网络管理
    三、Kvm虚拟化存储管理
    二、kvm虚拟机管理
    一、kvm虚拟化介绍
    九、Linux网络技术管理及进程管理
    园主的码云网站,可以在里面查看园主的练习代码哦
    万能Makefile,前戏做足项目做起来才顺畅。
  • 原文地址:https://www.cnblogs.com/cjoieryl/p/8206721.html
Copyright © 2020-2023  润新知