FFT

1、hihocoder1388（2016 acm 北京网络赛e题）

给出等长的A,B序列，求

input：

2
9
3 0 1 4 1 5 9 2 6
5 3 5 8 9 7 9 3 2
5
1 2 3 4 5
2 3 4 5 1

output：

80

0

解题思路：这题其实也是2个一维卷积的应用。公式可以化简为

此时我们可以把一个串变长一倍，再用另一个串反转求他们的卷积，如

1 2 3 4 5 1 2 3 4 5与2 3 4 5 1的卷积，然后就能看出有一串连续位置的fft值即为结果，求出其中最小值即可，又因为fft的精度问题，直接fft很可能WA，所以需要在

[cpp]  view plain  copy

 

1. for(int i = 0;i < len;i++)  
2.             sum[i]= (int)(x1[i].r+0.5);  

这种地方直接比较x1[i].r的大小，然后得出什么位置最小，重新按数组计算即可。

代码：

[cpp]  view plain  copy

 

1. #include <cstdio>  
2. #include <cstring>  
3. #include <iostream>  
4. #include <algorithm>  
5. #include <cmath>  
6. using namespace std;  
7.   
8. typedef long long ll;  
9.   
10. const double PI = acos(-1.0);  
11. //复数结构体  
12. struct complex  
13. {  
14.     double r,i;  
15.     complex(double _r = 0.0,double _i = 0.0)  
16.     {  
17.         r = _r; i = _i;  
18.     }  
19.     complex operator +(const complex &b)  
20.     {  
21.         return complex(r+b.r,i+b.i);  
22.     }  
23.     complex operator -(const complex &b)  
24.     {  
25.         return complex(r-b.r,i-b.i);  
26.     }  
27.     complex operator *(const complex &b)  
28.     {  
29.         return complex(r*b.r-i*b.i,r*b.i+i*b.r);  
30.     }  
31. };  
32. /* 
33.  * 进行FFT和IFFT前的反转变换。 
34.  * 位置i和 （i二进制反转后位置）互换 
35.  * len必须去2的幂 
36.  */  
37. void change(complex y[],int len)  
38. {  
39.     int i,j,k;  
40.     for(i = 1, j = len/2;i < len-1; i++)  
41.     {  
42.         if(i < j)swap(y[i],y[j]);  
43.         //交换互为小标反转的元素，i<j保证交换一次  
44.         //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的  
45.         k = len/2;  
46.         while( j >= k)  
47.         {  
48.             j -= k;  
49.             k /= 2;  
50.         }  
51.         if(j < k) j += k;  
52.     }  
53. }  
54. /* 
55.  * 做FFT 
56.  * len必须为2^k形式， 
57.  * on==1时是DFT，on==-1时是IDFT 
58.  */  
59. void fft(complex y[],int len,int on)  
60. {  
61.     change(y,len);  
62.     for(int h = 2; h <= len; h <<= 1)  
63.     {  
64.         complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));  
65.         for(int j = 0;j < len;j+=h)  
66.         {  
67.             complex w(1,0);  
68.             for(int k = j;k < j+h/2;k++)  
69.             {  
70.                 complex u = y[k];  
71.                 complex t = w*y[k+h/2];  
72.                 y[k] = u+t;  
73.                 y[k+h/2] = u-t;  
74.                 w = w*wn;  
75.             }  
76.         }  
77.     }  
78.     if(on == -1)  
79.         for(int i = 0;i < len;i++)  
80.             y[i].r /= len;  
81. }  
82. const int maxn = 400010;  
83. complex x1[maxn],x2[maxn];  
84. ll a[maxn],b[maxn],num[maxn],sum;  
85.   
86. int main()  
87. {  
88.     int T,i,j,len,n,len1,len2;  
89.     scanf("%d",&T);  
90.     while (T--)  
91.     {  
92.         scanf("%d",&n);  
93.         sum=0;  
94.         len1=2*n;  
95.         len2=n;  
96.         for (i=0;i<n;i++)  
97.         {  
98.             scanf("%lld",&a[n-i-1]);  
99.             a[2*n-i-1]=a[n-i-1];  
100.             sum+=a[n-i-1]*a[n-i-1];  
101.         }  
102.         for (i=0;i<n;i++)  
103.         {  
104.             scanf("%lld",&b[n-i-1]);  
105.             sum+=b[n-i-1]*b[n-i-1];  
106.         }  
107.   
108.   
109.         len=1;  
110.         while (len<len1*2) len<<=1;  
111.         for (i=0;i<len1;i++)  
112.         {  
113.             x1[i]=complex(a[i],0);  
114.         }  
115.         for (i=len1;i<len;i++)  
116.         {  
117.             x1[i]=complex(0,0);  
118.         }  
119.         for (i=0;i<len2;i++)  
120.         {  
121.             x2[i]=complex(b[len2-i-1],0);  
122.         }  
123.         for (i=len2;i<len;i++)  
124.         {  
125.             x2[i]=complex(0,0);  
126.         }  
127.         fft(x1,len,1);  
128.         fft(x2,len,1);  
129.         for(int i = 0;i < len;i++)  
130.             x1[i] = x1[i]*x2[i];  
131.         fft(x1,len,-1);  
132.         ll ans=0;  
133.         double temp=-1;  
134.         int pos;  
135.         for(i = n-1;i < 2*n;i++)  
136.         {  
137.             if (temp<x1[i].r)  
138.             {  
139.                 temp=x1[i].r;  
140.                 pos=i;  
141.             }  
142.         }  
143.         int t=0;  
144.         for (i=pos;i>pos-n;i--)  
145.         {  
146.             ans+=a[i]*b[--len2];  
147.         }  
148.         printf("%lld\n",sum-2*ans);  
149.   
150.   
151.     }  
152.     return 0;  
153. }