fork(3) download 
  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<R> operator *(const polynomial<R> &another) const{
  85.         if(is_same<R, complex<double>>::value){
  86.             int n = size()+another.size()-1;
  87.             if(!size() || !another.size() || n<=32){
  88.                 return slow_multiplication(another);
  89.             }
  90.             for(; (n&-n)!=n; n+=n&-n);
  91.             vector<complex<double>> x(n), y(n);
  92.             for(int i=0; i<(int)size(); i++){
  93.                 x[i] = a[i];
  94.             }
  95.             for(int i=0; i<(int)another.size(); i++){
  96.                 y[i] = another[i];
  97.             }
  98.             x = fft(x, 1);
  99.             y = fft(y, 1);
  100.             for(int i=0; i<n; i++){
  101.                 x[i] *= y[i];
  102.             }
  103.             polynomial<complex<double>> result(fft(x, -1));
  104.             result.resize(size()+another.size()-1);
  105.             return result;
  106.         }else if(is_same<R, double>::value || is_same<R, int>::value || is_same<R, long long>::value){
  107.             polynomial<complex<double>> f(size()), g(another.size());
  108.             for(int i=0; i<(int)size(); i++){
  109.                 f[i] = a[i];
  110.             }
  111.             for(int i=0; i<(int)another.size(); i++){
  112.                 g[i] = another[i];
  113.             }
  114.             polynomial<complex<double>> h = f * g;
  115.             polynomial<R> result(h.size());
  116.             if(is_same<R, double>::value){
  117.                 for(int i=0; i<(int)h.size(); i++){
  118.                     result[i] = h[i].real();
  119.                 }
  120.             }else{
  121.                 for(int i=0; i<(int)h.size(); i++){
  122.                     result[i] = (R)floor(h[i].real()+0.5);
  123.                 }
  124.             }
  125.             return result;
  126.         }else{
  127.             return slow_multiplication(another);
  128.         }
  129.     }
  130. };
  131.  
  132. int main(){
  133.     polynomial<int> f(2), g(2);
  134.     f[1] = 1, f[0] = 1;
  135.     g[1] = 1, g[0] = -1;
  136.     polynomial<int> h = f*g;
  137.     for(int i=(int)h.size()-1; i>=0; i--){
  138.         printf("%d%c", h[i], i? ' ': '\n');
  139.     }
  140.     return 0;
  141. }
  142.  
Compilation error #stdin compilation error #stdout 0s 0KB
comments (?)
stdin 
Standard input is empty
compilation info 
prog.cpp: In instantiation of 'polynomial<R> polynomial<R>::operator*(const polynomial<R>&) const [with R = int]':
prog.cpp:136:27:   required from here
prog.cpp:105:20: error: could not convert 'result' from 'polynomial<std::complex<double> >' to 'polynomial<int>'
             return result;
                    ^
stdout 
Standard output is empty
https://ideone.com/JbRjhU
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!