第一次写FFT。这里贴出我学习FFT的博客地址,希望大家也能有所收获。
http://blog.csdn.net/ljhandlwt/article/details/51999762
https://wenku.baidu.com/view/79162aa1bed5b9f3f90f1ca8.html
当这两篇资料都看完的时候就可以领略到这个算法的美了。
#include <cmath>
#include <queue>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iostream>
#include <algorithm>
#define MAX_N 500005
using namespace std;
const double PI = acos(-1.0);
///自定义复数类
struct Complex
{
double real, imag;
Complex(double r = 0.0, double i = 0.0) : real(r), imag(i) {}
Complex operator + (const Complex& x) const
{
return Complex(real + x.real, imag + x.imag);
}
Complex operator - (const Complex& x) const
{
return Complex(real - x.real, imag - x.imag);
}
Complex operator * (const Complex& x) const
{
return Complex(real * x.real - imag * x.imag, real * x.imag + imag * x.real);
}
};
///雷德算法:将数组a中的每个数字的二进制倒转
void Rader(Complex* a, int len)
{
int j = len >> 1;
for (int i = 1; i < len - 1; i++)
{
if (i < j) swap(a[i], a[j]);
int k = len >> 1;
while (j >= k)
{
j -= k;
k >>= 1;
}
if (j < k) j += k;
}
}
void FFT(Complex* a, int len, int flag)
{
Rader(a, len);
for (int cur_length = 2; cur_length <= len; cur_length <<= 1)
{
///初始化单位复根
double ang = flag * 2 * PI / cur_length;
Complex Wn(cos(ang), sin(ang));
for (int i = 0; i < len; i += cur_length)
{
Complex W(1.0, 0.0);
for (int j = i; j < i + cur_length / 2; j++)
{
Complex u = a[j];
Complex t = W * a[j + cur_length / 2];
a[j] = u + t;
a[j + cur_length / 2] = u - t;
W = W * Wn;
}
}
}
if (flag == -1)
{
for (int i = 0; i < len; i++)
a[i].real /= len;
}
}
Complex a[MAX_N], b[MAX_N];
int result[MAX_N];
///利用快速傅里叶变换实现的高精度乘法
void Fast_multi(char* sa, char* sb)
{
int lena = strlen(sa), lenb = strlen(sb);
///两个长度为n的数字相乘,结果最长为2n
///所以答案长度最大为max(2lena, 2lenb)
///要将答案长度拓展为2的次幂,便于分治
int len = 1;
while (len < 2 * lena || len < 2 * lenb) len <<= 1;
///将两个操作数长度拓展成答案的长度
int i = 0;
for (i = 0; i < lena; i++)
{
a[i].real = sa[lena - 1 - i] - '0';
a[i].imag = 0.0;
}
while (i < len)
{
a[i].real = a[i].imag = 0.0;
i++;
}
for (i = 0; i < lenb; i++)
{
b[i].real = sb[lenb - 1 - i] - '0';
b[i].imag = 0.0;
}
while (i < len)
{
b[i].real = b[i].imag = 0.0;
i++;
}
///求向量a和向量b的卷积
FFT(a, len, 1);
FFT(b, len, 1);
for (i = 0; i < len; i++)
a[i] = a[i] * b[i];
FFT(a, len, -1);
///将a和b卷积的每位数字四舍五入保存到result数组中
for (i = 0; i < len; i++)
result[i] = a[i].real + 0.5;
///将结果进位
for (i = 0; i < len; i++)
{
result[i + 1] += result[i] / 10;
result[i] %= 10;
}
///可能有前导0
int high;
for (i = len - 1; i >= 0; i--)
{
if (result[i])
{
high = i;
break;
}
}
for (i = high; i >= 0; i--)
printf("%d", result[i]);
putchar('\n');
}
char sa[MAX_N], sb[MAX_N];
int main()
{
//freopen("in1.txt", "r", stdin);
while (~scanf("%s%s", sa, sb))
Fast_multi(sa, sb);
return 0;
}