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;
}