#ifndef __FFT_H__
#define __FFT_H__

#include <vector>
#include <iterator>
#include <algorithm>

#ifdef _MSC_VER
  #define _USE_MATH_DEFINES
#endif

#include <cmath>

#include "param/complex.h"

class FFT{
  public:
    template <class T, class IndexType>
    static Complex<T> ft(const std::vector<T> &v, IndexType k){
      Complex<T> result(0);
      int n(0);
      for(typename std::vector<T>::const_iterator it(v.begin()); it != v.end(); ++it, n++){
        result += Complex<T>::exp(-2. * M_PI * n * k / v.size()) * *(it);
      }
      return result;
    }
    template <class T, class IndexType>
    static Complex<T> ft(const std::vector<Complex<T> > &v, IndexType k){
      Complex<T> result(0);
      int n(0);
      for(typename std::vector<Complex<T> >::const_iterator it(v.begin()); it != v.end(); ++it, n++){
        result += Complex<T>::exp(-2. * M_PI * n * k / v.size()) * *(it);
      }
      return result;
    }
    template <class T>
    static std::vector<Complex<T> > fft_ct(const std::vector<Complex<T> > &v){
      
      typedef std::vector<Complex<T> > V1;
      typedef std::vector<V1> V2;
      
      //まず何回2で割り切れるか数える
      int power2(0);
      typename V1::size_type size(v.size());
      if(size > 0){
        while(!(size & 0x1)){power2++; size >>= 1;}
      }

      //FFTをする前に並べ替え
      typename V1::size_type border(v.size());
      V2 sorted;
      sorted.push_back(v);
      for(int i = 0; i < power2; i++){
        border >>= 1;
        V2 temp_plus, temp_minus;
        for(typename V2::iterator element(sorted.begin()); element != sorted.end(); ++element){
          V1 sorted_plus, sorted_minus;
          for(typename V1::size_type j(0); j < border; j++){
            sorted_plus.push_back((*element)[j] + (*element)[j + border]);
            sorted_minus.push_back(((*element)[j] - (*element)[j + border]) * Complex<T>::exp(-2 * j * M_PI / element->size()));
          }
          temp_plus.push_back(sorted_plus);
          temp_minus.push_back(sorted_minus);
        }
        /*for(typename V2::iterator element(temp_plus.begin()); element != temp_plus.end(); ++element){
          for(typename V1::iterator element2(element->begin()); element2 != element->end(); ++element2){
            cout << "DEBUG:" << *element2 << endl;
          }
        }
        for(typename V2::iterator element(temp_minus.begin()); element != temp_minus.end(); ++element){
          for(typename V1::iterator element2(element->begin()); element2 != element->end(); ++element2){
            cout << "DEBUG:" << *element2 << endl;
          }
        }*/
        sorted.clear();
        for(typename V2::iterator it(temp_plus.begin()); it != temp_plus.end(); ++it){
          sorted.push_back(*it);
        }
        for(typename V2::iterator it(temp_minus.begin()); it != temp_minus.end(); ++it){
          sorted.push_back(*it);
        }
        
        /*for(typename V2::iterator element(sorted.begin()); element != sorted.end(); element++){
          for(typename V1::iterator element2(element->begin()); element2 != element->end(); element2++){
            cout << "DEBUG:" << *element2 << endl;
          }
        }*/
      }
      
      // FFTする
      V1 result(v.size());
      int base(0);
      for(typename V1::size_type j(0); j < v.size() / sorted.size(); j++){
        int i(0);
        int k(j * (1 << power2));
        for(typename V2::iterator element(sorted.begin()); element != sorted.end(); ++element, i++){    
          result[base + i] = ft(*element, k);
          //cout << "base + i: " << (base + i) << "; val = " << ft(*element, k) << endl;
          //cout << "i:" << i << "; j:" << j << "; k: " << k << ";" << endl;
        }
        base += (1 << (power2 * 2));
        base %= v.size();
      }
      
      return result;
    }
    template <class T>
    static std::vector<Complex<T> > fft_ct(const std::vector<T> &v){
      std::vector<Complex<T> > casted;
      for(typename std::vector<T>::const_iterator it(v.begin()); it != v.end(); ++it){
        casted.push_back(Complex<T>(*it));  
      }
      return fft_ct(casted);
    }
    template <class T>
    static Complex<T> ift_no_divide(const std::vector<Complex<T> > &v, int n){
      Complex<T> result(0);
      int k(0);
      for(typename std::vector<Complex<T> >::const_iterator it(v.begin()); it != v.end(); ++it, k++){
        result += Complex<T>::exp(2. * M_PI * n * k / v.size()) * *(it);
      }
      return result;
    }
    template <class T>
    static Complex<T> ift(const std::vector<Complex<T> > &v, int n){
      return ift_no_divide(v, n) / v.size();
    }
    template <class T>
    static std::vector<Complex<T> > ifft_ct(const std::vector<Complex<T> > &v){
      
      typedef std::vector<Complex<T> > V1;
      typedef std::vector<V1> V2;
      
      //まず何回2で割り切れるか数える
      int power2(0);
      typename V1::size_type size(v.size());
      if(size > 0){
        while(!(size & 0x1)){power2++; size >>= 1;}
      }

      //IFFTをする前に並べ替え
      typename V1::size_type border(v.size());
      V2 sorted;
      sorted.push_back(v);
      for(int i = 0; i < power2; i++){
        border >>= 1;
        V2 temp_plus, temp_minus;
        for(typename V2::iterator element(sorted.begin()); element != sorted.end(); element++){
          V1 sorted_plus, sorted_minus;
          for(typename V1::size_type j(0); j < border; j++){
            sorted_plus.push_back((*element)[j] + (*element)[j + border]);
            sorted_minus.push_back(((*element)[j] - (*element)[j + border]) * Complex<T>::exp(2 * j * M_PI / element->size()));
          }
          temp_plus.push_back(sorted_plus);
          temp_minus.push_back(sorted_minus);
        }
        sorted.clear();
        for(typename V2::size_type j(0); j < temp_plus.size(); j++){sorted.push_back(temp_plus[j]);}
        for(typename V2::size_type j(0); j < temp_minus.size(); j++){sorted.push_back(temp_minus[j]);}
      }
      
      // IFFTする
      V1 result(v.size());
      int base(0);
      for(typename V1::size_type j(0); j < v.size() / sorted.size(); j++){
        int i(0);
        int k(j * (1 << power2));
        for(typename V2::iterator element(sorted.begin()); element != sorted.end(); ++element, i++){    
          result[base + i] = ift_no_divide(*element, k) / v.size();
          //cout << "base + i: " << (base + i) << "; val = " << ft(*element, k) << endl;
          //cout << "i:" << i << "; j:" << j << "; k: " << k << ";" << endl;
        }
        base += (1 << (power2 * 2));
        base %= v.size();
      }
      
      return result;
    }
    
    template <class T>
    inline static std::vector<Complex<T> > fft(const std::vector<Complex<T> > &v){
      return fft_ct(v);
    }
    template <class T>
    inline static std::vector<Complex<T> > fft(const std::vector<T> &v){
      return fft_ct(v);
    }
    
    template <class T>
    inline static std::vector<Complex<T> > ifft(const std::vector<Complex<T> > &v){
      return ifft_ct(v);
    }
};

#endif /* __FFT_H__ */