fork download
copy 
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3.  
  4. const auto PI = acos(-1);
  5.  
  6. /**
  7.  * @brief Permutes an array of numbers by the inverse of its binary index.
  8.  * @tparam T The type of value contained in $a$. 
  9.  * @param a A reference to a vector that is to be sorted.
  10.  * @return Nothing.
  11.  * 
  12.  * Changes $a$ by the reverse of the binary representation of each index
  13.  * from $[0, \texttt{a}.size())$. 
  14.  * 
  15.  * For instance, let $a = [0, 1, 2, 3, 4, 5, 6, 7]$. Applying `
  16.  * reverse_bit_sort(a)` would change $a$ into $[0, 4, 2, 6, 1, 5, 3, 7]$. 
  17.  * 
  18.  * Note that the contents of $a$ are modified. 
  19.  */
  20. template<typename T>
  21. void reverse_bit_sort(vector<T> &a) {
  22. 	int n = a.size();
  23. 	for (int i = 1, j = 0; i < n; i++) {
  24. 		// i - current index
  25. 		// j - current reversed index
  26. 		// t - current bit
  27. 		// To transition from i to i + 1, flip
  28. 		// all prefix 1s in j and the first set bit.
  29.  
  30. 		int t = n >> 1;
  31. 		for (; t & j; t >>= 1)
  32. 			j ^= t;
  33. 		j ^= t;
  34. 		if (i < j)
  35. 			swap(a[i], a[j]);
  36. 	}
  37. }
  38.  
  39. /**
  40.  * @brief Computes the discrete fourier transform of $a$.
  41.  * @tparam T A floating point type, such that $a$ contains `std::complex<T>`.
  42.  * @param a A reference to a vector that is to have the DFT computed.
  43.  * @return Nothing.
  44.  * 
  45.  * Given a polynomial $a = a[0] + a[1]x + a[2]x^2 + \dots = \sum_{i = 0}^{n - 1} a[i]x^i$, 
  46.  * computes a vector $z$ such that $z[i] = a(w_n^i)$ for $i \in [0, n)$; that is, evaluates
  47.  * $a$ at each of the roots of unity for order $n$. 
  48.  * 
  49.  * Note that the contents of $a$ are modified. 
  50.  */
  51. template<typename T> 
  52. typename enable_if<is_floating_point<T>::value, void>::type
  53. fast_fourier_transform(vector<complex<T>> &a) {
  54. 	int n = a.size();
  55.  
  56. 	// eliminate need for recursion using butterfly transform
  57. 	reverse_bit_sort(a); 
  58.  
  59. 	// iterate over length of segment
  60. 	for (int l = 2; l <= n; l <<= 1) { 
  61. 		// l - length of segment
  62. 		// theta - angle len
  63. 		// dw - change in w
  64. 		// i - current block
  65. 		// j - current index
  66. 		// w - current value
  67.  
  68. 		T theta = 2 * PI / l;
  69. 		complex<T> dw(cos(theta), sin(theta));
  70. 		for (int i = 0; i < n; i += l) { 
  71. 			complex<T> w = 1; // trivial root
  72. 			for (int j = 0; j < l / 2; j++, w *= dw) {
  73. 				auto t1 = a[i + j], t2 = a[i + j + l / 2] * w;
  74. 				a[i + j] = t1 + t2, a[i + j + l / 2] = t1 - t2;
  75. 			}
  76. 		}
  77. 	}
  78. }
  79.  
  80. /**
  81.  * @brief Computes the inverse discrete fourier transform of $a$. 
  82.  * @tparam T A floating point type, such that $a$ contains `std::complex<T>`.
  83.  * @param a A reference to a vector that is to have the inverse DFT computed.
  84.  * @return Nothing.
  85.  * 
  86.  * Given a vector $a$ such that $a[i] = p(w_n^i)$ for some polynomial $p$ of
  87.  * degree $n$, computes the vector of coefficients that constitute $p$. In 
  88.  * other words, it interpolates a polynomial evaluated at the roots of unity
  89.  * of order $n$. 
  90.  * 
  91.  * Note that the contents of $a$ are modified. 
  92.  */
  93. template<typename T> 
  94. typename enable_if<is_floating_point<T>::value, void>::type
  95. inverse_fast_fourier_transform(vector<complex<T>> &a) {
  96. 	int n = a.size();
  97.  
  98. 	// eliminate need for recursion using butterfly transform
  99. 	reverse_bit_sort(a);
  100. 	for (int l = 2; l <= n; l <<= 1) {
  101. 		// l - length of segment
  102. 		// theta - angle len (note sign change)
  103. 		// dw - change in w
  104. 		// i - current block
  105. 		// j - current index
  106. 		// w - current value
  107.  
  108. 		T theta = -2 * PI / l;
  109. 		complex<T> dw(cos(theta), sin(theta));
  110. 		for (int i = 0; i < n; i += l) {
  111. 			complex<T> w = 1;
  112. 			for (int j = 0; j < l / 2; j++) {
  113. 				auto t1 = a[i + j], t2 = a[i + j + l / 2] * w;
  114. 				a[i + j] = t1 + t2;
  115. 				a[i + j + l / 2] = t1 - t2;
  116. 				w *= dw;
  117. 			}
  118. 		}
  119. 	}
  120.  
  121. 	// divide all coefficients by $n$ 
  122. 	// (derived from inverse of vandermonde matrix)
  123. 	for (int i = 0; i < n; i++)
  124. 		a[i] /= n; 
  125. }
  126.  
  127. /**
  128.  * @brief Computes the convolution of two vectors using fast fourier transform.
  129.  * @tparam T The type contained in the vectors.
  130.  * @tparam U A floating point type (defaulted to double) used in FFT.
  131.  * @param a The first vector
  132.  * @param b The second vector
  133.  * @return The convolution of $a$ and $b$. 
  134.  * 
  135.  * Given two vectors $a$ and $b$, computes $c$ such that 
  136.  * $c[i] = \sum_{j = 0}^i a[j]b[i - j]$ for $i \in [0, n + m - 1)$. 
  137.  * 
  138.  * Note that this method runs in $\mathcal{O}(n\log n)$ (as opposed to the $\mathcal{O}(n^2)$
  139.  * naive solution). However, there may be issues with floating-point arithmatics and overflow. 
  140.  * For best results, set `U` to `long double`. 
  141.  */
  142. template<typename T, typename U = double>
  143. vector<T> convolution(const vector<T> &a, const vector<T> &b) {
  144. 	vector<complex<U>> pa(a.begin(), a.end()), pb(b.begin(), b.end());
  145.  
  146. 	// scale $n$ up to a power of two for divide and conquer to work
  147. 	int n = 1;
  148. 	while (n < a.size() + b.size()) 
  149. 		n <<= 1;
  150. 	pa.resize(n), pb.resize(n);
  151.  
  152. 	// find discrete fourier transforms of a and b
  153. 	fast_fourier_transform(pa);
  154. 	fast_fourier_transform(pb);
  155.  
  156. 	// compute point product of a and b
  157. 	for (int i = 0; i < n; i++)
  158. 		pa[i] *= pb[i];
  159.  
  160. 	// compute inverse discrete fourier transform
  161. 	inverse_fast_fourier_transform(pa);
  162.  
  163. 	// return answer (assuming that T is a real type)
  164. 	n = a.size() + b.size() - 1;
  165. 	vector<T> ret(n);
  166. 	for (int i = 0; i < n; i++)
  167. 		ret[i] = round(pa[i].real());
  168. 	return ret;
  169. }
  170.  
  171. int main() {
  172. 	int N, M;
  173. 	cin >> N >> M;
  174. 	vector<int> A(N + 1), B(M + 1);
  175. 	for (int &a : A)
  176. 		cin >> a;
  177. 	for (int &b : B)
  178. 		cin >> b;
  179. 	auto C = convolution(A, B);
  180. 	for (int c : C)
  181. 		cout << c << ' ';
  182. 	cout << '\n';
  183. }
Runtime error #stdin #stdout 0.75s 1828684KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
https://ideone.com/fYB5iv
language:
C++14 (gcc 8.3)
created:
2 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!