#include <cmath>
#include <iostream>
#include <algorithm> // for copy
#include <iterator> // for ostream_iterator

#ifdef CBLAS
#ifdef Darwin
#include <Accelerate/Accelerate.h>
#elif Linux
#include <cblas.h>
#endif
#endif

using namespace std;

#include "laguerreTransform.h"
#include "lagsht_exceptions.h"

#define DEBUG 1

namespace LagSHT {


LaguerreTransform::LaguerreTransform(int N, r_8 R, int alpha) : 
  BaseLaguerreTransform(N,alpha), R_(R) {
  
#if DEBUG >= 1
  cout << "LaguerreTransform start...." <<endl;
#endif

  if(alpha<0)  throw LagSHTError("LagTransform call with alpha<0");
  LaguerreFuncQuad lag(N_,alpha_);
  lag.QuadWeightNodes(nodes_,weights_);


  int alphaFact = 1;
  for(int i=1;i<=alpha_; i++) alphaFact *= i;

 
  alphaFact_ = sqrt((r_8)alphaFact); //sqrt(alpha!)

}//Ctor

int LaguerreTransform::R2Index(r_8 r) const {
  using namespace std;
  r_8 rscaled = r*nodes_[N_-1]/R_;
  vector<r_8> dist(N_);
  transform(nodes_.begin(),nodes_.end(),dist.begin(), distance(rscaled));
  return min_element(dist.begin(),dist.end()) - dist.begin();
}//RadIndex


void LaguerreTransform::MultiAnalysis(const vector< complex<r_8> >& fi, 
				      vector< complex<r_8> >& fn, int stride) {

  if(fi.size() != (size_t)(N_*stride) ) throw LagSHTError("LagTransform::Analysis size error");
    fn.resize(N_*stride);

#if DEBUG >= 1
    cout << "Analysis with Function...." << endl;
#endif
#if DEBUG >= 2
    std::copy(fi.begin(), fi.end(), std::ostream_iterator<char>(std::cout, " "));
#endif                     
    r_8 invalphaFact = 1.0/alphaFact_; //1/alpha!                   

    vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
    facts[0] = invalphaFact;
    for(int n=1;n<N_;n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
    
    LaguerreFuncQuad lag(N_-1);
    vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1

    vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
    for (int k=0; k<N_; k++){
      r_8 rk = nodes_[k];
      lag.Values(rk,LnAll);
      for (int n = 0; n<N_; n++ ){
	//	LnkMtx[n*N_+k] = LnAll[n]*facts[n];
	LnkMtx[n+N_*k] = LnAll[n]*facts[n]*weights_[k];
      }
    }



#ifdef CBLAS
    cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (double*)(&fi[0]), 2*stride, &LnkMtx[0], N_, 0, (double *)(&fn[0]), 2*stride);
#else
    vector<complex<r_8> > vtmp(N_);
    for (int l=0; l<stride; l++) {
      vtmp.assign(N_,0.);
      for (int i = 0; i<N_; i++ ){
	complex<r_8> fli = fi[l+ i*stride];
	r_8 wi = weights_[i];
	for (int n=0; n<N_; n++){
	  //	  vtmp[n] +=  fli * wi * LnkMtx[n*N_+i];
	  vtmp[n] +=  fli * LnkMtx[n+N_*i];
	}//loop on k
      }//loop on n
      for (int n=0; n<N_; n++) {
	fn[l+n*stride] += vtmp[n];
      }
    }//loop on l
#endif

}//MultiAnalysis


void LaguerreTransform::MultiSynthesis(const vector< complex<r_8> >& fn, 
				       vector< complex<r_8> >& fi, int stride) {

  if(fn.size() != (size_t)(N_*stride) ) throw LagSHTError("LagTransform::Synthesis size error");
    fi.resize(N_*stride);

#if DEBUG >= 1
    cout << "Multi Synthesis with Function....:" << N_ << endl;
#endif
#if DEBUG >= 2
    std::copy(fn.begin(), fn.end(), std::ostream_iterator<char>(std::cout, " "));
#endif
    r_8 invalphaFact = 1.0/alphaFact_; //1/alpha! 

    vector<r_8> facts(N_); //sqrt(n!/(n+alpha)!) en iteratif
    facts[0] = invalphaFact;
    for(int n=1; n<N_; n++) facts[n] = facts[n-1]*sqrt(((r_8)n)/((r_8)(n+alpha_)) );
    

    LaguerreFuncQuad lag(N_-1);
    vector<r_8> LnAll(N_); //all the values Lag_n(r) n:0...N-1
    vector<r_8> LnkMtx(N_*N_); //all the values Lag_n(r) n:0...N-1 TRY r_16->r_8
    for (int k=0; k<N_; k++){
      r_8 rk = nodes_[k];
      lag.Values(rk,LnAll);
      for (int n = 0; n<N_; n++ ){
	LnkMtx[n*N_+k] = LnAll[n]*facts[n];
      }
    }
    
#if CBLAS
    cblas_dgemm (CblasColMajor, CblasNoTrans, CblasTrans, 2*stride, N_, N_, 1., (double*)(&fn[0]), 2*stride, &LnkMtx[0], N_, 0, (double *)(&fi[0]), 2*stride);
#else
    vector<complex<r_8> > vtmp(N_);
    for (int l=0; l<stride; l++) {
      vtmp.assign(N_,0.);
      for (int n = 0; n<N_; n++ ){
	complex<r_8> fln = fn[l+ n*stride];
	for (int k=0; k<N_; k++){
	  vtmp[k] +=  fln * LnkMtx[n*N_+k];
	}//loop on k
      }//loop on n
      for (int k=0; k<N_; k++) fi[l+k*stride] += vtmp[k];
    }//loop on l
#endif
}//MultiSynthesis



}//end namespace