zoukankan      html  css  js  c++  java
  • FFT

    W(nk,N)=e^(-j2pi/Nnk)//蝶形因子
    W(nk,N)=W(-nk,N)=W((N-n)k,N)=W(n(N-k),N) 相当于W(-nk,N)+e(-j*2npi)或+e(-j
    2kpi)对称性
    W(nk,N)=W((N+n)k,N)=W(n(N+k),N)周期性
    W(nk,N)=W(nmk,mN)=W(n/mk,N/m) = e(-j2pi/N*nk)=e(-j2mpi/Nmnk) 可约性
    W(0,N)=1 = e^0 = 1 W(N/2,N)=-1=e-jpi W(k+N/2,N)=-W(k,N) 特殊点
    X(k)=Y(k)+W(k,N)
    Z(k)//因为E(W((2n+1)k,N))=W(k,N)W(20+1)k,N)+W(k,N)W(21+1)k,N)+...=W(k,N)E(W(2nk,N)=W(k,N)E(W(nk,N/2))
    X(k+2/N)=Y(k)-W(k,N)Z(k)//这里Y(k)和Z(k)都有N/2的周期又W(k+2/N,N)=-W(k,N)
    一直进行奇偶分解到2点DFT,一共log(2,N)级,每级蝶形运算N/2,每次蝶形运算复数乘法1次,复数加法2次
    比直接DFT运算量N
    N降低至N/2*log(2,N)
    然后就是进行位交换了,比如第3个为011,位逆序为110交换到第6个处
    添加一个bin神链接:http://www.cnblogs.com/kuangbin/archive/2013/07/24/3210389.html
    其实这个3层还是没太懂

    #include <bits/stdc++.h>
    using namespace std;
    const int N = 200010;
    int a[N];
    map<int,int> p;
    const double pi = acos(-1.0);
    class complex {
    public:
    	double r, i;
    	complex(){}
    	complex(double r, double i) {
    		this->r = r;
    		this->i = i;
    	}
    	complex operator +(const complex &rhs) {
    		return complex(r + rhs.r, i + rhs.i);
    	}
    	complex operator -(const complex &rhs) {
    		return complex(r - rhs.r, i - rhs.i);
    	}
    	complex operator *(const complex &rhs) {
    		return complex(r*rhs.r - i*rhs.i, r*rhs.i + i*rhs.r);
    	}
    };
    //改变输入序列,通过2进制逆序值
    void change(complex y[], int len) {
    	int x = 0;
    	while ((1 << x) < len)x++;
    	for (int i = 0;i < len;i++) {
    		int j = i;
    		int k = 0;
    		while (k < x/2) {
    			if ((j&(1 << k))&&!(j&(1<<(x-k-1)))){
    				j = j - (1 << k) + (1 << (x - k - 1));
    			}
    			else if (!(j&(1 << k)) && (j&(1 << (x - k - 1)))) {
    				j = j + (1 << k) - (1 << (x - k - 1));
    			}
    			k++;
    		}
    		if (j > i)
    			swap(y[i], y[j]);
    	}
    }
    void fft(complex y[], int len, int on)
    {
    	change(y, len);
    	/*complex Y[1000];
    	for (int i = 0;i < len;i++)
    		Y[i].r = y[i].r;*/
    	/*for (int k = 0;k < len;k++)
    	{
    		complex sc=complex(0, 0);
    		for (int j = 0;j < len;j++) {
    			complex c = complex(Y[j].r * cos(on*2 * pi*k*j / len), -Y[j].r * sin(on*2 * pi*k*j / len));
    			sc = sc+c;
    		}
    		y[k] = sc;
    	}*/
    	for (int l = 2;l <= len;l <<= 1) {
    		complex wn = complex(cos(-on * 2 * pi / l), sin(-on * 2 * pi / l));
    		for (int j = 0;j < len;j += l) {
    			complex wnp = complex(1, 0);
    			for (int k = j;k < j + l / 2;k++) {
    				complex u = y[k];
    				complex t = wnp*y[k + l / 2];
    				y[k] = u + t;
    				y[k + l / 2] = u - t;
    				wnp = wnp*wn;
    			}
    		}
    	}
    	//for (int h = 2; h <= len; h <<= 1)
    	//{
    		/*complex wn(cos(-on * 2 * pi / len), sin(-on * 2 * pi / len));
    		for (int j = 0;j < len;j ++)
    		{
    			complex w(1, 0);
    			for (int k = 0;k < len;k++)
    			{
    				complex u = y[k];
    				complex t = w*y[k + len / 2];
    				y[k] = u + t;
    				y[k + len / 2] = u - t;
    				w = w*wn;
    			}
    		}*/
    	//}
    	if (on == -1)
    		for (int i = 0;i < len;i++)
    			y[i].r /= len;
    }
    char s1[N], s2[N];
    complex a1[N], a2[N];
    int ans[N],tmp[N];
    int main()
    {
    	while (~scanf("%s%s", s1, s2)) {
    		int len1 = strlen(s1);
    		int len2 = strlen(s2);
    		int tlen = len1 + len2;
    		int len = 1;
    		while (len < tlen)len *= 2;//长度>=len1+len2-1的圆周卷积才等于线性卷积,又用的是基2-fft所以长度为2^x
    		for (int i = 0;i < len;i++) {
    			if (i < len1)a1[i].r = s1[len1-1-i] - '0';
    			else a1[i].r = 0;
    			a1[i].i = 0;
    		}//补齐0
    		for (int i = 0;i < len;i++) {
    			if (i < len2)a2[i].r = s2[len2-1-i] - '0';
    			else a2[i].r = 0;
    			a2[i].i = 0;
    		}//补齐0
    		//由于输出序列为顺序的所以采取抽取时间的fft方法
    		fft(a1, len, 1);//fft
    		fft(a2, len, 1);//fft
    		for (int i = 0;i < len;i++)
    			a1[i] = a1[i]*a2[i];
    			//Y(k)=H(k)*X(k)
    		//ifft可以用fft子程序实现,先求一次共轭,再fft,序列集体除以len,再求一次共轭,由于虚数部分为0,基本不用求共轭
    		fft(a1, len, -1);//ifft
    		for (int i = 0;i < len;i++) {
    			ans[i] = (int)(a1[i].r + 0.5);
    		}
    		for (int i = 1;i <= len;i++) {
    			ans[i] += ans[i - 1] / 10;
    			ans[i - 1] %= 10;
    		}
    		int x = 0;
    		len = len1 + len2 - 1;
    		for (;len>0;len--)
    			if (ans[len] != 0)break;
    		for (;len>=0;len--)
    			printf("%d", ans[len]);
    		printf("
    ");
    	}
    	return 0;
    }
    
    
  • 相关阅读:
    flex + bison multiple parsers
    Educational Codeforces Round 95 (Rated for Div. 2)
    python学习笔记 day20 序列化模块(二)
    python学习笔记 day20 常用模块(六)
    python 学习笔记 常用模块(五)
    python学习笔记 day19 常用模块(四)
    python学习笔记 day19 常用模块(三)
    python学习笔记 day19 常用模块(二)
    python学习笔记 day19 作业讲解-使用正则表达式实现计算器
    python学习笔记 day19 常用模块
  • 原文地址:https://www.cnblogs.com/HaibaraAi/p/5043872.html
Copyright © 2011-2022 走看看