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

class f1 extends f_x
{
  public double func(double x)
  {double C1=3.74177489e8;
   double C2=1.438769e4;
   return C1/(C2*C2*C2*C2)*x*x*x/(Math.exp(x)-1.0);
  }
}

public class integ_gauss_laguerre 
{
    public static double [] lngamma(double x)
    {//logaritmik hata fonksiyonu
    //xx[0] hata fonksiyonunun logaritmasının değerini 
    //xx[1] hata fonksiyonunun kendisinin işaretini verir
    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_laguerre_coefficients(int n,double alpha)
{   //Gauss-Laguerre quadratürü katsayılarını hesaplar
	double sgngam=1;
	//===========================================
    //end of internal method 
    double a[][]=new double[2][n];
    int i;
    int j;
    double r=0;
    double r1;
    double p1;
    double p2;
    double p3;
    double dp3;
    double tsg;
    double eps=1.0e-15;
    for(i = 0; i <= n-1; i++)
    {
        if( i==0 )
        {
            r = (1.0+alpha)*(3.0+0.92*alpha)/(1+2.4*n+1.8*alpha);
        }
        else
        {
            if( i==1 )
            {
                r = r+(15.0+6.25*alpha)/(1.0+0.9*alpha+2.5*n);
            }
            else
            {
                r = r+((1.0+2.55*(i-1))/(1.9*(i-1))+1.26*(i-1)*alpha/(1+3.5*(i-1)))/(1.0+0.3*alpha)*(r-a[0][i-2]);
            }
        }
        do
        {
            p2 = 0;
            p3 = 1;
            for(j = 0; j <= n-1; j++)
            {
                p1 = p2;
                p2 = p3;
                p3 = ((-r+2*j+alpha+1)*p2-(j+alpha)*p1)/(j+1);
            }
            dp3 = (n*p3-(n+alpha)*p2)/r;
            r1 = r;
            r = r-p3/dp3;
        }
        while(Math.abs(r-r1)>=eps*(1+Math.abs(r))*100);
        a[0][i] = r;
        double lng[]=lngamma((double)n);
        double lng1[]=lngamma((alpha+n));
        a[1][i] = -Math.exp(lng1[0]-lng[0])/(dp3*n*p2);
        tsg=lng[1];
    }    
return a;
}

public static double gauss_laquerre_integral(f_x f_xnt,double xa,int n,double alpha)
{
//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_laguerre_coefficients(n,alpha);
double z=0;
for(int i=0;i<n;i++)
{ z+=a[1][i]*f_xnt.func(a[0][i]+xa)*Math.exp(a[0][i])*Math.pow(a[0][i],-alpha);}
return z;
}

public static double gauss_laquerre_integral(f_x f_xnt,double xa,int n)
{
double alpha=2.0;	
double a[][]=new double[2][n];
a=gauss_laguerre_coefficients(n,alpha);
//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]+xa)*Math.exp(a[0][i])/(a[0][i]*a[0][i]);}
return z;
}

      public static void main(String[] args)
      { 
	  f1 b1=new f1();
      System.out.println("integral  : "+gauss_laquerre_integral(b1,0,100)+" tam = "+5.67051e-8);
      }
}
