Advertisement
willy108

fft

Apr 9th, 2022
1,406
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 3.42 KB | None | 0 0
  1.  
  2. //misaka will carry me to master
  3. #include <iostream>
  4. #include <cstdio>
  5. #include <cstring>
  6. #include <cmath>
  7. #include <utility>
  8. #include <cassert>
  9. #include <algorithm>
  10. #include <vector>
  11. #include <functional>
  12. #include <numeric>
  13. #include <set>
  14. #include <map>
  15. #include <complex>
  16.  
  17. #define ll long long
  18. #define lb long double
  19. #define sz(vec) ((int)(vec.size()))
  20. #define all(x) x.begin(), x.end()
  21. #define pb push_back
  22. #define mp make_pair
  23.  
  24. const lb eps = 1e-9;
  25. const ll mod = 1e9 + 7, ll_max = 1e18;
  26. //const ll mod = (1 << (23)) * 119 +1;
  27. const int MX = 2e5 +10, int_max = 0x3f3f3f3f;
  28.  
  29. using namespace std;
  30.  
  31. #define lp std::complex<double>
  32. #define hsb(x) (31 - __builtin_clz(x))
  33.  
  34. const int D = 19;
  35. const int NN = (1 << D);
  36. const lb tau = 4.0 * acos(0);
  37.  
  38. //this function returns the reversed version of x with b bits
  39. int rev(unsigned int x, int b){ //lots of trust
  40.     x = (((x & 0xaaaaaaaa) >> 1) | ((x & 0x55555555) << 1));
  41.     x = (((x & 0xcccccccc) >> 2) | ((x & 0x33333333) << 2));
  42.     x = (((x & 0xf0f0f0f0) >> 4) | ((x & 0x0f0f0f0f) << 4));
  43.     x = (((x & 0xff00ff00) >> 8) | ((x & 0x00ff00ff) << 8));
  44.     //https://aticleworld.com/5-way-to-reverse-bits-of-an-integer/ method 4
  45.     return((x >> 16) | (x << 16)) >> (32 - b);
  46. }
  47.  
  48. //evaluates the roots of unity of poly
  49. //if op == 1, reverses it as well
  50. vector<lp> halfft(vector<lp> poly, int op = 0){
  51.     vector<lp> ans(sz(poly));  //answers will be stored here
  52.     static lp tmp[NN]; //auxilary memory
  53.     if(sz(poly) == 1){ //hardcode when its just a constant term
  54.         ans[0] = poly[0];
  55.         return ans;
  56.     }
  57.     int p2 = hsb(sz(poly)); //this denotes log2(sz(poly))
  58.     for(int i = 0; i<sz(poly); i++){
  59.         ans[rev(i, p2)] = poly[i]; //this is the base case of storing it into the bitwise reversal
  60.     }
  61.     for(int k = 1, n = 2; k<=p2; k++, n *= 2){ //find the answer for blocks size 2^k as k goes from 1 -> log2(sz(poly))
  62.         for(int pos = 0; pos < sz(poly); pos += n){ //iterate over the starting spot
  63.             for(int i = 0; i<n; i++){
  64.                 lp x_0;
  65.                 if(op == 0) x_0 = exp(lp(0, tau*(lb)i/(lb)(n))); // normal root
  66.                 else x_0 = exp(lp(0, tau*(lb)(n-i)/(lb)(n))); //inverted root
  67.                 tmp[pos + i] = ans[pos + i%(n/2)] + x_0 * ans[pos + n/2 + (i%(n/2))];
  68.                 //i am just really smart so this works, i write into ans[] smartly so that i dont overlap
  69.             }
  70.             for(int i = 0; i<n; i++){
  71.                 ans[pos + i] = tmp[pos + i]; //copy tmp into ans
  72.             }
  73.         }
  74.     }
  75.     return ans;
  76. }
  77.  
  78. vector<lp> fft(const vector<int>& a, const vector<int>& b){
  79.     int s;
  80.     for(s = 1; s<(sz(a) +sz(b)); s*=2);
  81.     vector<lp> A(s), B(s);
  82.     for(int i = 0; i<s; i++){
  83.         A[i] = B[i] = 0;
  84.         if(i < sz(a)) A[i] = a[i];
  85.         if(i < sz(b)) B[i] = b[i];
  86.     } //make the two size 2^k polynomials A and B
  87.     vector<lp> aa = halfft(A), bb = halfft(B); //halfft A and B
  88.     for(int i = 0; i<sz(aa); i++){
  89.         aa[i] *= bb[i]; //multiply them together
  90.     }  
  91.     vector<lp> res = halfft(aa, 1), ans(s);
  92.     for(int i = 0; i<s; i++){
  93.         ans[i] = (lb)(real(res[i]))/(lb)(s); //divide by s and write into the answer
  94.        
  95.     }
  96.     return ans;
  97. }
  98.  
  99. void solve(){
  100.     int n, m;
  101.     cin >> n >> m;
  102.     n++, m++;
  103.     vector<int> p1(n), p2(m);
  104.     for(int i = 0; i<n; i++){
  105.         cin >> p1[i];
  106.     }
  107.     for(int i =0; i<m; i++){
  108.         cin >> p2[i];
  109.     }
  110.     vector<lp> ans = fft(p1, p2);
  111.     for(int i = 0; i<n + m -1; i++){
  112.         cout << (int)round(real(ans[i])) << " ";
  113.     }
  114.     cout << "\n";
  115. }
  116.  
  117. int main(){
  118.   cin.tie(0) -> sync_with_stdio(0);
  119.  
  120.   int T = 1;
  121.   //cin >> T;
  122.   while(T--){
  123.         solve();
  124.   }
  125.   return 0;
  126. }
  127.  
  128.  
  129.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement