import java.util.*;
import java.awt.*;
import javax.swing.*;
import java.awt.event.*;

class f1 extends f_x
{
  public double func(double x)
  {
  return Math.exp(x);	  
  }
}

public class NA63
{
    public static double [] lngamma(double x)
    {// logaritmik gamma function
     // xx[0]logarithm of gamma function
     // xx[1] sign of gamma function
	double xx[]=new double[2];    
    double result;
    double a;
    double b;
    double c;
    double p;
    double q;
    double u;
    double w;
    double z;
    int i;
    double logpi;
    double ls2pi;
    double tmp;
    double sgngam = 1;
    logpi = 1.14472988584940017414;
    ls2pi = 0.91893853320467274178;
    if( x<-34.0 )
    {
        q = -x;
        double yy[]=lngamma(q);
        w = yy[0];
        tmp=yy[1];
        p = Math.floor(q);
        i = (int)Math.round(p);
        if( i%2==0 )
        {
            sgngam = -1;
        }
        else
        {
            sgngam = 1;
        }
        z = q-p;
        if( z>0.5 )
        {
            p = p+1;
            z = p-q;
        }
        z = q*Math.sin(Math.PI*z);
        result = logpi-Math.log(z)-w;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    if( x<13 )
    {
        z = 1;
        p = 0;
        u = x;
        while(u>=3)
        {
            p = p-1;
            u = x+p;
            z = z*u;
        }
        while(u<2)
        {
            z = z/u;
            p = p+1;
            u = x+p;
        }
        if( z<0 )
        {
            sgngam = -1;
            z = -z;
        }
        else
        {
            sgngam = 1;
        }
        if( u==2 )
        {
            result = Math.log(z);
            //return result;
            xx[0]=result;
            xx[1]=sgngam;
            return xx;
        }
        p = p-2;
        x = x+p;
        b = -1378.25152569120859100;
        b = -38801.6315134637840924+x*b;
        b = -331612.992738871184744+x*b;
        b = -1162370.97492762307383+x*b;
        b = -1721737.00820839662146+x*b;
        b = -853555.664245765465627+x*b;
        c = 1;
        c = -351.815701436523470549+x*c;
        c = -17064.2106651881159223+x*c;
        c = -220528.590553854454839+x*c;
        c = -1139334.44367982507207+x*c;
        c = -2532523.07177582951285+x*c;
        c = -2018891.41433532773231+x*c;
        p = x*b/c;
        result = Math.log(z)+p;
        //return result;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    q = (x-0.5)*Math.log(x)-x+ls2pi;
    if( x>100000000 )
    {
        result = q;
        //return result;
        xx[0]=result;
        xx[1]=sgngam;
        return xx;
    }
    p = 1/(x*x);
    if( x>=1000.0 )
    {
        q = q+((7.9365079365079365079365*0.0001*p-2.7777777777777777777778*0.001)*p+0.0833333333333333333333)/x;
    }
    else
    {
        a = 8.11614167470508450300*0.0001;
        a = -5.95061904284301438324*0.0001+p*a;
        a = 7.93650340457716943945*0.0001+p*a;
        a = -2.77777777730099687205*0.001+p*a;
        a = 8.33333333333331927722*0.01+p*a;
        q = q+a/x;
    }
    xx[0]= q;
    xx[1]=sgngam;
    return xx;
    }
//====================

    public static double[][] gauss_jacobi_coefficients(int n,double alpha,double beta)
    {
        double eps=1.0e-30;
	    int i = 0;
        int j = 0;
        double r = 0;
        double r1 = 0;
        double t1 = 0;
        double t2 = 0;
        double t3 = 0;
        double p1 = 0;
        double p2 = 0;
        double p3 = 0;
        double pp = 0;
        double an = 0;
        double bn = 0;
        double a = 0;
        double b = 0;
        double c = 0;
        double tmpsgn = 0;
        double alfbet = 0;
        double temp = 0;
        int its = 0;
        double A[][]=new double[2][n];
        for(i=0; i<=n-1; i++)
        {   if( i==0 )
            {
                an = alpha/n;
                bn = beta/n;
                t1 = (1+alpha)*(2.78/(4+n*n)+0.768*an/n);
                t2 = 1+1.48*an+0.96*bn+0.452*an*an+0.83*an*bn;
                r = (t2-t1)/t2;
            }
            else
            {
                if( i==1 )
                {
                    t1 = (4.1+alpha)/((1+alpha)*(1+0.156*alpha));
                    t2 = 1+0.06*(n-8)*(1+0.12*alpha)/n;
                    t3 = 1+0.012*beta*(1+0.25*Math.abs(alpha))/n;
                    r = r-t1*t2*t3*(1-r);
                }
                else
                {
                    if( i==2 )
                    {
                        t1 = (1.67+0.28*alpha)/(1+0.37*alpha);
                        t2 = 1+0.22*(n-8)/n;
                        t3 = 1+8*beta/((6.28+beta)*n*n);
                        r = r-t1*t2*t3*(A[0][0]-r);
                    }
                    else
                    {
                        if( i<n-2 )
                        {
                            r = 3*A[0][i-1]-3*A[0][i-2]+A[0][i-3];
                        }
                        else
                        {
                            if( i==n-2 )
                            {
                                t1 = (1+0.235*beta)/(0.766+0.119*beta);
                                t2 = 1/(1+0.639*(n-4)/(1+0.71*(n-4)));
                                t3 = 1/(1+20*alpha/((7.5+alpha)*n*n));
                                r = r+t1*t2*t3*(r-A[0][i-2]);
                            }
                            else
                            {
                                if( i==n-1 )
                                {
                                    t1 = (1+0.37*beta)/(1.67+0.28*beta);
                                    t2 = 1/(1+0.22*(n-8)/n);
                                    t3 = 1/(1+8*alpha/((6.28+alpha)*n*n));
                                    r = r+t1*t2*t3*(r-A[0][i-2]);
                                }
                            }
                        }
                    }
                }
            }
            alfbet = alpha+beta;
            do
            {
                temp = 2+alfbet;
                p1 = (alpha-beta+temp*r)*0.5;
                p2 = 1;
                for(j=2; j<=n; j++)
                {
                    p3 = p2;
                    p2 = p1;
                    temp = 2*j+alfbet;
                    a = 2*j*(j+alfbet)*(temp-2);
                    b = (temp-1)*(alpha*alpha-beta*beta+temp*(temp-2)*r);
                    c = 2*(j-1+alpha)*(j-1+beta)*temp;
                    p1 = (b*p2-c*p3)/a;
                }
                pp = (n*(alpha-beta-temp*r)*p1+2*(n+alpha)*(n+beta)*p2)/(temp*(1-r*r));
                r1 = r;
                r = r1-p1/pp;
            }
            while( Math.abs(r-r1)>=eps*(1+Math.abs(r))*100 );
            A[0][i] = r;
            double yy[]=lngamma(alpha+n);
            double w1 = yy[0];
            tmpsgn=yy[1];
            yy=lngamma(beta+n);
            double w2 = yy[0];
            tmpsgn=yy[1];
            yy=lngamma(n+1);
            double w3 = yy[0];
            tmpsgn=yy[1];
            yy=lngamma(n+alfbet+1);
            double w4 = yy[0];
            tmpsgn=yy[1];
            A[1][i] = Math.exp(w1+w2-w3-w4)*temp*Math.pow(2, alfbet)/(pp*p2);
        }
   return A; 
   }


public static double gauss_jacobi_integral(f_x f_xnt,int n,double alpha,double beta)
{
//n : number of integral coefficients
// this routine first generates gauss legendre coefficients
// for [x1,x2] band 
// then calculates gauss hermite integral
double a[][]=new double[2][n];
a=gauss_jacobi_coefficients(n,alpha,beta);
System.out.println(Matrix.toString(a));
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*f_xnt.func(a[0][i]);}
return z;
}

      public static void main(String[] args)
      { 
	  f1 b1=new f1();
      System.out.println(gauss_jacobi_integral(b1,5,-0.5,-0.5));
      }
}
