证明请参考算法概论的快速博立叶。
FFT递归:
#include <bits/stdc++.h>
using namespace std;
const int maxn = 5005;
const double PI = acos ( -1.0 );
struct node
{
double real, image;
node ( double a = 0, double b = 0 ) : real ( a ), image ( b ) { }
node operator + ( const node &b )
{
return node ( real+b.real, image+b.image );
}
node operator - ( const node &b )
{
return node ( real-b.real, image-b.image );
}
node operator * ( const node &b )
{
return node ( real*b.real-image*b.image,
real*b.image+image*b.real );
}
};
node a[maxn], b[maxn];
int re[maxn];
node * FFT ( node x[maxn], int n, int dft )
{
node * ans = ( node * ) malloc ( sizeof ( node )*n );
if ( n == 1 )
{
ans[0] = x[0];
return ans;
}
node lf[maxn], rt[maxn];
for ( int i = 0; i < n; i += 2 )
lf[i/2] = x[i];
for ( int i = 1; i < n; i += 2 )
rt[i/2] = x[i];
node * L = FFT ( lf, n/2, dft );
node * R = FFT ( rt, n/2, dft );
node nw = node ( cos ( dft*2*PI/n ), sin ( dft*2*PI/n ) );
node w = node ( 1, 0 );
for ( int i = 0; i < n/2; i ++ )
{
ans[i] = L[i]+w*R[i];
ans[i+n/2] = L[i]-w*R[i];
w = w*nw;
}
return ans;
}
void solve ( )
{
int n;
scanf ( "%d", &n );
for ( int i = 0; i < n; i ++ )
{
scanf ( "%lf", &a[i].real );
a[i].image = 0;
}
for ( int i = 0; i < n; i ++ )
{
scanf ( "%lf", &b[i].real );
b[i].image = 0;
}
int p = 1;
while ( p < 2*n )
p <<= 1;
for ( int i = n; i < p; i ++ )
a[i].real = a[i].image = b[i].real = b[i].image = 0;
node * pa = FFT ( a, p, 1 );
node * pb = FFT ( b, p, 1 );
for ( int i = 0; i < p; i ++ )
a[i] = pa[i]*pb[i];
pa = FFT ( a, p, -1 );
for ( int i = 0; i < p; i ++ )
re[i] = pa[i].real/p+0.5;
for ( int i = 0; i < p; i ++ )
printf ( "%d ", re[i] );
printf ( "\n" );
}
int main ( )
{
solve ( );
return 0;
}
FFT迭代代码(HDU1402):
#include <bits/stdc++.h>
using namespace std;
const int maxn = 500005;
const double PI = acos ( -1.0 );
struct node
{
double real, image;
node ( double a = 0, double b = 0 ) : real ( a ), image ( b ) { }
node operator + ( const node &b )
{
return node ( real+b.real, image+b.image );
}
node operator - ( const node &b )
{
return node ( real-b.real, image-b.image );
}
node operator * ( const node &b )
{
return node ( real*b.real-image*b.image,
real*b.image+image*b.real );
}
};
node a[maxn], b[maxn], W[maxn];
int re[maxn];
char x[maxn], y[maxn];
void change ( node a[], int n )
{
int p = 1;
while ( ( 1 << p ) < n )
p ++;
for ( int i = 0; i < n; i ++ )
{
int j = 0, k = i, t = p;
while ( t -- )
{
j <<= 1;
j |= k&1;
k >>= 1;
}
if ( j > i )
swap ( a[i], a[j] );
}
}
void FFT ( node a[], int n, int dft )
{
change ( a, n );
for ( int i = 1; i < n; i *= 2 )
{
node nw ( cos ( dft*PI*2/2/i ), sin ( dft*PI*2/2/i ) );
for ( int j = 0; j < n; j += 2*i )
{
node w ( 1, 0 );
for ( int k = 0; k < i; k ++ )
{
node tmp = a[j+k];
node res = w*a[j+k+i];
a[j+k] = tmp+res;
a[j+k+i] = tmp-res;
w = w*nw;
}
}
}
if ( dft == -1 )
{
for ( int i = 0; i < n; i ++ )
a[i].real /= n;
}
}
void solve ( )
{
int lenx, leny;
while ( ~ scanf ( "%s%s", x, y ) )
{
lenx = strlen ( x );
for ( int i = 0; i < lenx; i ++ )
{
a[i].real = x[lenx-1-i]-'0';
a[i].image = 0;
}
leny = strlen ( y );
for ( int i = 0; i < leny; i ++ )
{
b[i].real = y[leny-1-i]-'0';
b[i].image = 0;
}
int p = 1;
while ( p < 2*lenx || p < 2*leny )
p <<= 1;
for ( int i = lenx; i < p; i ++ )
a[i].real = a[i].image = 0;
for ( int i = leny; i < p; i ++ )
b[i].real = b[i].image = 0;
FFT ( a, p, 1 );
FFT ( b, p, 1 );
for ( int i = 0; i < p; i ++ )
a[i] = a[i]*b[i];
FFT ( a, p, -1 );
for ( int i = 0; i < p; i ++ )
re[i] = a[i].real+0.5;
for ( int i = 0; i < p; i ++ )
{
re[i+1] += re[i]/10;
re[i] %= 10;
}
int h = 0;
for ( int i = p-1; i >= 0; i -- )
if ( re[i] )
{
h = i;
break ;
}
for ( int i = h; i >= 0; i -- )
printf ( "%d", re[i] );
printf ( "\n" );
}
}
int main ( )
{
solve ( );
return 0;
}