fork(1) download
copy 
  1. #include <cstdio>
  2. #include <cmath>
  3. #include <cassert>
  4. #include <complex>
  5. #include <vector>
  6. #include <type_traits>
  7.  
  8. using namespace std;
  9.  
  10. #define FFT_LEN 65536
  11.  
  12. const double pi = acos(-1);
  13.  
  14. vector<complex<double>> fft(const vector<complex<double>> &x, const int &inv) {
  15.   static bool fft_ready = false;
  16.   static complex<double> loli[FFT_LEN];
  17.   if (!fft_ready) {
  18.     for (int k = 0; k<FFT_LEN; k++) {
  19.       loli[k] = exp(complex<double>(0, 2 * pi*k / FFT_LEN));
  20.     }
  21.     fft_ready = true;
  22.   }
  23.   int n = x.size();
  24.   assert((n&-n) == n && n <= FFT_LEN && abs(inv) == 1);
  25.   vector<complex<double>> X = x;
  26.   for (int i = 1, j = 0; i<n; i++) {
  27.     for (int k = n >> 1; !((j ^= k)&k); k >>= 1);
  28.     if (i < j) {
  29.       swap(X[i], X[j]);
  30.     }
  31.   }
  32.   for (int i = 2; i <= n; i *= 2) {
  33.     int d = (inv == 1) ? FFT_LEN - (FFT_LEN / i) : FFT_LEN / i;
  34.     for (int j = 0; j<n; j += i) {
  35.       for (int k = 0, a = 0; k<i / 2; k++, a = (a + d) % FFT_LEN) {
  36.         complex<double> s = X[j + k], t = loli[a] * X[j + k + i / 2];
  37.         X[j + k] = s + t;
  38.         X[j + k + i / 2] = s - t;
  39.       }
  40.     }
  41.   }
  42.   if (inv == -1) {
  43.     for (int i = 0; i<(int)X.size(); i++) {
  44.       X[i] /= n;
  45.     }
  46.   }
  47.   return X;
  48. }
  49.  
  50. template<class R> class polynomial {
  51. private:
  52.   vector<R> a;
  53.   polynomial<R> slow_multiplication(const polynomial<R> &another) const {
  54.     if (!size() || !another.size()) {
  55.       return polynomial<R>();
  56.     }
  57.     polynomial<R> result(size() + another.size() - 1);
  58.     for (int i = 0; i<(int)size(); i++) for (int j = 0; j<(int)another.size(); j++) {
  59.       result[i + j] += a[i] * another[j];
  60.     }
  61.     return result;
  62.   }
  63. public:
  64.   polynomial(const size_t &n = 0) {
  65.     a = vector<R>(n);
  66.   }
  67.   polynomial(const vector<R> &coef) {
  68.     a = coef;
  69.   }
  70.   size_t size() const {
  71.     return a.size();
  72.   }
  73.   void resize(const size_t &n) {
  74.     a.resize(n);
  75.   }
  76.   R& operator [](const int &i) {
  77.     assert(0 <= i && i<(int)a.size());
  78.     return a[i];
  79.   }
  80.   const R& operator [](const int &i) const {
  81.     assert(0 <= i && i<(int)a.size());
  82.     return a[i];
  83.   }
  84.   polynomial<complex<double>> operator*(const polynomial<complex<double>> &another) const {
  85.     int n = size() + another.size() - 1;
  86.     if (!size() || !another.size() || n <= 32) {
  87.       return slow_multiplication(another);
  88.     }
  89.     for (; (n&-n) != n; n += n&-n);
  90.     vector<complex<double>> x(n), y(n);
  91.     for (int i = 0; i<(int)size(); i++) {
  92.       x[i] = a[i];
  93.     }
  94.     for (int i = 0; i<(int)another.size(); i++) {
  95.       y[i] = another[i];
  96.     }
  97.     x = fft(x, 1);
  98.     y = fft(y, 1);
  99.     for (int i = 0; i<n; i++) {
  100.       x[i] *= y[i];
  101.     }
  102.     polynomial<complex<double>> result(fft(x, -1));
  103.     result.resize(size() + another.size() - 1);
  104.     return result;
  105.   }
  106.   template<class T>
  107.   polynomial<T> operator*(const polynomial<T> &another) const {
  108.     if (is_same<T, complex<double>>::value) {
  109.     } else if (is_same<T, double>::value || is_same<T, int>::value || is_same<T, long long>::value) {
  110.       polynomial<complex<double>> f(size()), g(another.size());
  111.       for (int i = 0; i<(int)size(); i++) {
  112.         f[i] = a[i];
  113.       }
  114.       for (int i = 0; i<(int)another.size(); i++) {
  115.         g[i] = another[i];
  116.       }
  117.       polynomial<complex<double>> h = f * g;
  118.       polynomial<T> result(h.size());
  119.       if (is_same<T, double>::value) {
  120.         for (int i = 0; i<(int)h.size(); i++) {
  121.           result[i] = h[i].real();
  122.         }
  123.       } else {
  124.         for (int i = 0; i<(int)h.size(); i++) {
  125.           result[i] = (T)floor(h[i].real() + 0.5);
  126.         }
  127.       }
  128.       return result;
  129.     } else {
  130.       return slow_multiplication(another);
  131.     }
  132.   }
  133. };
  134.  
  135. int main() {
  136.   polynomial<int> f(2), g(2);
  137.   f[1] = 1, f[0] = 1;
  138.   g[1] = 1, g[0] = -1;
  139.   polynomial<int> h = f*g;
  140.   for (int i = (int)h.size() - 1; i >= 0; i--) {
  141.     printf("%d%c", h[i], i ? ' ' : '\n');
  142.   }
  143.   return 0;
  144. }
  145.  
Success #stdin #stdout 0s 4500KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1 0 -1
https://ideone.com/mosfZW
language:
C++14 (gcc 8.3)
created:
8 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!