FFT乘法模板

思路:

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

一个知识点是：

  A(x)=a₀+a₁x+a₂x²+a₃x³+……+a_n-1x^n-1

 A^[0](x)=a₀+a₂x+a₄x²+……+a_n-2x^n/2-1

 A^[1](x)=a₁+a₃x+a₅x²+……+a_n-1x^n/2-1

 

 A^[0](x²)+x*A^[1](x²)=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 200000using 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[N], x2[N];char a[N / 2], b[N / 2];int sum[N];    //存储最后的结果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;}