多項式乗算（高速フーリエ変換）
4676 ワード
アルゴリズムの導論ではフーリエ変換の全体的な流れが大体分かりました。
後はどのように設計プログラムを計画しますか？プログラムが綺麗で効率的です。しばらく時間がなくて、先に記録を作って、後はゆっくり実現します。
http://blog.miskcoo.com/2015/04/polynomial-multiplication-and-fast-fourier-transform
https://www.zybuluo.com/397915842/note/37965
ネット上である部分のコードを抜粋して、参考にすることができます。
後はどのように設計プログラムを計画しますか？プログラムが綺麗で効率的です。しばらく時間がなくて、先に記録を作って、後はゆっくり実現します。
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('
');
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.}