fork download
copy 
  1. #include <bits/stdc++.h>
  2. using namespace std;
  3.  
  4. void fastIO(){
  5.     ios_base::sync_with_stdio(false);
  6.     cin.tie(0);
  7. }
  8.  
  9. int N, M;
  10.  
  11. typedef complex<long double> cp;
  12. const long double PI = acos(-1);
  13.  
  14. void fft(vector<cp> &a, bool invert = false){
  15.     int n = (int)a.size();
  16.     // if n == 1 then A(w^0) = A(0) = a0
  17.     if (n == 1){
  18.         return;
  19.     }
  20.     //Let A_even(x) = a0+a2*x+a4*x^2+a6*x^3+...
  21.     //A_odd(x) = a1+a3*x+a5*x^2+a7*x^3+...
  22.     vector<cp> even, odd;
  23.     //0 2 4 6 8..
  24.     for (int i = 0; i < n; i += 2){
  25.         even.push_back(a[i]);
  26.     }
  27.     //1 3 5 7 9..
  28.     for (int i = 1; i < n; i += 2){
  29.         odd.push_back(a[i]);
  30.     }
  31.     /**
  32.      * divide and conquer - calculate what A_even(x^2) and A_odd(x^2) are, the sizes of each are n/2
  33.      * Put n roots of unity into A in this function, put n/2 roots of unity into the recursive one
  34.      * */
  35.     fft(even, invert);
  36.     fft(odd, invert);
  37.     //in each iteration, w = e^(2*pi*i/n)
  38.     //initially w = e^0, then multiply by e^(2*pi/n)
  39.     cp w(1, 0), mult(cos(2*PI/n), sin(2*PI/n));
  40.     //in the inverted case everything is to the power of -1
  41.     if (invert) mult = {cos(-2*PI/n), sin(-2*PI/n)};
  42.     for (int i = 0; i < n/2; i++){
  43.         //A(x) = A_even(x^2)+x*A_odd(x^2), let x = w_n^i
  44.         //if i < n/2, A(w_n^i) = A_even(w_(n/2)^(2*i)) + w_n^i * A_odd(w_(n/2)^(2*i))
  45.         a[i] = even[i]+w*odd[i];
  46.         //if i >= n/2, the same equation holds but i > n/2 would be out of bounds of A_even and A_odd
  47.         //rearranging and substituting eq for i < n/2 yields the equation for i >= n/2
  48.         a[i+n/2] = even[i]-w*odd[i];
  49.         //in the inverted matrix, everything is divided by n at the end
  50.         //since n = 2^(layers), divide the inverted matrix by 2 each time does the same thing as dividing by n at the end
  51.         if (invert){
  52.             a[i] /= 2;
  53.             a[i+n/2] /= 2;
  54.         }
  55.         //update w_n^i
  56.         w *= mult;
  57.     }
  58. }
  59.  
  60. int main(){
  61.     fastIO();
  62.     cin >> N >> M;
  63.     N++, M++;
  64.     vector<cp> acoeff, bcoeff;
  65.     for (int i = 1; i <= N; i++){
  66.         int a; cin >> a;
  67.         acoeff.push_back(a);
  68.     }
  69.     for (int i = 1; i <= M; i++){
  70.         int b; cin >> b;
  71.         bcoeff.push_back(b);
  72.     }
  73.     //array sizes must be a power of 2 since the size is divided by 2 in each recursive layer
  74.     int ind = 0;
  75.     while (N+M > (1<<ind)){
  76.         ind++;
  77.     }
  78.     for (int i = N; i < (1<<ind); i++) acoeff.push_back(0);
  79.     for (int i = M; i < (1<<ind); i++) bcoeff.push_back(0);
  80.     //compute A(w_n^0), A(w_n^1), ...
  81.     fft(acoeff);
  82.     //compute B(w_n^0), B(w_n^1), ...
  83.     fft(bcoeff);
  84.     //convolution - C(w_n^i) = A(w_n^i)*B(w_n^i)
  85.     vector<cp> c((1<<ind));
  86.     for (int i = 0; i < (1<<ind); i++){
  87.         c[i] = acoeff[i]*bcoeff[i];
  88.     }
  89.     //inverted matrix for fft - calculate C(w_n^(-i))
  90.     fft(c, true);
  91.     vector<long long> result(1<<ind);
  92.     for (int i = 0; i < N+M-1; i++){
  93.         result[i] = round(c[i].real());
  94.         cout << result[i] << " ";
  95.     }
  96. }
Success #stdin #stdout 0.01s 5392KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
29822521
https://ideone.com/35iCAi
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!