1. 程式人生 > >[洛谷P1919 A*B Problem] [FFT板子題]

[洛谷P1919 A*B Problem] [FFT板子題]

FFT就不多說了吧,反正已經有寫得很好的blog了,強烈推薦去看menci的FFT學習筆記,寫得十分詳細,條理也很清晰,我就是看這個才懂FFT的!!(傳送門:menci FFT學習筆記

這裡只是想提幾個注意點。

[1] 蝴蝶操作需謹慎,就是將每個a[i]對映到對應位置的時候,由於二進位制的位置從0開始標號,容易搞錯;

[2] 自定義的operator運算最好加括號,防止優先順序的問題

[3] 最後a[i]/n應該改成(a[i]+0.5)/n,因為由於精度誤差,可能a[i]/n出現12.9999/13這種情況,導致整除的時候結果比原來小。

[4] 中途迭代的時候下標不能搞錯,我一開始就把a[i+l]寫成a[i]了。

FFT不容易除錯,能讀程式找錯最好,不能的話只好出幾個樣例來看了,要細心一點。

#include <cstdio>
#include <cmath>
#include <cstring>
#define db double
#define rep(i,j,k) for (i=j;i<=k;i++)
#define down(i,j,k) for (i=j;i>=k;i--)
using namespace std;
const int N=12e5+5;
const db PI=acos(-1.0);
char sa[N],sb[N];
struct cpl{
	db a,b;
	cpl () { a=b=0; }
	cpl (db aa,db bb) { a=aa; b=bb; }
	cpl operator + (const cpl &B) {	return cpl(a+B.a,b+B.b); }
	cpl operator - (const cpl &B) {	return cpl(a-B.a,b-B.b); }
	cpl operator * (const cpl &B) {	return cpl(a*B.a-b*B.b,a*B.b+b*B.a); }
}A[N],B[N],C[N];
struct large{
	int n,d[N];
	cpl a[N];
	int pos(int x,int K)
	{
		int i,ret=0;
		rep(i,0,K-1)
			if (x&(1<<i)) ret+=(1<<(K-i-1));
		return ret;
	}
	void FFT(cpl *A,int K,int sgn) //K-positions
	{
		int i,l,m,k,len=1<<K;
		cpl bsc,w,L,R;
		rep(i,0,len-1) a[pos(i,K)]=A[i]; 
		for (k=2;k<=len;k<<=1)
		{
			m=k/2; bsc=cpl(cos(2*PI/k),sgn*sin(2*PI/k));
			for (l=0;l<len;l+=k)
			{
				w=cpl(1,0);
				rep(i,0,m-1)
				{
					L=a[l+i]; R=a[l+m+i];
					a[l+i]=L+(w*R); //ÓÅÏȼ¶!! 
					a[l+i+m]=L-(w*R);
					w=w*bsc;
				}
			}
		}
		rep(i,0,len-1)
			if (sgn==1) A[i]=a[i];
			else A[i]=cpl((a[i].a+0.5)/len,0); //·ÀÖ¹³öÏÖ12.99999ÕâÑùµÄ¾«¶ÈÎó²î,Òª+0.5 
	}
	void mult (const large &Y,large &ret)
	{
		int i,k,len=n+Y.n-1;
		for (k=0;(1<<k)<len;k++);
		len=1<<k; //0~len-1
		rep(i,0,len-1) A[i]=cpl(d[i],0),B[i]=cpl(Y.d[i],0);
		FFT(A,k,1);
		FFT(B,k,1);
		rep(i,0,len-1) C[i]=A[i]*B[i];
		FFT(C,k,-1); ret.n=len;
		rep(i,0,len-1) ret.d[i]=C[i].a;
		rep(i,0,ret.n-1)
		{
			if (ret.d[i]>9)
			{
				if (i==ret.n-1) ret.n++;
				ret.d[i+1]+=ret.d[i]/10;
				ret.d[i]%=10;
			}
		}
		for (;ret.n>1 && !ret.d[ret.n-1];ret.n--);
	}
}numA,numB,ans;
int n,i;
int main()
{
	scanf("%d",&n); numA.n=numB.n=n;
	scanf("%s%s",sa,sb);
	rep(i,0,n-1) numA.d[i]=sa[n-i-1]-'0',numB.d[i]=sb[n-i-1]-'0';
	numA.mult(numB,ans);
	down(i,ans.n-1,0) printf("%d",ans.d[i]);
	return 0;
}