// 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   

public class nderiv  // for function  deriv
{
  public nderiv()
  {
      // no initialization needed
  }

  public static 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 


  public static 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, term, b;
    int k = 0;
    int jorder, ipower, isum, iat, jterm, r;
    int debug = 0;
  
    if(npoints<=order)
    {
      System.out.println("ERROR in call to deriv, npoints="+npoints+
                         " < order="+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 coefficients 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 
    {
      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(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++) // 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 
    bb[0] = b;
  } // end deriv 


  public static double pown(double x, int pwr)
  {
    int i;
    double val = 1.0;

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


  public static void check(int order, int npoints, int point, int a[], int b)
  {
    int pwr; // simple test using x^pwr 
    double sum, x, h, term, err, xpoint, deriv;
    double hh[]={0.125, 0.0625, 0.03125, 0.015625};
    int i, k, m, degree;
  
    m = 2;
    for(degree=npoints-2; degree<npoints+2; degree++)
    {
      if(order>3 && degree==npoints+1) m=3; // extra check on error 
      for(k=0; k<m; k++)
      {
        if(k>0 && degree<npoints) break;
        h = hh[k];
        if(order>3) h=h/2.0; // cumulative, for reasonable error 
        if(order>4) h=h/2.0;
        if(order>5) h=h/2.0;
        pwr = degree;
        sum = 0.0;
        x = 1.0; // let x = 1.0, compute derivative using this x 
        for(i=0; i<npoints; i++)
        {
          term = (double)a[i]*pown(x, pwr); // x^pwr + x^pwr-1 + ... 1 
          sum = sum + term;
          x = x + h;
        }
        sum = sum / (double)b;
        sum = sum / pown(h, order); // approximate derivative at 'point' 
        xpoint = 1.0+(double)point*h; // xpoint is x offset by 'point' 
        deriv = 1.0;
        for(i=0; i<order; i++)
        {  
          deriv = deriv * pwr; // analytic derivative 
          pwr = pwr -1;
        }
        deriv = deriv * pown(xpoint, pwr);
        err = deriv - sum;
        System.out.println("                err="+err+", h="+h+", deriv="+
                           deriv+", x="+xpoint+", pwr="+degree);
      } 
    }
  } // end check 

  public static void main (String[] args) // test deriv
  {
    nderiv N = new nderiv();

    int norder = 6;
    int npoints = 12;
    int order, points, term, jterm;
    int a[]=new int[50];         // coefficients of f(x0), f(x1), ... 
    int bb[]=new int[1];         // dumb way to get out parameter
    int b;
  
    System.out.println("Formulas for numerical computation of derivatives  ");
    System.out.println("  f^(3) means third derivative  ");
    System.out.println("  f^(2)(x[1]) means second derivative computed at point x[1]  ");
    System.out.println("  f(x[2]) means function evaluated at point x[2]  ");
    System.out.println("  Points are x[0], x[1]=x[0]+h, x[2]=x0+2h, ...  ");
    System.out.println("  The error term gives the power of h and  ");
    System.out.println("  derivative order with z chosen for maximum value. ");
    System.out.println("  The error terms, err=, are for test polynomial sum(x^pwr)  ");
    System.out.println("  with h = 0.125 or smaller, and order = 'pwr'   ");
  
    for(order=1; order<=norder; order++)
    {
      npoints = npoints -2; // fewer for larger derivatives 
      if(npoints<=order+1) npoints = order + 2;
      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, bb);
          b = bb[0];
        
          for(jterm=0; jterm<points; jterm++)
          {
            if(jterm==0)
	      System.out.println("f^("+order+")(x["+term+"])=(1/"+b+"h^"+
                                 order+")( "+a[jterm]+" * f(x["+jterm+"]) +");
            else if(jterm==points-1)
	      System.out.println("                      "+a[jterm]+" * f(x["+
                                 jterm+"] ) + O(h^"+(points-order)+")f^("+
                                 points+")(z)");
            else
	      System.out.println("                      "+a[jterm]+
                                 " * f(x["+jterm+"] +");
                     
          } // end jterm loop 
          check(order, points, term, a, b);
          System.out.println(" ");
        } // end term loop 
        System.out.println(" ");
      }  // end points loop 
      System.out.println(" ");
    } // end order loop 
    System.out.println("many formulas checked against Abramowitz and Stegun, ");
    System.out.println("Handbook of Mathematical Functions Table 25.2");
  } // end main of nderiv
} // end class nderiv.java