/* nderiv.c  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 computer                 */
/*           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   */

#include <stdio.h>
#undef abs
#define abs(x) ((x)<0)?(-(x)):(x)

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

  if(a==0 || b==0) return 1;
  if(abs(a)>abs(b)) { a1 = abs(a); b1 = abs(b); }
  else              { a1 = abs(b); b1 = 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[100];          /* coefficients of h */
  int x[100];          /* coefficients of x, variable for differentiating */
  int numer[100][100]; /* numerator of a term numer[term][pos] */
  int denom[100];      /* denominator of coefficient */
  int i, j, k, term, b;
  int jorder, ipower, isum, iat, jterm, r;
  int debug = 0;
  
  if(npoints<=order)
  {
    printf("ERROR in call to deriv, npoints=%d < order=%d+1=%d \n\n",
           npoints, order);
    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++)
    {
      printf("denom[%d]=%d \n", i, denom[i]);
      for(j=0; j<npoints-1; j++)
        printf("numer[%d][%d]=%d \n", 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)printf("p(x)= %d x^%d at term=%d \n",
                        numer[term][j], j, 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++)
        printf("p^(%d)(x)= %d x^%d at term=%d \n",
               order, numer[term][i], i, 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)
      printf("f^(%d)(x[%d]) = (1/h^%d) (%d/%d f(x[%d]) + \n",
             order, iat, order, a[jterm], denom[jterm], jterm);
  }
  if(debug>0) printf("\n");
  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)
      printf("f^(%d)(x[%d])=(1/%dh^%d)(%d f(x[%d]) + \n",
             order, iat, b, order, a[jterm], jterm);
  } /* end computing terms of coefficients */
  *bb = b;
} /* end deriv */


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


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[4]={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;
      printf("                err=%g, h=%g, deriv=%g, x=%g, pwr=%d \n",
             err, h, deriv, xpoint, degree);
    } 
  }
} /* end check */


int main(int argc, char * argv[])
{
  int norder = 6;
  int npoints = 12;
  int order, points, term, jterm;
  int a[50];         /* coefficients of f(x0), f(x1), ... */
  int b;
  
  printf("Formulas for numerical computation of derivatives \n");
  printf("  f^(3) means third derivative \n");
  printf("  f^(2)(x[1]) means second derivative computed at point x[1] \n");
  printf("  f(x[2]) means function evaluated at point x[2] \n");
  printf("  Points are x[0], x[1]=x[0]+h, x[2]=x0+2h, ... \n");
  printf("  The error term gives the power of h and \n");
  printf("  derivative order with z chosen for maximum value.\n");
  printf("  The error terms, err=, are for test polynomial sum(x^pwr) \n");
  printf("  with h = 0.125 or smaller, and order = 'pwr' \n\n");
  
  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++)
      {
        printf("computing order=%d, npoints=%d, at term=%d \n",
                order, points, term);
        deriv(order, points, term, a, &b);
        
        for(jterm=0; jterm<points; jterm++)
        {
          if(jterm==0)
            printf("f^(%d)(x[%d])=(1/%dh^%d)(%6d * f(x[%d]) + \n",
                     order, term, b, order, a[jterm], jterm);
          else if(jterm==points-1)
            printf("                      %6d * f(x[%d] ) + O(h^%d)f^(%d)(z) \n",
                    a[jterm], jterm, points-order, points);
          else
            printf("                      %6d * f(x[%d]) + \n",
                    a[jterm], jterm);
        } /* end jterm loop */
        check(order, points, term, a, b);
        printf("\n");
      } /* end term loop */
      printf("\n");
    }  /* end points loop */
    printf("\n");
  } /* end order loop */
  printf("many formulas checked against Abramowitz and Stegun,\n");
  printf("Handbook of Mathematical Functions Table 25.2\n");
  return 0;
} /* end main of nderiv.c */