// Big_nuderiv.java  non uniformly spaced derivative coefficients 
//    given  x0, x1, x2, x3  non uniformly spaced
//    find the coefficients of y0, y1, y2, y3
//    in order to compute the derivatives:
//      y'(x0) = c00*y0 + c01*y1 + c02*y2 + c03*y3
//      y'(x1) = c10*y0 + c11*y1 + c12*y2 + c13*y3
//      y'(x2) = c20*y0 + c21*y1 + c22*y2 + c23*y3
//      y'(x3) = c30*y0 + c31*y1 + c32*y2 + c33*y3
//
//  Method:  fit data to  y(x)    = a + b*x + c*x^2 + d*x^3
//                        y'(x)   =     b   + 2*c*x + 3*d*x^2
//                        y''(x)  =           2*c   + 6*d*x
//                        y'''(x) =                   6*d
//
//           with y0, y1, y2 and y3 variables:
//
//                              y0  y1  y2  y3 
//  form:  | 1  x0  x0^2  x0^3   1   0   0   0 | = a
//         | 1  x1  x1^2  x1^3   0   1   0   0 | = b
//         | 1  x2  x2^2  x2^3   0   0   1   0 | = c
//         | 1  x3  x3^2  x3^3   0   0   0   1 | = d
//
// reduce  | 1   0     0     0  a0  a1  a2  a3 | = a
//         | 0   1     0     0  b0  b1  b2  b3 | = b
//         | 0   0     1     0  c0  c1  c2  c3 | = c
//         | 0   0     0     1  d0  d1  d2  d3 | = d
//
//                              this is just the inverse
//
//   y'(x0) =  b0*y0        + b1*y1        + b2*y2        + b3*y3        +
//             2*c0*y0*x0   + 2*c1*y1*x0   + 2*c2*y2*x0   + 2*c3*y3*x0   +
//             3*d0*y0*x0^2 + 3*d1*y1*x0^2 + 3*d2*y2*x0^2 + 3*d3*y3*x0^2
//
//   or  c00 = b0 + 2*c0*x0 + 3*d0*x0^2
//       c01 = b1 + 2*c1*x0 + 3*d1*x0^2
//       c02 = b2 + 2*c2*x0 + 3*d2*x0^2
//       c03 = b3 + 2*c3*x0 + 3*d3*x0^2
//
//   y'(x0) = c00*y0 + c01*y1 + c02*y2 +c03*y3
//
//   y'(x1) =  b0*y0        + b1*y1        + b2*y2        + b3*y3        +
//             2*c0*y0*x1   + 2*c1*y1*x1   + 2*c2*y2*x1   + 2*c3*y3*x1   +
//             3*d0*y0*x1^2 + 3*d1*y1*x1^2 + 3*d2*y2*x1^2 + 3*d3*y3*x1^2
//
//   or  c10 = b0 + 2*c0*x1 + 3*d0*x1^2
//       c11 = b1 + 2*c1*x1 + 3*d1*x1^2
//       c12 = b2 + 2*c2*x1 + 3*d2*x1^2
//       c13 = b3 + 2*c3*x1 + 3*d3*x1^2
//
//       cij = bj + 2*cj*xi + 3*dj*xi^2
//
//
//   y'(x1) = c10*y0 + c11*y1 + c12*y2 +c13*y3
//
//   y''(x0) = 2*c0*y0    + 2*c1*y1    + 2*c2*y2    + 2*c3*y3    +
//             6*d0*y0*x0 + 6*d1*y1*x0 + 6*d2*y2*x0 + 6*d3*y3*x0
//
//   or  c00 = 2*c0 + 6*d0*x0
//       c01 = 2*c1 + 6*d1*x0
//       c02 = 2*c2 + 6*d2*x0
//       c03 = 2*c3 + 6*d3*x0
//
//   y'''(x0) = 6*d0*y0 + 6*d1*y1 + 6*d2*y2 + 6*d3*y3
//
//   or  c00 = 6*d0    independent of x
//       c01 = 6*d1
//       c02 = 6*d2
//       c03 = 6*d3

import java.math.BigDecimal;

class Big_nuderiv
{
  int nbits = 332; // precision in bits, 100 digits

  Big_nuderiv(int order, int npoint, int point, double x[], double c[])
  {
    int n = npoint;
    BigDecimal A[][] = new BigDecimal[n][n];
    BigDecimal fct[] = new BigDecimal[n];
    BigDecimal pwr;

    for(int i=0; i<n; i++) // build matrix that will be inverted 
    {
      pwr=new BigDecimal(1.0);
      for(int j=0; j<n; j++)
      {
        A[i][j] = pwr;
        pwr = pwr.multiply(new BigDecimal(x[i]));
      }
    }
    Big_invert(A); // invert the matrix 
    // form derivative for order and point 
    for(int i=0; i<n; i++) // compute factors for order 
    {
      fct[i] = new BigDecimal(1.0);
    }
    for(int j=0; j<order; j++)
    {
      for(int i=0; i<n; i++)
      {
	fct[i] = fct[i].multiply(new BigDecimal(i-j));
      }
    }
    for(int i=0; i<n; i++)
    {
      c[i] = 0.0;
      pwr = new BigDecimal(1.0);
      for(int j=order; j<n; j++)
      { 
	  c[i] = c[i] + (fct[j].multiply(A[j][i].multiply(pwr))).doubleValue();
	pwr = pwr.multiply(new BigDecimal(x[point]));
      }
    }
  } // end Big_nuderiv constructor (double)cx

  Big_nuderiv(int order, int npoint, int point, double x[], BigDecimal c[])
  {
    int n = npoint;
    BigDecimal A[][] = new BigDecimal[n][n];
    BigDecimal fct[] = new BigDecimal[n];
    BigDecimal pwr;

    for(int i=0; i<n; i++) // build matrix that will be inverted 
    {
      pwr=new BigDecimal(1.0);
      for(int j=0; j<n; j++)
      {
        A[i][j] = pwr;
        pwr = pwr.multiply(new BigDecimal(x[i]));
      }
    }
    Big_invert(A); // invert the matrix 
    // form derivative for order and point 
    for(int i=0; i<n; i++) // compute factors for order 
    {
      fct[i] = new BigDecimal(1.0);
    }
    for(int j=0; j<order; j++)
    {
      for(int i=0; i<n; i++)
      {
	  fct[i] = fct[i].multiply(new BigDecimal(i-j));
      }
    }
    for(int i=0; i<n; i++)
    {
      c[i] = new BigDecimal(0);
      pwr = new BigDecimal(1.0);
      for(int j=order; j<n; j++)
      { 
	  c[i] = c[i].add((fct[j].multiply(A[j][i].multiply(pwr))));
	  pwr = pwr.multiply(new BigDecimal(x[point]));
      }
    }
  } // end Big_nuderiv constructor (bigDecimal)cx

  void Big_invert(BigDecimal AI[][]) // in place
  {
    int n = AI.length;
    int row[] = new int[n];
    int col[] = new int[n];
    BigDecimal temp[] = new BigDecimal[n];
    int hold , I_pivot , J_pivot;
    BigDecimal pivot, abs_pivot;
    BigDecimal one = new BigDecimal("1.0");
    BigDecimal zero = new BigDecimal("0.0");

    one = one.setScale(nbits,BigDecimal.ROUND_DOWN);

    // set up row and column interchange vectors
    for(int k=0; k<n; k++)
    {
      row[k] = k ;
      col[k] = k ;
    }
    // begin main reduction loop
    for(int k=0; k<n; k++)
    {
      // find largest element for pivot
      pivot = AI[row[k]][col[k]] ;
      I_pivot = k;
      J_pivot = k;
      for(int i=k; i<n; i++)
      {
        for(int j=k; j<n; j++)
        {
          abs_pivot = pivot.abs();
          if((AI[row[i]][col[j]].abs()).compareTo(abs_pivot)>0)
          {
            I_pivot = i ;
            J_pivot = j ;
            pivot = AI[row[i]][col[j]] ;
          }
        }
      }
      if((pivot.abs()).compareTo(new BigDecimal(1.0E-300)) <0)
      {
        System.out.println("Matrix is singular !");
        return;
      }
      hold = row[k];
      row[k]= row[I_pivot];
      row[I_pivot] = hold ;
      hold = col[k];
      col[k]= col[J_pivot];
      col[J_pivot] = hold ;
       // reduce about pivot
      AI[row[k]][col[k]] = one.divide(pivot, BigDecimal.ROUND_DOWN);
      AI[row[k]][col[k]] = AI[row[k]][col[k]].
                           setScale(nbits,BigDecimal.ROUND_DOWN);
      for(int j=0; j<n; j++)
      {
        if(j != k)
        {
	  AI[row[k]][col[j]] = AI[row[k]][col[j]].multiply(AI[row[k]][col[k]]);
          AI[row[k]][col[j]] = AI[row[k]][col[j]].
                               setScale(nbits,BigDecimal.ROUND_DOWN);
        }
      }
      // inner reduction loop
      for(int i=0; i<n; i++)
      {
        if(k != i)
        {
          for(int j=0; j<n; j++)
          {
            if( k != j )
            {
	      AI[row[i]][col[j]] = AI[row[i]][col[j]].subtract(
				   AI[row[i]][col[k]].multiply(
                                   AI[row[k]][col[j]])) ;
              AI[row[i]][col[j]] = AI[row[i]][col[j]].
                                   setScale(nbits,BigDecimal.ROUND_DOWN);
            }
          }
          AI[row[i]][col[k]] = (AI[row[i]][col[k]].multiply(
                                AI[row[k]][col[k]])).negate();
          AI[row[i]][col[k]] = AI[row[i]][col[k]].
                               setScale(nbits,BigDecimal.ROUND_DOWN);
        }
      }
    }
    // end main reduction loop

    // unscramble rows
    for(int j=0; j<n; j++)
    {
      for(int i=0; i<n; i++)
      {
        temp[col[i]] = AI[row[i]][j];
      }
      for(int i=0; i<n; i++)
      {
        AI[i][j] = temp[i] ;
      }
    }
    // unscramble columns
    for(int i=0; i<n; i++)
    {
      for(int j=0; j<n; j++)
      {
        temp[row[j]] = AI[i][col[j]] ;
      }
      for(int j=0; j<n; j++)
      {
        AI[i][j] = temp[j] ;
      }
    }
  } // end Big_invert

} // end class Big_nuderiv