Time Limit: 2000/1000 MS (Java/Others)    Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 9413    Accepted Submission(s): 1468


Problem Description
Calculate A * B.
 

 

Input
Each line will contain two integers A and B. Process to end of file.

Note: the length of each integer will not exceed 50000.
 

 

Output
For each case, output A * B in one line.
 

 

Sample Input
1 2 1000 2
 

 

Sample Output
2 2000
 

 

Author
DOOM III
 

 

Recommend
DOOM III
 

就一个高精度乘法 FFT加速。

最近正好要捡起fft,就顺便整理了模板。

FFT的原理还是算法导论靠谱,没有那么艰深难懂,就涉及怎么进行FFT和FFT需要的原理和定理。

看看算法导论里FFT的部分,一定要读到迭代实现那部分!!

看了好久求和引理,才发觉他是为了保证$w_n^k$与$w_n^{k+2/h}$的对称性(即$w_n^{k+2/h}=-w_n^k$)的,这个引理是必要的。

对于多项式序列,我们可以用两个O(nlgn)(n>max(len1,len2)*2)的FFT将其系数表示转化为点值表示(DFT),然后用O(n) 相乘,接着用FFT把结果的点值表示变为系数表示(IDFT),总体算起来是3O(nlgn)+O(n),即O(nlgn)的时间复杂度。比O(n^2)好多了。

以下是学习的两个版本。

  1 #include<bits/stdc++.h>
  2 #define clr(x) memset(x,0,sizeof(x))
  3 #define clr_1(x) memset(x,-1,sizeof(x))
  4 #define clrmax(x) memset(x,0x3f3f3f3f,sizeof(x))
  5 #define LL long long
  6 #define mod 1000000007
  7 #define PI 3.1415926535
  8 using namespace std;
  9 char s1[200010],s2[200010];
 10 int a[200010],b[200010];
 11 //复数序列结构体
 12 struct complexed
 13 {
 14         double r,i;
 15         complexed(double _r=0.0,double _i=0.0)
 16         {
 17             r=_r;
 18             i=_i;
 19         }
 20         complexed operator +(complexed b)
 21         {
 22             return complexed(r+b.r,i+b.i);
 23         }
 24         complexed operator -(complexed b)
 25         {
 26             return complexed(r-b.r,i-b.i);
 27         }
 28         complexed operator *(complexed b)
 29         {
 30             return complexed(r*b.r-i*b.i,i*b.r+r*b.i);
 31         }
 32 }num1,num2;
 33 vector<complexed> multi1,multi2;
 34 inline int max(int a,int b)
 35 {
 36     return a>b?a:b;
 37 }
 38 //并将长度变为2…^(k+1)
 39 void changelen(int &len)
 40 {
 41     int mul=1;
 42     while(mul<len)
 43         mul<<=1;
 44     mul<<=1;
 45     len=mul;
 46     return ;
 47 }
 48 //将整数序列复制到复数序列中
 49 void copyed(int *a,vector<complexed> &multi,int len)
 50 {
 51     multi.resize(len);
 52     for(int i=0;i<len;i++)
 53         multi[i]=(complexed){a[i],0};
 54     return;
 55 }
 56 //DFT的话on=1,IDFT on=-1;
 57 void fft(vector<complexed> &multi,int len,int on)
 58 {
 59     complexed wn,w,u,t;
 60     //wn,w,u,t如算法导论中所示
 61     vector<complexed> ans;
 62     ans.resize(len);
 63     //ans存每次操作计算后的y,最后再作为下次的multi。
 64     for(int h=len/2;h>=1;h>>=1)
 65     {
 66         wn=(complexed){cos(2*on*PI/(len/h)),sin(2*on*PI/(len/h))};
 67         for(int i=0;i<h;i++)
 68         {
 69             w=(complexed){1,0};
 70             for(int j=0;j<len/h/2;j++)
 71             {
 72                 //蝴蝶操作
 73                 u=multi[i+2*h*j];
 74                 t=multi[i+2*h*j+h]*w;
 75                 ans[i+h*j]=u+t;
 76                 ans[i+h*j+len/2]=u-t;
 77                 w=w*wn;
 78             }
 79         }
 80         //ans作为下次计算的multi
 81         multi=ans;
 82     }
 83     //IDFT每个元素都得除以n
 84     if(on==-1)
 85         for(int i=0;i<len;i++)
 86             multi[i].r/=len;
 87     return ;
 88 }
 89 int main()
 90 {
 91     int len1,len2,len;
 92     while(scanf("%s%s",s1,s2)!=EOF)
 93     {
 94         len1=strlen(s1);
 95         len2=strlen(s2);
 96         clr(a);
 97         clr(b);
 98         for(int i=0;i<len1;i++)
 99         {
100             a[len1-i-1]=s1[i]-'0';
101         }
102         for(int i=0;i<len2;i++)
103         {
104             b[len2-i-1]=s2[i]-'0';
105         }
106         len=max(len1,len2);
107         //取长度较长者作为长度,并将长度变为2…^(k+1)
108         changelen(len);
109         //将两个整数序列复制到复数序列中
110         copyed(a,multi1,len);
111         copyed(b,multi2,len);
112         //对两个复数序列进行DFT,变为点值表示
113         fft(multi1,len,1);
114         fft(multi2,len,1);
115         //对应点点值相乘
116         for(int i=0;i<len;i++)
117             multi1[i]=multi1[i]*multi2[i];
118         //将的出来的点值表示进行IDFT变为系数表示
119         fft(multi1,len,-1);
120         //四舍五入减小损失精度
121         for(int i=0;i<len;i++)
122         {
123             a[i]=(int)(multi1[i].r+0.5);
124         }
125         //进位
126         for(int i=0;i<len;i++)
127         {
128             a[i+1]=a[i+1]+a[i]/10;
129             a[i]%=10;
130         }
131         len=len1+len2-1;
132         //去掉前导0
133         while(a[len]<=0 && len>0) len--;
134         for(int i=len;i>=0;i--)
135             printf("%d",a[i]);
136         printf("\n");
137     }
138     return 0;
139 }
无位逆序置换的步长实现

相关文章:

  • 2021-11-04
  • 2021-06-23
  • 2021-07-30
  • 2021-05-18
  • 2021-09-28
  • 2021-06-17
  • 2022-12-23
  • 2022-12-23
猜你喜欢
  • 2022-12-23
  • 2021-06-18
  • 2021-10-02
  • 2021-09-14
  • 2022-12-23
  • 2022-03-07
  • 2022-12-23
相关资源
相似解决方案