2179: FFT快速傅立叶 - #include <Seter> - using namespace Orz;

1861: [Zjoi2006]Book 书架

2179: FFT快速傅立叶

Seter posted @ 2011年11月22日 19:30 in BZOJ with tags HighPrecision NumberTheory Divide , 3270 阅读

http://www.zybbs.org/JudgeOnline/problem.php?id=2179

RunID	User	Problem	Result	Memory	Time	Language	Code Length	Submit Time
175346	Seter	2179	Accepted	1756 kb	1120 ms	C/Edit	1597 B	2011-11-22 19:17:39

于是花了3个小时写了个分治版的大数乘法……压9位后常数果然巨小，60000位时比FFT还快一点！膜拜AC大神208MS稳居榜首！

用FFT过掉后看到TIM是用分治的……于是回去想了好久……终于想出来了……

多项式(Ax+B)(Cx+D)=ACxx+BD+(ADx+BCx)=ACxx+BD+[(A+B)(C+D)-AC-BD]x

令x=10^(n/2)，于是每次算出AC,BD,(A+B)(C+D)即可。T(n)=3T(n/2)+O(n)

当n较小时直接暴力出解……测了下n<=10时效果最好，压3位的话60效果最好。

2179
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38	`#include <stdio.h>` `#include <string.h>` `#define logn 10` `#define maxn 6670` `#define base 1000000000` `int` `re[maxn<<1],a[maxn],b[maxn],ABCD[logn][maxn],A[logn][maxn],B[logn][maxn],N;` `char` `str[120003];` `multi(re,a,b,n,d)int` `re,a,b,n,d;{` `int` `i,j=0,ln=n>>1,rn=n-ln,l;long` `long` `t;` `memset(re,0,n<<3);` `if(n<=10){` `for(i=-1;++i!=n;)for(j=-1;++j!=n;re[i+j]=t-lbase)re[i+j+1]+=l=(t=(long` `long)a[i]b[j]+re[i+j])/base;` `for(l=n<<1;--l&&!re[l];);return` `l;` `}memset(A[d],0,rn+1<<2);memset(B[d],0,rn+1<<2);` `for(i=-1;++i!=ln;){` `A[d][i+1]+=(A[d][i]+=a[i]+a[i+ln])/base;A[d][i]%=base;` `B[d][i+1]+=(B[d][i]+=b[i]+b[i+ln])/base;B[d][i]%=base;` `}if(rn==ln+1){` `A[d][rn]+=(A[d][i]+=a[i<<1])/base;A[d][i]%=base;` `B[d][rn]+=(B[d][i]+=b[i<<1])/base;B[d][i]%=base;` `}l=multi(ABCD[d],A[d],B[d],rn+1,d+1);` `for(i=multi(B[d],a,b,ln,d+1)+1;i--;re[i]=B[d][i])ABCD[d][i]-=B[d][i];` `for(i=multi(A[d],a+ln,b+ln,rn,d+1)+1;i--;re[i+(ln<<1)]+=A[d][i])ABCD[d][i]-=A[d][i];` `for(i=-1;i++!=l+1;){` `ABCD[d][i+1]+=ABCD[d][i]/base;` `if((ABCD[d][i]%=base)<0)ABCD[d][i]+=base,ABCD[d][i+1]--;` `if(ABCD[d][i])re[i+ln]+=ABCD[d][i];` `}for(i=-1;++i!=n<<1;(re[i]%=base)&&(j=i))re[i+1]+=re[i]/base;` `return` `j;` `}main(){` `int` `i,p,k,d;` `scanf("%d",&N);d=N%9?9-N%9:0;fread(str,1,(N<<1)+2,stdin);` `for(k=(N+d)/9-1,p=d,i=N;i--;(++p==9)&&(p=0,k--)){` `a[k]=a[k]10+str[N-i]-48;` `b[k]=b[k]*10+str[(N<<1)-i+1]-48;` `}for(N=(N+d)/9,printf("%d",re[i=multi(re,a,b,N,0)]);i--;printf("%09d",re[i]));putchar('\n');` `return` `0;` `}`