/* mpf_nuderiv.c  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:
 *           (Remember the y values are unknown at this time.
 *            The y values are the unknowns in the PDE.
 *            We are finding the c coefficients in order to
 *            build a system of simultaneous equations.)
 *
 *                              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
 *
 */

#include "mpf_inverse.h"
/* must be before "digits" is defined */
#include "mpf_nuderiv.h"
#include <stdlib.h>

void mpf_nuderiv(int order, int npoint, int point, mpf_t x[], mpf_t c[])
{
  mpf_t *A;
  mpf_t *AI;
  mpf_t pwr, tmp;
  mpf_t *fct;
  int i, j, k, n;

  mpf_set_default_prec(digits*3.32);

  mpf_init(tmp);
  mpf_init(pwr); /* others are init when first used */

  n = npoint;
  A = (mpf_t *)malloc(n*n*sizeof(mpf_t));
  AI = (mpf_t *)malloc(n*n*sizeof(mpf_t));
  fct = (mpf_t *)malloc(n*sizeof(mpf_t));
  for(i=0; i<n; i++)
    for(j=0; j<n; j++)
      mpf_init(AI[j*n+i]);

  for(i=0; i<n; i++) /* build matrix that will be inverted */
  {
    mpf_set_si(pwr, 1);
    for(j=0; j<n; j++)
    {
      mpf_init_set(A[i*n+j], pwr);
      mpf_mul(pwr, pwr, x[i]);
    }
  }
  mpf_inverse(n, digits, A, AI); /* invert the matrix */
  /* form derivative for order and point */
  for(i=0; i<n; i++) /* compute factors for order */
  {
    mpf_init_set_si(fct[i], 1);
  }
  for(j=0; j<order; j++)
  {
    for(i=0; i<n; i++)
    {
      mpf_set_si(tmp, i-j);
      mpf_mul(fct[i], fct[i], tmp);
    }
  }
  for(i=0; i<n; i++)
  {
    mpf_set_si(c[i], 0);
    mpf_set_si(pwr, 1);
    for(j=order; j<n; j++)
    { 
      mpf_mul(tmp, AI[j*n+i], pwr);
      mpf_mul(tmp, tmp, fct[j]);
      mpf_add(c[i], c[i], tmp);
      mpf_mul(pwr, pwr, x[point]);
    }
  }
  free(A);
  free(AI);
  free(fct);
} /* end mpf_nuderiv */