fork download
copy 
  1. //https://l...content-available-to-author-only...j.ac/p/108
  2. #include <bits/stdc++.h>
  3. #include <ext/pb_ds/assoc_container.hpp>
  4. #include <ext/pb_ds/tree_policy.hpp>
  5. using namespace __gnu_pbds;
  6. using namespace std;
  7.  
  8. #define pb push_back
  9. #define ff first
  10. #define ss second
  11.  
  12. typedef long long ll;
  13. typedef long double ld;
  14. typedef pair<int, int> pii;
  15. typedef pair<ll, ll> pll;
  16. typedef pair<ld, ld> pld;
  17.  
  18. const int INF = 1e9;
  19. const ll LLINF = 1e18;
  20. const int MOD = 1e9 + 7;
  21.  
  22. template<class K> using sset =  tree<K, null_type, less<K>, rb_tree_tag, tree_order_statistics_node_update>;
  23.  
  24. inline ll ceil0(ll a, ll b) {
  25.     return a / b + ((a ^ b) > 0 && a % b);
  26. }
  27.  
  28. void setIO() {
  29.     ios_base::sync_with_stdio(0); cin.tie(0);
  30. }
  31.  
  32. const double PI = acos(-1);
  33. typedef complex<ld> cint;
  34.  
  35. //Returns a complex number in the form of e^2*PI*k*i/n (i is an imaginary number)
  36. cint root(int n, int k){
  37.     return exp(cint(0, 2*PI*k/n));
  38. }
  39.  
  40. //Returns the bitwise reverse of x using only the first mx bits
  41. int rev(int x, int mx){
  42.     int ret = 0;
  43.     for(int i = 0; i <= mx; i++){
  44.         if(x & (1 << i)) ret |= 1 << (mx - i);
  45.     }
  46.     return ret;
  47. }
  48.  
  49. /**
  50.  * Finds the Discrete Fourier Transform (DFT) of a polynomial using the
  51.  * Fast Fourier Transform (FTT) method in O(n log n) time complexity.
  52.  * 
  53.  * Definitions:
  54.  * 
  55.  *  root(n, k) = e^(2*PI*i*k/n), where i is an imaginary number
  56.  *  X = {x_1, x_2, x_3, ..., x_n}
  57.  *  Y = {y_1, y_2, y_3, ..., y_n}
  58.  *  F(x) = {x_1, x_2 * x, x_3 * x^2, ..., x_n * x^(n - 1)}
  59.  * 
  60.  *  DFT(X) = Y, defined as follows
  61.  * 
  62.  *      y_1 = F(root(n, 0))
  63.  *      y_2 = F(root(n, 1))
  64.  *      y_3 = F(root(n, 2))
  65.  *         ...
  66.  *      y_n = F(root(n, n - 1))
  67.  * 
  68.  * But how can we calculate Y, when DFT(X) = Y?
  69.  *
  70.  * Solution: 
  71.  *  First, we see that we can split F(x) into even and odd parts, lets call them E(x), and O(x)
  72.  * 
  73.  *      O(x) = x_1 + x_3 * x + x_5 * x^2 ... x_(n - 1) * x^(n/2 - 1)
  74.  *      E(x) = x_2 + x_4 * x + x_6 * x^2 ... x_n * x^(n/2 - 1)
  75.  * 
  76.  *  Then, F(x) becomes
  77.  * 
  78.  *      F(x) = O(x^2) + x * E(x^2)
  79.  * 
  80.  *  Since y_i = F(root(n, i - 1)), y_i becomes
  81.  * 
  82.  *      y_i = O(root(n, i - 1)^2) + root(n, i - 1) * E(root(n, i - 1)^2)
  83.  * 
  84.  *  Simplifying further, this becomes
  85.  *  
  86.  *      y0_i = O(root(n/2, i - 1)^2) = O(root(n, i - 1))
  87.  *      y1_i = E(root(n/2, i - 1)^2) = E(root(n, i - 1))
  88.  *      y_i = y0_i + root(n, i)*y1_i, i < n/2
  89.  *      y_(n/2 + i) = y0_i + root(n, n/2 + i)*y1_i, n/2 <= i
  90.  * 
  91.  *  From here, we can find Y recursively. Since every split halves the n, and every merge takes O(n) complexity,
  92.  *  amortized time complexity is O(n log n) 
  93.  *  
  94.  *  To calculate invDFT, it is the same as calculating DFT with only a few changes. Simply negate k when computing
  95.  *  root(n, k), as well as dividing the final answer by n. This can be proven by using inverse matrices.
  96.  * 
  97.  * @param v The polynomial to perfom FFT on. The size of v should be a power of 2
  98.  * @param sign Tndicates whether or not the function computes the 
  99.  *             DFT if sign is 1 or invDFT if sign is -1.
  100.  */
  101. void fft(vector<cint> &v, int sign){
  102.     int n = v.size(), lg = log2(n) - 1;
  103.     cint tmp[n];  
  104.     /**
  105.      * swaps each position with its bitwise inverse. 
  106.      * This step is needed to ensure correct ordering for the iterative implementation,
  107.      * as some of the even and odd indices are swapped every iteration.
  108.      */
  109.     for(int i = 0; i < n; i++){
  110.         if(i < rev(i, lg)) swap(v[i], v[rev(i, lg)]);
  111.     }
  112.     for(int i = 2; i <= n; i *= 2){
  113.         //compute each segment of length i
  114.         for(int j = 0; j < n; j += i){
  115.             for(int k = 0; k < i/2; k++){
  116.                 //computes the left half of the segment using the formula described above
  117.                 tmp[k] = v[j + k] + root(i, sign*k)*v[j + i/2 + k];
  118.                 //computes the right half of the segment using the formula described above
  119.                 tmp[i/2 + k] = v[j + k] + root(i, sign*(i/2 + k))*v[j + i/2 + k];
  120.             }
  121.             //updates the actual values to the new values
  122.             for(int k = 0; k < i; k++) v[j + k] = tmp[k];
  123.         }
  124.     }
  125.     //If we are computing invDFT, we must divide every value by n
  126.     if(sign == -1) for(cint &i : v) i /= n;
  127. }
  128.  
  129. /**
  130.  * Multiplies two polynomials using Discrete Fourier Transform (DFT) in O(n log n) time complexity.
  131.  * 
  132.  * Definitions:
  133.  *  X = {x_1, x_2, x_3, ..., x_n}
  134.  *  Y = {y_1, y_2, y_3, ..., y_n}
  135.  * 
  136.  *  DFT(X) = Y
  137.  *  invDFT(Y) = X
  138.  *  invDFT(DFT(X)) = X
  139.  * 
  140.  *  X * Y = the product of polynomials X and Y
  141.  *  DFT(X) * DFT(Y) = {x_1 * y_1, x_2 * y_2, x_3 * y_3, ..., x_n * y_n}
  142.  *  DFT(X * Y) = DFT(X) * DFT(Y)
  143.  * 
  144.  * Solution:
  145.  *  Given two vectors A and B, lets define their product vector as ANS = A * B
  146.  *  Then, using the distributive property, we can see that
  147.  * 
  148.  *      DFT(A * B) = DFT(A) * DFT(B)
  149.  *      DFT(ANS) = DFT(A) * DFT(B)
  150.  * 
  151.  *  Applying invDFT to both sides, we get
  152.  * 
  153.  *      ANS = invDFT(DFT(A) * DFT(B))
  154.  * 
  155.  *  Performing DFT(A) * DFT(B) is much simpler than A * B, and can be done in O(n) time complexity
  156.  *  We can calculate DFT and invDFT in O(n log n) time complexity, giving us a total complexity of O(n log n)
  157.  * 
  158.  * @param a the first polynomial
  159.  * @param b the second polynomial
  160.  * @return returns the polynomial a * b
  161.  */
  162. vector<ll> mult(const vector<ll> &a, const vector<ll> &b){
  163.     //resizes the answer to a power of 2
  164.     int n = 1;
  165.     while(n < a.size() + b.size()) n *= 2;
  166.     //constructs two copies of polynomials a and b, but using complex numbers instead
  167.     vector<cint> aa(a.begin(), a.end()), bb(b.begin(), b.end());
  168.     //resize the polynomials to have a size with a power of 2
  169.     aa.resize(n), bb.resize(n);
  170.     //computes DFT(a) and DFT(b)
  171.     fft(aa, 1), fft(bb, 1);
  172.     //sets a = DFT(a) * DFT(b) as described above
  173.     for(int i = 0; i < n; i++) aa[i] *= bb[i];
  174.     //returns a to a polynomial by setting a = invDFT(a)
  175.     fft(aa, -1);
  176.     //extracts the answer from the final polynomial. The values must be rounded due to double precision
  177.     vector<ll> ret(n);
  178.     for(int i = 0; i < n; i++){
  179.         ret[i] = round(aa[i].real());
  180.     }
  181.     return ret;
  182. }
  183.  
  184. int main(){
  185.     setIO();
  186.     int n, m;
  187.     cin >> n >> m;
  188.     n++, m++;
  189.     vector<ll> v1(n), v2(m);
  190.     for(int i = 0; i < n; i++) cin >> v1[i];
  191.     for(int i = 0; i < m; i++) cin >> v2[i];
  192.     vector<ll> ans = mult(v1, v2);
  193.     //for(auto i : ans) cout << i << endl;
  194.     for(int i = 0; i < n + m - 1; i++) cout << ans[i] << " ";
  195.     cout << endl;
  196. }
Runtime error #stdin #stdout #stderr 0.01s 5528KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Standard output is empty
stderr
copy 
terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc
https://ideone.com/cbyTP4
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!