Advertisement
_no0B

Untitled

Nov 14th, 2021
903
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 2.57 KB | None | 0 0
  1. /**
  2. Iterative Implementation of Discrete Fast Fourier Transform
  3. Complexity: O(N log N)
  4.  
  5. Possible Optimizations:
  6. 1. Remove leading zeros (No, seriously!!!!)
  7. 2. Writing a custom complex class reduces runtime.
  8. 3. If there are multiple testcases of similar sizes, it
  9.    might be faster to precalculate the roots of unity.
  10. **/
  11.  
  12. #include<bits/stdc++.h>
  13. using namespace std;
  14.  
  15. typedef complex<double> CD;
  16. typedef long long LL;
  17. const double PI = acos(-1);
  18.  
  19. struct FFT
  20. {
  21.     int N;
  22.     vector<int> perm;
  23.  
  24.     void precalculate() {
  25.         perm.resize(N);
  26.         perm[0] = 0;
  27.  
  28.         for (int k=1; k<N; k<<=1) {
  29.             for (int i=0; i<k; i++) {
  30.                 perm[i] <<= 1;
  31.                 perm[i+k] = 1 + perm[i];
  32.             }
  33.         }
  34.     }
  35.  
  36.     void fft(vector<CD> &v, bool invert = false) {
  37.         if (v.size() != perm.size()) {
  38.             N = v.size();
  39.             assert(N && (N&(N-1)) == 0);
  40.             precalculate();
  41.         }
  42.  
  43.         for (int i=0; i<N; i++)
  44.             if (i < perm[i])
  45.                 swap(v[i], v[perm[i]]);
  46.  
  47.         for (int len = 2; len <= N; len <<= 1) {
  48.             double angle = 2 * PI / len;
  49.             if (invert) angle = -angle;
  50.             CD factor = polar(1.0, angle);
  51.  
  52.             for (int i=0; i<N; i+=len) {
  53.                 CD w(1);
  54.                 for (int j=0; j<len/2; j++) {
  55.                     CD x = v[i+j], y = w * v[i+j+len/2];
  56.                     v[i+j] = x+y;
  57.                     v[i+j+len/2] = x-y;
  58.                     w *= factor;
  59.                 }
  60.             }
  61.         }
  62.         if (invert)
  63.             for (CD &x : v) x/=N;
  64.     }
  65.  
  66.     vector<LL> multiply(const vector<LL> &a, const vector<LL> &b) {
  67.         vector<CD> fa(a.begin(), a.end()), fb(b.begin(), b.end());
  68.  
  69.         int n = 1;
  70.         while (n < a.size()+ b.size())  n<<=1;
  71.         fa.resize(n);
  72.         fb.resize(n);
  73.  
  74.         fft(fa);
  75.         fft(fb);
  76.         for (int i=0; i<n; i++) fa[i] *= fb[i];
  77.         fft(fa, true);
  78.  
  79.         vector<LL> ans(n);
  80.         for (int i=0; i<n; i++)
  81.             ans[i] = round(fa[i].real());
  82.         return ans;
  83.     }
  84. };
  85.  
  86. ///Solves SPOJ POLYMUL
  87. ///Given two polynomials, find their product
  88.  
  89. int main()
  90. {
  91.     int t;
  92.     cin>>t;
  93.     FFT fft;
  94.  
  95.     while (t--)
  96.     {
  97.         int n;
  98.         cin>>n;
  99.  
  100.         vector<LL> a(n+1), b(n+1), ans;
  101.         for (int i=0; i<=n; i++)    cin>>a[n-i];
  102.         for (int i=0; i<=n; i++)    cin>>b[n-i];
  103.  
  104.         ans = fft.multiply(a, b);
  105.  
  106.         for (int i=2*n; i>=0; i--)
  107.             cout<<ans[i]<<" ";
  108.         cout<<endl;
  109.     }
  110. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement