思路:
算法导论第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[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; }
FFT乘法模板
最新推荐文章于 2024-03-22 15:43:29 发布