fork download
copy 
  1. #include <iostream>
  2. #include <cmath>
  3. #include <cassert>
  4. #include <algorithm>
  5. #include <vector>
  6. #include <complex>
  7.  
  8. #define lb long double
  9. #define sz(vec) ((int)(vec.size()))
  10. #define all(x) x.begin(), x.end()
  11.  
  12. using namespace std;
  13.  
  14. #define lp std::complex<double>
  15. #define hsb(x) (31 - __builtin_clz(x))
  16.  
  17. const int D = 19;
  18. const int NN = (1 << D);
  19. const lb tau = 4.0 * acos(0);
  20.  
  21. //this function returns the reversed version of x with b bits
  22. int rev(unsigned int x, int b){ //lots of trust
  23. 	x = (((x & 0xaaaaaaaa) >> 1) | ((x & 0x55555555) << 1));
  24. 	x = (((x & 0xcccccccc) >> 2) | ((x & 0x33333333) << 2));
  25. 	x = (((x & 0xf0f0f0f0) >> 4) | ((x & 0x0f0f0f0f) << 4));
  26. 	x = (((x & 0xff00ff00) >> 8) | ((x & 0x00ff00ff) << 8));
  27. 	//https://a...content-available-to-author-only...d.com/5-way-to-reverse-bits-of-an-integer/ method 4
  28. 	return((x >> 16) | (x << 16)) >> (32 - b);
  29. }
  30.  
  31. //evaluates the roots of unity of poly
  32. //if op == 1, reverses it as well
  33. vector<lp> halfft(const vector<lp>& poly, const vector<lp>& rt, int op = 0){
  34. 	vector<lp> ans(sz(poly));
  35. 	if(sz(poly) == 1){ //hardcode size 1 polynomials
  36. 		ans[0] = poly[0];
  37. 		return ans;
  38. 	}
  39. 	int p2 = hsb(sz(poly)); //2^p2 is the size of the polynomial
  40. 	for(int i = 0; i<sz(poly); i++){
  41. 		ans[rev(i, p2)] = poly[i]; 
  42. 	} //this is the base case of storing each coeff into the reverse bitwise form
  43.  
  44. 	for(int k = 1, n = 2; k<=p2; k++, n *= 2){ //iterate n from small to large and find the answers for the n-th roots at each step
  45. 		for(int pos = 0; pos < sz(poly); pos += n){ //solve for some block on the polynomial
  46. 			for(int i = 0; i<n/2; i++){
  47. 				int i1 = i, i2 = i + n/2;
  48. 				lp r1 = ans[pos + i] + rt[i1 * (1 << (p2 - k))] * ans[pos + n/2 + i];
  49. 				lp r2 = ans[pos + i] + rt[i2 * (1 << (p2 - k))] * ans[pos + n/2 + i];
  50. 				//this is the square sum lemma which states
  51. 				//f(x) = f+(x^2) + x * f-(x^2)
  52. 				//ans already stores the answers for f+(x^2) and f-(x^2)
  53. 				//x is acessed by figuring out which position in rt the root actually corrosponds to
  54. 				ans[pos + i1] = r1;
  55. 				ans[pos + i2] = r2;
  56. 				//by storying r1 and r2 i dont need auxilary memory
  57. 				//ans is always written into in a way such that ans[pos + i] has the answer
  58. 				//for the index i of the block of [pos, pos + 2^k)
  59. 			}
  60. 		}
  61. 	}
  62. 	if(op) //reverse the roots
  63. 		reverse(1 + all(ans));
  64. 	return ans;
  65. }
  66.  
  67. vector<lp> fft(const vector<int>& a, const vector<int>& b){
  68. 	int s;
  69. 	for(s = 1; s<(sz(a) +sz(b)); s*=2); //find the smallest power of 2 the answer polynomial can fit into
  70. 	vector<lp> A(s, 0), B(s, 0), rt(s);
  71. 	for (int i = 0; i < sz(a); i++) {
  72. 		A[i] = a[i];
  73. 	}
  74. 	for (int i = 0; i < sz(b); i++) {
  75. 		B[i] = b[i];
  76. 	}
  77. 	for(int i = 0; i<s; i++){
  78. 		rt[i] = exp(lp(0, tau*(lb)i/(lb)(s))); //precompute the roots
  79. 	}
  80. 	//make the two size 2^k polynomials A and B
  81. 	vector<lp> aa = halfft(A, rt), bb = halfft(B, rt); //halfft A and B
  82. 	for(int i = 0; i<sz(aa); i++){
  83. 		aa[i] *= bb[i]; //multiply them together
  84. 	}	
  85. 	vector<lp> ans = halfft(aa, rt, 1); //do the reverse fft on the new product
  86. 	for(int i = 0; i<s; i++){
  87. 		ans[i] /= (lb)s;	//divide by s and write into the answer	
  88. 	}
  89. 	ans.resize(sz(a) + sz(b) - 1);
  90. 	return ans; 
  91. }
  92.  
  93.  
  94. int main(){
  95.   	int n, m;
  96. 	cin >> n >> m;
  97. 	n++, m++;
  98. 	vector<int> p1(n), p2(m);
  99. 	for(int i = 0; i<n; i++){
  100. 		cin >> p1[i];
  101. 	}
  102. 	for(int i =0; i<m; i++){
  103. 		cin >> p2[i];
  104. 	}
  105. 	vector<lp> ans = fft(p1, p2);
  106. 	for(int i = 0; i<sz(ans); i++){
  107. 		cout << (int)round(real(ans[i])) << " ";
  108. 	}
  109. 	cout << "\n";
  110. 	return 0;
  111. }
  112.  
  113.  
  114.  
  115.  
Success #stdin #stdout 0.01s 5548KB
comments (?)
stdin
copy 
1 2
1 2
1 3 1
stdout
copy 
1 5 7 2
https://ideone.com/yEkBY2
language:
C++14 (gcc 8.3)
created:
2 years ago
visibility:
secret

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!