Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- #include <bits/stdc++.h>
- using namespace std;
- // define pi
- const double PI = 3.14159265359;
- // compute the DFT of x or the IDFT if invert
- void dft(vector<complex<double>>& x, bool invert = false) {
- int N = x.size();
- if (N <= 1) return;
- // spilt x into its even-indexed and odd-indexed terms
- vector<complex<double>> E(N / 2), O(N / 2);
- for (int i = 0; i < N / 2; i++) E[i] = x[2*i], O[i] = x[2*i+1];
- // compute the DFT/IDFT of each half
- dft(E, invert);
- dft(O, invert);
- // w_k are the Nth roots of unity
- complex<double> w(1);
- double theta = 2 * PI / N * (invert ? -1 : 1);
- complex<double> rot(cos(theta), sin(theta));
- // combine the halves
- for (int k = 0; k < N / 2; k++) {
- x[k] = E[k] + w * O[k]; // x[k] for k = 0...N/2-1
- x[k + N/2] = E[k] - w * O[k]; // x[k+N/2] by periodicity
- if (invert) x[k] /= 2, x[k + N / 2] /= 2; // needed if inverting
- w *= rot; // rotate w to the next root of unity
- }
- }
- // compute the coefficients of A(x)*B(x) where a and b are the coefficients of A(x) and B(x) respectively
- vector<int> poly_mul(vector<int>& a, vector<int>& b)
- {
- // convert to complex data type
- vector<complex<double>> ca(a.begin(), a.end());
- vector<complex<double>> cb(b.begin(), b.end());
- vector<int> res;
- // zero pad a and b until length is a power of 2
- int N = 1;
- while(N < ca.size() + cb.size()) N <<= 1;
- ca.resize(N);
- cb.resize(N);
- // compute the DFT of A and B
- dft(ca);
- dft(cb);
- // multiply the results and compute the IDFT
- for(int i=0;i<N;i++) ca[i] *= cb[i];
- dft(ca, true);
- // extract the coefficients and return
- for(int i=0;i<N;i++) res.push_back(round(ca[i].real()));
- while(res[res.size()-1]==0) res.pop_back();
- return res;
- }
- int main()
- {
- vector<int> a = {-1, 2, 5};
- vector<int> b = {1, 2};
- vector<int> ab = poly_mul(a, b);
- for(int e: ab) cout << e << ' '; cout << '\n';
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement