#include<math.h>

/*

   Joel Heinrich
   February 10 2005

Returns Gauss-Laguerre quadrature abscissas and log(weights) which can
be used to approximate

      integral u=0 to infinity pow(u,alpha)*exp(-u)*f(u) du
as
      sum k=0 to n-1  exp(lw[k])*f(x[k])

or equivalently

      sum k=0 to n-1  exp(lw[k]+log(f(x[k])))

The quadrature is exact for polynomial f of degree 2n-1 or less.

*/




void gausslaguerre(double x[],double lw[],int n,double alpha){
  const int nshift = 20;
  const double shift = 1<<nshift, rshift=1/shift;
  int i;
  double z=0;
  
  for(i=0;i<n;++i) {
    int j=0, k=2, nscale=0;
    double dz=0.0, p1=0, p2=0;
    if(i==0) {
      z=(1.0+alpha)*(3.0+0.92*alpha)/(1.0+2.4*n+1.8*alpha);
    } else if(i==1) {
      z += (15.0+6.25*alpha)/(1.0+2.5*n+0.9*alpha);
    } else if(i==2) {
      const double ai=i-1;
      z += ( (1.0+2.55*ai)/(1.9*ai) + 1.26*ai*alpha/(1.0+3.5*ai) )*
	(z-x[i-2])/(1.0+0.3*alpha);
    } else if(i==3) {
      z = 3.0*(x[2]-x[1])+x[0];
    } else if(i==4) {
      z = 4.0*x[3] - 6.0*x[2] + 4.0*x[1] - x[0];
    } else if(i==5) {
      z = 5.0*x[4] - 10.0*x[3] + 10.0*x[2] - 5.0*x[1] + x[0];
    } else {
      z = 6.0*x[i-1] - 15.0*x[i-2] + 20.0*x[i-3] -
	15.0*x[i-4] + 6.0*x[i-5] - x[i-6];
    }
    while(k>0) {
      p1=1;
      p2=0;
      nscale=0;
      z -= dz;
      for(j=1;j<=n;++j){
	const double p3=p2;
	p2=p1;
	p1=((2*j-1+alpha-z)*p2 - (j-1+alpha)*p3)/j;
	if(fabs(p2)>shift) {
	  ++nscale;
	  p1 *= rshift;
	  p2 *= rshift;
	}
      }
      dz = p1*z/(n*p1-(n+alpha)*p2);
      if(fabs(dz)<1.0e-10*z)--k;
    }
    x[i]=z;
    lw[i] = log(z/(p2*p2)) - 2*nshift*nscale*M_LN2 ;
  }
  
  {
    double t = 0.0;
    for(i=n-1;i>=0;--i)
      t += exp(lw[i]);
    t = lgamma(alpha+1)-log(t);
    for(i=0;i<n;++i)
      lw[i] += t;
  }

  return;
}