A*B problem(FFT)
设两个多项式\(A(x)\)和\(B(x)\),它们的系数镜像反转一下,得到的多项式是\(A'(x)\)和\(B'(x)\)。那么\(C(x)=A(x)*B(x)\)和\(C'(x)=A'(x)*B'(x)\)的系数也是镜像反转的。这个,,感性理解一下吧。
于是倒过来搞会很方便。由于是大整数乘法,算完从个位到最高位进一下位就行了。注意数组要开四倍!
#include <cmath>
#include <cctype>
#include <cstdio>
using namespace std;
const int maxn=6e4+5;
const double pi=3.1415926535898;
struct Cpx{
double x, y;
Cpx (double t1=0, double t2=0){ x=t1, y=t2; }
}A[maxn*4], B[maxn*4], C[maxn*4];
Cpx operator +(Cpx &a, Cpx &b){ return Cpx(a.x+b.x, a.y+b.y); }
Cpx operator -(Cpx &a, Cpx &b){ return Cpx(a.x-b.x, a.y-b.y); }
Cpx operator *(Cpx &a, Cpx &b){ return Cpx(a.x*b.x-a.y*b.y, a.x*b.y+a.y*b.x); }
void swap(Cpx &x, Cpx &y){ Cpx t=x; x=y; y=t; }
int n, r[maxn*4], limit=1, l, ans[maxn*4]; //这里也要乘4!
void fdft(Cpx *a, int n, int flag){
for (int i=0; i<n; ++i) if (i<r[i]) swap(a[i], a[r[i]]);
for (int mlen=1; mlen<n; mlen<<=1){
Cpx w1(cos(pi/mlen), flag*sin(pi/mlen)), x, y;
for (int i=0; i<n; i+=(mlen<<1)){ //起点
Cpx w(1, 0);
for (int j=i; j<i+mlen; ++j, w=w*w1){ //遍历区间系数转点值
x=a[j], y=w*a[j+mlen];
a[j]=x+y; a[j+mlen]=x-y; }
}
}
}
int main(){
scanf("%d", &n); //倒着存方便补零 两个n-1次的多项式
char c; while (!isdigit(c=getchar()));
for (int i=0; i<n; ++i, c=getchar()) A[n-i-1].x=c-48;
while (!isdigit(c=getchar()));
for (int i=0; i<n; ++i, c=getchar()) B[n-i-1].x=c-48;
while (limit<=2*n-2) limit<<=1, ++l; //需要2×n-2个点,但是有2×n-1个系数
for (int i=1; i<limit; ++i)
r[i]=(r[i>>1]>>1)+((i&1)<<(l-1)); //1<<r就是倒数第r+1位
fdft(A, limit, 1); fdft(B, limit, 1);
for (int i=0; i<limit; ++i) C[i]=A[i]*B[i];
fdft(C, limit, -1);
for (int i=0; i<limit; ++i) ans[i]+=+C[i].x/limit+0.5;
for (int i=0; i<limit; ++i)
ans[i+1]+=ans[i]/10, ans[i]%=10;
bool flag=false;
for (int i=limit; i>=0; --i){
if (ans[i]>0) flag=true;
if (flag) printf("%d", ans[i]);
}
return 0;
}