#ifndef LAGUERREBUILDER_SEEN
#define LAGUERREBUILDER_SEEN
/*
 *  This file is part of LagSHT.
 *
 *  LagSHT is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation; either version 2 of the License, or
 *  (at your option) any later version.
 *
 *  LagSHT is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with libsharp; if not, write to the Free Software
 *  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
 */

/*
 *  LagSHT is being developed at the Linear Accelerateur Laboratory (LAL)
 *  91898 ORSAY CEDEX - FRANCE
 *  main author: J.E Campagne
 */

#include <cmath> //very important for abs definition
#include <vector>
#include <iostream>
#include <iomanip>      // std::setprecision

//JEC 12/1/16: alpha becomes double instead of int
//JEC (N+alpha)! -> Gamma[N+alpha+1]. Use Boost::tgamma function
#include "boost/math/special_functions/gamma.hpp"

#include "lagsht_numbers.h"

using namespace std;


#define DEBUG 0

namespace LagSHT {


/////////////////////////////////////////////////////////////////////////////////////////////
//// - -Class generalized Laguerre Function and Gauss-Laguerre Quadrature --------------   //
////   The Laguerre Function is scaled by exp(-r/2) compared to the Laguerre Polynomials   //
////   the exp(r_i) scales the weights compared to the Polynomial Gauss-Laguerre Quadrature//
//// (nb. an public heritance from LaguerrePolyQuad would made confusion)                  //
/////////////////////////////////////////////////////////////////////////////////////////////

class LaguerreFuncQuad   {
public:
  //! Creator
  LaguerreFuncQuad(int N, r_8 alpha=2): N_(N), alpha_(alpha),
    ln2(0.693147180559945309417232121458176),
    inv_ln2(1.4426950408889634073599246810018921),
    fsmall(std::ldexp(1.,-large_exponent2)),
    fbig(std::ldexp(1.,large_exponent2)),
    xi(N+2), cf(maxscale+1-minscale) {

    using namespace std;
    xi[0]=0.;
    for (size_t i=1; i<xi.size(); ++i)
      xi[i]=1./i;
    for (size_t i=0; i<cf.size(); ++i)
      cf[i] = ldexp(1.,(int(i)+minscale)*large_exponent2);
  }

  //! Destructor
  virtual ~LaguerreFuncQuad() {
  }
   
  //! Getters
  int GetOrder() const { return N_; } 
  r_8 GetAlpha() const { return alpha_; }

  //! Compute Weights and Nodes of the Gauss-Laguerre quadrature
  void QuadWeightNodes(vector<r_8>& nodes, vector<r_8>& weights) const;

  //! Compute  L_n^alpha(x) * exp(-x/2) * scale with  n = N_ 
  r_8 Value(r_8 x, r_8 scale=1.0) const;
  
  //! Compute  L_n^alpha(x) * exp(-x/2) * scale with  n=0,...,N_ (included)
  void Values(r_8 x, vector<r_8>& vals, r_8 scale=1.0);  

protected:

  /*! returns the amount by which x must be shifted to go closer to a root of
    L_N^alpha(x)
  */
  r_8 newtoncorr (r_8 x, int N) const {
    using namespace std;
    if(N<1) return 0.;
    
    r_8 pnm1=0., pn=1.;
    r_8 c1=alpha_-1., c2=alpha_-1.;
    for(int i=1;i<=N;i++)
      {
        r_8 pnm2 = pnm1;
        pnm1 = pn;
        c1+=2.;
        c2+=1.;
        pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
        // if necessary, scale down
        // the function result is independent of the absolute scale 
        while (abs(pn)>=fbig)
          { pn*=fsmall; pnm1*=fsmall; }
      }
    r_8 ppn=(N*pn-(N+alpha_)*pnm1)/x;
    return -pn/ppn;
   }//newtoncorr

  //JEC 12/1/16 add a normalization factor 
  /*! returns L_N^alpha(x)*exp(-x/2)*norm
   */
  r_8 value(r_8 x, int N, r_8 norm = 1.0) const {
    using namespace std;
    if (N<0) return 0;

    r_8 log2expfac = -0.5*x*inv_ln2;
    int scale=int(log2expfac/large_exponent2);
    log2expfac-=scale*large_exponent2;

    r_8 pnm1=0., pn=exp(ln2*log2expfac)*norm;
    r_8 c1=alpha_-1., c2=alpha_-1.;
    for(int i=1;i<=N;i++)
      {
        r_8 pnm2 = pnm1;
        pnm1 = pn;
        c1+=2.;
        c2+=1.;
        pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
        while (abs(pn)>=fbig) // if necessary, scale down and count rescalings
          { pn*=fsmall; pnm1*=fsmall; ++scale; }
      }
    return (scale<minscale) ? 0 : pn*cf[scale-minscale];
  }//value

  //JEC 12/1/16 add a normalization factor 
  /*! returns L_n^alpha(x)*exp(-x/2)*norm [n=0..N] in res.
   */
  void values (r_8 x, int N, std::vector<r_8> &res, r_8 norm=1.0) {
    using namespace std;
    if (N<0) { res.resize(0); return; }
    res.resize(N+1);
    res[0]=exp(-0.5*x)*norm;

    r_8 log2expfac = -0.5*x*inv_ln2;
    int scale=int(log2expfac/large_exponent2);
    log2expfac-=scale*large_exponent2;

    r_8 pnm1=0., pn=exp(ln2*log2expfac)*norm;
    r_8 c1=alpha_-1., c2=alpha_-1.;
    for(int i=1;i<=N;i++)
      {
        r_8 pnm2 = pnm1;
        pnm1 = pn;
        c1+=2.;
        c2+=1.;
        pn = ((c1-x)*pnm1 - c2*pnm2)*xi[i];
        while (abs(pn)>=fbig) // if necessary, scale down and count rescalings
          { pn*=fsmall; pnm1*=fsmall; ++scale; }
        res[i] = (scale<minscale) ? 0 : pn*cf[scale-minscale];
      }
  }//values

  /*! computes nodes and weights for Gauss-Laguerre integration of order N
    NOTE: the weights are correct for modified Laguerre polynomials, i.e.
    L_N^alpha(x)*exp(-x/2)
    The nodes/weights are indexed i:0,...,N-1
  */
    void nodes_and_weights(int N, std::vector<r_8> &nodes,
			   std::vector<r_8> &weights) const {
    using namespace std;
    nodes.resize(N);
    weights.resize(N);

    /* JEC 12/1/16 alpha is r_8. tmp1 = (N+alpha)!/N!/(N+1)^2 = Gamma(N+alpha+1)/Gamma(N+1)/(N+1)^2
      r_8 tmp1 = N+alpha_;
      for (int i=1;i<alpha_-1;i++) tmp1 *= (tmp1-1.);
      tmp1 /= (N+1.);

      Use the ratio of Gamma functions as in Boost library: Todo see if one should tune the Policy
    */

    r_8 tmp1 = boost::math::tgamma_ratio((r_8)N+alpha_+1.0, r_8(N)+1.0);
    tmp1 /= (N+1.)*(N+1.);

    r_8 x=0.;
    // Approximation by Stroud & Secrest, Gaussian quadrature formulas, 1966
    for(int k=0; k<N; k++) { 
      if (k==0)
	x = (1.0+alpha_)*(3.+0.92*alpha_)/(1.+2.4*N+1.8*alpha_);
      else if (k==1)
	x += (15.+6.25*alpha_)/(1.+0.9*alpha_+2.5*N);
      else
	{
	  int km1 = k-1;
	  r_8 c0 = (1.+2.55*km1)/(1.9*km1)+1.26*km1*alpha_/(1.+3.5*km1);
	  r_8 c1 = 1.+0.3*alpha_;
	  x += c0*(x-nodes[k-2])/c1;
	}
      
      // Refinement using Newton's method
      // Once we notice that we are close, do one more step and then exit
      bool done=false;
      while (true)
	{
	  r_8 xold=x;
	  x+=newtoncorr (x, N);
	  
	  if (done) break;
	  if (abs(x-xold)<max(1e-12,1e-12*abs(x))) done=true;
	}
      
      nodes[k] = x;
      
      r_8 tmp0 = value(x, N+1);
      weights[k] = x*tmp1/(tmp0*tmp0);
    }//end for

#if DEBUG >= 1  
    //verification 
    r_16 sumR = 0.;
    r_16 sumW = 0.;
    r_16 wip;
    for (int i=0;i<N;i++){
      
      cout << "root["<<i<<"]: " << setprecision(16) << nodes[i] << " LagN(ri) = " << Value(nodes[i]) << endl; 

      sumR += nodes[i];
      wip = exp(log(abs(weights[i]))-nodes[i]); if(weights[i]<0) wip *= -1.0;
      sumW += wip;
      //    sumW += weights[i]*exp(-nodes[i]/tau);
    }

    /* JEC 12/1/16 alpha ! = alphaFact = Gamma(alpha+1)
      int alphaFact = 1;
      for(int i=1;i<=alpha_; i++) alphaFact *= i;
      r_16 sumRTh = (r_16)(N*(N+alpha_));
      r_16 sumWTh = (r_16)(alphaFact);
    */
    
    r_16 sumRTh = (r_16)(N*(N+alpha_));
    r_16 sumWTh = boost::math::tgamma((r_16)(alpha_+1.0));


    cout << "Sum roots = " << sumR << " -> diff with theory: " << sumR - sumRTh << endl;
    cout << "Sum weights = " << sumW << " -> diff with theory: " << sumW - sumWTh << endl;
    
    if (fabs(sumR-sumRTh)>(0.000001*sumRTh)) {
      cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
	   << "WARNING!!!! verify LaguerrePolyQuad sum of nodes not accurate: "
	   << " current value " << sumR << " theoretical value : " << sumRTh 
	   << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
	   << endl;
    }
    
    if (fabs(sumW-sumWTh)>(0.000001*sumWTh)) {
      cout << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
	   << "WARNING!!!! verify LaguerrePolyQuad sum of weights not accurate: "
	   << " current value " << sumW << " theoretical value : " << sumWTh 
	   << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!\n" 
	   << endl;
    }
#endif
    
  }//nodes_and_weights

protected:
  int N_; //!< Order of the polynomial
  r_8 alpha_; //!< second parameter of the generalized Laguerre polynomial 

private:
  enum { large_exponent2=90, minscale=-14, maxscale=14 };
  r_8 ln2, inv_ln2, fsmall, fbig;
  std::vector<r_8> xi; 
  std::vector<r_8> cf; 


}; //LaguerreFuncQuad



} // Fin du namespace

#undef DEBUG

#endif //LAGUERREBUILDER_SEEN