BZOJ 2179: FFT快速傅立叶

2179: FFT快速傅立叶

Time Limit: 10 Sec  Memory Limit: 259 MB
Submit: 3425  Solved: 1795
[Submit][Status][Discuss]

Description

给出两个n位10进制整数x和y，你需要计算x*y。

Input

第一行一个正整数n。 第二行描述一个位数为n的正整数x。 第三行描述一个位数为n的正整数y。

Output

输出一行，即x*y的结果。

Sample Input

1
3
4

Sample Output

12

数据范围：
n<=60000

FFT模板题。

直接进行高精度乘法是 O(n2) 的，于是我们采用FFT来 O(nlogn) 实现： 
（注明：以下着重介绍算法流程和算法思想，具体细节参考《算法导论》）

1.我们把乘数的每一位看作多项式的系数，得到多项式 A(x) （因为高精度乘法的本质就是多项式乘法）

2.首先求出 A(ωkn) ，其中 k∈[0,n−1] ， ωn 是n次单位复根。 
由于n次单位复根的一些奇妙性质： 
相消引理 
 
折半引理 
 
我们可以采用分治 O(nlogn) 的时间求出这 n 项的值，但是递归实现常数较大，我们采用蝴蝶算法来迭代实现。 
 
如图，把原来顺次排列的数列变成叶子中的顺序就可以迭代了~ 
（叶子中的顺序就是原序列的二进制逆序）

3.（这一步叫插值） 
如下图， 
ai 表示多项式的系数； 
ωkn 就是我们带入多项式的值： 
第一个多项式 x=ω0n ，第二个是 x=ω1n ，第三个是 x=ω2n ….知道第n个是 x=ωn−1n ； 
yi 就是带入不同的 x 的求出的多项式的值。 

那么，通过上一步已经求出了两个多项式（两个乘数）的 yi ，现在我们把对应的 yi 乘起来就是乘积的多项式的 Yi ，我们要求的是这个多项式的系数即 ai ，因此我们只要求出 ωkn 矩阵的逆矩阵即可。

最后推出 ai 的表达式 
 
和第二步要求的式子几乎一样~再来一次FFT即可解决~ 
时间复杂度依然是 O(nlogn)

转载出处：http://blog.csdn.net/regina8023/article/details/44908753

附代码模板

#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<complex>
#include<iostream>
#include<algorithm>
#include<iomanip>
#include<vector>
#include<string>
#include<queue>
#include<set>
#include<map>
#define pi acos(-1)
using namespace std;
inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
	return f*x;
}
const int N=1<<17;
complex<double> a[N],b[N],t[N];
char s[N];int c[N],r[N];
void FFT(complex<double> x[],int n,int p)
{
	for(int i=0;i<n;i++)if(i<r[i])swap(x[i],x[r[i]]);
	for(int m=2;m<=n;m<<=1)
	{
		complex<double> wn(cos(2*p*pi/m),sin(2*p*pi/m));
		for(int i=0;i<n;i+=m)
		{
			complex<double> w(1,0),temp;int k=m>>1;//不知道为什么这个在这里定义跑得快
			for(int j=0;j<k;j++,w*=wn)
			{
				temp=x[i+j+k]*w;
				x[i+j+k]=x[i+j]-temp;
				x[i+j]=x[i+j]+temp;
			}
		}
	}
}
int main()
{
	int n=read();
	scanf("%s",s+1);
	for(int i=0;i<n;i++)a[i]=s[n-i]-'0';
	scanf("%s",s+1);
	for(int i=0;i<n;i++)b[i]=s[n-i]-'0';
	for(int i=1,j=n;(i>>2)<j;i<<=1)n=i;
	for(int i=0;i<n;i++) r[i]=(i&1)*(n>>1)+(r[i>>1]>>1); 
	FFT(a,n,1);FFT(b,n,1);
	for(int i=0;i<n;i++)t[i]=a[i]*b[i];
	FFT(t,n,-1);
	for(int i=0;i<n;i++)c[i]=t[i].real()/n+0.1;
	int len=0;
	for (int i=0;i<n;i++)if (c[i])len=i,c[i+1]+=c[i]/10,c[i]%=10;
    for (int i=len;i>=0;i--)printf("%d",c[i]);puts("");
    return 0;
}
/*
1
3
4

*/

2017.9.3更新模板

采用重载复数，快很很多

#include<cmath>
#include<ctime>
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<complex>
#include<iostream>
#include<algorithm>
#include<iomanip>
#include<vector>
#include<string>
#include<bitset>
#include<queue>
#include<map>
#include<set>

using namespace std;

typedef long long ll;
typedef double db;

inline int read()
{
	int x=0,f=1;char ch=getchar();
	while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
	while(ch<='9'&&ch>='0'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
	return x*f;
}
void print(int x)
{if(x<0)putchar('-'),x=-x;if(x>=10)print(x/10);putchar(x%10+'0');}

#define pi acos(-1)

const int N=200100;

struct cp
{
	db real,image;
	friend cp operator -(const cp &x,const cp &y){return (cp){x.real-y.real,x.image-y.image};}
	friend cp operator +(const cp &x,const cp &y){return (cp){x.real+y.real,x.image+y.image};}
	friend cp operator *(const cp &x,const cp &y){return (cp){x.real*y.real-x.image*y.image,x.real*y.image+x.image*y.real};}
}a[N],b[N],c[N];

int r[N],ans[N];

void fft(cp x[],int n,int p)
{
	register int i,j,k,m;
	for(i=0;i<n;++i)if(r[i]<i)swap(x[i],x[r[i]]);
	for(m=2;m<=n;m<<=1)
	{
		cp wn=(cp){cos(2*pi*p/m),sin(2*pi*p/m)};
		for(i=0;i<n;i+=m)
		{
			cp w=(cp){1,0},tmp;k=m>>1;
			for(j=0;j<k;++j,w=w*wn)
			{
				tmp=x[i+j+k]*w;
				x[i+j+k]=x[i+j]-tmp;
				x[i+j]=x[i+j]+tmp;
			}
		}
	}
}

char s[N];

int main()
{
	int n=read();
	register int i,j,len=0;
	scanf("%s",s);
	for(i=0;i<n;++i)a[i].real=s[n-i-1]-'0';
	scanf("%s",s);
	for(i=0;i<n;++i)b[i].real=s[n-i-1]-'0';
	for(i=1,j=n;(i>>2)<j;i<<=1)n=i;
	for(i=0;i<n;++i)r[i]=((i&1)*(n>>1))+(r[i>>1]>>1);
	fft(a,n,1);fft(b,n,1);
	for(i=0;i<n;++i)c[i]=a[i]*b[i];
	fft(c,n,-1);
	for(i=0;i<n;++i)ans[i]=c[i].real/n+0.1;
	for(i=0;i<n;++i)if(ans[i]){len=i;ans[i+1]+=ans[i]/10;ans[i]%=10;}
	for(i=len;~i;i--)putchar(ans[i]+'0');puts("");
	return 0;
}