#include <iostream>
#include <cmath>
#include <vector>
#include <complex>
#include <iterator>
#include <algorithm>
#include <fstream>
#include <valarray>
#include "gsl_fft.hpp"
#include <iomanip>

using namespace std;
typedef complex<double> dtype ;
void valarray_out(ostream& stream, valarray<dtype> v){
  cout<<v.size()<<endl;
  stream<<setprecision(8);
  stream<<fixed;// stream<<scientific;
  for(unsigned i=0; i<v.size(); i++){
    stream<<i<<" "<<setw(12)<< v[i].real()<<" "<< setw(12)<<v[i].imag()<<endl;
  }
}

int main (void){
  const int n=128;
  const int nones= 10; // Keep less than n!
  valarray<dtype> data(n);
  ofstream myfile("result");
  //
  // Symmetric pulse:  1111111111100...001111111111
  //                   11 ones            10 ones
  data = 0.0;
  data[slice(0,nones+1,1)] =+ 1.0;
  data[slice(n-nones,nones,1)] =+ 1.0;

  cout<<"before fft :\n";
  valarray_out(cout,data);
  valarray_out(myfile,data);
  myfile<<endl;
  // ----------------------
  // forward FFT
  int status = fft(data,1);
  // ----------------------
  if(status!=0){
    cout << "fft fails\n";
    return 1;
  }
  cout<<"after forward fft (divided by sqrt(n)) :\n";
  data /= sqrt(data.size());
  valarray_out(cout,data);
  valarray_out(myfile,data);

  // ----------------------
  // backward FFT
  status = fft(data,-1);
  // ----------------------
  if(status!=0){
    cout << "fft fails\n";
    return 1;
  }

  cout<<"after backward fft (divided by sqrt(n)) :\n";
  data /= sqrt(data.size());
  valarray_out(cout,data);

  myfile.close();
  return 0;
}