HDU 1402 FFT 求 大数乘法
这题的数据量是5w, 也就是传统意义上的n^2算法是不可取的。这里就用到了FFT FFT一般的作用就是使得多项式乘法的复杂度降到nlogn。利用FFT可以快速求出循环卷积。 那么卷积又是什么样一个东西。 ----------------------------------------以下内容转自http://blog.sina.com.cn/s/blog_6733026501019ubf.html-------------------- 信号处理中的一个重要运算是卷积.初学卷积的时候,往往是在连续的情形, ----------------------------------------------------------完毕------------------------------ 然后我们就知道卷积大概的作用了。 那么FFT本来是信号里面的东西,而我没学过信号。 所以看的也不怎么懂。 大概就是对离散的信号,先将其转变为一些正弦函数,然后这些正弦函数叠加能构成这个离散信号,但是这些正弦函数易于处理。处理完之后就可以再转变回来。 两个过程叫做DFT和IDFT。 对于本道题。意义就很明显了。 可以把两个大整数相乘看做是多项式乘法。 最后求出各系数后再进位即可 代码如下、 #include <iostream> #include <cstdio> #include <algorithm> #include <cstring> #include <cmath> #include <map> #include <queue> #include <set> #include <vector> using namespace std; #define L(x) (1 << (x)) const double PI = acos(-1.0); const int Maxn = 133015; double ax[Maxn],ay[Maxn],bx[Maxn],by[Maxn]; char sa[Maxn/2],sb[Maxn/2]; int sum[Maxn]; int x1[Maxn],x2[Maxn]; int revv(int x,int bits) { int ret = 0; for (int i = 0; i < bits; i++) { ret <<= 1; ret |= x & 1; x >>= 1; } return ret; } void fft(double * a,double * b,int n,bool rev) { int bits = 0; while (1 << bits < n) ++bits; for (int i = 0; i < n; i++) { int j = revv(i,bits); if (i < j) swap(a[i],a[j]),swap(b[i],b[j]); } for (int len = 2; len <= n; len <<= 1) { int half = len >> 1; double wmx = cos(2 * PI / len),wmy = sin(2 * PI / len); if (rev) wmy = -wmy; for (int i = 0; i < n; i += len) { double wx = 1,wy = 0; for (int j = 0; j < half; j++) { double cx = a[i + j],cy = b[i + j]; double dx = a[i + j + half],dy = b[i + j + half]; double ex = dx * wx - dy * wy,ey = dx * wy + dy * wx; a[i + j] = cx + ex,b[i + j] = cy + ey; a[i + j + half] = cx - ex,b[i + j + half] = cy - ey; double wnx = wx * wmx - wy * wmy,wny = wx * wmy + wy * wmx; wx = wnx,wy = wny; } } } if (rev) { for (int i = 0; i < n; i++) a[i] /= n,b[i] /= n; } } int solve(int a[],int na,int b[],int nb,int ans[]) { int len = max(na,nb),ln; for(ln=0; L(ln)<len; ++ln); len=L(++ln); for (int i = 0; i < len ; ++i) { if (i >= na) ax[i] = 0,ay[i] =0; else ax[i] = a[i],ay[i] = 0; } fft(ax,ay,len,0); for (int i = 0; i < len; ++i) { if (i >= nb) bx[i] = 0,by[i] = 0; else bx[i] = b[i],by[i] = 0; } fft(bx,by,0); for (int i = 0; i < len; ++i) { double cx = ax[i] * bx[i] - ay[i] * by[i]; double cy = ax[i] * by[i] + ay[i] * bx[i]; ax[i] = cx,ay[i] = cy; } fft(ax,1); for (int i = 0; i < len; ++i) ans[i] = (int)(ax[i] + 0.5); return len; } int main() { int l1,l2,l; int i; while(gets(sa)) { gets(sb); memset(sum,sizeof(sum)); l1 = strlen(sa); l2 = strlen(sb); for(i = 0; i < l1; i++) x1[i] = sa[l1 - i - 1]-'0'; for(i = 0; i < l2; i++) x2[i] = sb[l2-i-1]-'0'; l = solve(x1,l1,x2,sum); for(i = 0; i<l || sum[i] >= 10; i++) // 进位 { sum[i + 1] += sum[i] / 10; sum[i] %= 10; } l = i; while(sum[l] <= 0 && l>0) l--; // 检索最高位 for(i = l; i >= 0; i--) putchar(sum[i] + '0'); // 倒序输出 putchar('n'); } return 0; } 然后模板来一发。 #define L(x) (1 << (x)) const double PI = acos(-1.0); const int Maxn = 133015; double ax[Maxn],by[Maxn]; int revv(int x,int ans[]) //两个数组求卷积,有时ans数组要开成long long { int len = max(na,1); for (int i = 0; i < len; ++i) ans[i] = (int)(ax[i] + 0.5); return len; } int solve(long long a[],int ans[]) //自己跟自己求卷积,有时候ans数组要开成long long { int len = na,ln; for(ln = 0; L(ln) < na; ++ln); len=L(++ln); for(int i = 0; i < len; ++i) { if (i >= na) ax[i] = 0,ay[i] = 0; else ax[i] = a[i],0); for(int i=0; i<len; ++i) { double cx = ax[i] * ax[i] - ay[i] * ay[i]; double cy = 2 * ax[i] * ay[i]; ax[i] = cx,1); for(int i=0; i<len; ++i) ans[i] = ax[i] + 0.5; return len; } (编辑:李大同) 【声明】本站内容均来自网络,其相关言论仅代表作者个人观点,不代表本站立场。若无意侵犯到您的权利,请及时与联系站长删除相关内容! |