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 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.     template<class S>
  64.     typename enable_if<is_same<R, S>::value && is_same<S, complex<double>>::value, polynomial<S>>::type
  65.     multiplication_impl(const polynomial<S>& another) const {
  66.             int n = size()+another.size()-1;
  67.             if(!size() || !another.size() || n<=32){
  68.                 return slow_multiplication(another);
  69.             }
  70.             for(; (n&-n)!=n; n+=n&-n);
  71.             vector<complex<double>> x(n), y(n);
  72.             for(int i=0; i<(int)size(); i++){
  73.                 x[i] = a[i];
  74.             }
  75.             for(int i=0; i<(int)another.size(); i++){
  76.                 y[i] = another[i];
  77.             }
  78.             x = fft(x, 1);
  79.             y = fft(y, 1);
  80.             for(int i=0; i<n; i++){
  81.                 x[i] *= y[i];
  82.             }
  83.             polynomial<complex<double>> result(fft(x, -1));
  84.             result.resize(size()+another.size()-1);
  85.             return result;
  86.     }
  87.     template<class S>
  88.     typename enable_if<
  89.              is_same<R, S>::value &&
  90.              (is_same<S, double>::value || is_same<S, int>::value || is_same<S, long long>::value)
  91.              ,polynomial<S>>::type
  92.     multiplication_impl(const polynomial<S>& another) const {
  93.             polynomial<complex<double>> f(size()), g(another.size());
  94.             for(int i=0; i<(int)size(); i++){
  95.                 f[i] = a[i];
  96.             }
  97.             for(int i=0; i<(int)another.size(); i++){
  98.                 g[i] = another[i];
  99.             }
  100.             polynomial<complex<double>> h = f * g;
  101.             polynomial<R> result(h.size());
  102.             if(is_same<R, double>::value){
  103.                 for(int i=0; i<(int)h.size(); i++){
  104.                     result[i] = h[i].real();
  105.                 }
  106.             }else{
  107.                 for(int i=0; i<(int)h.size(); i++){
  108.                     result[i] = (R)floor(h[i].real()+0.5);
  109.                 }
  110.             }
  111.             return result;
  112.     }
  113.     template<class S>
  114.     typename enable_if<
  115.     	      is_same<R, S>::value &&
  116.              !is_same<S, complex<double>>::value &&
  117.              !is_same<S, double>::value &&
  118.              !is_same<S, int>::value &&
  119.              !is_same<S, long long>::value
  120.              ,polynomial<S>>::type
  121.     multiplication_impl(const polynomial<S>& another) const {
  122.         return slow_multiplication(another);
  123.     }
  124. public:
  125.     polynomial(const size_t &n = 0){
  126.         a = vector<R>(n);
  127.     }
  128.     polynomial(const vector<R> &coef){
  129.         a = coef;
  130.     }
  131.     size_t size() const{
  132.         return a.size();
  133.     }
  134.     void resize(const size_t &n){
  135.         a.resize(n);
  136.     }
  137.     R& operator [](const int &i){
  138.         assert(0<=i && i<(int)a.size());
  139.         return a[i];
  140.     }
  141.     const R& operator [](const int &i) const{
  142.         assert(0<=i && i<(int)a.size());
  143.         return a[i];
  144.     }
  145.     polynomial<R> operator *(const polynomial<R> &another) const{
  146.         return multiplication_impl(another);
  147.     }
  148. };
  149.  
  150. int main(){
  151.     polynomial<int> f(2), g(2);
  152.     f[1] = 1, f[0] = 1;
  153.     g[1] = 1, g[0] = -1;
  154.     polynomial<int> h = f*g;
  155.     for(int i=(int)h.size()-1; i>=0; i--){
  156.         printf("%d%c", h[i], i? ' ': '\n');
  157.     }
  158.     return 0;
  159. }
  160.  
Success #stdin #stdout 0s 4500KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
1 0 -1
https://ideone.com/PH4Sej
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!