HDU1402 A * B Problem Plus 快速傅里叶变换FFT

看了算法导论（第2版）第30章的《多项式与快速傅里叶变换》

多项式有系数表示法和点值表示法。

两个次数界为n的多项式A（x）和B（x）相乘，输入输出均采用系数表示法。（假定n为2的幂）

1）使次数界增加一倍：A（x）和B（x）扩充为次数界为2n的多项式，并构造起系数表示

2）求值：两次应用2n阶FFT，计算出A（x）和B（x）的长度为2n的点值表示

3）点乘：计算多项式C(x)=A（x）*B（x）的点值表示

4）插值：对2n个点值对应用一次FFT计算出其逆DFT，就可以构造出多项式C（x）的系数表示

第1步和第3步需要时间为O(n)，第2步和第4步运用FFT需要时间为O(nlgn)。

第2步求取n个点的值所需时间为O(n^2)，如果我们巧妙地选取x（k）则其运行时间变为O(nlgn)，FFT主要利用了单位复根（w^n=1的w就是n次单位复根，他们是e^(2PI*k/n),k=0,1,……n-1）的特殊性质。

用FFT做SPOJ235竟然TLE，传说要用NTT（数论变换number theoretic transform），有待提高ing

Accepted 1402 375MS 6740K 3099 B

#include < iostream >
#include < cmath >
using namespace std;
const double PI = acos( - 1.0 );
struct vir{
double re,im; // 实部和虚部
vir( double a = 0 , double b = 0 ){
re = a; im = b;
}
vir operator + ( const vir & b){
return vir(re + b.re,im + b.im);
}
vir operator - ( const vir & b){
return vir(re - b.re, im - b.im);
}
vir operator * ( const vir & b){
return vir(re * b.re - im * b.im , re * b.im + im * b.re);
}
};
vir x1[ 200005 ],x2[ 200005 ];
void change(vir * x, int len, int loglen)
{
int i,j,k,t;
for (i = 0 ;i < len;i ++ )
{
t = i;
for (j = k = 0 ; j < loglen; j ++ ,t >>= 1 )
k = (k << 1 ) | (t & 1 );
if (k < i){
vir wt = x[k];
x[k] = x[i];
x[i] = wt;
}
}
}
void fft(vir * x, int len, int loglen)
{
int i,j,t,s,e;
change(x,len,loglen);
t = 1 ;
for (i = 0 ;i < loglen;i ++ ,t <<= 1 )
{
s = 0 ; e = s + t;
while (s < len)
{
vir a,b,wo(cos(PI / t),sin(PI / t)),wn( 1 , 0 );
for (j = s;j < s + t;j ++ )
{
a = x[j]; b = x[j + t] * wn;
x[j] = a + b; x[j + t] = a - b;
wn = wn * wo;
}
s = e + t; e = s + t;
}
}
}
void dit_fft(vir * x, int len, int loglen)
{
int i,j,s,e,t = 1 << loglen;
for (i = 0 ;i < loglen;i ++ )
{
t >>= 1 ;
s = 0 ; e = s + t;
while (s < len)
{
vir a,b,wn( 1 , 0 ),wo(cos(PI / t), - sin(PI / t));
for (j = s;j < s + t;j ++ )
{
a = x[j] + x[j + t]; b = (x[j] - x[j + t]) * wn;
x[j] = a; x[j + t] = b;
wn = wn * wo;
}
s = e + t; e = s + t;
}
}
change(x,len,loglen);
for (i = 0 ;i < len;i ++ )
x[i].re /= len;
}
int main()
{
char a[ 100005 ],b[ 100005 ];
int i,len1,len2,len,loglen;
int t,over;
while (scanf( " %s%s " ,a,b) != EOF)
{
len1 = strlen(a) << 1 ;
len2 = strlen(b) << 1 ;
len = 1 ;
loglen = 0 ;
while (len < len1)
len <<= 1 , loglen ++ ;
while (len < len2)
len <<= 1 , loglen ++ ;
// init
for (i = 0 ;a[i];i ++ )
x1[i].re = a[i] - ' 0 ' , x1[i].im = 0 ;
for (;i < len;i ++ )
x1[i].re = x1[i].im = 0 ;
for (i = 0 ;b[i];i ++ )
x2[i].re = b[i] - ' 0 ' , x2[i].im = 0 ;
for (;i < len;i ++ )
x2[i].re = x2[i].im = 0 ;
fft(x1,len,loglen);
fft(x2,len,loglen);
for (i = 0 ;i < len;i ++ )
x1[i] = x1[i] * x2[i];
dit_fft(x1,len,loglen);

// 将x1.re从浮点数转化为十进制整型存入a
for (i = (len1 + len2) / 2 - 2 ,over = len = 0 ;i >= 0 ;i -- )
{
t = x1[i].re + over + 0.5 ;
a[len ++ ] = t % 10 ;
over = t / 10 ;
}
while (over){
a[len ++ ] = over % 10 ;
over /= 10 ;
}

// output
for (len -- ;len >= 0 &&! a[len];len -- );
if (len < 0 )
putchar( ' 0 ' );
else
for (;len >= 0 ;len -- )
putchar(a[len] + ' 0 ' );
putchar( ' \n ' );
}
return 0 ;
}

 

转载于:https://www.cnblogs.com/DreamUp/archive/2010/09/05/1818576.html