zoukankan      html  css  js  c++  java
  • FFT 的一些技巧

    三次变两次 FFT

    我们发现:

    [(F(x)+iG(x))^2=F(x)^2-G(x)^2+2iF(x)G(x) ]

    也就是说,我们把 (F(x)) 作为实部,(G(x)) 作为虚部,那么它的平方的虚部的 (1/2) 就是 (F(x)G(x))

    可惜精度比较低。

    四次 FFT 求任意模数多项式乘法

    假设我们要求 (M(x) imes N(x)pmod{p}),因为如果我们直接 FFT 就会爆 double ,所以我们可以把 (M(x)) 拆成 (kA(x)+B(x))(N(x)) 拆成 (kC(x)+D(x)),其中 (kapprox sqrt{p}),那么,你的值域就大约变为 (p imes n) 的了,但是你就需要 (7) 次 FFT 了。

    我们假设有 (Q(x)=A(x)+iB(x),E(x)=C(x)+iD(x)),设 (q(x_1)) 表示 (Q(x))(x_1) 的取值,(A_j) 表示 (A(x))(j) 项系数,那么我们就有:

    [q(omega_n^x)=sum_{j=0} A_jomega_n^{xj}+iB_jomega_n^{xj} ]

    [=sum_{j=0} A_j(cos{frac{2pi xj}{n}}+isin{frac{2pi xj}{m}})+B_j(icos{frac{2pi xj}{n}}-sin{frac{2pi xj}{n}}) ]

    [q(omega_n^{-x})=sum_{j=0} A_j(cos{frac{2pi xj}{n}}-isin{frac{2pi xj}{n}})+B_j(icos{frac{2pi xj}{n}}+sin{frac{2pi xj}{m}}) ]

    [a(omega_n^{x})=sum_{j=0} A_j(cos{frac{2pi xj}{n}}+isin{frac{2pi xj}{n}}) ]

    [b(omega_n^{x})=sum_{j=0} B_j(cos{frac{2pi xj}{n}}+isin{frac{2pi xj}{n}}) ]

    我们假设 (q(x).r) 表示它的实部,(q(x).f) 表示它的虚部,那么我们就可以得到:

    [a(omega_n^{x}).r=frac{q(omega_n^x).r+q(omega_n^{-x}).r}{2},a(omega_n^{x}).f=frac{q(omega_n^x).f-q(omega_n^{-x}).f}{2} ]

    [b(omega_n^{x}).r=frac{q(omega_n^{-x}).f+q(omega_n^x).f}{2},a(omega_n^{x}).f=frac{q(omega_n^{-x}).r-q(omega_n^x).r}{2} ]

    然后,我们求出 (A(x)E(x))(B(x)E(x)),我们就可以得到 (A(x)C(x),A(x)D(x),B(x)C(x),B(x)D(x)),然后算就好了。

    暂时没有懂 (3.5) 次 FFT 的做法,所以不写了。

    Code

    #include <bits/stdc++.h>
    using namespace std;
    
    #define double long double
    #define Int register int
    #define MAXN 270005
    
    template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
    template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
    template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
    
    struct Complex{
    	double x,y;
    	Complex(){}
    	Complex (double _x,double _y){x = _x,y = _y;}
    	Complex operator / (const double &p)const{return Complex{x / p,y / p};}
    	Complex operator + (const Complex &p)const{return Complex{x + p.x,y + p.y};}
    	Complex operator - (const Complex &p)const{return Complex{x - p.x,y - p.y};}
    	Complex operator * (const Complex &p)const{return Complex{x * p.x - y * p.y,x * p.y + p.x * y};}
    };
    
    #define pi (double)acos(-1)
    int l,lim,rev[MAXN];
    void fft (Complex *a,int type){
    	for (Int i = 0;i < lim;++ i) if (i < rev[i]) swap (a[i],a[rev[i]]);
    	for (Int i = 1;i < lim;i <<= 1){
    		Complex Wn(cos(pi / i),type * sin(pi / i));
    		for (Int j = 0,r = i << 1;j < lim;j += r){
    			Complex w(1,0);
    			for (Int k = 0;k < i;++ k,w = w * Wn){
    				Complex x = a[j + k],y = w * a[i + j + k];
    				a[j + k] = x + y,a[i + j + k] = x - y;
    			}
    		}
    	}  
    	if (type == -1) for (Int i = 0;i < lim;++ i) a[i] = a[i] / lim;
    }
    
    int n,m,mod,a[MAXN],b[MAXN],ans[MAXN];
    Complex Q[MAXN],E[MAXN],C[MAXN],D[MAXN];
    
    #define ll long long
    signed main(){
    	read (n,m,mod),lim = 1;
    	while (lim < n + m) lim <<= 1,++ l;
    	for (Int i = 0;i < lim;++ i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l - 1);int up = (1 << 15) - 1;
    	for (Int i = 0;i <= n;++ i) read (a[i]),Q[i] = Complex (a[i] >> 15,a[i] & up);
    	for (Int i = 0;i <= m;++ i) read (b[i]),E[i] = Complex (b[i] >> 15,b[i] & up);
    	fft (Q,1),fft (E,1);
    	for (Int i = 0;i < lim;++ i){
    		int re = (lim - 1) & (lim - i);
    		C[i] = Complex((Q[i].x + Q[re].x) / 2,(Q[i].y - Q[re].y) / 2) * E[i];
    		D[i] = Complex((Q[re].y + Q[i].y) / 2,(Q[re].x - Q[i].x) / 2) * E[i];
    	}
    	fft (C,-1),fft (D,-1);
    	for (Int i = 0;i < lim;++ i){
    		ll v1 = (ll)(C[i].x + 0.5) % mod,v2 = (ll)(C[i].y + D[i].x + 0.5) % mod,v3 = (ll)(D[i].y + 0.5) % mod;
    		ans[i] = ((v1 << 30) + (v2 << 15) + v3) % mod;
    	}
    	for (Int i = 0;i <= n + m;++ i) write ((ans[i] % mod + mod) % mod),putchar (' ');
    	putchar ('
    ');
    	return 0;
    }
    
  • 相关阅读:
    Layui表格之动态添加数据(表格多选解决方案)
    Ztree勾选节点后取消勾选其父子节点
    如何判断某个元素在页面上是否存在
    Linux环境搭建SVN服务
    Linux环境下验证码不显示F12报500
    Ztree节点前加上两个自定义按钮
    解决window.location.href参数太长
    Ztree加载完成默认选中根节点右侧生成表格
    输入框点击下滑Ztree菜单
    Mybatis中and和or的细节处理
  • 原文地址:https://www.cnblogs.com/Dark-Romance/p/14532043.html
Copyright © 2011-2022 走看看