HDU 1402 A * B Problem Plus FFT
大数乘法,暴力逐项乘后相加 O(n^2) 或者 用JAVA的BigInteger都超时
这时候就要用到FFT了...多项式乘法,大数乘法都能在O(nlogn)时间内解决
因为之前没有复变和数字信号课程的基础...所以今天看FFT异常的吃力...
一个上午都在研究复数的概念看FFT的定义,下午的时候才开始看FFT的算法...晚上就学了一个别人的板...过了这个水题
发现FFT完全就是模板,直接会用就行了...
但是这里还是讲解一下原理:
首先这个课件应该是最能解决概念性问题的: http://wenku.baidu.com/view/14457119de80d4d8d15a4f34.html
然后就是kuangbin神推荐的这个课件: http://wenku.baidu.com/view/8bfb0bd476a20029bd642d85.html
这题对于FFT,我们一共就要做三件事 -> 对a,b求值(做DFT) , 两个值点乘 , 对结果插值(IDFT)
其中值得一提的也就是求值插值的过程了...(插值是求值的逆过程,也就是对Vandermonde矩阵求逆后的DFT,具体的算导上写的很详细...结论在在515页上面的30.11式,和30.08式很像对吧 :>
递归法的求值在算导513页,迭代法是在517页,其中递归法我暂时还没去敲这个板...我先讲解一下楼下要贴的迭代的板
首先一个难点是反转置换,这个太神奇了...我无法解释为什么反转下标的二进制能够达到这个效果...举个图中的例子 4是100 在图中的位置是001 , 6是110 在图中的位置是011
通过这个反转操作, 我们就可以自底向上进行算式合并了....下面这个图很直观的说明了反转的重要性
这个板很逆天的给了一个我看不懂的过程...具体实现如下
1 void rev(complex *y,int len) // logn 2 { 3 // 对向量反转置换 见算导517页 4 register int i,j,k; 5 for(i = 1 , j = len / 2 ; i < len - 1 ; i++) 6 { 7 if( i < j ) swap(y[i] , y[j]); 8 k = len / 2; 9 while(j >= k) 10 { 11 j -= k; 12 k /= 2; 13 } 14 if(j < k) j += k; 15 } 16 }
然后就是求值的过程了...
也就两个算式进行一下蝴蝶操作 ... (名字很屌...其实很好理解就是交叉加减
前n/2个项的算式:
后n/2个项的算式:
你可以发现前后两个算式除了符号不同以外其他都一样...
我们将其中共同的算式提取出来,这也就是俗称的蝴蝶操作...
关于单位复根和Vandermonde矩阵我就不做多余的解释了...可以参见算导或者上面的课件
以下是代码
1 void fft(complex *y , int len , double op) // nlogn 2 { 3 // a为向量, h为向量长度, OP=1时求值 OP=-1时插值 4 register int i,j,k,h; 5 complex u,t; 6 rev(y,len); // 对向量反转置换 7 for(h = 2;h <= len; h <<= 1) // 控制层数 8 { 9 complex wn(cos(op*2*PI/h),sin(op*2*PI/h));//该层单位复根,h为层数 10 for(j = 0 ; j < len ; j += h) // 控制起始下标 11 { 12 complex w(1,0); // 初始化螺旋因子 13 for(k = j ; k < j + h / 2 ; k++) // 配对 14 { 15 u = y[k]; 16 t = w * y[k+h/2]; 17 y[k] = u + t; 18 y[k + h / 2] = u - t; 19 w = w * wn; // 更新螺旋因子 20 } 21 } 22 } 23 if(op == -1) 24 for(i = 0 ; i < len ; i++) 25 y[i].real /= len; // IDFT 26 }
不是自己的代码看起来还是不爽,今天再推掉重写一遍好了...
HDU1402:
1 struct complex 2 { 3 double real,imag; 4 bool vis; 5 complex(double a = 0.,double b = 0.){ 6 real = a; 7 imag = b; 8 } 9 complex operator + (complex e){ 10 return complex(real+e.real,imag+e.imag); 11 } 12 complex operator - (complex e){ 13 return complex(real-e.real,imag-e.imag); 14 } 15 complex operator * (complex e){ 16 return complex(real*e.real-imag*e.imag,real*e.imag+imag*e.real); 17 } 18 }x1[MAXN],x2[MAXN]; 19 char a[MAXN],b[MAXN]; 20 int sum[MAXN]; 21 void rev(complex *y,int len) // logn 22 { 23 // 对向量反转置换 见算导517页 24 register int i,j,k; 25 for(i = 1 , j = len / 2 ; i < len - 1 ; i++) 26 { 27 if( i < j ) swap(y[i] , y[j]); 28 k = len / 2; 29 while(j >= k) 30 { 31 j -= k; 32 k /= 2; 33 } 34 if(j < k) j += k; 35 } 36 } 37 void fft(complex *y , int len , double op) // nlogn 38 { 39 // a为向量, h为向量长度, OP=1时求值 OP=-1时插值 40 register int i,j,k,h; 41 complex u,t; 42 rev(y,len); // 对向量反转置换 43 for(h = 2;h <= len; h <<= 1) // 控制层数 44 { 45 complex wn(cos(op*2*PI/h),sin(op*2*PI/h));//该层单位复根,h为层数 46 for(j = 0 ; j < len ; j += h) // 控制起始下标 47 { 48 complex w(1,0); // 初始化螺旋因子 49 for(k = j ; k < j + h / 2 ; k++) // 配对 50 { 51 u = y[k]; 52 t = w * y[k+h/2]; 53 y[k] = u + t; 54 y[k + h / 2] = u - t; 55 w = w * wn; // 更新螺旋因子 56 } 57 } 58 } 59 if(op == -1) 60 for(i = 0 ; i < len ; i++) 61 y[i].real /= len; // IDFT 62 } 63 64 int main() 65 { 66 //freopen("in.txt","r",stdin); 67 //CHEAT; 68 while(~scanf("%s%s",a,b)) 69 { 70 memset(sum,0,sizeof(sum)); 71 int i ; 72 int len1 = strlen(a); 73 int len2 = strlen(b); 74 int len = 1 ; 75 while(len < len1 * 2 || len < len2 * 2) len <<= 1; //求最大次界 76 //cout<<len<<endl; 77 //倒置 将多余次数界初始化为0 78 for(i = 0 ; i < len1 ; i++) 79 { 80 x1[i].real = a[len1-i-1] - '0'; 81 x1[i].imag = 0.0; 82 } 83 for(; i < len ; i++) 84 { 85 x1[i].real = 0.0; 86 x1[i].imag = 0.0; 87 } 88 for(i = 0 ; i < len2 ; i++) 89 { 90 x2[i].real = b[len2-i-1] - '0'; 91 x2[i].imag = 0.0; 92 } 93 for(; i < len ; i++) 94 { 95 x2[i].real = 0.0; 96 x2[i].imag = 0.0; 97 } 98 //对x1,x2求值 99 fft(x1,len,1); 100 fft(x2,len,1); 101 for(i = 0 ; i < len ; i++) x1[i] = x1[i] * x2[i]; // 点乘 102 //对乘积插值 103 fft(x1,len,-1); 104 for(i = 0 ; i < len ; i++) 105 sum[i] = x1[i].real + 0.5; // 四舍五入 106 for(i = 0 ; i < len ; i++) // 进位 107 { 108 sum[i+1] += sum[i] / 10; 109 sum[i] %= 10; 110 } 111 len = len1 + len2 - 1; 112 while(sum[len] <= 0 && len > 0) len--; // 除前导0 113 for(i = len ; i >= 0 ; i--) printf("%c",sum[i]+'0'); 114 puts(""); 115 } 116 return 0; 117 }