FFT bzoj2179

2179: FFT快速傅立叶

Time Limit: 10 Sec  Memory Limit: 259 MB
Submit: 724  Solved: 304
[Submit][Status]

Description

给出两个n位10进制整数x和y，你需要计算x*y。

Input

第一行一个正整数n。 第二行描述一个位数为n的正整数x。 第三行描述一个位数为n的正整数y。

Output

输出一行，即x*y的结果。

Sample Input

1
3
4

Sample Output

12

数据范围：
n<=60000

FFT

快速计算两个多项式的乘积

对于多项式A(x)=a0+a1 x+a2 x^2 + ... + an x^n，若an≠0，则称之为n次多项式

考虑两个n次多项式A(x)和B(x)

记他们的乘积为C(x)

两个十进制数可以认为是A(10),B(10)

他们的积就是C(10)

所以求两个数的乘积可以通过求C(x)的多项式系数来求

朴素实现是O(n^2)的，很简单不赘述

假设我们知道多项式在n个不同点处的值，那这个多项式的系数就确定了

因为我们可以将其当做一个n阶非齐次方程组

Y = X A

（Y是值向量，X是在n个不同点x1..xn处的范德蒙德矩阵，形如

1 x1 x1^2 x1^3 ... x1^n

1 x2 x2^2 x2^3 ... x2^n

.......................................

1 xn xn^2 xn^3 ... xn^n

A是要求的系数向量，求A可以通过给Y左乘一个X的逆来实现）

A，B是n阶的，则C最高位2n阶，为方便起见，令N为大于等于2n的2的整数次幂（处理比较方便，后面就明白）

把A，B补为N阶的（补零）

这时候引入另外一个东西，i为复数单位根

e^ix=cosx+i sinx（常识）

记e^2π/n ix为wn

有wn^0 = 1，wn^n = 1

所以wn^x=wn^(x%n)

考虑n为偶数的情况

wn^2=e^4π/n ix=e^2π/(n/2) ix=w(n/2) ①

另外wn^(n/2)=e^iπ=-1，则wn^x= - wn^(x+n/2) ②

那么(wn^x)^2 = (wn^(x+n/2))^2 ③

这个就是实现FFT的关键，因为wn的特殊性，我们可以快速求出wn^0,wn^1,..,wn^(n-1)处的点值

考虑N-1阶多项式A(x)，因为假定N是2的整数次幂，N=1就不说了

其余情况把多项式拆为两个N/2-1阶多项式

A1(x)=a0+a2x+a4x^2+...+aN-2x^(N/2-1)

A2(x)=a1+a3x+a5x^2+..+aN-1x^(N/2-1)

（就是奇数项和偶数项分开）

A(x)=A1(x^2)+xA2(x^2)

这时候我们求A(x)在wN^0,wN^1,..wN^(N-1)处的点值就可以转化为

求A1(x)和A2(x)在w(N/2)^0,w(N/2)^1,..w(N/2)^(N/2-1)的点值（由①得出）

由③可得A1((wN^k)^2)=A1(wN^(k+N/2)^2)，A2((wN^k)^2)=A2(wN^(k+N/2)^2)

那么由②可得，对于任意k<N/2，有

A(wN^k)=A1(w(N/2)^k)+w(N/2)^k A2(w(N/2)^k)
A(wN^(k+N/2))=A1(w(N/2)^k)-w(N/2)^k A2(w(N/2)^k)

可以发现，如果知道了A1(x)和A2(x)在w(N/2)^0,w(N/2)^1,...,w(N/2)^(N/2-1)这N/2个点的值，就可以在线性时间求出A(x)在wN^0,wN^1,...,wN^N-1处的点值

这个过程我们称之为离散傅里叶变换DFT，复杂度递推式为T(N)=2T(N/2)+O(N)，明显T(N)=O(NlogN)

这时候已经求出了N个点值，接下来就是用点值求出系数表达

观察

Y = X A

X已知，Y已知，X可逆（参照范德蒙德行列式的计算），很明显A=inv(X) Y，即给Y左乘一个X的逆矩阵就能求得A

X为

1 wn^0 wn^0^2 wn^0^3 ... wn^0^n

1 wn^1 wn^1^2 wn^1^3 ... wn^1^n

.......................................

1 wn^n-1 wn^n-1^2 wn^n-1^3 ... wn^n-1^n

可以发现DFT的过程实际上是利用取值的特殊性快速计算向量A左乘X的向量Y

X的逆矩阵为

1 wn^0 wn^0^2 wn^0^3 ... wn^0^n

1 wn^-1 wn^-1^2 wn^-1^3 ... wn^-1^n

.......................................

1 wn^-(n-1) wn^-(n-1)^2 wn^-(n-1)^3 ... wn^-(n-1)^n

再除以n

可以发现不考虑除以n的话，逆矩阵和原矩阵的差别只有指数的正负之差，验证一下发现指数取负性质还是相同的

所以只要将基础元素wn变成wn^-1再做一遍DFT，之后元素除以n就能得出多项式的系数表达了

#include <cstdio>#include <cstring>#include <complex>#include <cmath>#include <algorithm>#define Maxn 280000using namespace std;const double eps = 1e-4;const double pi = (double)3.1415926535897932384626;complex<double> A[Maxn],B[Maxn],C[Maxn],e[Maxn];int N, Len, Ans[Maxn];char ch;void DFT(complex<double> A[], complex<double> B[], int S, int L, int dt){if (L == 1){B[S] = A[S];return ;}int lim = L >> 1, mid = S + lim;for (int i = 0; i < lim; i++){B[S + i] = A[S + (i << 1)];B[mid + i] = A[S + (i << 1) + 1];}DFT(B, A, S, lim, dt << 1);DFT(B, A, mid, lim, dt << 1);for (int i = 0; i < lim; i++){B[S + i] = A[S + i] + e[i * dt] * A[mid + i];B[mid + i] = A[S + i] - e[i * dt] * A[mid + i];}}int main(){//Initscanf("%d\n", &N);for (int i = N - 1; i >=0; i--){ch = getchar();A[i] = ch - '0';}scanf("\n");for (int i = N - 1; i >=0; i--){ch = getchar();B[i] = ch - '0';}//Prepareint k = 1;while (k < N)k <<= 1;N = k << 1;for (int i = 0; i <= N; i++)e[i] = complex<double>(cos(2 * pi * i / N), sin(2 * pi * i / N));DFT(A, C, 0, N, 1);memcpy(A, C, sizeof(C));DFT(B, C, 0, N, 1);memcpy(B, C, sizeof(C));for (int i = 0; i < N; i++)C[i] = A[i] * B[i];for (int i = 0; i <= N; i++)e[i] = complex<double>(cos( - 2 * pi * i / N), sin( - 2 * pi * i / N));DFT(C, A, 0, N, 1);Len = N - 1;for (int i = 0; i <= Len || Ans[i] != 0; i++){Ans[i] += A[i].real() / N + eps;Ans[i + 1] += Ans[i] / 10;Ans[i] %= 10;Len = max(Len, i);}while (Len > 0 && Ans[Len] == 0)Len--;for (int i = Len; i >= 0; i--)printf("%d", Ans[i]);return 0;}