多项式乘法 （快速傅里叶变换）
算法导论里基本看懂了傅里叶变换的整个运行流程。
后面剩下的就是如何规划设计程序，使得程序漂亮又高效，暂时没时间，先做个记录，后面慢慢实现。
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
https://www.zybuluo.com/397915842/note/37965
网上摘录某段代码，可供参考：
【多项式乘法 （快速傅里叶变换）】
#include
2.#include
3.#include
4.const int MAXN = 2e5 + 5;
5.const double PI = acos(-1.0);
6.#define max(a, b) (a) > (b) ? (a) : (b)
7.
8.class Complex
9.{
10.public:
11.double real, imag;
12.
13.Complex(double real = 0.0, double imag = 0.0)
14.{
15.this->real = real, this->imag = imag;
16.}
17.
18.Complex operator - (const Complex &elem) const
19.{
20.return Complex(this->real - elem.real, this->imag - elem.imag);
21.}
22.
23.Complex operator + (const Complex &elem) const
24.{
25.return Complex(this->real + elem.real, this->imag + elem.imag);
26.}
27.
28.Complex operator * (const Complex &elem) const
29.{
30.return Complex(this->real * elem.real - this->imag * elem.imag, this->real * elem.imag + this->imag * elem.real);
31.}
32.
33.void setValue(double real = 0.0, double imag = 0.0)
34.{
35.this->real = real, this->imag = imag;
36.}
37.};
38.
39.Complex A[MAXN], B[MAXN];
40.int res[MAXN], len, mlen, len1, len2;
41.char str1[MAXN >> 1], str2[MAXN >> 1];
42.
43.void Swap(Complex &a, Complex &b)
44.{
45.Complex tmp = a;
46.a = b;
47.b = tmp;
48.}
49.
50.void Prepare()
51.{
52.len1 = strlen(str1), len2 = strlen(str2);
53.mlen = max(len1, len2);
54.len = 1;
55.
56.// 将 len 扩大到 2 的整数幂
57.while(len < (mlen << 1))
58.len <<= 1;
59.
60.//初始化多项式的系数
61.for(int i = 0;
i < len1;
++ i)
62.A[i].setValue(str1[len1 - i - 1] - '0', 0);
63.
64.for(int i = 0;
i < len2;
++ i)
65.B[i].setValue(str2[len2 - i - 1] - '0', 0);
66.
67.// 补 0
68.for(int i = len1;
i < len;
++ i)
69.A[i].setValue();
70.
71.for(int i = len2;
i < len;
++ i)
72.B[i].setValue();
73.}
74.
75.//雷德算法 位逆序置换
76.void Rader(Complex y[])
77.{
78.for(int i = 1, j = len >> 1, k;
i < len - 1;
++ i)
79.{
80.if(i < j)
81.Swap(y[i], y[j]);
82.
83.k = len >> 1;
84.
85.while(j >= k)
86.{
87.j -= k;
88.k >>= 1;
89.}
90.
91.if(j < k)
92.j += k;
93.}
94.}
95.
96.//DFT : op == 1
97.//IDFT : op == -1
98.void FFT(Complex y[], int op)
99.{
100.//先位逆序置换
101.Rader(y);
102.
103.// h 为每次要处理的长度， h = 1 时不需处理
104.for(int h = 2;
h <= len;
h <<= 1)
105.{
106.// Wn = e^(2 * PI / n)，如果是插值，那么 Wn = e^(-2 * PI / n)
107.Complex Wn(cos(op * 2 * PI / h), sin(op * 2 * PI / h));
108.
109.for(int i = 0;
i < len;
i += h)
110.{
111.//旋转因子，初始化为 e^0
112.Complex W(1, 0);
113.
114.for(int j = i;
j < i + h / 2;
++ j)
115.{
116.Complex u = y[j];
117.Complex t = W * y[j + h / 2];
118.//蝴蝶操作
119.y[j] = u + t;
120.y[j + h / 2] = u - t;
121.//每次更新旋转因子
122.W = W * Wn;
123.}
124.}
125.}
126.
127.// 插值的时候要除以 len
128.if(op == -1)
129.for(int i = 0;
i < len;
++ i)
130.y[i].real /= len;
131.}
132.
133.//DFT 后将 A 和 B 相应点值相乘，将结果放到 res 里面
134.void Convolution(Complex *A, Complex *B)
135.{
136.//evaluation
137.FFT(A, 1), FFT(B, 1);
138.
139.for(int i = 0;
i < len;
++ i)
140.A[i] = A[i] * B[i];
141.
142.//interpolation
143.FFT(A, -1);
144.
145.for(int i = 0;
i < len;
++ i)
146.res[i] = (int)(A[i].real + 0.5);
147.}
148.
149.void Adjustment(int *arr)
150.{
151.//次数界为 len，所以不用担心进位不会进到第 len 位
152.for(int i = 0;
i < len;
++ i)
153.{
154.res[i + 1] += res[i] / 10;
155.res[i] %= 10;
156.}
157.
158.//去除多余的 0
159.while(-- len && res[len] == 0);
160.}
161.
162.void Display(int *arr)
163.{
164.for(int i = len;
i >= 0;
-- i)
165.putchar(arr[i] + '0');
166.
167.putchar('\n');
168.}
169.
170.int main()
171.{
172.while(gets(str1) && gets(str2))
173.{
174.Prepare();
175.Convolution(A, B);
176.Adjustment(res);
177.Display(res);
178.}
179.
180.return 0;
181.}
推荐阅读
- neo4j|neo4j cql语句 快速查询手册
- 通过复盘快速成长（附模板）
- jar|springboot项目打成jar包和war包，并部署（快速打包部署）
- 快速阅读作业【2/21】《阅读（革命性新定义》）
- Spring注解05|Spring注解05 @Import 给容器快速导入一个组件
- JS/JavaScript|JS/JavaScript CRC8多项式 16进制
- 蒙哥马利乘法原理
- 虐虐快速通关秘籍放这，高度概括就是听话照做！
- 想要用正确的方法快速进阶，这三个方法推荐给你。
- vue快速上手-3