// rderiv.java
//       compute formulas for numerical derivatives on hx equal spaces points
//       order is order of derivative, 1 = first derivative, 2 = second
//       points is number of points where value of function is known
//       f(x0), f(x1), f(x2) ... f(x points-1)
//       term is the point where derivative is computed
//       f'(x0) =  a(0)*f(x0) + a(1)*f(x1)
//                    + ... + a(points-1)*f(x points-1)
//
// algorithm: use divided differences to get polynomial p(x) that           
//            approximates  f(x).  f(x)=p(x)+error term                     
//            f'(x) = p'(x) + error term'                                   
//            substitute xj = x0 + j*h                                      
//            substitute x = x0 to get  p'(x0) etc   

public class rderiv  // for function  deriv
{
  rderiv(int order, int npoints, int point, double hx, double c[])
  {
    // compute the exact coefficients to numerically compute a derivative 
    // of order 'order' using 'npoints' at ordinate point,                
    // where order>=1, npoints>=order+1, 0 <= point < npoints,            
    // 'a' array returned with numerator of coefficients,                 
    // 'bb' returned with denominator of h^order                          
  
    int a[]=new int[npoints];  // integer coefficients of solution 
    int h[]=new int[npoints+1];  // coefficients of h 
    int x[]=new int[npoints+1];  // coefficients of x, variable for differentiating 
    int numer[][]=new int[npoints][npoints]; // numerator of a term numer[term][pos] 
    int denom[]=new int[npoints];        // denominator of coefficient 
    int jj, b;
    int k = 0;
    int ipower, isum, iat, r;
    int debug = 0;
    double bh;

    if(npoints<=order)
    {
      System.out.println("ERROR in call to deriv, npoints="+npoints+
                         " < order="+order+"+1");
      return;
    }
    for(int term=0; term<npoints; term++)
    {
      denom[term] = 1;
      jj = 0;
      for(int i=0; i<npoints-1; i++)
      {
        if(jj==term) jj++; // no index j in jth term 
        numer[term][i] = -jj; // from (x-x[j]) 
        denom[term] = denom[term]*(term-jj);
        jj++;
      }
    } // end term setting up numerator and denominator 
    if(debug>0)
      for(int i=0; i<npoints; i++)
      {
        System.out.println("denom["+i+"]="+denom[i]);
        for(int j=0; j<npoints-1; j++)
          System.out.println("numer["+i+"]["+j+"]="+numer[i][j]);
      }
        
    // substitute xj = x0 + j*h, just keep j value in numer[term][i] 
    // don't need to count x0 because same as x's 
    // done by above setup 
      
    // multiply out polynomial in x with coefficients of h 
    // starting from (x+j1*h)*(x+j2*h) * ... 
    // then differentiate d/dx 
    for(int term=0; term<npoints; term++)
    {
      x[0] = 1; // initialize 
      x[1] = 0;
      h[1] = 0;
      k = 1;
      for(int i=0; i<npoints-1; i++)
      {
        for(int j=0; j<k; j++)            // multiply by (x + j h) 
        {
          h[j] = x[j]*numer[term][i]; // multiply be constant x[j] 
        }
        for(int j=0; j<k; j++)         // multiply by x 
        {
          h[j+1] = h[j+1]+x[j]; // multiply by x and add 
        }
        k++;
        for(int j=0; j<k; j++) x[j] = h[j];
        x[k] = 0;
        h[k] = 0;
      }
      for(int j=0; j<k; j++)
      {
        numer[term][j] = x[j];
        if(debug>0)System.out.println("p(x)= "+numer[term][j]+
                                      " x^"+j+" at term="+term);
                          
      } // have p(x) for this 'term' 

      // differentiate 'order' number of times 
      for(int jorder=0; jorder<order; jorder++)
      {
        for(int i=1; i<k; i++) // differentiate p'(x) 
        { 
          numer[term][i-1] = i*numer[term][i];
        }
        k--;
      } // have p^(order) for this term 
      if(debug>0)
        for(int i=0; i<k; i++)
	    System.out.println("p^("+order+")(x)= "+numer[term][i]+
                               " x^"+i+" at term="+term);
                 
    } // end 'term' loop  for one 'order', for one 'npoints' 

    // substitute each point for x, evaluate, get coefficients of f(x[j]) 
    iat = point; // evaluate at x[iat], substitute 
    for(int jterm=0; jterm<npoints; jterm++) // get each term 
    {
      ipower = iat;
      isum = numer[jterm][0];
      for(int j=1; j<k; j++)
      {
        isum = isum + ipower * numer[jterm][j];
        ipower = ipower * iat;
      }
      a[jterm] = isum;
      if(debug>0)
	  System.out.println("f^("+order+")(x["+iat+"]) = (1/h^"+order+
                             ") ("+a[jterm]+"/"+denom[jterm]+" f(x["+
                             jterm+"]) +");
    }
    if(debug>0) System.out.println(" ");
    b = 0;
    for(int jterm=0; jterm<npoints; jterm++) // clean up each term 
    {
      if(a[jterm]==0) continue;
      r = gcd(a[jterm],denom[jterm]);
      a[jterm] = a[jterm] / r;
      denom[jterm] = denom[jterm] /r;
      if(denom[jterm]<0)
      {
        denom[jterm] = -denom[jterm];
        a[jterm] = - a[jterm];
      }
      if(denom[jterm]>b) b=denom[jterm]; // largest denominator 
    }
    for(int jterm=0; jterm<npoints; jterm++) // check for divisor 
    {
      r = (a[jterm] * b) / denom[jterm];
      if(r * denom[jterm] != a[jterm] * b)
      {
        r = gcd(b, denom[jterm]);
        b = b * (denom[jterm] / r);
      }
    }    
    for(int jterm=0; jterm<npoints; jterm++) // adjust for divisor 
    {
      a[jterm] = (a[jterm] * b) / denom[jterm];
      if(debug>0)
	  System.out.println("f^("+order+")(x["+iat+"])=(1/"+b+"h^"+order+
                             ")("+a[jterm]+" f(x["+jterm+"]) +");
    } // end computing terms of coefficients 
    bh = 1.0/(double)b;
    for(int i=0; i<order; i++) bh = bh/hx;
    for(int i=0; i<npoints; i++) c[i] = bh*(double)a[i];
  } // end rderiv 

  int gcd(int a, int b)  // needed by deriv  
  {
    int a1, b1, q, r;

    if(a==0 || b==0) return 1;
    if(Math.abs(a)>Math.abs(b)) { a1 = Math.abs(a); b1 = Math.abs(b); }
    else                        { a1 = Math.abs(b); b1 = Math.abs(a); }
    r=1;
    while(r!=0)
    {
      q = a1 / b1;
      r = a1 - q * b1;
      a1 = b1;
      b1 = r;
    }
    return a1;
  } // end gcd 

} // end class rderiv.java