fork 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 T>
  51. class polynomial_base {
  52. public:
  53.   polynomial_base(size_t n = 0) : a(n) {
  54.   }
  55.  
  56.   polynomial_base(const std::vector<R> &coef) : a(coef) {
  57.   }
  58.  
  59.   size_t size() const {
  60.     return a.size();
  61.   }
  62.  
  63.   void resize(size_t n) {
  64.     a.resize(n);
  65.   }
  66.  
  67.   R &operator[](int i) {
  68.     assert((0 <= i) && (i < static_cast<int>(a.size())));
  69.     return a[i];
  70.   }
  71.  
  72.   const R &operator[](int i) const {
  73.     return (*const_cast<polynomial_base *>(this))[i];
  74.   }
  75.  
  76.   T operator*(const polynomial_base &another) const {
  77.     if (!size() || !another.size()) {
  78.       return T();
  79.     }
  80.     T result(size() + another.size() - 1);
  81.     for (int i = 0; i<(int)size(); i++) for (int j = 0; j<(int)another.size(); j++) {
  82.       result[i + j] += a[i] * another[j];
  83.     }
  84.     return result;
  85.   }
  86.  
  87. private:
  88.   std::vector<R> a;
  89. };
  90.  
  91. template<class R>
  92. class polynomial : public polynomial_base<R, polynomial<R>> {
  93. public:
  94.   polynomial(size_t n = 0) : polynomial_base<R, polynomial<R>>(n) {
  95.   }
  96.  
  97.   polynomial(const std::vector<R> &coef) : polynomial_base<R, polynomial<R>>(coef) {
  98.   }
  99. };
  100.  
  101. template<>
  102. class polynomial<std::complex<double>> : public polynomial_base<std::complex<double>, polynomial<std::complex<double>>> {
  103. public:
  104.   polynomial(size_t n = 0) : polynomial_base<std::complex<double>, polynomial<std::complex<double>>>(n) {
  105.   }
  106.  
  107.   polynomial(const std::vector<std::complex<double>> &coef) : polynomial_base<std::complex<double>, polynomial<std::complex<double>>>(coef) {
  108.   }
  109.  
  110.   polynomial<complex<double>> operator*(const polynomial<complex<double>> &another) const {
  111.     int n = size() + another.size() - 1;
  112.     if (!size() || !another.size() || n <= 32) {
  113.       return polynomial_base<std::complex<double>, polynomial<std::complex<double>>>::operator*(another);
  114.     }
  115.     for (; (n&-n) != n; n += n&-n);
  116.     vector<complex<double>> x(n), y(n);
  117.     for (int i = 0; i<(int)size(); i++) {
  118.       x[i] = (*this)[i];
  119.     }
  120.     for (int i = 0; i<(int)another.size(); i++) {
  121.       y[i] = another[i];
  122.     }
  123.     x = fft(x, 1);
  124.     y = fft(y, 1);
  125.     for (int i = 0; i<n; i++) {
  126.       x[i] *= y[i];
  127.     }
  128.     polynomial<complex<double>> result(fft(x, -1));
  129.     result.resize(this->size() + another.size() - 1);
  130.     return result;
  131.   }
  132. };
  133.  
  134. template<class R, class T>
  135. class polynomial_primitive : public polynomial_base<R, T> {
  136. public:
  137.   polynomial_primitive(size_t n = 0) : polynomial_base<R, T>(n) {
  138.   }
  139.  
  140.   polynomial_primitive(const std::vector<R> &coef) : polynomial_base<R, T>(coef) {
  141.   }
  142.  
  143.   T operator*(const polynomial_primitive &another) const {
  144.     polynomial<complex<double>> f(this->size()), g(another.size());
  145.     for (int i = 0; i<(int)this->size(); i++) {
  146.       f[i] = (*this)[i];
  147.     }
  148.     for (int i = 0; i<(int)another.size(); i++) {
  149.       g[i] = another[i];
  150.     }
  151.     polynomial<complex<double>> h = f * g;
  152.     T result(h.size());
  153.     for (int i = 0; i<(int)h.size(); i++) {
  154.       result[i] = (R)T::round(h[i].real());
  155.     }
  156.     return result;
  157.   }
  158. };
  159.  
  160. template<class R, class T>
  161. class polynomial_real : public polynomial_primitive<R, T> {
  162. public:
  163.   static double round(double x) {
  164.     return x;
  165.   }
  166.  
  167.   polynomial_real(size_t n = 0) : polynomial_primitive<R, T>(n) {
  168.   }
  169.  
  170.   polynomial_real(const std::vector<R> &coef) : polynomial_primitive<R, T>(coef) {
  171.   }
  172. };
  173.  
  174. template<class R, class T>
  175. class polynomial_integer : public polynomial_primitive<R, T> {
  176. public:
  177.   static double round(double x) {
  178.     return floor(x + 0.5);
  179.   }
  180.  
  181.   polynomial_integer(size_t n = 0) : polynomial_primitive<R, T>(n) {
  182.   }
  183.  
  184.   polynomial_integer(const std::vector<R> &coef) : polynomial_primitive<R, T>(coef) {
  185.   }
  186. };
  187.  
  188. template<>
  189. class polynomial<double> : public polynomial_real<double, polynomial<double>> {
  190. public:
  191.   polynomial(size_t n = 0) : polynomial_real<double, polynomial<double>>(n) {
  192.   }
  193.  
  194.   polynomial(const std::vector<double> &coef) : polynomial_real<double, polynomial<double>>(coef) {
  195.   }
  196. };
  197.  
  198. template<>
  199. class polynomial<int> : public polynomial_integer<int, polynomial<int>> {
  200. public:
  201.   polynomial(size_t n = 0) : polynomial_integer<int, polynomial<int>>(n) {
  202.   }
  203.  
  204.   polynomial(const std::vector<int> &coef) : polynomial_integer<int, polynomial<int>>(coef) {
  205.   }
  206. };
  207.  
  208. template<>
  209. class polynomial<long long> : public polynomial_integer<long long, polynomial<long long>> {
  210. public:
  211.   polynomial(size_t n = 0) : polynomial_integer<long long, polynomial<long long>>(n) {
  212.   }
  213.  
  214.   polynomial(const std::vector<long long> &coef) : polynomial_integer<long long, polynomial<long long>>(coef) {
  215.   }
  216. };
  217.  
  218. int main() {
  219.   polynomial<int> f(2), g(2);
  220.   f[1] = 1, f[0] = 1;
  221.   g[1] = 1, g[0] = -1;
  222.   polynomial<int> h = f*g;
  223.   for (int i = (int)h.size() - 1; i >= 0; i--) {
  224.     printf("%d%c", h[i], i ? ' ' : '\n');
  225.   }
  226.   return 0;
  227. }
  228.  
Success #stdin #stdout 0s 4496KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1 0 -1
https://ideone.com/Ecs0iN
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!