· 博主语文体育老师教的
· 本文年龄限定 16+
· 吐槽上面 2 条的都是⑨
老实说我一点也不想写这东西……
但Delayyy君都写掉了那你说我能不写么> <~
需要深入了解的话推荐:
1、《算法合集 —— 多项式乘法》( 张嘉琳 )
2、《算法导论》关于快速傅里叶变换的章节.
3、《一切皆整数的 FFT》关于点值 X 取整数的方法 (百度空间 AekdyCoin)
自己就把关于复数的证明写了……纠结了半天的东西.
设 W ( k, n ) 表示 1 的 N 次单位根 ( 注意在复数情况下是正好有 n 个解的 ). 则 W (k, n) 具有以下性质:
W (k, n) = cos (2 * π * k / n) + i * sin(2 * π * k / n)
折半性 : W (k, n) = W (2k, 2n) (代入上式 系数2 直接约掉)
周期性 : W (k, n) = W ( k + n, n) (三角函数的周期是 n )
半反性 : W (k, n) = - W(k + n / 2, n) (同上, 自己YY个三角函数就懂了)
则对于 FFT 主算法中的式子证明 :
S(x) = G(x) + x · F(x)
S(w(k, n)) = G(w(k, n)²) + w(k, n) · F(w(k, n)²)
S(w(k, n)) = G(w(k, n / 2)) + w(k, n) · F(w(k, n / 2))
并且:
S(w(k + n / 2, n)) = G(w(k, n / 2)) + w(k + n / 2, n) · F(w(k, n / 2))
S(w(k + n / 2, n)) = G(w(k, n / 2)) - w(k + n / 2, n) · F(w(k, n / 2))
另外, 弄完主算法后的插值的证明就不放了……如此巨坑还是请自行Google
以上证明完毕. FFT的分治还是很强大的......
参考 : CSDN evil live 《FFT高精乘法》
CSDN delayyy 《FFT快速傅里叶变换》
以下是很挫的Code ( 用 FFT 优化高精度乘法 )
#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <algorithm>
const double pi = acos((double)-1);
using namespace std;
char s[10010];
long long c[10010];
int len_a, len_b, n, t;
struct line
{
double x, y;
} a[10010], b[10010], w[10010], y[10010];
line operator - (line &a, line &b) { line c; c.x = a.x - b.x; c.y = a.y - b.y; return c; }
line operator + (line &a, line &b) { line c; c.x = a.x + b.x; c.y = a.y + b.y; return c; }
line operator * (line &a, line &b) { line c; c.x = a.x * b.x - a.y * b.y; c.y = a.x * b.y + a.y * b.x; return c; }
void fft(line *a, int s, int t)
{
int len = n >> t;
if (len == 1) return; line wk;
fft(a, s, t + 1); fft(a, s + (1 << t), t + 1);
for (int k = 0; k < (len >> 1); ++k)
{
int p = s + (k << (t + 1));
wk = w[k << t] * a[p + (1 << t)];
y[k] = a[p] + wk,
y[k + (len >> 1)] = a[p] - wk;
}
for (int k = 0; k < len; ++k) a[s + (k << t)] = y[k];
}
void set(char *s, line *a, int &len) //读入, 四位一压
{
int lenk, sta, i;
gets(s + 1);
lenk = strlen(s + 1); sta = lenk % 4; len = (lenk + 4) / 4 - 1;
for (i = 1; i <= sta; ++i) a[len].x = a[len].x * 10 + s[i] - '0';
for (; i <= lenk; i += 4)
a[--len].x = (s[i] - '0') * 1000 + (s[i + 1] - '0') * 100 + (s[i + 2] - '0') * 10 + s[i + 3] - '0';
len = (lenk + 3) / 4;
}
int main()
{
freopen("input.txt", "r", stdin);
freopen("output.txt", "w", stdout);
set(s, a, len_a); set(s, b, len_b);
n = len_a + len_b;
for (t = 1; t < n; t <<= 1); n = t;
for (int i = 0; i < n; ++i) w[i].x = cos(pi * i * 2 / n), w[i].y = sin(pi * 2 * i / n);
fft(a, 0, 0); //a 化点值过程
fft(b, 0, 0); //b 化点值过程
for (int i = 0; i < n; ++i) a[i] = a[i] * b[i], w[i].y = -w[i].y; //逆矩阵处理
fft(a, 0, 0); //插值过程
for (int i = 0; i < n; ++i) c[i] += (long long) round(a[i].x / n), c[i + 1] = c[i] / 10000, c[i] %= 10000;
int top = n;
while (!c[top] && top > 0) --top;
printf("%d", c[top]);
for (int i = top - 1; i >= 0; --i) printf("%04d", c[i]);
}