FFT乘法模板

思路:

   算法导论第30章有详细说明。此处只是简略说明其主要的步骤。

一个知识点是：

  A(x)=a0+a1x+a2x2+a3x3+……+an-1xn-1

 A[0](x)=a0+a2x+a4x2+……+an-2xn/2-1

 A[1](x)=a1+a3x+a5x2+……+an-1xn/2-1

 

 A[0](x2)+x*A[1](x2)=A(x)  

以上是 二进制平摊反转置换跟求和的主要式子。

多项式有两种表示形式：点值表示，系数表示。

快速FFT主要有以下四点：

   1. 使次数界（上界）增加一倍。A(x)、B(x)的长度扩充到2*n

   2. 求值。主要是求点值表示A(x)、B(x)的点值表示

   3. 点乘。C(x)=A(x)*B(x)

   4. 插值。对C(x)进行插值，求出其系数表示。

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#define Pi acos(-1.0)//定义Pi的值
#define N 200000
using namespace std;
struct complex //定义复数结构体
{
double re, im;
complex ( double r = 0.0, double i = 0.0 )
{  re = r, im = i; } //初始化
//定义三种运算
complex operator + ( complex o )
{ return complex ( re + o.re, im + o.im );}
complex operator - ( complex o )
{ return complex ( re - o.re, im - o.im );}
complex operator * ( complex o )
{ return complex ( re * o.re - im * o.im, re * o.im + im * o.re );}
} x1
, x2
;
char a[N / 2], b[N / 2];
int sum
;    //存储最后的结果

void BRC ( complex *y, int len ) //二进制反转倒置
{
int i, j, k;
for ( i = 1, j = len / 2; i < len - 1; i++ )
{
if ( i < j ) { swap ( y[i], y[j] ); } //i<j保证只交换一次
k = len / 2;
while ( j >= k )
{
j -= k; k = k / 2;
}
if ( j < k ) { j += k; }
}
}
void FFT ( complex *y, int len , double on ) //on=1表示顺，-1表示逆
{
int i, j, k, h;
complex u, t;
BRC ( y, len );
for ( h = 2; h <= len; h <<= 1 ) //控制层数
{
//初始化单位复根
complex wn ( cos ( on * 2 * Pi / h ), sin ( on * 2 * Pi / h ) );
for ( j = 0; j < len; j += h ) //控制起始下标
{
//初始化螺旋因子
complex w ( 1, 0 );
for ( k = j; k < j + h / 2; k++ )
{
u = y[k];
t = w * y[k + h / 2];
y[k] = u + t;
y[k + h / 2] = u - t;
w = w * wn; //更新螺旋因子
}
}
}
if ( on == -1 )
for ( i = 0; i < len; i++ ) //逆FFT(IDFT)
{
y[i].re /= len;
}

}
int main()
{
int len1, len2, len, i;
while ( scanf ( "%s%s", a, b ) != EOF )
{
len1 = strlen ( a );
len2 = strlen ( b );
len = 1;
//扩充次数界至2*n
while ( len < 2 * len1 || len < 2 * len2 ) { len <<= 1; } //右移相当于len=len*2
//倒置存储
for ( i = 0; i < len1; i++ )
{ x1[i].re = a[len1 - 1 - i] - '0'; x1[i].im = 0.0;}
for ( ; i < len1; i++ ) //多余次数界初始化为0
{x1[i].re = x1[i].im = 0.0;}
for ( i = 0; i < len2; i++ )
{ x2[i].re = b[len2 - 1 - i] - '0'; x2[i].im = 0.0;}
for ( ; i < len2; i++ ) //多余次数界初始化为0
{x2[i].re = x2[i].im = 0.0;}
//FFT求值
FFT ( x1, len, 1 ); //FFT(a) 1表示顺 -1表示逆
FFT ( x2, len, 1 ); //FFT(b)
//点乘，结果存入x1
for ( i = 0; i < len; i++ )
{
x1[i] = x1[i] * x2[i];
}
//插值，逆FFT（IDTF）
FFT ( x1, len, -1 );

//细节处理
for ( i = 0; i < len; i++ )
{
sum[i] = x1[i].re + 0.5;    //四舍五入
}
for ( i = 0; i < len; i++ ) //进位
{
sum[i + 1] += sum[i] / 10;
sum[i] %= 10;
}
//输出
len = len1 + len2 - 1;
while ( sum[len] <= 0 && len > 0 ) { len--; } //检索最高位
for ( i = len; i >= 0; i-- ) //倒序输出
{
cout << sum[i];
}
cout << endl;
}
return 0;
}