BZOJ-2179 FFT快速傅立叶

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <cctype>
#include <complex>
#define rep(i, l, r) for(int i=l; i<=r; i++)
#define down(i, l, r) for(int i=l; i>=r; i--)
#define cd complex <double>
#define PI acos(0.0)*2.0
#define ll long long
#define base 10000
#define maxn 50009
using namespace std;
inline int read()
{
	int x=0, f=1; char ch=getchar();
	while (!isdigit(ch)) {if (ch==‘-‘) f=-1; ch=getchar();}
	while (isdigit(ch)) x=x*10+ch-‘0‘, ch=getchar();
	return x*f;
}
cd a[maxn], b[maxn], c[maxn], A[maxn];
ll ans[maxn];
int n, m, len, n1[maxn], n2[maxn], s[maxn], ansn=0;
char num[maxn];
void fft(cd *a, bool flag)
{
	rep(i, 0, n-1) s[i]=0;
	for(int i=1, j=n; i<n; i*=2, j/=2) rep(h, j/2, j-1) s[h]+=i;
	for(int i=1; i<n; i*=2) rep(j, 0, i-1) s[j+i]+=s[j];
	rep(i, 0, n-1) A[i]=a[s[i]];
	double pi=flag?PI:-PI;
	for(int step=1; step<n; step*=2)
	{
		cd e=exp(cd(0, 2.0*pi/double(step*2))), w=cd(1, 0);
		for(int pos=0; pos<step; ++pos, w*=e) 
			for(int i=pos; i<n; i+=step*2)
			{
				cd ret=A[i], rec=w*A[i+step];
				A[i]=ret+rec, A[i+step]=ret-rec;
			}
	}
	if (!flag) rep(i, 0, n-1) A[i]/=n;
	rep(i, 0, n-1) a[i]=A[i];
}
int main()
{
	m=read(); len=(m-1)/4+1;
	scanf("%s", num);
	rep(i, 0, m-1) n1[(m-1-i)/4]=n1[(m-1-i)/4]*10+num[i]-‘0‘;
	scanf("%s", num);
	rep(i, 0, m-1) n2[(m-1-i)/4]=n2[(m-1-i)/4]*10+num[i]-‘0‘;
	n=1; m=len*2; while (n<m) n*=2;
	rep(i, 0, n-1) a[i]=cd(n1[i], 0); fft(a, true);
	rep(i, 0, n-1) b[i]=cd(n2[i], 0); fft(b, true);
	rep(i, 0, n-1) c[i]=a[i]*b[i]; fft(c, false);
	rep(i, 0, n-1) ans[i]=(ll)(c[i].real()+0.5);
	rep(i, 0, n-1) ans[i+1]+=ans[i]/base, ans[i]%=base;
	ansn=n; while (ansn>0 && !ans[ansn]) ansn--;
	printf("%lld", ans[ansn]);
	down(i, ansn-1, 0)
		if (ans[i]>=1000) printf("%lld", ans[i]);
		else if (ans[i]>=100) printf("0%lld", ans[i]);
		else if (ans[i]>=10) printf("00%lld", ans[i]);
		else printf("000%lld", ans[i]);
	printf("\n");
	return 0;
}