//Evan Warner
	//Class computes and stores point and weight arrays for Gauss-Jacobi quadrature 
	//Routine from Numerical Recipes in C, Second Edition
   //Initial approximations for roots given by Stroud and Secrest
   //Program continues with Newton's Method, using relation given by a polynomial and the polynomial
   //of one lower order to derive the derivative
	//Only provides weights/points for range [-1,1]
	
	//See MathWorld for more info including error term
   
    public class GaussJacobiWeights
   {
      private static double myAlpha, myBeta;
      private static double[] myPoints, myWeights;
      private static int myN;
      private static final int MAXIT=10;
      private static final double EPS=3.0*Math.pow(10,-14); //Floating-point tolerance of system
       public GaussJacobiWeights(double al, double be, int n)
      {
         myN=n;
         myAlpha=al;
         myBeta=be;
         myPoints = new double[myN];
         myWeights = new double[myN];
         double r1,r2,r3,z=0;
         for(int i=0;i<myN;i++) //Computes initial guesses for the first and last roots
         {  
            if(i==0)
            {
               double an=myAlpha/((double)myN);
               double bn=myBeta/((double)myN);
               r1=(1.0+myAlpha)*(2.78/(4.0+myN*myN)+0.768*an/((double)myN));
               r2=1.0+1.48*an+0.96*bn+0.452*an*an+0.83*an*bn;
               z=1.0-r1/r2;
            }
            else if(i==1)
            {
               r1=(4.1+myAlpha)/((1.0+myAlpha)*(1.0+0.156*myAlpha));
               r2=1.0+0.06*(myN-8.0)*(1.0+0.12*myAlpha)/((double)myN);
               r3=1.0+0.012*myBeta*(1.0+0.25*Math.abs(myAlpha))/((double)myN);
               z-=(1.0-z)*r1*r2*r3;
            }
            else if(i==2)
            {
               r1=(1.67+0.28*myAlpha)/(1.0+0.37*myAlpha);
               r2=1.0+0.22*(myN-8.0)/((double)myN);
               r3=1.0+8.0*myBeta/((6.28+myBeta)*myN*myN);
               z-=(myPoints[0]-z)*r1*r2*r3;
            }
            else if(i==(myN-2))
            {
               r1=(1.0+0.235*myBeta)/(0.766+0.119*myBeta);
               r2=1.0/(1.0+0.639*(myN-4.0)/(1.0+0.71*(myN-4.0)));
               r3=1.0/(1.0+20.0*myAlpha/((7.5+myAlpha)*myN*myN));
               z+=(z-myPoints[myN-4])*r1*r2*r3;
            }
            else if(i==(myN-1))
            {
               r1=(1.0+0.37*myBeta)/(1.67+0.28*myBeta);
               r2=1.0/(1.0+0.22*(myN-8.0)/myN);
               r3=1.0/(1.0+8.0*myAlpha/((6.28+myAlpha)*myN*myN));
               z+=(z-myPoints[myN-3])*r1*r2*r3;
            }
            else
               z=3.0*myPoints[i-1]-3.0*myPoints[i-2]+myPoints[i-3];
            double ab=myAlpha+myBeta;
            int its;
            double temp=0,pp=0,p2=0;
            for(its=1;its<=MAXIT;its++) //Iterates Newton's method to find roots of Jacobi polynomials
            {
               temp=2.0+ab;
               double p1=(myAlpha-myBeta+temp*z)/2.0;
               p2=1.0;
               double p3;
               for(int j=2;j<=myN;j++)
               {
                  p3=p2;
                  p2=p1;
                  temp=2*j+ab;
                  double a=2*j*(j+ab)*(temp-2.0);
                  double b=(temp-1.0)*(myAlpha*myAlpha-myBeta*myBeta+temp*(temp-2.0)*z);
                  double c=2.0*(j-1.0+myAlpha)*(j-1.0+myBeta)*temp;
                  p1=(b*p2-c*p3)/a;
               }
               pp=(myN*(myAlpha-myBeta-temp*z)*p1+2.0*(myN+myAlpha)*(myN+myBeta)*p2)/(temp*(1.0-z*z));
               double z1=z;
               z=z1-p1/pp;
               if(Math.abs(z-z1)<=EPS)
                  break;
            }
            if(its>MAXIT)
               System.out.println("SOMETHING'S FISHY"); //Should never get here under problem assumptions
            myPoints[i]=z;
            myWeights[i]=Math.exp(gammln(myAlpha+myN)+gammln(myBeta+myN)-gammln(myN+1.0)-gammln(myN+ab+1.0))*temp*Math.pow(2.0,ab)/(pp*p2);
         }
      }
      
   	//Returns the natural logarithm of the gamma function applied to a real number x
      //Routine from Numerical Recipes in C, Second Edition, due to formula derived by Lanczos
      //Error is bounded by |e|<2*10e-10, not including roundoff error
       public double gammln(double x)
      {
         final double[] cof={76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.001208650973866179,-0.000005395239384953};
         double y=x;
         double temp=x+5.5;
         temp-=(x+0.5)*Math.log(temp);
         double ser=1.000000000190015;
         for(int j=0;j<=5;j++)
         {
            y++;
            ser+=cof[j]/y;
         }
         return -temp+Math.log(2.5066282746310005*ser/x);
      }
      
   	//Returns weight array
       public double[] getWeights()
      {
         return myWeights;
      }
      
   	//Returns point array
       public double[] getPoints()
      {
         return myPoints;
      }
   }