/*
 * General Math package that includes
 * functions like Bessel, error, gamma, Legendre, Chebychev etc.
 */
/**
 *
 * @author turhan_coban
 */
import static java.lang.Math.*;
import java.io.*;
import java.util.*;
import javax.swing.*;
import javax.swing.table.*;
import java.awt.*;
    class data
	{
	public double x[];	
	public data(double[] xi)
	{
	int n=xi.length;
	x=new double[n];
	for(int i=0;i<n;i++) {x[i]=xi[i];}
	}	
	public String toString() {return ""+x[0]+" "+x[1];}
	}

	
public class Mathd 
{   //some of the turkish equivalents are also included
   //The double value that is closer than any other to pi, the ratio of the circumference of a circle to its diameter. 
	public static double PI=Math.PI;
   
	//The double value that is closer than any other to e, the base of the natural logarithms.
   public static double E=Math.E;
   
   //C Euler constant
   public static double C=0.577215664901532860606512090082402431;
   
   //K Catalan's constant
   public static double K=0.915798391644461;

   public static double[] setX(double xmin,double dx,double xmax)
   {
    int n=(int)((xmax-xmin)/dx+0.00000000000001);
    double x[]=new double[n+1];
    for(int i=0;i<=n;i++)
    {x[i]=xmin+dx*i;}
    return x;	   
   }
   public static double[] setX(int xmin,int xmax)
   {
    int n=xmax-xmin;
    double x[]=new double[n+1];
    for(int i=0;i<=n;i++)
    {x[i]=xmin+i;}
    return x;	   
   }
   //Returns the absolute value of a double value.
   public static double abs(double a)
   {return Math.abs(a);}
   public static double[] abs(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.abs(a[i]);}
	   return b;
   }
   //Returns the absolute value of a float value.
   public static float abs(float a)
   {return Math.abs(a);}
   public static float[] abs(float a[])
   {  int n=a.length;
      float b[]=new float[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.abs(a[i]);}
	  return b;
   }
   //Returns the absolute value of a double value, similar to c function.
   public static double fabs(double a)
   {return Math.abs(a);}
   public static double[] fabs(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.abs(a[i]);}
	   return b;
   }
   //Returns the absolute value of a float value, similar to c function.
   public static float fabs(float a)
   {return Math.abs(a);}
   public static float[] fabs(float a[])
   {  int n=a.length;
      float b[]=new float[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.abs(a[i]);}
	  return b;
   }
   //Returns the absolute value of a long value.
   public static long abs(long a)
   {return Math.abs(a);}
   public static long[] abs(long a[])
   {  int n=a.length;
      long b[]=new long[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.abs(a[i]);}
	  return b;
   }
   //Returns the absolute value of an integer value
   public static int abs(int a)
   {return Math.abs(a);}
   public static int[] abs(int a[])
   {  int n=a.length;
      int b[]=new int[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.abs(a[i]);}
	  return b;
   }
   // sine of a radian
   public static double sin(double a)
   {return Math.sin(a);}
   
   public static double[] sin(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.sin(a[i]);}
	   return b;
   }
   //sine hyperbolic function
   public static double sinh(double a)
   {return Math.sinh(a);}
   public static double[] sinh(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.sinh(a[i]);}
	   return b;
   }
   //cosine hyperbolic function
   public static double cosh(double a)
   {return Math.cosh(a);}
   public static double[] cosh(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.cosh(a[i]);}
	   return b;
   }
   //tangent hyperbolic function
   public static double tanh(double a)
   {return Math.tanh(a);}
   public static double[] tanh(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.tanh(a[i]);}
	   return b;
   }
   //cosine function
   public static double cos(double a)
   {return Math.cos(a);}
   public static double[] cos(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.cos(a[i]);}
	   return b;
   }
     //tangent function
   public static double tan(double a)
   {return Math.tan(a);}
   public static double[]tan(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.tan(a[i]);}
	   return b;
   }
   //inverse sin function
   public static double asin(double a)
   {return Math.asin(a);}
   public static double[] asin(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.asin(a[i]);}
	   return b;
   }
   //inverse cosine function
   public static double acos(double a)
   {return Math.asin(a);}
   public static double[] acos(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.acos(a[i]);}
	   return b;
   }
   //inverse tangent function
   public static double atan(double a)
   {return Math.atan(a);}
   public static double[] atan(double a[])
   {  int n=a.length;
      double b[]=new double[n];
      for(int i=0;i<n;i++)
      {b[i]=Math.atan(a[i]);}
	   return b;
   }
   //inverse tangent function with x and y cartesian input
   public static double atan2(double a,double b)
   {return Math.atan2(a,b);}
   public static double[] atan2(double a[],double b[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.atan2(a[i],b[i]);}
	   return b;
   }

   //square root
   public static double sqrt(double a)
   {return Math.sqrt(a);}
   
   public static double[] sqrt(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.sqrt(a[i]);}
	   return c;
   }
   public static double karekok(double a)
   {return Math.sqrt(a);}
   
   public static double[] karekok(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.sqrt(a[i]);}
	   return c;
   }
   //cube root    
   public static double cbrt(double a)
   {return Math.cbrt(a);}
      public static double[] cbrt(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.cbrt(a[i]);}
	   return c;
   }

   public static double kubkok(double a)
   {return Math.cbrt(a);}
   public static double[] kubkok(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.cbrt(a[i]);}
	   return c;
   }
   
   // natural logarithm
   public static double ln(double a)
   {return Math.log(a);}
   
   public static double[] ln(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.log(a[i]);}
	   return c;
   }
   // natural logarithm
   public static double log(double a)
   {return Math.log(a);}
   public static double[] log(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.log(a[i]);}
	   return c;
   }
   // 10 based logarithm
   public static double log10(double a)
   {return Math.log10(a);}
   public static double[] log10(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.log10(a[i]);}
	   return c;
   }
   // N based logarithm
   public static double logN(double a,double N)
   {return Math.log(a)/Math.log(N);}
    public static double[] logN(double a[],double N)
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.log(a[i])/Math.log(N);}
	   return c;
   }  
   // Returns the natural logarithm of the sum of the argument and 1.
   public static double log1p(double a)
   {return Math.log1p(a);}
   
   //Returns the unbiased exponent used in the representation of a double.
   public static int getExponent(double d)
   {return Math.getExponent(d);}
   
   // Returns the unbiased exponent used in the representation of a float.
   public static int getExponent(float d)
   {return Math.getExponent(d);}     
   
   // Returns the greater of two double values.
   public static double max(double a,double b)
   {return Math.max(a,b);}
   public static double max(double a[])
   {  int n=a.length;
      double max1=a[0];
      for(int i=1;i<n;i++)
      {max1=Math.max(max1,a[i]);}
	   return max1;
   }
   // Returns the greater of two int values.
   public static int max(int a,int b)
   {return Math.max(a,b);}
   public static int max(int a[])
   {  int n=a.length;
      int max1=a[0];
      for(int i=1;i<n;i++)
      {max1=Math.max(max1,a[i]);}
	   return max1;
   }
   // Returns the greater of two long values.
   public static long max(long a,long b)
   {return Math.max(a,b);}
   
   public static long min(long a,long b)
   {return Math.min(a,b);}
   
   public static long max(long a[])
   {  int n=a.length;
      long max1=a[0];
      for(int i=1;i<n;i++)
      {max1=Math.max(max1,a[i]);}
	   return max1;
   }
   
   public static double max(double xmin,double dx,double xmax)
   {  double a[]=setX(xmin,dx,xmax); 
	   int n=a.length;
      double max1=a[0];
      for(int i=1;i<n;i++)
      {max1=Math.max(max1,a[i]);}
	   return max1;
   }
   
   public static double min(double xmin,double dx,double xmax)
   {  double a[]=setX(xmin,dx,xmax); 
	   int n=a.length;
      double max1=a[0];
      for(int i=1;i<n;i++)
      {max1=Math.min(max1,a[i]);}
	   return max1;
   }
   
   public static long min(long a[])
   {  int n=a.length;
      long max1=a[0];
      for(int i=1;i<n;i++)
      {max1=Math.min(max1,a[i]);}
	   return max1;
   }
   // Returns Euler's number e raised to the power of a double value.
   public static double exp(double a)
   {return Math.exp(a);}
   public static double[] exp(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.exp(a[i]);}
	   return c;
   }
   //same as exp (Turkish equivalent)
   public static double us(double a)
   {return Math.exp(a);}
   public static double[] us(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.exp(a[i]);}
	   return c;
   }
   // e^x-1
   public static double expm1(double a)
   {return Math.expm1(a);}
      public static double[] expm1(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.expm1(a[i]);}
	   return c;
   }
   //x^y
   public static double pow(double a,double b)
   {return Math.pow(a,b);}
   public static double[] pow(double a[],double b)
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.pow(a[i],b);}
	   return c;
   }
   //same as pow x^y (Turkish equivalent)
   public static double us(double a,double b)
   {return Math.pow(a,b);}
   public static double[] us(double a[],double b)
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.pow(a[i],b);}
	   return c;
   }
   
   // Returns sqrt(x2 +y2) without intermediate overflow or underflow.  
   public static double hypot(double a,double b)
   {return Math.hypot(a,b);}
   public static double[] hypot(double a[],double b[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.hypot(a[i],b[i]);}
	   return c;
   }
   // Returns the signum function of the argument; zero if the argument is zero, 
   // 1.0 if the argument is greater than zero, -1.0 if the argument is less than 
   // zero.
   public static double signum(double d)
   {return Math.signum(d);}
   
   //Converts an angle measured in degrees to an approximately equivalent 
   //angle measured in radians.
   public static double toRadians(double a)
   {return Math.toRadians(a);}
   public static double[] toRadians(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.toRadians(a[i]);}
	   return c;
   }
   //Converts an angle measured in radians to an approximately equivalent 
   //angle measured in degrees.
   public static double toDegrees(double a)
   {return Math.toDegrees(a);}
   public static double[] toDegrees(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.toDegrees(a[i]);}
	   return c;
   }
   //same as toRadians (Turkish equivalent)
   public static double radyanacevir(double a)
   {return Math.toRadians(a);}
   public static double[] radyanacevir(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.toRadians(a[i]);}
	   return c;
   }
   //same as toDegree (Turkish equivalent)
   public static double dereceyecevir(double a)
   {return Math.toDegrees(a);}
   //degree input sine function
   
   public static double[] dereceyecevir(double a[])
   {  int n=a.length;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.toDegrees(a[i]);}
	   return c;
   }
   //sine of a degree
   public static double sind(double a) 
   { 
	double r=(a/180.0)*PI;
	return Math.sin(r);  
   } 
   public static double[] sind(double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {  r=(a[i]/180.0)*PI;
	      c[i]=Math.sin(r);}
	   return c;
   }
   //cosine of a degree
   public static double cosd(double a) 
   { 
	double r=(a/180.0)*PI;
	return Math.cos(r);  
   }
   public static double[] cosd(double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {  r=(a[i]/180.0)*PI;
	      c[i]=Math.cos(r);}
	   return c;
   }
  //tangent o a degree
   public static double tand(double a) 
   { 
	double r=(a/180.0)*PI;
	return Math.tan(r);  
   }
   public static double[] tand(double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      {  r=(a[i]/180.0)*PI;
	      c[i]=Math.tan(r);}
	   return c;
   }
   //Returns a double value with a positive sign, 
   //greater than or equal to 0.0 and less than 1.0.
   public static double random() 
   { 
	return Math.random();  
   }
   public static double[] random(int n) 
   {  double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=Math.random();}
	   return c;  
   }
   public static double[] random(double a,double b,int n) 
   {  double c[]=new double[n];
      for(int i=0;i<n;i++)
      {c[i]=a+(b-a)*Math.random();}
	   return c; 
   }
   // inverse sine , converts to degree value
   public static double asind(double b) 
   { 
	return Math.asin(b)*180.0/PI;  
   }
   public static double[] asind(double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=Math.asin(a[i])*180.0/PI;}
	   return c;
   }
   // inverse cosine , converts to degree value 
   public static double acosd(double b) 
   { 
	return Math.acos(b)*180/PI;  
   }
   public static double[] acosd(double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=Math.acos(a[i])*180.0/PI;}
	   return c;
   }
      // inverse tangent , converts to degree value
  public static double atand(double b) 
   { 
	return Math.atan(b)*180/PI;  
   }
   public static double[] atand(double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=Math.atan(a[i])*180.0/PI;}
	   return c;
   }
   // inverse tangent , converts to degree value use x and y cartesian input
   public static double atan2d(double a,double b) 
   { 
	return Math.atan2(a,b)*180.0/PI;  
   }
   public static double[] atan2d(double a[],double b[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=Math.atan2(a[i],b[i])*180.0/PI;}
	   return c;
   }
//gamma function   
public static double gamma(double xx)
{   //Gamma fonksiyonu
	double a[]=lngamma(xx);
	return a[1]*Math.exp(a[0]);
} 
   public static double[] gamma(double a[])
   {  int n=a.length;
     double aa[][]=new double[2][n];
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { aa[i]=lngamma(a[i]);
	      c[i]=aa[i][1]*Math.exp(aa[i][0]);}
	   return c;
   } 
//natural logarithm of gamma function 
	public static double [] lngamma(double x)
    {//logaithm of gamma function
    //xx[0] value of logarithm of gamma function 
    //xx[1] sign of gamma functon
	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;
    }
 // C Euler constant   
public static double C(long n)
{
double cc=0;
for(double i=1;i<=n;i++)
{cc+=1.0/i;}
cc-=Math.log(n);
return cc;
}

// factorial
public static double factrl(int n)
{
	// n! of an integer number
	int ntop=5;
	double a[]=new double[50];
	a[0]=1.0;
	a[1]=1.0;
	a[2]=2.0;
	a[3]=6.0;
	a[5]=24.0;
	int j;
	if (n < 0) System.out.println("Negative factorial in routine FACTRL");
	if (n > 50) gamma(n+1.0);
	while (ntop<=n) {
		j=ntop-1;
		a[ntop]=a[j]*ntop;
		ntop++;
	}
	return a[n];
}
//natural logarithm of gamma function
public static double factln(int n)
{
	//Math.log(factorial) of an integer number
	double a[]={0.0,0.0,0.6931471805599453,1.0986122886681098,1.3862943611198906,
	1.6094379124341003,1.791759469228055,1.9459101490553132,2.0794415416798357,
	2.1972245773362196,2.302585092994046,2.3978952727983707,2.4849066497880004,
	2.5649493574615367,2.6390573296152584,2.70805020110221,2.772588722239781,
	2.833213344056216,2.8903717578961645,2.9444389791664403,2.995732273553991,
	3.044522437723423,3.091042453358316,3.1354942159291497,3.1780538303479458,
	3.2188758248682006,3.258096538021482,3.295836866004329,3.332204510175204,
	3.367295829986474,3.4011973816621555,3.4339872044851463,3.4657359027997265,
	3.4965075614664802,3.5263605246161616,3.5553480614894135,3.58351893845611,
	3.6109179126442243,3.6375861597263857,3.6635616461296463,3.6888794541139363,
	3.713572066704308,3.7376696182833684,3.7612001156935624,3.784189633918261,
	3.8066624897703196,3.828641396489095,3.8501476017100584,3.871201010907891,
	3.8918202981106265,3.912023005428146,3.9318256327243257,3.9512437185814275,
	3.970291913552122,3.9889840465642745,4.007333185232471,4.02535169073515,
	4.04305126783455,4.060443010546419,4.07753744390572,4.0943445622221,4.110873864173311,
	4.127134385045092,4.143134726391533,4.1588830833596715,4.174387269895637,4.189654742026425,
	4.204692619390966,4.219507705176107,4.23410650459726,4.248495242049359,4.2626798770413155,
	4.276666119016055,4.290459441148391,4.30406509320417,4.31748811353631,4.330733340286331,
	4.343805421853684,4.356708826689592,4.3694478524670215,4.382026634673881,4.394449154672439,
	4.406719247264253,4.418840607796598,4.430816798843313,4.442651256490317,4.454347296253507,
	4.465908118654584,4.477336814478207,4.48863636973214,4.499809670330265,4.51085950651685,
	4.5217885770490405,4.532599493153256,4.543294782270004,4.553876891600541,4.564348191467836,
	4.574710978503383,4.584967478670572,4.59511985013459,4.605170185988092};
	if (n < 0) {System.out.println("Negative factorial in routine FACTLN"); return 0.0;}
	else if(n > 100) {double b[]=lngamma(n+1.0);return b[0];}
	else {return a[n];}
}
public static double PI()
{
double a=1.4142135623730951;
double b=0;
double x=2+a;
double t=0;
for(int k=1;k<=5;k++)
{t=Math.sqrt(a);
 b=t*(1+b)/(a+b);
 a=0.5*(t+1.0/t);
 x=x*b*(1+a)/(1+b);
}
return x;
}

/*
//Ramanujan PI formula
public static double PI()
{  //Ramanujan formula
   double p=0,pi=0.0;
   double c3=2.0*Math.sqrt(2.0)/9801.0;
   double c,c1,c2;
   for(int n=0;n<20;n++)
   {c1=factln(n);c=factln(4*n)-c1*c1*c1*c1+Math.log(1103.0+26390.0*n)-4.*n*Math.log(396);c2=c3*Math.exp(c);  
   p+=c2;
   pi=1.0/p;
   }  
   return pi;   
}
*/

//binomial coefficients
public static double  bico(int n,int k)
{
	return Math.floor(0.5+Math.exp(factln(n)-factln(k)-factln(n-k)));
}

//beta function
public static double beta(double z,double w)
{
	//gamma(z)*gamma(w)/gamma(z+w);
	return Math.exp(lngamma(z)[0]+lngamma(w)[0]-lngamma(z+w)[0]);
}
	//incomplete gamma function P(a,x)
public static double gammp(double a,double x)
{
	double gamser,gammcf;
	if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine GAMMP");
	if (x < (a+1.0)) {
		gamser=gser(a,x);
		return gamser;
	} else {
		gammcf=gcf(a,x);
		return 1.0-gammcf;
	}
}

//incomplete gamma function Q(a,x)=1-P(a,x);
public static double gammq(double a,double x)
{
	double gamser,gammcf;
	if (x < 0.0 || a <= 0.0) System.out.println("Invalid arguments in routine GAMMQ");
	if (x < (a+1.0)) {
		gamser=gser(a,x);
		return 1.0-gamser;
	} else {
		gammcf=gcf(a,x);
		return gammcf;
	}
}

public static double gser(double a,double x)
{
	int ITMAX=200;
	double EPS=3.0e-7;
	int n;
	double sum,del,ap;
	double gln=lngamma(a)[0];
	double gamser=0;
	if (x <= 0.0) {
		if (x < 0.0) System.out.println("x less than 0 in routine GSER");
		gamser=0.0;
		return gamser;
	} else {
		ap=a;
		del=sum=1.0/a;
		for (n=1;n<=ITMAX;n++) {
			ap += 1.0;
			del *= x/ap;
			sum += del;
			if (Math.abs(del) < Math.abs(sum)*EPS) {
				gamser=sum*Math.exp(-x+a*Math.log(x)-(gln));
				return gamser;
			}
		}
		//System.out.println("a too large, ITMAX too small in routine GSER");
		return gamser;
	}
}

public static double gcf(double a,double x)
{
	int ITMAX=200;
	double EPS=3.0e-7;
	int n;
	double gln;
	double gammcf=0;
	double gold=0.0,g,fac=1.0,b1=1.0;
	double b0=0.0,anf,ana,an,a1,a0=1.0;
	gln=lngamma(a)[0];
	a1=x;
	for (n=1;n<=ITMAX;n++) {
		an=(double) n;
		ana=an-a;
		a0=(a1+a0*ana)*fac;
		b0=(b1+b0*ana)*fac;
		anf=an*fac;
		a1=x*a0+anf*a1;
		b1=x*b0+anf*b1;
		if (a1!=0) {
			fac=1.0/a1;
			g=b1*fac;
			if (Math.abs((g-gold)/g) < EPS) {
				gammcf=Math.exp(-x+a*Math.log(x)-gln)*g;
				return gammcf;
			}
			gold=g;
		}
	}
	//System.out.println("a too large, ITMAX too small in routine GCF");
    return gammcf;
}
//error function
public static double  erf(double x)
{
   //error function
   return x < 0.0 ? -gammp(0.5,x*x) : gammp(0.5,x*x);
}
   public static double[] erf(double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=erf(a[i]);}
	   return c;
   }
//complimentary error function
public static double  erfc(double x)
{
	return x < 0.0 ? 1.0+gammp(0.5,x*x) : gammq(0.5,x*x);
}
   public static double[] erfc(double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=erfc(a[i]);}
	   return c;
   }
public static double  erfcc(double x)
{
	double t,z,ans;
	z=Math.abs(x);
	t=1.0/(1.0+0.5*z);
	ans=t*Math.exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
		t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
		t*(-0.82215223+t*0.17087277)))))))));
	return  x >= 0.0 ? ans : 2.0-ans;
}

//beta function
public static double betacf(double a,double b,double x)
{
	int ITMAX=100;
	double EPS=3.0e-7;
	double qap,qam,qab,em,tem,d;
	double bz,bm=1.0,bp,bpp;
	double az=1.0,am=1.0,ap,app,aold;
	int m;
	qab=a+b;
	qap=a+1.0;
	qam=a-1.0;
	bz=1.0-qab*x/qap;
	for (m=1;m<=ITMAX;m++) {
		em=(double) m;
		tem=em+em;
		d=em*(b-em)*x/((qam+tem)*(a+tem));
		ap=az+d*am;
		bp=bz+d*bm;
		d = -(a+em)*(qab+em)*x/((qap+tem)*(a+tem));
		app=ap+d*az;
		bpp=bp+d*bz;
		aold=az;
		am=ap/bpp;
		bm=bp/bpp;
		az=app/bpp;
		bz=1.0;
		if (Math.abs(az-aold) < (EPS*Math.abs(az))) return az;
	}
	System.out.println("a or b too big, or ITMAX too small in BETACF");
    return az;
}

//incomplete beta function
public static double betai(double a,double b,double x)
{
	double bt;
	if (x < 0.0 || x > 1.0) System.out.println("Bad x in routine BETAI");
	if (x == 0.0 || x == 1.0) bt=0.0;
	else
		bt=Math.exp(lngamma(a+b)[0]-lngamma(a)[0]-lngamma(b)[0]+a*Math.log(x)+b*Math.log(1.0-x));
	if (x < (a+1.0)/(a+b+2.0))
		return bt*betacf(a,b,x)/a;
	else
		return 1.0-bt*betacf(b,a,1.0-x)/b;
}
//J Bessel function of double order v
   public static double[] J(double v,double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=J(v,a[i]);}
	    return c;
   }
   
public static double J(double v,double x)
{ 
double JJ=0;double a=1.0;int fact=1;double b;
if( x <=20)  
{
double x1,x2,x3,x4;
long x5=1;
double lna=0.0;
x1=v*Math.log(0.5*x);
double lna1=Math.log(0.25*x*x);
for(int k=0;k<100;k++)
{   
     x2=lna;//Math.log(a);
     x3=lngamma(k+1)[0];
     x4=lngamma(v+k+1)[0];
	 b=x1+x2-x3-x4;
	 b=x5*Math.exp(b);
	 lna+=lna1;
	 x5*=-1;
	 JJ+=b;	
}
if(v==0 && x==0) JJ=1;
 //if(v==1 && x==0) JJ=0;
}
else //x>30
{JJ=Math.sqrt(2.0/Math.PI/x)*Math.cos(x-0.5*v*Math.PI-0.25*Math.PI);}
return JJ;
}
//Y Bessel function of double order v
public static double[] Y(double v,double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=Y(v,a[i]);}
	    return c;
   }
public static double Y(double v,double x)
{
	double YY=0;
int vtest=(int)(v*10)/10;	
if(v==0 && v==(double)vtest) YY=bessy0(x);
else if(v==1 && v==(double)vtest) YY=bessy1(x); 
else if (v==(double)vtest) YY=bessy((int)v,x);
else {YY=(J(v,x)*Math.cos(v*Math.PI)-J(-v,x))/Math.sin(v*Math.PI);}
return YY;
}
//K Bessel function of double order v	
public static double[] K(double v,double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=K(v,a[i]);}
	    return c;
   }
public static double K(double v,double x)
{
double YY=0;
int vtest=(int)(v*10)/10;	
if(v==0 && v==(double)vtest) YY=bessk0(x);
else if(v==1 && v==(double)vtest) YY=bessk1(x); 
else if (v==(double)vtest) YY=bessk((int)v,x);
else {YY=Math.PI/2.0*(I(-v,x)-I(v,x))/Math.sin(v*Math.PI);}
return YY;
}		
//I Bessel function of double order v	
public static double[] I(double v,double a[])
   {  int n=a.length;
   double r=0;
      double c[]=new double[n];
      for(int i=0;i<n;i++)
      { c[i]=I(v,a[i]);}
	    return c;
   }
public static double I(double v,double x)
{ //I Bessel fonksiyonu
double JJ=0;double a=1.0;int fact=1;double b; 
{
double x1,x2,x3,x4;
long x5=1;
double lna=0.0;
x1=v*Math.log(0.5*x);
double lna1=Math.log(0.25*x*x);
for(int k=0;k<100;k++)
{   
     x2=lna;//Math.log(a);
     x3=lngamma(k+1)[0];
     x4=lngamma(v+k+1)[0];
	 b=x1+x2-x3-x4;
	 b=Math.exp(b);
	 lna+=lna1;
	 JJ+=b;	
}
if(v==0 && x==0) JJ=1;
 //if(v==1 && x==0) JJ=0;
}
return JJ;
}
//J Bessel function of integer order 0
public static double bessj0(double x)
{
	double ax,z;
	double xx,y,ans,ans1,ans2;

	if ((ax=Math.abs(x)) < 8.0) {
		y=x*x;
		ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7
			+y*(-11214424.18+y*(77392.33017+y*(-184.9052456)))));
		ans2=57568490411.0+y*(1029532985.0+y*(9494680.718
			+y*(59272.64853+y*(267.8532712+y*1.0))));
		ans=ans1/ans2;
	} else {
		z=8.0/ax;
		y=z*z;
		xx=ax-0.785398164;
		ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
			+y*(-0.2073370639e-5+y*0.2093887211e-6)));
		ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			+y*(-0.6911147651e-5+y*(0.7621095161e-6
			-y*0.934935152e-7)));
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
	}
	return ans;
}
//J Bessel function of integer order 1
public static double bessj1(double x)
{
	double ax,z;
	double xx,y,ans,ans1,ans2;

	if ((ax=Math.abs(x)) < 8.0) {
		y=x*x;
		ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
			+y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))));
		ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
			+y*(99447.43394+y*(376.9991397+y*1.0))));
		ans=ans1/ans2;
	} else {
		z=8.0/ax;
		y=z*z;
		xx=ax-2.356194491;
		ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
			+y*(0.2457520174e-5+y*(-0.240337019e-6))));
		ans2=0.04687499995+y*(-0.2002690873e-3
			+y*(0.8449199096e-5+y*(-0.88228987e-6
			+y*0.105787412e-6)));
		ans=sqrt(0.636619772/ax)*(cos(xx)*ans1-z*sin(xx)*ans2);
		if (x < 0.0) ans = -ans;
	}
	return ans;
}
//J Bessel function of integer order n
public static double bessj(int n, double x)
{
	double ACC=40.0;
	double BIGNO=1.0e10;
	double BIGNI=1.0e-10;
	int j,m;
	boolean jsum=true;
	double ax,bj,bjm,bjp,sum,tox,ans;
	if (n < 2) System.out.println("Index n less than 2 in bessj");
	ax=Math.abs(x);
	if (ax == 0.0)
		return 0.0;
	else if (ax > (double) n) {
		tox=2.0/ax;
		bjm=bessj0(ax);
		bj=bessj1(ax);
		for (j=1;j<n;j++) {
			bjp=j*tox*bj-bjm;
			bjm=bj;
			bj=bjp;
		}
		ans=bj;
	} else {
		tox=2.0/ax;
		m=2*((n+(int) sqrt(ACC*n))/2);
		jsum=false;
		bjp=ans=sum=0.0;
		bj=1.0;
		for (j=m;j>0;j--) {
			bjm=j*tox*bj-bjp;
			bjp=bj;
			bj=bjm;
			if (Math.abs(bj) > BIGNO) {
				bj *= BIGNI;
				bjp *= BIGNI;
				ans *= BIGNI;
				sum *= BIGNI;
			}
			if (jsum) sum += bj;
			jsum=!jsum;
			if (j == n) ans=bjp;
		}
		sum=2.0*sum-bj;
		ans /= sum;
	}
	double JJ=0;
	if(x < 0.0 ) {JJ=-ans;}
	else {JJ=ans;}
	return  ans;
}
//I Bessel function of integer order 0
public static double bessi0(double x)
{
	double ax,ans;
	double y;

	if ((ax=Math.abs(x)) < 3.75) {
		y=x/3.75;
		y*=y;
		ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
			+y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
	} else {
		y=3.75/ax;
		ans=(exp(ax)/sqrt(ax))*(0.39894228+y*(0.1328592e-1
			+y*(0.225319e-2+y*(-0.157565e-2+y*(0.916281e-2
			+y*(-0.2057706e-1+y*(0.2635537e-1+y*(-0.1647633e-1
			+y*0.392377e-2))))))));
	}
	return ans;
}
//I Bessel function of integer order 1
public static double bessi1(double x)
{
	double ax,ans;
	double y;

	if ((ax=Math.abs(x)) < 3.75) {
		y=x/3.75;
		y*=y;
		ans=ax*(0.5+y*(0.87890594+y*(0.51498869+y*(0.15084934
			+y*(0.2658733e-1+y*(0.301532e-2+y*0.32411e-3))))));
	} else {
		y=3.75/ax;
		ans=0.2282967e-1+y*(-0.2895312e-1+y*(0.1787654e-1
			-y*0.420059e-2));
		ans=0.39894228+y*(-0.3988024e-1+y*(-0.362018e-2
			+y*(0.163801e-2+y*(-0.1031555e-1+y*ans))));
		ans *= (exp(ax)/sqrt(ax));
	}
	return x < 0.0 ? -ans : ans;
}
//I Bessel function of integer order n
public static double bessi(int n, double x)
{
    double ACC=40.0;
	double BIGNO=1.0e10;
	double BIGNI=1.0e-10;
	int j;
	double bi,bim,bip,tox,ans;

	if (n < 2) System.out.println("Index n less than 2 in bessi");
	if (x == 0.0)
		return 0.0;
	else {
		tox=2.0/Math.abs(x);
		bip=ans=0.0;
		bi=1.0;
		for (j=2*(n+(int) sqrt(ACC*n));j>0;j--) {
			bim=bip+j*tox*bi;
			bip=bi;
			bi=bim;
			if (Math.abs(bi) > BIGNO) {
				ans *= BIGNI;
				bi *= BIGNI;
				bip *= BIGNI;
			}
			if (j == n) ans=bip;
		}
		ans *= bessi0(x)/bi;
		if(x<0) return -ans;
		else return ans;
	}
}
//Y Bessel function of integer order 0
public static double bessy0(double x)
{
	double z;
	double xx,y,ans,ans1,ans2;

	if (x < 8.0) {
		y=x*x;
		ans1 = -2957821389.0+y*(7062834065.0+y*(-512359803.6
			+y*(10879881.29+y*(-86327.92757+y*228.4622733))));
		ans2=40076544269.0+y*(745249964.8+y*(7189466.438
			+y*(47447.26470+y*(226.1030244+y*1.0))));
		ans=(ans1/ans2)+0.636619772*bessj0(x)*log(x);
	} else {
		z=8.0/x;
		y=z*z;
		xx=x-0.785398164;
		ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4
			+y*(-0.2073370639e-5+y*0.2093887211e-6)));
		ans2 = -0.1562499995e-1+y*(0.1430488765e-3
			+y*(-0.6911147651e-5+y*(0.7621095161e-6
			+y*(-0.934945152e-7))));
		ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
	}
	return ans;
}
//Y Bessel function of integer order 1
public static double bessy1(double x)
{
	double z;
	double xx,y,ans,ans1,ans2;
	if (x < 8.0) {
		y=x*x;
		ans1=x*(-0.4900604943e13+y*(0.1275274390e13
			+y*(-0.5153438139e11+y*(0.7349264551e9
			+y*(-0.4237922726e7+y*0.8511937935e4)))));
		ans2=0.2499580570e14+y*(0.4244419664e12
			+y*(0.3733650367e10+y*(0.2245904002e8
			+y*(0.1020426050e6+y*(0.3549632885e3+y)))));
		ans=(ans1/ans2)+0.636619772*(bessj1(x)*log(x)-1.0/x);
	} else {
		z=8.0/x;
		y=z*z;
		xx=x-2.356194491;
		ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
			+y*(0.2457520174e-5+y*(-0.240337019e-6))));
		ans2=0.04687499995+y*(-0.2002690873e-3
			+y*(0.8449199096e-5+y*(-0.88228987e-6
			+y*0.105787412e-6)));
		ans=sqrt(0.636619772/x)*(sin(xx)*ans1+z*cos(xx)*ans2);
	}
	return ans;
}
//Y Bessel function of integer order n
public static double bessy(int n, double x)
{
	int j;
	double by,bym,byp,tox;

	if (n < 2) System.out.println("Index n less than 2 in bessy");
	tox=2.0/x;
	by=bessy1(x);
	bym=bessy0(x);
	for (j=1;j<n;j++) {
		byp=j*tox*by-bym;
		bym=by;
		by=byp;
	}
	return by;
}
public static double bessk0(double x)
{
	double y,ans;

	if (x <= 2.0) {
		y=x*x/4.0;
		ans=(-log(x/2.0)*bessi0(x))+(-0.57721566+y*(0.42278420
			+y*(0.23069756+y*(0.3488590e-1+y*(0.262698e-2
			+y*(0.10750e-3+y*0.74e-5))))));
	} else {
		y=2.0/x;
		ans=(exp(-x)/sqrt(x))*(1.25331414+y*(-0.7832358e-1
			+y*(0.2189568e-1+y*(-0.1062446e-1+y*(0.587872e-2
			+y*(-0.251540e-2+y*0.53208e-3))))));
	}
	return ans;
}
//K Bessel function of integer order 1
public static double bessk1(double x)
{
	double y,ans;

	if (x <= 2.0) {
		y=x*x/4.0;
		ans=(log(x/2.0)*bessi1(x))+(1.0/x)*(1.0+y*(0.15443144
			+y*(-0.67278579+y*(-0.18156897+y*(-0.1919402e-1
			+y*(-0.110404e-2+y*(-0.4686e-4)))))));
	} else {
		y=2.0/x;
		ans=(exp(-x)/sqrt(x))*(1.25331414+y*(0.23498619
			+y*(-0.3655620e-1+y*(0.1504268e-1+y*(-0.780353e-2
			+y*(0.325614e-2+y*(-0.68245e-3)))))));
	}
	return ans;
}
//K Bessel function of integer order n
public static double bessk(int n, double x)
{
	int j;
	double bk,bkm,bkp,tox;

	if (n < 2) System.out.println("Index n less than 2 in bessk");
	tox=2.0/x;
	bkm=bessk0(x);
	bk=bessk1(x);
	for (j=1;j<n;j++) {
		bkp=bkm+j*tox*bk;
		bkm=bk;
		bk=bkp;
	}
	return bk;
}

//Legendre polynomials of integer order n
// -1<=x<=1
public static double P(int n, double x)
{
if(n==0) return 1.0;
else if(n==1) return x;
else
{return (P((n-1),x)*x*(double)(2*n-1)-P((n-2),x)*(double)(n-1))/(double)n;}
// end legendre polynomials
}
//Derivative of Legendre Polynomial
public static double dP(int n, double x)
{
{return (double)n/(x*x-1)*(x*P(n,x)-P((n-1),x));}
// end legendre polynomials
}

//Chebychev polynomials of integer order n
// -1<=x<=1
public static double T(int n, double x)
{
if(n==0) return 1.0;
else if(n==1) return x;
else
{return (2.0*T((n-1),x)*x-T((n-2),x));}
//end Chebychev polynomials
}

public static double fin_efficiency(String type,double m,double L)
{   double eta=1.0;
	if(type.equals("rectangular") || type.equals("Rectangular"))
	{eta=tanh(m*L)/(m*L);}
    else if(type.equals("rectangular pin") || type.equals("Rectangular pin"))
	{eta=tanh(m*L)/(m*L);}
	else if(type.equals("triangualar") || type.equals("Triangualar"))
	{eta=1.0/(m*L)*bessi1(2.0*m*L)/bessi0(2.0*m*L);}
    else if(type.equals("triangualar pin") || type.equals("Triangualar pin"))
	{eta=2.0/(m*L)*bessi(2,2.0*m*L)/bessi1(2.0*m*L);}
	else if(type.equals("parabolic") || type.equals("Parabolic"))
	{eta=2.0/(sqrt(4.0*(m*L)*(m*L)+1.0)+1.0);}
	else if(type.equals("parabolic pin") || type.equals("Parabolic pin"))
	{eta=2.0/(sqrt(4.0/9.0*(m*L)*(m*L)+1.0)+1.0);}
	return eta;
}
public static double fin_efficiency(String type,double m,double r1,double r2)
{   double eta=1.0;
	if(type.equals("circular") || type.equals("Circular"))
	{double C2=(2*r1/m)/(r2*r2-r1*r1);
	 eta=C2*(bessk1(m*r1)*bessi1(m*r2)-bessi1(m*r1)*bessk1(m*r2))/(bessi0(m*r1)*bessk1(m*r2)+bessk0(m*r1)*bessi1(m*r2));
	 return eta;
	}
	return eta;
}
//Radiation shape function
  public static double F(double x)
  {
  double C2=1.438769e4;
  double y=C2/x;
  double total=0;
  for(int n=1;n<40;n++)
  {total+=Math.exp(-n*y)/n*(y*y*y+3.0*y*y/n+6.0*y/(n*n)+6.0/(n*n*n));}  
  double pi=Math.PI;
  total=15.0/(pi*pi*pi*pi)*total;
  return total;
  }

     public static double[][] read(String filename) throws IOException
     {
	 //reading double values from one line of data    
	 BufferedReader fin=new BufferedReader(new FileReader(filename));
     double a[][]=read(fin);
     return a;
     }
     
    public static double[][] read(BufferedReader fin) throws IOException
    {
	 ArrayList<data> a=new ArrayList<data>();
	 int mmax=0;
	 double b[][];
	 int n=0,m=0;
	 String s1="";	 
	 try{
       while(fin!=null)
       { s1=fin.readLine();
       if(s1==null)break;
         StringTokenizer token=new StringTokenizer(s1);
         m=token.countTokens();
         if(m>mmax) mmax=m;
         double aa[]=new double[m];
         int j=0;           
         while(token.hasMoreTokens())
         { Double ax=new Double(token.nextToken());
           aa[j++]=ax.doubleValue();
         }
         data vx=new data(aa);
         a.add(vx);
       n++;
	   }
       } catch(EOFException e_eof)
          {
          //close ffile
          try{fin.close();
      } 
      catch(NullPointerException e) {}   
      catch(IOException e)
             {
             System.err.println("Error Closing File\n"+e.toString());
             System.exit(1);
             }
          } //End of EOFException
     b=new double[n][mmax];
     int k=0;
     Iterator<data> i=a.iterator();
     while(i.hasNext()) {data vx=i.next();b[k++]=vx.x;}
     return b;
     }     
          
     public static String toString(double a[][])
     {
	 String s="";
     for(int i=0;i<a.length;i++)
     {
	 for(int j=0;j<a[0].length;j++)    
     {s+="a["+i+"]["+j+"] = "+a[i][j]+" ";}
     s+="\n";
     }
     return s;
	 }
}