//=============================================================//
//     Numerical Analysis Package in Java                      //
//     Numerical analysis calculation class Numeric            //
//                                                             //
//     Author : Dr. M. Turhan Çoban                            //
//     EGE Üniversitesi Mühendislik Fakultesi                  //
//     Makina Mühendisliği Bölümü                              //
//     email : turhan.coban@ege.edu.tr                         //
//     web   : www.axtelsoft.com/turhan.coban                  //
//=============================================================//
// Numeric packages utilises functions in abstract classes
// Actual functions should be derived from this abstract classes
// by extends term. example:  class f1 extends f_x

import java.io.*;
import java.text.*;
import java.util.Locale;

abstract class f_x
{
  //single function single independent variable
  // example f=x*x
  abstract double func(double x);
}

abstract class f_xj
{
  // single function  multi independent variable
  // a single value is returned  (
  // example f=x[0]+sin(x[1])
  // func(x) returns the value of f
  abstract double func(double x[]);
}

abstract class f_xa
{
  // single function  multi independent variable
  // a single value is returned 
  // example f=t*x[0]+t^2*sin(x[1])
  // func(x) returns the value of f
  abstract double func(double t,double x[]);
}

abstract class f_xi
{
  // multifunction multi independent variable
  // a single value is returned indiced to equation_ref
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[0]-x[1]
  // func(x,1) returns the value of f[1]
  // func(x,0) returns the value of f[0]

  abstract double func(double x[],int equation_ref);
}

abstract class fi_xi
{
  // multifunction multi independent variable
  // vector of dependent variables are returned
  // example f[0]=x[0]+sin(x[1])
  //         f[1]=x[0]*x[0]-x[1]
  // func(x) returns the value of f[0] and f[1]
  // as a two dimensional vector

  abstract double[] func(double x[]);
}

abstract class fij_xi
{
  // multifunction multi independent variable
  // for n independent variable n*n matrix of
  // dependent variables are returned
  // example f[0][0]=x[0]+Math.sin(x[1])   f[0][1]=x[0]-x[1]
  //         f[1][0]=x[0]*x[0]-x[1]        f[1][1]=Math.exp(x[0]+x[1]*x[1]
  // func(x) returns the value of f[0][0], f[0][1]
  //                              f[1][0], f[1][1]
  // as a two dimensional matrix
  abstract double[][] func(double x[]);
}


public class Numeric
{
//This is a library of numerical methods
//=======Part 1 Minimisation (Optimisation) of a function
public static double[] grads(f_xj Fn,fi_xi Gn,double P0[])
{
   return grads(Fn,Gn,P0,100,1.0e-10,1.0e-10,0);
}

//  derivative(f_xj f_deriv,double x[])


public static double[] grads(f_xj Fn,fi_xi Gn,double P0[],double max1,double delta,double epsilon,int print)
{
    //GRADS   Gradient search for a minimum of non-linear multivariable function
    // Inputs
    //   Fn        name of the vector function
    //   Gn        gradient for Fn
    //   P0        starting point
    //   max1      maximum number of iterations
    //   delta     convergence tolerance for the independent variables
    //   epsilon   convergence tolerance for the dependent variable
    //   show      if show==1 the iterations are displayed
    //
    //   P0        point V0 for the minimum
    //   y0        function value  Fn(V0)
    //   h         minimum step size
    //   err       error  estimate
    //---------------------------------------------------------------------------

    int n = P0.length;
    int maxj = 20;
    double big = 1e8;
    double h = 1.0;
    double len = Matrix.norm(P0);
    double y0 = Fn.func(P0);
    if (len > 1e4) {
      h = len / 1e4;
    }
    double err = 1;
    double cnt = 0;
    double cond = 0;
    double P[] = new double[n];
    double P1[] = new double[n];
    double P2[] = new double[n];
    double Pmin[] = new double[n];
    double Y, y1, y2, ymin, hmin;
    double S[]=new double[n];
    double h0,h1,h2;
    double e0,e1,e2,d;
    int i;
    for(i=0;i<n;i++)
    { P[i]=P0[i]; }
    Y= y0;
    while (cnt<max1 && cond!=5 && (h>delta || err>epsilon))
    {
      double z[] = Gn.func(P0);
      for (i = 0; i < n; i++) {
            S[i] = z[i];
          }
      P1 = Matrix.add(P0, Matrix.multiply(S, h));
      P2 = Matrix.add(P0, Matrix.multiply(S, (2.0 * h)));
      y1 = Fn.func(P1);
      y2 = Fn.func(P2);
      cond = 0;
      int j = 0;
      while (j < maxj && cond == 0) {
        len = Matrix.norm(P0);
        if (y0 < y1) {
          for (i = 0; i < n; i++) {
            P2[i] = P1[i];
          }
          y2 = y1;
          h = h / 2.0;
          P1 = Matrix.add(P0, Matrix.multiply(S, h));
          y1 = Fn.func(P1);
        }
        else {
          if (y2 < y1) {
            for (i = 0; i < n; i++) {
              P1[i] = P2[i];
            }
            y1 = y2;
            h = 2 * h;
            P2 = Matrix.add(P0, Matrix.multiply(S, (2.0 * h)));
            y2 = Fn.func(P2);
          }
          else cond = -1;
        }
        j++;
        if (h < delta) {
          cond = 1;
        }
        if (Math.abs(h) > big || len > big) {cond = 5;}
      }
      if (cond == 5) {
        for (i = 0; i < n; i++) {
          Pmin[i] = P1[i];
        }
        ymin = y1;
      }
      else {
        d = 4 * y1 - 2 * y0 - 2 * y2; // Start of a long block:
        if (d < 0) {hmin = h * (4 * y1 - 3 * y0 - y2) / d;}
        else {
          cond = 4;
          hmin = h / 3.0;
        }
        Pmin = Matrix.add(P0, Matrix.multiply(S, hmin));
        ymin = Fn.func(Pmin);
        h0 = Math.abs(hmin);
        h1 = Math.abs(hmin - h);
        h2 = Math.abs(hmin - 2.0 * h);
        if (h0 < h) {
          h = h0;
        }
        if (h1 < h) {
          h = h1;
        }
        if (h2 < h) {
          h = h2;
        }
        if (h == 0) {
          h = hmin;
        }
        if (h < delta) {
          cond = 1;
        }
        e0 = Math.abs(y0 - ymin);
        e1 = Math.abs(y1 - ymin);
        e2 = Math.abs(y2 - ymin);
        if (e0 != 0 && e0 < err) {
          err = e0;
        }
        if (e1 != 0 && e1 < err) {
          err = e1;
        }
        if (e2 != 0 && e2 < err) {
          err = e2;
        }
        if (e0 == 0 && e1 == 0 && e2 == 0) {
          err = 0;
        }
        if (err < epsilon) {
          cond = 2;
        }
        if (cond == 2 && h < delta) {
          cond = 3;
        }
      } //end of if(cond==5)
      cnt++;
      for (i = 0; i < n; i++) {
        P0[i] = Pmin[i];
      }
      y0 = ymin;
      if(print==1)
       {System.out.println("x min = "+Matrix.toString(P0)+"/nymin = "+y0+"/nerror estimate ="+err+"/nminimum step size "+h); }
    } //end of while
    return P0;
}

public static double[] grads(f_xj Fn,double P0[])
{
   return grads(Fn,P0,200,1.0e-20,1.0e-20,0);
}

public static double[] grads(f_xj Fn,double P0[],double max1,double delta,double epsilon,int print)
{
    //GRADS   Gradient search for a minimum of non-linear multivariable function
    // Inputs
    //   Fn        name of the vector function
    //   Gn        gradient for Fn
    //   P0        starting point
    //   max1      maximum number of iterations
    //   delta     convergence tolerance for the independent variables
    //   epsilon   convergence tolerance for the dependent variable
    //   show      if show==1 the iterations are displayed
    //
    //   P0        point V0 for the minimum
    //   y0        function value  Fn(V0)
    //   h         minimum step size
    //   err       error  estimate
    //---------------------------------------------------------------------------

    int n = P0.length;
    int maxj = 20;
    double big = 1e8;
    double h = 1.0;
    double len = Matrix.norm(P0);
    double y0 = Fn.func(P0);
    if (len > 1e4) {
      h = len / 1e4;
    }
    double err = 1;
    double cnt = 0;
    double cond = 0;
    double PP[] = new double[n];
    double P[] = new double[n];
    double P1[] = new double[n];
    double P2[] = new double[n];
    double Pmin[] = new double[n];
    double Y, y1, y2, ymin, hmin;
    double S[]=new double[n];
    double h0,h1,h2;
    double e0,e1,e2,d;
    int i;
    for(i=0;i<n;i++)
    { P[i]=P0[i];PP[i]=P0[i]; }
    Y= y0;
    while (cnt<max1 && cond!=5 && (h>delta || err>epsilon))
    {
      double z[] = derivative(Fn,P0);
      for (i = 0; i < n; i++) {
            S[i] = z[i];
          }
      P1 = Matrix.add(P0, Matrix.multiply(S, h));
      P2 = Matrix.add(P0, Matrix.multiply(S, (2.0 * h)));
      y1 = Fn.func(P1);
      y2 = Fn.func(P2);
      cond = 0;
      int j = 0;
      while (j < maxj && cond == 0) {
        len = Matrix.norm(P0);
        if (y0 < y1) {
          for (i = 0; i < n; i++) {
            P2[i] = P1[i];
          }
          y2 = y1;
          h = h / 2.0;
          P1 = Matrix.add(P0, Matrix.multiply(S, h));
          y1 = Fn.func(P1);
        }
        else {
          if (y2 < y1) {
            for (i = 0; i < n; i++) {
              P1[i] = P2[i];
            }
            y1 = y2;
            h = 2 * h;
            P2 = Matrix.add(P0, Matrix.multiply(S, (2.0 * h)));
            y2 = Fn.func(P2);
          }
          else cond = -1;
        }
        j++;
        if (h < delta) {
          cond = 1;
        }
        if (Math.abs(h) > big || len > big) {cond = 5;}
      }
      if (cond == 5) {
        for (i = 0; i < n; i++) {
          Pmin[i] = P1[i];
        }
        ymin = y1;
      }
      else {
        d = 4 * y1 - 2 * y0 - 2 * y2; // Start of a long block:
        if (d < 0) {hmin = h * (4 * y1 - 3 * y0 - y2) / d;}
        else {
          cond = 4;
          hmin = h / 3.0;
        }
        Pmin = Matrix.add(P0, Matrix.multiply(S, hmin));
        ymin = Fn.func(Pmin);
        h0 = Math.abs(hmin);
        h1 = Math.abs(hmin - h);
        h2 = Math.abs(hmin - 2.0 * h);
        if (h0 < h) {
          h = h0;
        }
        if (h1 < h) {
          h = h1;
        }
        if (h2 < h) {
          h = h2;
        }
        if (h == 0) {
          h = hmin;
        }
        if (h < delta) {
          cond = 1;
        }
        e0 = Math.abs(y0 - ymin);
        e1 = Math.abs(y1 - ymin);
        e2 = Math.abs(y2 - ymin);
        if (e0 != 0 && e0 < err) {
          err = e0;
        }
        if (e1 != 0 && e1 < err) {
          err = e1;
        }
        if (e2 != 0 && e2 < err) {
          err = e2;
        }
        if (e0 == 0 && e1 == 0 && e2 == 0) {
          err = 0;
        }
        if (err < epsilon) {
          cond = 2;
        }
        if (cond == 2 && h < delta) {
          cond = 3;
        }
      } //end of if(cond==5)
      cnt++;
      for (i = 0; i < n; i++) {
        P0[i] = Pmin[i];
      }
      y0 = ymin;
      if(print==1)
       {System.out.println("x min = "+Matrix.toString(P0)+"/nymin = "+y0+"/nerror estimate ="+err+"/nminimum step size "+h); }
    } //end of while
    return P0;
}

  public static double quadmin(f_x f,double a,double b,double epsilon)
  {return quadmin(f,a,b,epsilon,1.0e-10,0);}

  public static double quadmin(f_x f,double a,double b)
  {return quadmin(f,a,b,1.0e-10,1.0e-10,0);}

  public static double quadmin(f_x f,double a,double b,double epsilon,double delta,int print)
  {
      // Minimum of a single variable function by using polynomial quadratures
      double p0 = a;
      int maxj = 20;
      int maxk = 30;
      double big = 1e6;
      double err = 1;
      int k = 1;
      int cond = 0;
      double  h = 1;
      double f1,p1,p2,pmin,y0,y1,y2,ymin,d,hmin,dh,h0,h1,h2;
      double e0,e1,e2;
      double p,dp,dy,yp;
      if (Math.abs(p0)>1e4)
        { h = Math.abs(p0)/1e4;}
      while (k<maxk && err>epsilon && cond !=5 )
      {
        f1 = (f.func(p0 + 0.00001) - f.func(p0 - 0.00001)) / 0.00002;
        if (f1 > 0) {h = -Math.abs(h);}
        p1 = p0 + h;
        p2 = p0 + 2 * h;
        pmin = p0;
        y0 = f.func(p0);
        y1 = f.func(p1);
        y2 = f.func(p2);
        ymin = y0;
        cond = 0;
        int j = 0;
        while (j < maxj && Math.abs(h) > delta && cond == 0)
        {
          if (y0 <= y1)
          {
            p2 = p1;
            y2 = y1;
            h = h / 2;
            p1 = p0 + h;
            y1 = f.func(p1);
          }
          else
          {
            if (y2 < y1)
            {
              p1 = p2;
              y1 = y2;
              h = 2 * h;
              p2 = p0 + 2 * h;
              y2 = f.func(p2);
            }
            else
            {  cond = -1;}
          } //end of if(y0 <= y1)
          j = j + 1;
        if (Math.abs(h) > big || Math.abs(p0) > big) {cond = 5;}
        }//end of while(j<maxj
        if (cond == 5)
        {
              pmin = p1;
              ymin = f.func(p1);
        }
        else
        {
          d = 4 * y1 - 2 * y0 - 2 * y2;
          //Start of a long block:
          if (d < 0) {
            hmin = h * (4 * y1 - 3 * y0 - y2) / d;
          }
          else
          {
            hmin = h / 3;
            cond = 4;
          }
          pmin = p0 + hmin;
          ymin = f.func(pmin);
          dh = Math.abs(h);
          h0 = Math.abs(hmin);
          h1 = Math.abs(hmin - h);
          h2 = Math.abs(hmin - 2 * h);
          if (h0 < dh) {h = h0;}
          if (h1 < dh) {h = h1;}
          if (h2 < dh) {h = h2;}
          if (h == 0) {h = Math.abs(hmin);}
          if (h < delta) {cond = 1;}
          if (Math.abs(h) > big || Math.abs(pmin) > big) {cond = 5;}
          e0 = Math.abs(y0 - ymin);
          e1 = Math.abs(y1 - ymin);
          e2 = Math.abs(y2 - ymin);
          if (e0 != 0 && e0 < err) {err = e0;}
          if (e1 != 0 && e1 < err) {err = e1;}
          if (e2 != 0 && e2 < err)  {err = e2;}
          if (e0 != 0 && e1 == 0 && e2 == 0) {err = 0;}
          if (err < epsilon) {cond = 2;}
          p0 = pmin;
          k = k + 1;
      } // End of the long else block
      if (cond == 2 && h < delta) {cond = 3;}
  }//end of the while
  p = p0;
  dp = h;
  yp = f.func(p);
  dy = err;
  if(print==1)
  {System.out.println("x min = "+p+"ymin = "+yp+"error of x ="+dp+"error of y"+dy); }
  return p;
  }


public static double golden(f_x f,double a,double b,double epsilon)
{return golden(f,a,b,epsilon,1.0e-10,0);}
public static double golden(f_x f,double a,double b)
{return golden(f,a,b,1.0e-10,1.0e-10,0);}
public static double golden(f_x f,double a,double b,double epsilon,double delta,int print)
{
// f abstract function of class f_x (f(x) type of single function single variable)
// fibonacchi golden ratio search
    double r1 = (Math.sqrt(5.0)-1.0)/2.0; // golden ratio
    double r2 = r1*r1;
    double h = b - a;
    double ya = f.func(a);
    double yb = f.func(b);
    double c = a + r2*h;
    double d = a + r1*h;
    double yc = f.func(c);
    double yd = f.func(d);
    int k = 1;
    double dp,dy,p,yp;
    while ((Math.abs(yb-ya)>epsilon) || (h>delta))
    {
      k++;
      if (yc<yd)
        {
        b = d;
        yb = yd;
        d = c;
        yd = yc;
        h = b - a;
        c = a + r2 * h;
        yc = f.func(c);
       }
      else
       {
         a = c;
         ya = yc;
         c = d;
         yc = yd;
         h = b - a;
         d = a + r1 * h;
         yd = f.func(d);
       }//end of if
    }//end of while
    dp = Math.abs(b-a);
    dy = Math.abs(yb-ya);
    p = a;
    yp = ya;
    if (yb<ya)
    {
      p = b;
      yp = yb;
    }
    if(print==1)
    {System.out.println("x min = "+p+"ymin = "+yp+"error of x ="+dp+"error of y"+dy); }
    return p;
}


public static double[] nelder(f_xj fnelder,double a[],double da[],int maxiteration,double tolerance,int printlist)
{
double x[][]=new double[a.length+1][a.length];
for(int i=0;i<x.length;i++)
   {for(int j=0;j<x[0].length;j++)
     {if(i==j){x[i][j]=a[j]+da[j];}
      else    {x[i][j]=a[j]; }
     }
   }
   // Nelder mead multidimensional simplex minimisation method
   //  Nelder & Mead 1965 Computer J, v.7, 308-313.

   // Parameter identification
   // fnelder : abstract multidimensional function f(x)
   // x : for n dimension n+1 simlex points
   // maxiteration : maximum number of iteration points
   // tolerance : convergence tolerance
         int NDIMS = x.length-1;
         int NPTS = x.length;
         int FUNC = NDIMS;
         int ncalls = 0;
            ////// set up the starting simplex //////////////////
            double p[][]=new double[NPTS][NPTS]; // [row][col] = [whichvx][coord,FUNC]
            double z[]=new double[NDIMS];
            double best = 1E99;
            //////////////// initialize the funcvals ////////////////
            for (int i=0; i<NPTS; i++)
            {
              for (int j=0; j<NDIMS; j++)
                  {p[i][j] = x[i][j];}
              p[i][NDIMS] = fnelder.func(p[i]);
            }
            int iter=0;
            for (iter=1; iter<maxiteration; iter++)
            {
                /////////// identify lo, nhi, hi points //////////////
                int  ilo=0, ihi=0, inhi = -1; // -1 means missing
                double flo = p[0][FUNC];
                double fhi = flo;
                double pavg,sterr;
                for (int i=1; i<NPTS; i++)
                {
                   if (p[i][FUNC] < flo)
                     {flo=p[i][FUNC]; ilo=i;}
                   if (p[i][FUNC] > fhi)
                     {fhi=p[i][FUNC]; ihi=i;}
                }
                double fnhi = flo;
                inhi = ilo;
                for (int i=0; i<NPTS; i++)
                  if ((i != ihi) && (p[i][FUNC] > fnhi))
                    {fnhi=p[i][FUNC]; inhi=i;}
                ////////// exit criterion //////////////
                if ((iter % 4*NDIMS) == 0)
                {
                    //standart error of yi  should be less than set criteria (tolerance)
                    //calculate average (including highest point)
                    pavg=0;
                    for(int i=0;i<NPTS;i++)
                      pavg+=p[i][FUNC];
                    pavg/=NPTS;
                    double tot=0;
                    if(printlist!=0)
                    {  System.out.print(toString(iter,6));
                       for (int j=0; j<=NDIMS; j++)
                          { System.out.print(toString(p[ilo][j],25,12));}
                       System.out.println("");
                    }
                    for(int i=0;i<NPTS;i++)
                    { tot=(p[i][FUNC]-pavg)*(p[i][FUNC]-pavg);}
                    sterr=Math.sqrt(tot/NPTS);
                    //if(sterr < tolerance)
                      { for (int j=0; j<NDIMS; j++)
                           { z[j]=p[ilo][j];}
                        //break;
                      }
                    best = p[ilo][FUNC];
                }

                ///// compute ave[] vector excluding highest vertex //////

                double ave[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  ave[j] = 0;
                for (int i=0; i<NPTS; i++)
                  if (i != ihi)
                    for (int j=0; j<NDIMS; j++)
                       ave[j] += p[i][j];
                for (int j=0; j<NDIMS; j++)
                   ave[j] /= (NPTS-1);


                ///////// try reflect ////////////////

                double r[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  r[j] = 2*ave[j] - p[ihi][j];
                double fr = fnelder.func(r);

                if ((flo <= fr) && (fr < fnhi))  // in zone: accept
                {
                    for (int j=0; j<NDIMS; j++)
                      p[ihi][j] = r[j];
                    p[ihi][FUNC] = fr;
                    continue;
                }

                if (fr < flo)  //// below zone; try expand, else accept
                {
                    double e[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      e[j] = 3*ave[j] - 2*p[ihi][j];
                    double fe = fnelder.func(e);
                    if (fe < fr)
                    {
                       for (int j=0; j<NDIMS; j++)
                         p[ihi][j] = e[j];
                       p[ihi][FUNC] = fe;
                       continue;
                    }
                    else
                    {
                       for (int j=0; j<NDIMS; j++)
                         p[ihi][j] = r[j];
                       p[ihi][FUNC] = fr;
                       continue;
                    }
                }

                ///////////// above midzone, try contractions:

                if (fr < fhi)  /// try outside contraction
                {
                    double c[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      c[j] = 1.5*ave[j] - 0.5*p[ihi][j];
                    double fc = fnelder.func(c);
                    if (fc <= fr)
                    {
                        for (int j=0; j<NDIMS; j++)
                          p[ihi][j] = c[j];
                        p[ihi][FUNC] = fc;
                        continue;
                    }
                    else   /////// contract
                    {
                        for (int i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (int j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                        continue;
                    }
                }

                if (fr >= fhi)   /// over the top; try inside contraction
                {
                    double cc[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      cc[j] = 0.5*ave[j] + 0.5*p[ihi][j];
                    double fcc = fnelder.func(cc);
                    if (fcc < fhi)
                    {
                        for (int j=0; j<NDIMS; j++)
                          p[ihi][j] = cc[j];
                        p[ihi][FUNC] = fcc;
                        continue;
                    }
                    else    ///////// contract
                    {
                        for (int i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (int j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                    }
                }
            }
       return z;
}
      public static double[] nelder(f_xj fnelder,double a[],double da[],double tolerance)
      {return nelder(fnelder,a,da,300,tolerance,0);}

      public static double[] nelder(f_xj fnelder,double a[],double da[])
      {return nelder(fnelder,a,da,300,1.0e-15,0);}

      public static double[] nelder(f_xj fnelder,double a[],double da[],int maxiteration,double tolerance)
      {return nelder(fnelder,a,da,maxiteration,tolerance,0);}
      
public static double[] nelder(f_xj fnelder,double x[][],int maxiteration,double tolerance,int printlist)
{
   // Nelder mead multidimensional simplex minimisation method
   //  Nelder & Mead 1965 Computer J, v.7, 308-313.

   // Parameter identification
   // fnelder : abstract multidimensional function f(x)
   // x : for n dimension n+1 simlex points
   // maxiteration : maximum number of iteration points
   // tolerance : convergence tolerance
         int NDIMS = x.length-1;
         int NPTS = x.length;
         int FUNC = NDIMS;
         int ncalls = 0;
            ////// set up the starting simplex //////////////////
            double p[][]=new double[NPTS][NPTS]; // [row][col] = [whichvx][coord,FUNC]
            double z[]=new double[NDIMS];
            double best = 1E99;
            //////////////// initialize the funcvals ////////////////
            for (int i=0; i<NPTS; i++)
            {
              for (int j=0; j<NDIMS; j++)
                  {p[i][j] = x[i][j];}
              p[i][NDIMS] = fnelder.func(p[i]);
            }
            int iter=0;
            for (iter=1; iter<maxiteration; iter++)
            {
                /////////// identify lo, nhi, hi points //////////////
                int  ilo=0, ihi=0, inhi = -1; // -1 means missing
                double flo = p[0][FUNC];
                double fhi = flo;
                double pavg,sterr;
                for (int i=1; i<NPTS; i++)
                {
                   if (p[i][FUNC] < flo)
                     {flo=p[i][FUNC]; ilo=i;}
                   if (p[i][FUNC] > fhi)
                     {fhi=p[i][FUNC]; ihi=i;}
                }
                double fnhi = flo;
                inhi = ilo;
                for (int i=0; i<NPTS; i++)
                  if ((i != ihi) && (p[i][FUNC] > fnhi))
                    {fnhi=p[i][FUNC]; inhi=i;}
                ////////// exit criterion //////////////
                if ((iter % 4*NDIMS) == 0)
                {
                    //standart error of yi  should be less than set criteria (tolerance)
                    //calculate average (including highest point)
                    pavg=0;
                    for(int i=0;i<NPTS;i++)
                      pavg+=p[i][FUNC];
                    pavg/=NPTS;
                    double tot=0;
                    if(printlist!=0)
                    {  System.out.print(toString(iter,6));
                       for (int j=0; j<=NDIMS; j++)
                          { System.out.print(toString(p[ilo][j],25,12));}
                       System.out.println("");
                    }
                    for(int i=0;i<NPTS;i++)
                    { tot=(p[i][FUNC]-pavg)*(p[i][FUNC]-pavg);}
                    sterr=Math.sqrt(tot/NPTS);
                    //if(sterr < tolerance)
                      { for (int j=0; j<NDIMS; j++)
                           { z[j]=p[ilo][j];}
                        //break;
                      }
                    best = p[ilo][FUNC];
                }

                ///// compute ave[] vector excluding highest vertex //////

                double ave[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  ave[j] = 0;
                for (int i=0; i<NPTS; i++)
                  if (i != ihi)
                    for (int j=0; j<NDIMS; j++)
                       ave[j] += p[i][j];
                for (int j=0; j<NDIMS; j++)
                   ave[j] /= (NPTS-1);


                ///////// try reflect ////////////////

                double r[] = new double[NDIMS];
                for (int j=0; j<NDIMS; j++)
                  r[j] = 2*ave[j] - p[ihi][j];
                double fr = fnelder.func(r);

                if ((flo <= fr) && (fr < fnhi))  // in zone: accept
                {
                    for (int j=0; j<NDIMS; j++)
                      p[ihi][j] = r[j];
                    p[ihi][FUNC] = fr;
                    continue;
                }

                if (fr < flo)  //// below zone; try expand, else accept
                {
                    double e[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      e[j] = 3*ave[j] - 2*p[ihi][j];
                    double fe = fnelder.func(e);
                    if (fe < fr)
                    {
                       for (int j=0; j<NDIMS; j++)
                         p[ihi][j] = e[j];
                       p[ihi][FUNC] = fe;
                       continue;
                    }
                    else
                    {
                       for (int j=0; j<NDIMS; j++)
                         p[ihi][j] = r[j];
                       p[ihi][FUNC] = fr;
                       continue;
                    }
                }

                ///////////// above midzone, try contractions:

                if (fr < fhi)  /// try outside contraction
                {
                    double c[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      c[j] = 1.5*ave[j] - 0.5*p[ihi][j];
                    double fc = fnelder.func(c);
                    if (fc <= fr)
                    {
                        for (int j=0; j<NDIMS; j++)
                          p[ihi][j] = c[j];
                        p[ihi][FUNC] = fc;
                        continue;
                    }
                    else   /////// contract
                    {
                        for (int i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (int j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                        continue;
                    }
                }

                if (fr >= fhi)   /// over the top; try inside contraction
                {
                    double cc[] = new double[NDIMS];
                    for (int j=0; j<NDIMS; j++)
                      cc[j] = 0.5*ave[j] + 0.5*p[ihi][j];
                    double fcc = fnelder.func(cc);
                    if (fcc < fhi)
                    {
                        for (int j=0; j<NDIMS; j++)
                          p[ihi][j] = cc[j];
                        p[ihi][FUNC] = fcc;
                        continue;
                    }
                    else    ///////// contract
                    {
                        for (int i=0; i<NPTS; i++)
                          if (i != ilo)
                          {
                              for (int j=0; j<NDIMS; j++)
                                p[i][j] = 0.5*p[ilo][j] + 0.5*p[i][j];
                              p[i][FUNC] = fnelder.func(p[i]);
                          }
                    }
                }
            }
        return z;
}
      public static double[] nelder(f_xj fnelder,double x[][],double tolerance)
      {return nelder(fnelder,x,300,tolerance,0);}

      public static double[] nelder(f_xj fnelder,double x[][])
      {return nelder(fnelder,x,300,1.0e-15,0);}

      public static double[] nelder(f_xj fnelder,double x[][],int maxiteration,double tolerance)
      {return nelder(fnelder,x,maxiteration,tolerance,0);}

//derivatives of a function

public static double derivative(f_x f_deriv,double x)
{
// df/dx
// This method calculates derivatives of a simple function f(x)
// accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_x
double h0=0.0256;
int i,m;
int n=7;
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
T[i][0]=( f_deriv.func(x + h[i]) - f_deriv.func( x - h[i]))/(2.0*h[i]);
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/
          (h[i]*h[i]- h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}

public static double derivative(f_xi f_deriv,int equation_ref,
double x[],int x_ref)
{
// dfi/dxj
// This method calculates derivative of a function selected from a set of
// functions. Accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_xi
// df(equation_ref)/dx(x_ref)
double h0=0.0256;
int i,m;
int n=7;
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
T[i][0]=( f_deriv.func(x1,equation_ref) -
          f_deriv.func(x2,equation_ref))/(2.0*h[i]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] -
  h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i] - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}

public static double[] derivative(f_xj f_deriv,double x[])
{
// df/dxj  j=0...x.length
// This method calculates derivative of a function with more than one independent variable. 
// Accuracy of method can be adjusted by changing variables h0 and n
// function input should be in the form given in abstract class
// f_xj,j=0...x.length = df/dx(x_ref)
double a[]=new double[x.length];
for(int x_ref=0;x_ref<x.length;x_ref++)
   {
   a[x_ref]=derivative(f_deriv,x,x_ref);
   }
return a;   	   
}

public static double derivative(f_xj f_deriv,double x[],int x_ref)
{
// df/dxj
// This method calculates derivative of a function with more than one independent variable. 
// Accuracy of method can be adjusted by changing variables h0 and n
// function input should be in the form given in abstract class
// f_xj = df/dx(x_ref)
double h0=0.256808;
int i,m;
int n=4;
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
//h[i]=h0*Math.pow(r,i);
h0/=2.0;
h[i]=h0;
}
//first derivative (difference formula)
for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
T[i][0]=(f_deriv.func(x1)-f_deriv.func(x2))/(2.0*h[i]);
//System.out.print("T["+i+","+0+"] ="+T[i][0]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
//System.out.println("\n");

for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] -
  h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i] - h[i+m]*h[i+m]);
  //System.out.print("T["+i+","+m+"] ="+T[i][m]);
  }
  //System.out.println("\n");
}
double xx=T[0][n-1];
return xx;
}



public static double derivative(fi_xi f_deriv,
int equation_ref,double x[],int x_ref)
{
// This method calculates derivative of a function selected from a set of
// functions with respect to x[x_ref].
// Accuracy of method can be adjusted by changing
// variables h0 and n
// function input should be in the form given in abstract class
// f_xi
double h0=0.0256;
int i,m;
int n=7;
double f1[];
f1=new double[x.length];
double f2[];
f2=new double[x.length];
double x1[];
x1=new double[x.length];
double x2[];
x2=new double[x.length];
for(i=0;i<x.length;i++)
{
x1[i]=x[i];
x2[i]=x[i];
}
//derivative of a simple function
double T[][];
T=new double[n][n];
double h[];
h=new double[n];
//vector<double> h(n,0);
for(i=0;i<n;i++)
{
  h[i]=0;
  for(int j=0;j<n;j++)
    T[i][j]=0;
}
h[0]=h0;
double r=0.5;
for( i=1;i<n;i++)
{
h[i]=h0*Math.pow(r,i);
}

for(i=0;i<n;i++)
{
x1[x_ref]+=h[i];
x2[x_ref]-=h[i];
f1=f_deriv.func(x1);
f2=f_deriv.func(x2);
T[i][0]=( f1[equation_ref] - f2[equation_ref])/(2.0*h[i]);
x1[x_ref]=x[x_ref];
x2[x_ref]=x[x_ref];
}
for(m=1;m<n;m++)
{
  for(i=0;i<n-m;i++)
  {
  T[i][m]=(h[i]*h[i]*T[i+1][m-1] - h[i+m]*h[i+m]*T[i][m-1])/(h[i]*h[i]
  - h[i+m]*h[i+m]);
  }
}
double xx=T[0][n-1];
return xx;
}

public static double[][] gauss_legendre_coefficients(double x1,double x2,int n)
{
	//calculates legendre gauss-coefficients as coefficients of the integral
	//for n terms
	double EPS=3.0e-15;
	int m,j,i;
	double z1,z,xm,xl,pp,p3,p2,p1;
    //double x[]=new double[n];
    //double w[]=new double[n];
    double a[][]=new double[2][n];//a[0][i]=x[i]   a[1][i]=w[i]
	m=(n+1)/2;
	xm=0.5*(x2+x1);
	xl=0.5*(x2-x1);
	for (i=1;i<=m;i++)  {
		z=Math.cos(Math.PI*((i-0.25)/(n+0.5)));
		do {
			p1=1.0;
			p2=0.0;
			for (j=1;j<=n;j++) {
				p3=p2;
				p2=p1;
				p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j;
			}
			pp=n*(z*p1-p2)/(z*z-1.0);
			z1=z;
			z=z1-p1/pp;
		} while (Math.abs(z-z1) > EPS);
		a[0][i-1]=xm-xl*z;
		a[0][n-i]=xm+xl*z;
		a[1][i-1]=2.0*xl/((1.0-z*z)*pp*pp);
		a[1][n-i]=a[1][i-1];
			}
return a;
}

public static double integral(f_x f_xnt,double x1,double x2,int n)
{
//n : number of integral coefficients
// this routine first generates gauss legendre coefficients
// for [x1,x2] band 
// then calculates gauss legendre integral
double a[][]=new double[2][n];
double s[]=new double[n];
double w[]=new double[n];
int i;
a=gauss_legendre_coefficients(x1,x2,n);
for(i=0;i<n;i++)
  { s[i]=a[0][i];
    w[i]=a[1][i];
  }
double z=0;
for(i=0;i<n;i++)
{ z+=w[i]*f_xnt.fonksiyon(s[i]);}
return z;
}

public static double integral(f_x f_xnt,double a,double b)
{
//integral f(x)dx
//integral of a function by using gauss-legendre quadrature
//coefficients are pre-calculated for 60 terms for [-1,1]
//band then utilises variable transform
 double s[],w[];
 int i;
 s=new double[30];
 w=new double[30];
 s[ 0] = .15532579626752470000E-02;
 s[ 1] = .81659383601264120000E-02;
 s[ 2] = .19989067515846230000E-01;
 s[ 3] = .36899976285362850000E-01;
 s[ 4] = .58719732103973630000E-01;
 s[ 5] = .85217118808615820000E-01;
 s[ 6] = .11611128394758690000E+00;
 s[ 7] = .15107475260334210000E+00;
 s[ 8] = .18973690850537860000E+00;
 s[ 9] = .23168792592899010000E+00;
 s[10] = .27648311523095540000E+00;
 s[11] = .32364763723456090000E+00;
 s[12] = .37268153691605510000E+00;
 s[13] = .42306504319570830000E+00;
 s[14] = .47426407872234120000E+00;
 s[15] = .52573592127765890000E+00;
 s[16] = .57693495680429170000E+00;
 s[17] = .62731846308394490000E+00;
 s[18] = .67635236276543910000E+00;
 s[19] = .72351688476904450000E+00;
 s[20] = .76831207407100990000E+00;
 s[21] = .81026309149462140000E+00;
 s[22] = .84892524739665800000E+00;
 s[23] = .88388871605241310000E+00;
 s[24] = .91478288119138420000E+00;
 s[25] = .94128026789602640000E+00;
 s[26] = .96310002371463720000E+00;
 s[27] = .98001093248415370000E+00;
 s[28] = .99183406163987350000E+00;
 s[29] = .99844674203732480000E+00;
w[ 0] = .39840962480827790000E-02;
w[ 1] = .92332341555455000000E-02;
w[ 2] = .14392353941661670000E-01;
w[ 3] = .19399596284813530000E-01;
w[ 4] = .24201336415292590000E-01;
w[ 5] = .28746578108808720000E-01;
w[ 6] = .32987114941090080000E-01;
w[ 7] = .36877987368852570000E-01;
w[ 8] = .40377947614710090000E-01;
w[ 9] = .43449893600541500000E-01;
w[10] = .46061261118893050000E-01;
w[11] = .48184368587322120000E-01;
w[12] = .49796710293397640000E-01;
w[13] = .50881194874202750000E-01;
w[14] = .51426326446779420000E-01;
w[15] = .51426326446779420000E-01;
w[16] = .50881194874202750000E-01;
w[17] = .49796710293397640000E-01;
w[18] = .48184368587322120000E-01;
w[19] = .46061261118893050000E-01;
w[20] = .43449893600541500000E-01;
w[21] = .40377947614710090000E-01;
w[22] = .36877987368852570000E-01;
w[23] = .32987114941090080000E-01;
w[24] = .28746578108808720000E-01;
w[25] = .24201336415292590000E-01;
w[26] = .19399596284813530000E-01;
w[27] = .14392353941661670000E-01;
w[28] = .92332341555455000000E-02;
w[29] = .39840962480827790000E-02;
int n=30;
double z=0;
double x,y;
for(i=0;i<n;i++)
{
x=(b+a)/2.0+(b-a)/2.0*s[i];
y=f_xnt.func(x);
z+=(b-a)/2*w[i]*y;
}
for(i=0;i<n;i++)
{
x=(b+a)/2.0+(b-a)/2.0*(-s[i]);
y=f_xnt.func(x);
z+=(b-a)/2.0*w[i]*y;
}
return z;
}

public static void trap(f_x ff,double a,double h, int j,int m,double R[][])
{
    //  successive trapezoidal integration rule
    //  this program will be utilised in romberg integration
        double sum=0;
        int p;
        for(p=1;p<=m;p++)
          { sum+=ff.func(a+h*(2*p-1)); }
        R[j][0]=R[j-1][0]/2.0+h*sum;
}

public static double integral_romberg(f_x ff,double a,double b)
{
//integral f(x)dx
//romberg integration
int n=8;//increase n to increase accuracy
double R[][];
R=new double[n+1][n+1];
int m=1;
double h=b-a;
double close=1;
double tol=1e-40;
int j=0,k=0;
double ret=0;
R[0][0]=h/2.0*(ff.func(a)+ff.func(b));
do
{
  j++;
  h/=2.0;
  Numeric.trap(ff,a,h,j,m,R);
  m*=2;
  for(k=1;k<=j;k++)
  {
    R[j][k]=R[j][k-1]+(R[j][k-1]-R[j-1][k-1])/(Math.pow(4,k)-1.0);
    close=Math.abs(R[j-1][j-1]-R[j][j]);
    ret=R[j][k];
  }
}while(close>tol && j<n);
return ret;
}
//==================================================
//Finding Non-linear Roots of Functions

//
public static double SIGN(double a,double b)
{
return Math.abs(a)*b/Math.abs(b);
}

public static double[] zbrac(double x1,double x2,f_x y)
{
	//returns a bracketed root 
	//if a root is not bracketed, enlarge region outwards each time
	int NTRY=50;
	double a[]=new double[2];
	double FACTOR=1.6;
	int j;
	double f1,f2;
	if (x1 == x2) System.out.println("Bad initial range in ZBRAC");
	f1=y.func(x1);
	f2=y.func(x2);
	for (j=1;j<=NTRY;j++) {
		if (f1*f2 < 0.0) {a[0]=x1;a[1]=x2;return a;}
		if (Math.abs(f1) < Math.abs(f2))
			f1=y.func(x1 += FACTOR*(x1-x2));
		else
			f2=y.func(x2 += FACTOR*(x2-x1));
	}
	{a[0]=x1;a[1]=x2;return a;}
}

public static double[][] zbrak(double x1,double x2,f_x y,int n,int nb)
{
	// n : number of sundivision required
	// x1,x2 : limit value to region
	// nb=number of sun region with roots
	// f_x y function value
	int nbb,i;
	double x,fp,fc,dx;
	double xb[][]=new double[2][nb];
	nbb=nb;
	nb=0;
	dx=(x2-x1)/n;
	x=x1;
	fp=y.func(x1);
	for (i=1;i<=n;i++) {
		x+=dx;
		fc=y.func(x);
		if (fc*fp < 0.0) {
			xb[0][nb]=x-dx;
			xb[1][nb]=x;
			nb++;
		}
		fp=fc;
		if (nbb == nb) return xb;
	}
System.out.println("full search is compleded results may not be valid");	
return xb;
}


public static double ridder(double x1,double x2,f_x y)
{
  // root of a non-linear equation	
  //reference : Numerical Recipes in C, second edition, William H. Press, 
  //Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
  //cambridge university press	
	int MAXIT=30;
    double xacc=1.0e-5;int j;
	double fl,fh,xl,xh,swap,dx,del,f;
	double rtf=x1;
	fl=y.func(x1);
	fh=y.func(x2);
	if (fl*fh > 0.0) System.out.println("Root must be bracketed in RTFLSP");
	if (fl < 0.0) {
		xl=x1;
		xh=x2;
	} else {
		xl=x2;
		xh=x1;
		swap=fl;
		fl=fh;
		fh=swap;
	}
	dx=xh-xl;
	for (j=1;j<=MAXIT;j++) {
		rtf=xl+dx*fl/(fl-fh);
		f=y.func(rtf);
		if (f < 0.0) {
			del=xl-rtf;
			xl=rtf;
			fl=f;
		} else {
			del=xh-rtf;
			xh=rtf;
			fh=f;
		}
		dx=xh-xl;
		if (Math.abs(del) < xacc || f == 0.0) return rtf;
	}
	System.out.println("Maximum number of iterations exceeded in RTFLSP");
    return rtf; 
}


public static double brent(double x1,double x2,f_x y)
{
    //Root of a non-linear equation
	//reference : Numerical Recipes in C, second edition, William H. Press, 
    //Saul A. Teukolsky, William T. Vetterling, Brian P. Flannery
    //cambridge university press	
	double ITMAX=100;
    double EPS=3.0e-8;
	int iter;
	double a=x1,b=x2,c,d,e,min1,min2;
	double tol=1.0e-5;
	c=0;
	d=0;
	e=0;
	double fa=y.func(a),fb=y.func(b),fc,p,q,r,s,tol1,xm;
	if (fb*fa > 0.0) System.out.println("Root must be bracketed in ZBRENT");
	fc=fb;
	for (iter=1;iter<=ITMAX;iter++) {
		if (fb*fc > 0.0) {
			c=a;
			fc=fa;
			e=d=b-a;
		}
		if (Math.abs(fc) < Math.abs(fb)) {
			a=b;
			b=c;
			c=a;
			fa=fb;
			fb=fc;
			fc=fa;
		}
		tol1=2.0*EPS*Math.abs(b)+0.5*tol;
		xm=0.5*(c-b);
		if (Math.abs(xm) <= tol1 || fb == 0.0) return b;
		if (Math.abs(e) >= tol1 && Math.abs(fa) > Math.abs(fb)) {
			s=fb/fa;
			if (a == c) {
				p=2.0*xm*s;
				q=1.0-s;
			} else {
				q=fa/fc;
				r=fb/fc;
				p=s*(2.0*xm*q*(q-r)-(b-a)*(r-1.0));
				q=(q-1.0)*(r-1.0)*(s-1.0);
			}
			if (p > 0.0)  q = -q;
			p=Math.abs(p);
			min1=3.0*xm*q-Math.abs(tol1*q);
			min2=Math.abs(e*q);
			if (2.0*p < (min1 < min2 ? min1 : min2)) {
				e=d;
				d=p/q;
			} else {
				d=xm;
				e=d;
			}
		} else {
			d=xm;
			e=d;
		}
		a=b;
		fa=fb;
		if (Math.abs(d) > tol1)
			b += d;
		else
			b += (xm > 0.0 ? Math.abs(tol1) : -Math.abs(tol1));
		fb=y.func(b);
	}
	System.out.println("Maximum number of iterations exceeded in BRENT");
return b;
}

public static double newton(double x,f_x y,f_x dy)
{
//root of a non-linear equation	
// x in f(x)=0 given also df(x)/dx
// Newton - Raphson method  for single equation
// with single variable
// required function and its derivative

int nmax=500;
double tolerance=1.0e-15;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
dfx=dy.func(x);
x-=fx/dfx;
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}


public static double newton(double x,f_x y)
{
//root of a non-linear equation
// x in f(x)=0 , df(x)/dx is calcuated by derivative
// Newton - Raphson method  for single equation
// required function only derivative is calculated
// numerically by using method derivative

int nmax=500;
double tolerance=1.0e-15;
double fx,dfx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
dfx=Numeric.derivative(y,x);
x-=fx/dfx;
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}


public static double newton_2nd_derivative(double x,f_x y,f_x dy)
{
//root of a non-linear equation	
// x in f(x)=0 given also df(x)/dx d2f(x)/d2x is calculated by method
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
//function and first order derivative is required
int nmax=500;
double dx=1e-15;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
dfx=dy.func(x);
dfxm=dy.func(x-dx);
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}

public static double newton_2nd_derivative(double x,f_x y)
{
//root of a non-linear equation	
// x in f(x)=0 given, df(x)/dx is calculated by derivative method
// d2f(x)/d2x is calculated by method
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
// function is required, first order derivative calculated
// numerically by derivative method.
int nmax=500;
double dx=1e-3;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
dfx=Numeric.derivative(y,x);
dfxm=Numeric.derivative(y,(x-dx));
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}

public static double secant_2nd_derivative(double x,f_x y)
{
//root of a non-linear equation	
// x in f(x)=0
// Newton - Raphson type method for single equation
// includes 2nd order derivative calculations
// function should be supplied
int nmax=500;
double dx=1.0e-3;
double dx1=2.0e-3;
double x1m;
double tolerance=1.0e-15;
double fx,fxm,fxm1,dfx,dfxm,d2fx;
for(int i=0;i<nmax;i++)
{
fx=y.func(x);
fxm=y.func(x-dx);
fxm1=y.func(x-dx1);
dfx=(fx-fxm)/dx;
dfxm=(fx-fxm1)/dx1;
d2fx=-6.0/dx/dx*(fx-fxm)+2.0/dx*(2.0*dfx+dfxm);
x-=(fx/dfx+.5*fx*fx/(dfx*dfx*dfx)*d2fx);
if(Math.abs(fx)<tolerance)
{
return x;
}
}
return x;
}


public static double bisection(double xl,double xu,f_x y)
{
//root of a non-linear equation	
// x in f(x)=0 bracket xl and xu is given
//bisection method to find roots of a function y.func(x)
//function should be supplied
// defination of variables :
// xl : lower guess
// xu : upper guess
// xr : root estimate
// es : stopping criterion
// ea :approximate error
// maxit : maximum iterations
// iter  : number of iteration
double test;
double xr=0;
double es,ea;
double fxl,fxr;
int maxit=500,iter=0;
es=0.001;
ea=1.1*es;
while((ea>es)&&(iter<maxit))
{
xr=(xl+xu)/2.0;
iter++;
if((xl+xu)!=0)
  { ea=Math.abs((xu-xl)/(xl+xu))*100;}
fxl= y.func(xl);
fxr= y.func(xr);
test= fxl*fxr;
if(test==0.0)       ea=0;
else if(test<0.0)  xu=xr;
else
  {
  xl=xr;
  }
} //end of while
return xr;
}

public static double bisectionsearch(double xl,double xu,f_x y)
{
//root of a non-linear equation
// x in f(x)=0 bracket xl and xu is given
//bisection search method to find roots of a function y.func(x)
//function should be supplied
// defination of variables :
// xl : lower guess
// xu : upper guess
// xr : root estimate
// es : stopping criterion
// ea :approximate error
// maxit : maximum iterations
// iter  : number of iteration
int maxit=100,iter=0;
double es=0.000001;
double dx=Math.abs((xu-xl)/maxit);
double dxp=dx;
double xr=xl;
double fx=y.func(xr);
double fxr=fx;
for(int i=0;i<maxit;i++)
{
xr+=dxp;
if(xr>xu) break;
fx=y.func(xr);
  if(fx*fxr<0)
  {
    if(dxp<es) return xr;
    else
    {
    xr-=dxp;
    dxp/=2.0;
    }
  }
}
return xr;
}

public static double[] newton( double x[],fi_xi y,fij_xi dy)
{
//root of a non-linear system of equations 	
// xi solution matrix in fj(xi)=0 set of non-linear equation
// derivative matrix dfj/dxi is required
//solution of nonlinear system of equations
//by using newton raphson method.
//ti :weigth factor
//x independent value vector:
//y dependent function vector
//dy derivative of dependent function matrix
double ti=1.0;
int i;
int nmax=500;
double tolerance=1.0e-30;
int n=x.length;
double b[];
b=new double[n];
for(i=0;i<n;i++)
{
b[i]=1.0;
}
i=0;
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
    b=Matrix.multiply(Matrix.divide(y.func(x),dy.func(x)),-ti);
    x=Matrix.add(x,b);
}
if(i >= nmax) System.out.println("warning maximum number"+
              " of iterations results may not be valid");
return x;
}

public static double[] newton( double x[],fi_xi y)
{
// xi solution matrix in fj(xi)=0 set of non-linear equation
// derivative matrix dfj/dxi is calculated by derivative method
//solution of nonlinear system of equations
//by using newton raphson method.
//this function does not require derivatives
//it is utilised method derivative to calculate derivatives
double ti=1.0;
int i,ii,jj;
int nmax=500;
double tolerance=1.0e-15;
int n=x.length;
double b[];
b=new double[n];
double dy[][];
dy=new double[n][n];
i=0;
for(i=0;i<n;i++)
{
b[i]=1.0;
}
while( i++ < nmax && Matrix.abs(b) > tolerance )
{
  for(ii=0;ii<n;ii++)
  {
    for(jj=0;jj<n;jj++)
    {
        dy[ii][jj]=Numeric.derivative(y,ii,x,jj);
    }
  }
    b=Matrix.multiply(Matrix.divide(y.func(x),dy),-ti);
    x=Matrix.add(x,b);
}
if(i >= nmax) System.out.println("warning maximum number of iterations"+
                                 " results may not be valid"          );
return x;
}


//========= least square curve fitting methods ==========

//-------- polynomial least square curve fitting --------

public static double[] spline(double x[],double y[],double yp1,double ypn)
{
	//yp1 end condition
	//ypn end condition
	// x,y given data set
	int n=x.length;
	double y2[]=new double[n];
	int i,k;
	double p,qn,sig,un,u[];
    u=new double[n];   
	if (yp1 > 0.99e30)
		y2[0]=u[0]=0.0;
	else {
		y2[0] = -0.5;
		u[0]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
	}
	for (i=1;i<(n-1);i++) {
		sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
		p=sig*y2[i-1]+2.0;
		y2[i]=(sig-1.0)/p;
		u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
		u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
	}
	if (ypn > 0.99e30)
		qn=un=0.0;
	else {
		qn=0.5;
		un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n-1]-y[n-2])/(x[n-1]-x[n-2]));
	}
	y2[n-1]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
	for (k=n-1;k>=0;k--)
		y2[k]=y2[k]*y2[k+1]+u[k];
return y2;
}

public static double splint(double xa[],double ya[],double y2a[],double x)
{
	int n=xa.length;
	int klo,khi,k;
	double h,b,a;
	double y;
	klo=1;
	khi=n;
	while (khi-klo > 1) 
	{
		k=(khi+klo) >> 1;
		if (xa[k] > x) khi=k;
		else klo=k;
	}
	h=xa[khi]-xa[klo];
	if (h == 0.0) System.out.println("Bad XA input to routine SPLINT");
	a=(xa[khi]-x)/h;
	b=(x-xa[klo])/h;
	y=a*ya[klo]+b*ya[khi]+((a*a*a-a)*y2a[klo]+(b*b*b-b)*y2a[khi])*(h*h)/6.0;
    return y;
}


public static double[] poly_least_square(double xi[],double yi[],int n)
// data will be fit into
// y=f(x)=a0+a1*x+a2*x2+a3*x^3+....+an*x^n
// polynomial least square curve fitting
// variables xi,yi vector of data points
// degree of curve fitting
{
// data input variables:
//xi vector of independent variable
//yi vector of dependent variable
//n : degree of curve fitting
int l=xi.length;
int i,j,k;
int np1=n+1;
double A[][];
A=new double[np1][np1];
double B[];
B=new double[np1];
double X[];
X=new double[np1];

for(i=0;i<n+1;i++)
  {
  for(j=0;j<n+1;j++)
    {
    if(i==0 && j==0)
      A[i][j]=l;
    else
      for(k=0;k<l;k++)
        A[i][j] += Math.pow(xi[k],(i+j));
    }
  for(k=0;k<l;k++)
        {
        if(i==0)
        B[i]+= yi[k];
        else
        B[i] += Math.pow(xi[k],i)*yi[k];
        }
  }
  X=Matrix.divide(B,A);
//X=B/A;

double max=0;
for(i=0;i<n+1;i++)
 if(Math.abs(X[i]) > max) max = Math.abs(X[i]);
for(i=0;i<n+1;i++)
 if((Math.abs(X[i]/max) > 0) && (Math.abs(X[i]/max) < 1.0e-100)) X[i]=0;
return X;
}

public static double f_least_square(double e[],double x)
{
// this function calculates the value of
// least square curve fitting function
int n=e.length;
double ff;
if(n!=0.0)
  { ff=e[n-1];
  for(int i=n-2;i>=0;i--)
   { ff=ff*x+e[i]; }
  }
else
  ff=0;
return ff;
}

public static double error(double x[],double y[],double e[])
{
//calculates absolute square root error of a least square approach
double n=x.length;
 int k;
 double total=0;
 for(k=0;k<n;k++)
 {
 total+=(y[k]-f_least_square(e,x[k]))*(y[k]-f_least_square(e,x[k]));
 }
total=Math.sqrt(total);
return total;
}

//-------End least square methods --------------------
//-------Orthogonal polynomial least square curve fitting ---

public static double[][] orthogonal_polynomial_least_square
(double xi[],double fi[],int m)
//orthogonal polynomial least square curve fitting
//this method do not require any matrix solution
//variables xi,fi vector of data points wi weight functions
//m degree of curve fitting
{
int i,j,k;
int n=xi.length;
int mp2=n+2;
int mp1=n+1;
//matrix<double> p(m+2,n,0);
double p[][];
p=new double[mp2][n];
//vector<double> gamma(m+1,0);
double gamma[];
gamma=new double[mp1];
//vector<double> beta(m+1,0);
double beta[];
beta=new double[mp1];
//vector<double> omega(m+1,0);
double omega[];
omega=new double[mp1];
//vector<double> alpha(m+1,0);
double alpha[];
alpha=new double[mp1];
//vector<double> b(m+1,0);
double b[];
b=new double[mp1];
//vector<double> wi(n,1.);
double wi[];
wi=new double[n];
//matrix<double> a(3,m+1);
double a[][];
a=new double[3][mp1];
for(i=0;i<n;i++)
  {
   p[1][i]=1.0;
   p[0][i]=0.0;
  }
gamma[0]=0;
for(i=0;i<n;i++)
  {
  gamma[0]+=wi[i];
  }
beta[0]=0.0;
for(j=0;j<m+1;j++)
  {
  omega[j]=0;
  for(i=0;i<n;i++)
    {
    omega[j]+=wi[i]*fi[i]*p[j+1][i];
    }
  b[j]=omega[j]/gamma[j];
if( j != m)
    {
    alpha[j+1]=0;
    for(i=0;i<n;i++)
      {
      alpha[j+1]+=wi[i]*xi[i]*p[j+1][i]*p[j+1][i]/gamma[j];
      }
    for(i=0;i<n;i++)
      {
      p[j+2][i]=(xi[i]-alpha[j+1])*p[j+1][i]-beta[j]*p[j][i];
      }
      gamma[j+1]=0;
    for(i=0;i<n;i++)
      {
      gamma[j+1]+=wi[i]*p[j+2][i]*p[j+2][i];
      }
    beta[j+1]=gamma[j+1]/gamma[j];
    }
}//end of j
for(j=0;j<m+1;j++)
  {
  a[0][j]=b[j];
  a[1][j]=alpha[j];
  a[2][j]=beta[j];
  }
return a;
}

public static double f_orhogonal_polynomial(double a[][],double x)
{

double yy=0;
int k;
int m=a[0].length-1;
int mp2=m+2;
double q[];
q=new double[mp2];
//vector<double> q(m+2,0.0);
for(k=m-1;k>=0;k--)
  {
  q[k]=a[0][k]+(x-a[1][k+1])*q[k+1]-a[2][k+1]*q[k+2];
  yy=q[k];
  }
return yy;
}


// differential equation solutions

public static double[] RK4(f_xi fp,double x0,double xn,double f0,int N)
{
//4th order Runge Kutta Method
//fp :given derivative function df/dx(f,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
double fi[]=new double[N+1];
double xi[]=new double[2];
double k[];
k=new double[4];
int i;
double x;
fi[0]=f0;
for(x=x0,i=0;x<xn;x+=h,i++)
  {
  xi[0]=x;
  xi[1]=fi[i];
   k[0]=h*fp.func(xi,0);
  xi[0]=x+h/2.0;
  xi[1]=fi[i]+k[0]/2;
   k[1]=h*fp.func(xi,0);
  xi[1]=fi[i]+k[1]/2.0;
   k[2]=h*fp.func(xi,0);
  xi[0]=x+h;
  xi[1]=fi[i]+k[2];
   k[3]=h*fp.func(xi,0);
  fi[i+1]=fi[i]+k[0]/6.0+k[1]/3.0+k[2]/3.0+k[3]/6.0;
  }
return fi;
}

public static double[][] RK4(f_xi fp,double x0,double xn,double f0[],int N)
{
//4th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[4];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[j]=fi[j-1][i]+k[1]/2.0;
     k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]+k[2];
     k[3]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/6.0+k[1]/3.0+k[2]/3.0+k[3]/6.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {
  a[j][i]=fi[j-1][i];
  }
}

return a;
}


public static double[][] RK6(f_xi fp,double x0,double xn,double f0[],int N)
{
//6th order Runge Kutta Method
//fp  : given set of derivative functions dfj/dxi(fj,x)
// xo : initial value of the independent variable
// xn : final value of the independent variable
// f0 : initial value of the dependent variable
// N  : number of dependent variable to be calculated
// fi : dependent variable
double h=(xn-x0)/N;
int M=f0.length;
double fi[][];
fi=new double[M][N+1];
double xi[]=new double[M+1];
double k[]=new double[6];
int i,j;
double x;
for(j=0;j<M;j++)
  {
  fi[j][0]=f0[j];
  xi[j+1]=f0[j];
  }
for(x=x0,i=0;i<N;x+=h,i++)
  {
  for(j=1;j<=M;j++)
    {
    xi[0]=x;
    xi[j]=fi[j-1][i];
    k[0]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/2;
    k[1]=h*fp.func(xi,j-1);
    xi[0]=x+h/2.0;
    xi[j]=fi[j-1][i]+k[0]/4.0+k[1]/4.0;
    k[2]=h*fp.func(xi,j-1);
    xi[0]=x+h;
    xi[j]=fi[j-1][i]-k[1]+2.0*k[2];
    k[3]=h*fp.func(xi,j-1);
    xi[0]=x+2.0/3.0*h;
    xi[j]=fi[j-1][i]+7.0/27.0*k[0]+10.0/27.0*k[1]+1.0/27.0*k[3];
    k[4]=h*fp.func(xi,j-1);
    xi[0]=x+1.0/5.0*h;
    xi[j]=fi[j-1][i]+28.0/625.0*k[0]-1.0/5.0*k[1]+546.0/625.0*k[2]+54.0/625.0*k[3]-378/625.0*k[4];
    k[5]=h*fp.func(xi,j-1);
    fi[j-1][i+1]=fi[j-1][i]+k[0]/24.0+5.0*k[3]/48.0+27.0*k[4]/56.0+125.0*k[5]/336.0;
    xi[j]=fi[j-1][i];
    }
  }
double a[][]=new double[M+1][N+1];
for(x=x0,i=0;i<=N;x+=h,i++)
{
  a[0][i]=x;
  for(j=1;j<=M;j++)
  {
  a[j][i]=fi[j-1][i];
  }
}

return a;
}


public static double[][]
RKF45(f_xi fp,double x0,double xn,double f0,int N) throws IOException
{
// Runge Kutta - Fehlberg Method
// fp  : derivative function df/dx(f,x)  (defined as a class)
// Tol :Tolerance
//Hmax : maximum possible step size
//Hmin : minimum possible step size
//k    : Runge kutta coefficients
//Err  : error
//x[i] : input variable to fp = df/dx(f,x) = df/dx(x[i])
//j    : actual step size
//fi[][]:solution matrix
//fj[][]:solution matrix in exact size(j) (return function)
double h=(xn-x0)/N;
double fi[][]=new double[2][1000];
double Tol=2.0e-3;
double Hmin=h/64.0;
double Hmax=h*64.0;
double k[]=new double[6];
double Err;
double s;
double xi[]=new double[2];
int j=0;
fi[0][j]=x0;
fi[1][j]=f0;
while(fi[0][j] < xn )
  {
  if(fi[0][j]+h > xn) h=xn-fi[0][j];
  xi[0]=fi[0][j];
  xi[1]=fi[1][j];
   k[0]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+h/4.0;
  xi[1]=fi[1][j]+k[0]/4.0;
   k[1]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+3.0/8.0*h;
  xi[1]=fi[1][j]+3.0/32.0*k[0]+9.0/32.0*k[1];
   k[2]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+12.0/13.0*h;
  xi[1]=fi[1][j]+1932.0/2197.0*k[0]-7200.0/2197.0*k[1]+7296.0/2197.0*k[2];
   k[3]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+h;
  xi[1]=fi[1][j]+439.0/216.0*k[0]-8.0*k[1]+3680.0/513.0*k[2]-845.0/4104.0*k[3];
   k[4]=h*fp.func(xi,0);
  xi[0]=fi[0][j]+0.5*h;
  xi[1]=fi[1][j]-8.0/27.0*k[0]+2.0*k[1]
  -3544/2565*k[2]+1859.0/4104.0*k[3]-11.0/40.0*k[4];
   k[5]=h*fp.func(xi,0);
   Err=Math.abs(1.0/360.0*k[0]-128/4275*k[2]
   -2197.0/75240*k[3]+1.0/5.0*k[4]+2.0/55.0*k[5]);
  if(Err<Tol || h<2*Hmin)
     {
     fi[0][j+1]=fi[0][j]+h;
     fi[1][j+1]=fi[1][j]+16.0/135.0*k[0]+6656.0/12825.0*k[2]+
     28561.0/56430.0*k[3]-9.0/50.0*k[4]+2.0/55.0*k[5];
     j++;
     }
  else
     {
     if(Err==0) s=0;
     else s=0.84*(Math.pow(Tol*h/Err,0.25));
     if(s<0.75 && h>(2*Hmin)  ) h/=2.0;
     if(s>1.5  && Hmax >(2*h) ) h*=2.0;
     }
  }
double fj[][]=new double[2][j+1];
for(int i=0;i<=j;i++)
  {
  fj[0][i]=fi[0][i];
  fj[1][i]=fi[1][i];
  }
return fj;
}

//special functions section

public static double gammln(double xx)
{
	//gamma function
	double x,tmp,ser;
	double cof[]={76.18009173,-86.50532033,24.01409822,
		-1.231739516,0.120858003e-2,-0.536382e-5};
	int j;

	x=xx-1.0;
	tmp=x+5.5;
	tmp -= (x+0.5)*Math.log(tmp);
	ser=1.0;
	for (j=0;j<=5;j++) {
		x += 1.0;
		ser += cof[j]/x;
	}
	return -tmp+Math.log(2.50662827465*ser);
}

public static double factrl(int n)
{
	// n! of an integer number
	int ntop=4;
	double a[]=new double[33];
	a[0]=1.0;
	a[1]=1.0;
	a[2]=2.0;
	a[3]=6.0;
	int j;
	if (n < 0) System.out.println("Negative factorial in routine FACTRL");
	if (n > 32) return Math.exp(gammln(n+1.0));
	while (ntop<n) {
		j=ntop++;
		a[ntop]=a[j]*ntop;
	}
	return a[n];
}

public static double factln(int n)
{
	//Math.og(factorial) of an integer number
	double a[]=new double[33];
	int ntop=4;
	a[0]=1.0;
	a[1]=1.0;
	a[2]=2.0;
	a[3]=6.0;
	int j;
	if (n < 0) System.out.println("Negative factorial in routine FACTLN");
	if (n <= 1) return 0.0;
	if (n > 32) return gammln(n+1.0);
	while (ntop<n) {
		j=ntop++;
		a[ntop]=a[j]*ntop;
	}
	return Math.log(a[n]);
}

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

public static double beta(double z,double w)
{
	//gamma(z)*gamma(w)/gamma(z+w);
	return Math.exp(gammln(z)+gammln(w)-gammln(z+w));
}

public static double gammp(double a,double x)
{
	//incomplete gamma function P(a,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;
	}
}

public static double gammq(double a,double x)
{
	//incomplete gamma function Q(a,x)=1-P(a,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=100;
	double EPS=3.0e-7;
	int n;
	double sum,del,ap;
	double gln=gammln(a);
	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=100;
	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=gammln(a);
	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;
}

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  erfc(double x)
{
	//complimentary error function
	return x < 0.0 ? 1.0+gammp(0.5,x*x) : gammq(0.5,x*x);
}

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;
}

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;
}

double betai(double a,double b,double x)
{
	//incomplete beta function
	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(gammln(a+b)-gammln(a)-gammln(b)+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;
}


static String toString(int n, int w)
// converts an int to a string with given width.
{
    String s = Integer.toString(n);
    while (s.length() < w)
       s = " " + s;
    return s;
}

static String toString(int left)
// converts an int to a string with a constant width.
{
    return toString(left,6);
}

public static String toString(double left, int w, int d)
// converts a double to a string with given width and decimals.
{
	NumberFormat df=NumberFormat.getInstance(Locale.US);
	df.setMaximumFractionDigits(d);
    df.setMinimumFractionDigits(d);
    df.setGroupingUsed(false);
    String s = df.format(left);
    while (s.length() < w)
      s = " " + s;
    if (s.length() > w)
    {
        s = "";
        for (int i=0; i<w; i++)
          s = s + "-";
    }
    return s;
}

public static String toString(double left)
{// converts a double to a string with a constant width and constant decimals.
  return toString(left,15,10);}

//End of Class Numeric
}