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 280000
using 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(){

//Init
	scanf("%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';
	}
	
//Prepare
	int 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;
}