// Nderiv.java
//           compute formulas for numerical derivatives                     
//           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) = (1/bh)*( 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   

import java.text.*;

public class Nderiv
{
  Nderiv()
  {
    System.out.println("Nderiv");
    int norder = 3;
    int npoints = 9;
    int order, points, term, jterm;
    int a[] = new int[50];  // coefficients of f(x0), f(x1), ... 
    int b[] = new int[1] ;  // use only b[0], output parameter
  
    System.out.println("compute formulas for numerical derivatives");  
    for(order=1; order<=norder; order++)
    {
      for(points=order+1; points<=npoints; points++)
      {
        for(term=0; term<points; term++)
        {
          System.out.println("computing order="+order+", npoints="+
                              points+", at term="+term);
          deriv(order, points, term, a, b); // retuens a and b
        
          for(jterm=0; jterm<points; jterm++)
          {
            if(jterm==0)
              System.out.println("f^("+order+")(x["+term+"])=(1/"+
                                 b[0]+"h^"+order+")( "+a[jterm]+
                                 " * f(x["+jterm+"] + ");
            else if(jterm==points-1)
              System.out.println("                      "+a[jterm]+
                                 " * f(x["+jterm+"] ) + O(1/h^"+
                                 points+") ");
            else
              System.out.println("                      "+a[jterm]+
                                 " * f(x["+jterm+"] + ");
          }
          check(order, points, term, a, b[0]);
          System.out.println();
        }
        System.out.println();
      }  // end points loop 
      System.out.println();
    } // end order loop 
  }

  int gcd(int a, int b)
  {
    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 


  void deriv(int order, int npoints, int point, int a[], int bb[])
  {
    // 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 h[]=new int[100];    // coefficients of h 
    int x[]=new int[100];    // coefficients of x, variable for differentiating 
    int numer[][]=new int[100][100]; // numerator of a term numer[term][pos] 
    int denom[]=new int[100];      // denominator of coefficient 
    int i, j, k=0, term, b;
    int jorder, ipower, isum, iat, jterm, r;
    int debug = 0;
  
    if(npoints<=order)
    {
      System.out.println("ERROR in call to deriv, npoints < order+1");
      return;
    }
    for(term=0; term<npoints; term++)
    {
      denom[term] = 1;
      j = 0;
      for(i=0; i<npoints-1; i++)
      {
        if(j==term) j++; // no index j in jth term 
        numer[term][i] = -j; // from (x-x[j]) 
        denom[term] = denom[term]*(term-j);
        j++;
      }
    } // end term setting up numerator and denominator 
    if(debug>0)
      for(i=0; i<npoints; i++)
      {
        System.out.println("denom["+i+"]="+denom[i]);
        for(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(term=0; term<npoints; term++)
    {
      x[0] = 1; // initialize 
      x[1] = 0;
      h[1] = 0;
      k = 1;
      for(i=0; i<npoints-1; i++)
      {
        for(j=0; j<k; j++)            // multiply by (x + j h) 
        {
          h[j] = x[j]*numer[term][i]; // multiply be constant x[j] 
        }
        for(j=0; j<k; j++)            // multiply by x 
        {
          h[j+1] = h[j+1]+x[j];       // multiply by x and add 
        }
        k++;
        for(j=0; j<k; j++) x[j] = h[j];
        x[k] = 0;
        h[k] = 0;
      }
      for(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(jorder=0; jorder<order; jorder++)
      {
        for(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(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 coeficients of f(x[j]) 
    iat = point; // evaluate at x[iat], substitute 
    for(jterm=0; jterm<npoints; jterm++) // get each term 
    {
      ipower = iat;
      isum = numer[jterm][0];
      for(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(jterm=0; jterm<npoints; jterm++) // clean up each term 
    {
      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(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(jterm=0; jterm<npoints; jterm++) // check 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 
    bb[0] = b;
  } // end deriv 

  double pow(double x, int pwr)
  {
    int i;
    double val = 1.0;

    for(i=0; i<pwr; i++) val = val * x;
    return val;
  } // end pow


  void check(int order, int npoints, int point, int a[], int b)
  {
    int pwr; // simple test using x^pwr 
    double sum, x, term, err, xpoint, deriv;
    double h = 0.125;
    int i, degree;
    DecimalFormat f1=new DecimalFormat("0.000E00");
  
    for(degree=npoints-2; degree<npoints+2; degree++)
    {
      pwr = degree;
      sum = 0.0;
      x = 1.0;
      for(i=0; i<npoints; i++)
      {
        term = (double)a[i]*pow(x, pwr);
        sum = sum + term;
        x = x + h;
      }
      sum = sum / (double)b;
      sum = sum / pow(h, order);
      xpoint = 1.0+(double)point*h;
      deriv = 1.0;
      for(i=0; i<order; i++)
      {  
        deriv = deriv * pwr;
        pwr = pwr -1;
      }
      deriv = deriv * pow(xpoint, pwr);
      err = deriv - sum;
      System.out.println("                err="+f1.format(err)+", sum="+
                         f1.format(sum)+", deriv="+f1.format(deriv)+
                         ", x="+xpoint+", pwr="+degree); 
    }
  } // end check

  public static void main (String[] args)
  {
    new Nderiv();
  } // end main
} // end class Nderiv