/* nuderiv2dg.c  non uniformly spaced derivative coefficients */
/*    given  x0,y0, x1,y1, x2,y2, ... xn,yn  non uniformly spaced
 *    find the coefficients  c???  of u0, u1, u2, ... un
 *    in order to compute the derivatives:
 *      Ux(x0,y0)  =  cx00*u0 +  cx01*u1 +  cx02*u2 + ...  cx0n*un
 *      Ux(x1,y1)  =  cx10*u0 +  cx11*u1 +  cx12*u2 + ...  cx1n*un
 *      Uy(xk,yk)  =  cyk0*u0 +  cyk1*u1 +  cyk2*u2 + ...  cykn*un
 *      Uxx(x3,y3) = cxx30*u0 + cxx31*u1 + cxx32*u2 + ... cxx3n*un
 *      ...
 *
 *           with u0, u1, u2, ...  un variables:
 *           (Remember the  values are unknown at this time.
 *            The u values are the unknowns in the PDE.
 *            We are finding the c coefficients in order to
 *            build a system of simultaneous equations, that
 *            when solved, provide the solution to the PDE.)
 *
 *  Method:   fit data   (the coefficients are lined up in four columns)
 *
 *  U(x,y) = c0           +  c1*x         +  c2*y          +  c3*x^2    +
 *           c4*xy        +  c5*y^2       +  c6*x^3        +  c7*x^2y   +
 *           c8*xy^2      +  c9*y^3       + c10*x^4        + c11*x^3y   +
 *          c12*x^2y^2    + c13*xy^3      + c14*y^4        + c15*x^5    +
 *          c16*x^4y      + c17*x^3y^2    + c18*x^2y^3     + c19*xy^4   +
 *          c20*y^5       + c21*x^6       + c22*x^5y       + c23*x^4y^2 +
 *          c24*x^3y^3    + c25*x^2y^4    + c26*x^y^5      + c27*y^6    +
 *          ...
 *
 *  then, taking the derivative of U(x,y) with respect to x :
 *  Ux(x,y) =                c1                          +  2*c3*x       +
 *           c4*y                       +  3*c6*x^2      +  2*c7*x*y     +
 *           c8*y^2                     + 4*c10*x^3      + 3*c11*x^2*y   +
 *        2*c12*x*y^2   +   c13*y^3                      + 5*c15*x^4     +
 *        4*c16*x^3*y   + 3*c17*x^2*y^2 + 2*c18*x*y^3    +   c19*y^4     +
 *                        6*c21*x^5     + 5*c22*x^4*y    + 4*c23*x^3*y^2 +
 *        3*c24*x^2*y^3 + 2*c25*x*y^4   +   c26*y^5      +
 *        ...
 *  
 *  then, taking the derivative of U(x,y) with respect to y:
 *  Uy(x,y)  =                               c2                          +
 *          c4*x       +  2*c5*y                         +    c7*x^2     +
 *        2*c8*x*y     +  3*c9*y^2                       +   c11*x^3     +
 *       2*c12*x^2*y   + 3*c13*x*y^2    + 4*c14*y^3                      +
 *         c16*x^4     + 2*c17*x^3*y    + 3*c18*x^2*y^2  + 4*c19*x*y^3   +
 *       5*c20*y^4                      +   c22*x^5      + 2*c23*x^4*y   +
 *       3*c24*x^3*y^2  +4*c25*x^2*y^3  + 5*c26*x*y^4    + 6*c27*y^5     +
 *       ...
 *
 *  ...
 *
 *  then, taking the forth derivative of U(x,y) with respect to x:
 *  Uxxxx(x,y) =  ...
 *           ...
 *
 *  Algorithm:
 *                                               u0 u1 u2 u3 u4 u5 ... 
 *  form:  | 1  x0  y0  x0^2  x0*y0 y0^2 x0^3 ... 1  0  0  0  0  0 ... | = c0
 *         | 1  x1  y1  x1^2  x1*y1 y1^2 x1^3 ... 0  1  0  0  0  0 ... | = c1
 *         | 1  x2  y2  x2^2  x2*y2 y2^2 x2^3 ... 0  0  1  0  0  0 ... | = c2
 *         | 1  x3  y3  x3^2  x3*y3 y3^2 x3^3 ... 0  0  0  1  0  0 ... | = c3
 *         | 1  x4  y4  x4^2  x4*y4 y4^2 x4^3 ... 0  0  0  0  1  0 ... | = c4
 *         | 1  x5  y5  x5^2  x5*y5 y5^2 x5^3 ... 0  0  0  0  0  1 ... | = c5
 *         | 1  x6  y6  x6^2  x6*y6 y6^2 x5^3 ... 0  0  0  0  0  0 ... | = c6
 *         |    ...
 *         | 1  x_(npoints-1)
 *
 *  reduce
 *     c0 c1 c2 c3 c4 c5 c6 ...
 *   | 1  0  0  0  0  0  0 ... c00  c01  c02  c03  c04  c05  c06 ... | = u0
 *   | 0  1  0  0  0  0  0 ... c10  c11  c12  c13  c14  c15  c16 ... | = u1
 *   | 0  0  1  0  0  0  0 ... c20  c21  c22  c23  c24  c25  c26 ... | = u2
 *   | 0  0  0  1  0  0  0 ... c30  c31  c32  c33  c34  c35  c36 ... | = u3
 *   | 0  0  0  0  1  0  0 ... c40  c41  c42  c43  c44  c45  c46 ... | = u4
 *   | 0  0  0  0  0  1  0 ... c50  c51  c52  c53  c54  c55  c56 ... | = u5
 *   | 0  0  0  0  0  0  1 ... c60  c61  c62  c63  c64  c65  c66 ... | = u6
 *     ...
 *                            |      this is just the inverse        |
 *                            | add points that do not make singular |
 * 
 * thus, the derivative of U(x,y) with respect to x at point x0,y0 is:
 * Ux(x0,y0) =
 *       c10*u0       +   c11*u1       +   c12*u2       +   c13*u3       + ...
 *     2*c30*u0*x0    + 2*c31*u1*x0    + 2*c32*u2*x0    + 2*c33*u3*x0    + ...
 *       c40*u0*y0    +   c41*u1*y0    +   c42*u2*y0    +   c43*u3*y0    + ...
 *     3*c60*u0*x0^2  + 3*c61*u1*x0^2  + 3*c62*u2*x0^2  + 3*c63*u3*x0^2  + ...
 *     2*c70*u0*x0*y0 + 2*c71*u1*x0*y0 + 2*c72*u2*x0*y0 + 2*c73*u2*x0*y0 + ...
 *           ...
 *
 * rewriting column coefficients as values to be computed:
 *       cx00 = c10 + 2*c30*x0 + c40*y0 + 3*c60*x0^2 + 2*c70*u0*x0*y0  + ...
 *       cx01 = c11 + 2*c31*x0 + c41*y0 + 3*c61*x0^2 + 2*c71*u0*x0*y0  + ...
 *       cx02 = c12 + 2*c32*x0 + c42*y0 + 3*c62*x0^2 + 2*c72*u0*x0*y0  + ...
 *       cx03 = c13 + 2*c33*x0 + c43*y0 + 3*c63*x0^2 + 2*c73*u0*x0*y0  + ...
 *       ...
 *
 *
 * To compute the first derivative with respect to x at point x0,y0,
 * given cx00, cx01, cx02, cx03, ... and knowing u0 at x0,y0, u1 ...
 *
 * Ux(x0,y0) = cx00*u0 + cx01*u1 + cx02*u2 +cx03*u3 ...
 *             hint: the cx00, cx01, cx03, .. cx0n are the values returned 
 *                   for point zero
 *
 *
 * For point one, x1,y1 compute coefficients:
 *
 *       cx10 = c10 + 2*c30*x1 + c40*y1 + 3*c60*x1^2 + 2*c70*u0*x1*y1  + ...
 *       cx11 = c11 + 2*c31*x1 + c41*y1 + 3*c61*x1^2 + 2*c71*u0*x1*y1  + ...
 *       cx12 = c12 + 2*c32*x1 + c42*y1 + 3*c62*x1^2 + 2*c72*u0*x1*y1  + ...
 *       cx13 = c13 + 2*c33*x1 + c43*y1 + 3*c63*x1^2 + 2*c73*u0*x1*y1  + ...
 *       ...
 *
 * Ux(x1,y1) = cx10*u0 + cx11*u1 + cx12*u2 +cx13*u3 + ...
 *                hint: the cx10, cx11, cx13, .. cx1n are the values returned 
 *                      for point one
 *
 * 
 * the first derivative of U(x,y) with respect to y at point x1,y1 is:
 * Uy(x1,y1) = c20*u0        + c21*u1      +   c22*u2      +   c23*u3      ...
 *             c40*u0*x1     + c41*u1*x1   +   c42*u2*x1   +   c43*u3*x1   ...
 *           2*c50*u0*y1   + 2*c51*u1*y1   + 2*c52*u2*y1   + 2*c53*u3*y1   ...
 *             c70*u0*x1^2 +   c71*u1*x1^2 +   c72*u2*x1^2 +   c73*u3*x1^2 ...
 *
 * rewriting column coefficients as values to be computed:
 *
 *       cy10 = c20 + c40*x1 + 2*c50*y1 + c70*x1^2 + ...
 *       cy11 = c21 + c41*x1 + 2*c51*y1 + c71*x1^2 + ...
 *       cy12 = c22 + c42*x1 + 2*c52*y1 + c72*x1^2 + ...
 *       cy13 = c23 + c43*x1 + 2*c53*y1 + c73*x1^2 + ...
 *
 * To compute the first derivative with respect to y at point x1,y1,
 * given cy00, cy01, cy02, cy03, ... and knowing u0 at x0,y0, u1 ...
 *
 * Uy(x1,y1) = cy10*u0 + cy11*u1 + cy12*u2 + cy13*u3 ...
 *
 *
 * ... Uxy, Uyy, Uxxx, Uxxy, Uxyy, Uyyy, Uxxxx, Uxxxy, Uxxyy, Uxyyy, Uyyyy
 *     all partial derivatives up to fourth
 *
 * Now, we really get clever and use xpow and ypow to compute
 * cx, cy, cxx, cxy, cyy, ... terms
 *
 * Then we are very careful to find the set of points that do
 * not create a singular matrix, and return indices of the useful points
 * Always make point the first.
 */

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

static int debug = 0;
static int stat = 0;
static double *P;  /* point matrix */
static double *AI; /* inverse */
static int ordpt[] = {0, 3, 6, 10, 15, 21, 28, 36};
static int n;      /* matrix size is n by n */
static int norder; /* actual maximum order */
static int maxnpwr = 36;
                     /* max of 36 points used now, thus max of 7th order */
static int xpow[] = {0, 1, 0, 2, 1, 0, 3, 2, 1, 0, 4, 3, 2, 1, 0,
                     5, 4, 3, 2, 1, 0, 6, 5, 4, 3, 2, 1, 0,
                     7, 6, 5, 4, 3, 2, 1, 0};
static int ypow[] = {0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4,
                     0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 6,
                     0, 1, 2, 3, 4, 5, 6, 7};
static double xpwr[8], ypwr[8];

static void build(int order, int npoint, int point,
                  double x[], double y[], int ip[]);

static int test_build(int ni, double xa[], double ya[]);

static int inverse(int n, double A[], double AA[]);

/* the number of x,y  npoints  must be at least sum i=1,order+1 of i
 * here, order means maximum order of terms used in computing cx, cy, ...
 * order 1  3 npoints   only dx, dy
 * order 2  6 npoints   only dxx, dxy, dyy and lower orders
 * order 3 10 npoints   only dxxx, dxxy, dxyy, dyyy and lower orders
 * order 4 15 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy and lower orders
 * order 5 21 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy and lower orders
 * order 6 28 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy and lower orders
 * order 7 36 npoints   only dxxxx, dxxxy, dxxyy, dxyyy, dyyyy and lower orders
 * point is 0..npoint-1 the point about which the derivative is computed
 */

/* stat is actual number of entries in c[] and ip[] */
/* xpow[] and ypow[] are correct for the entries in c[] */

int nu2dx(int order, int npoint, int point, double x[],
          double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(xpow[i]<1) continue;
      ci = (double)xpow[i];
      xp = xpow[i]-1; /* taking first derivative */
      yp = ypow[i];
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dx */

int nu2dy(int order, int npoint, int point,
          double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(ypow[i]<1) continue;
      ci = (double)ypow[i];
      xp = xpow[i];
      yp = ypow[i]-1;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dy */

int nu2dxx(int order, int npoint, int point,
           double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cxx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(xpow[i]<2) continue;
      ci = (double)xpow[i]*(double)(xpow[i]-1);
      xp = xpow[i]-2; /* taking second derivative */
      yp = ypow[i];
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxx */

int nu2dxy(int order, int npoint, int point,
           double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cxy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(xpow[i]<1 || ypow[i]<1) continue;
      ci = (double)xpow[i]*(double)(ypow[i]);
      xp = xpow[i]-1;
      yp = ypow[i]-1;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxy */

int nu2dyy(int order, int npoint, int point,
           double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(ypow[i]<2) continue;
      ci = (double)ypow[i]*(double)(ypow[i]-1);
      xp = xpow[i];
      yp = ypow[i]-2;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dyy */

int nu2dxxx(int order, int npoint, int point,
            double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cxxx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(xpow[i]<3) continue;
      ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(xpow[i]-2);
      xp = xpow[i]-3;
      yp = ypow[i];
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxxx */

int nu2dxxy(int order, int npoint, int point,
            double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<n; j++) /* for cxxx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<maxnpwr; i++)
    {
      if(i>=n) continue;
      if(xpow[i]<2 || ypow[i]<1) continue;
      ci = (double)xpow[i]*(double)(xpow[i]-1)*(double)(ypow[i]);
      xp = xpow[i]-2;
      yp = ypow[i]-1;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxxy */

int nu2dxyy(int order, int npoint, int point,
            double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<n; j++) /* for cxxx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<maxnpwr; i++)
    {
      if(i>=n) continue;
      if(xpow[i]<1 || ypow[i]<2) continue;
      ci = (double)xpow[i]*(double)(ypow[i])*(double)(ypow[i]-1);
      xp = xpow[i]-1;
      yp = ypow[i]-2;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxyy */

int nu2dyyy(int order, int npoint, int point,
            double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cyyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(ypow[i]<3) continue;
      ci = (double)ypow[i]*(double)(ypow[i]-1)*(double)(ypow[i]-2);
      xp = xpow[i];
      yp = ypow[i]-3;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dyyy */

int nu2dxxxx(int order, int npoint, int point,
             double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cxxxx-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(xpow[i]<4) continue;
      ci = (double)xpow[i]*(double)(xpow[i]-1)*
           (double)(xpow[i]-2)*(double)(xpow[i]-3);
      xp = xpow[i]-4;
      yp = ypow[i];
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxxxx */

int nu2dxxxy(int order, int npoint, int point,
             double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cxxxy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(xpow[i]<3 || ypow[i]<1) continue;
      ci = (double)xpow[i]*(double)(xpow[i]-1)*
           (double)(xpow[i]-2)*(double)(ypow[i]);
      xp = xpow[i]-3;
      yp = ypow[i]-1;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxxxy */

int nu2dxxyy(int order, int npoint, int point,
             double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cxxyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(xpow[i]<2 || ypow[i]<2) continue;
      ci = (double)xpow[i]*(double)(xpow[i]-1)*
           (double)(ypow[i])*(double)(ypow[i]-1);
      xp = xpow[i]-2;
      yp = ypow[i]-2;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxxyy */

int nu2dxyyy(int order, int npoint, int point,
             double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cxyyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(xpow[i]<1 || ypow[i]<3) continue;
      ci = (double)xpow[i]*(double)(ypow[i])*
           (double)(ypow[i]-1)*(double)(ypow[i]-2);
      xp = xpow[i]-1;
      yp = ypow[i]-3;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dxyyy */

int nu2dyyyy(int order, int npoint, int point,
             double x[], double y[], double c[], int ip[])
{
  int i, j;
  double ci;
  int xp, yp;

  build(order, npoint, point, x, y, ip);
  for(j=0; j<stat; j++) /* for cyyyy-point[u] */
  {
    c[j] = 0.0;
    for(i=0; i<stat; i++)
    {
      if(ypow[i]<4) continue;
      ci = (double)ypow[i]*(double)(ypow[i]-1)*
	   (double)(ypow[i]-2)*(double)(ypow[i]-3);
      xp = xpow[i];
      yp = ypow[i]-4;
      c[j] += ci*AI[i*n+j]*xpwr[xp]*ypwr[yp];
    } 
  }
  return stat;
} /* end nu2dyyyy */


static void build(int order, int npoint, int point,
                  double x[], double y[], int ip[])
{
  int i, j, k, sing;
  int n2i, nd, ni, nj;
  double xa[36], ya[36], tmp;

  /* checks  */
  if(point<0 || point>=npoint)
     printf("ERROR nuderiv2dg point outside range\n"); 
  n = npoint; /* global */
  norder = order; /* global */
  P = (double *)calloc(n*n, sizeof(double));
  if(P==NULL) printf("**** nuderiv2dg ERROR P not allocated \n");
  AI = (double *)calloc(n*n, sizeof(double));
  if(AI==NULL) printf("**** nuderiv2dg ERROR AI not allocated \n");

  n2i = 1; /* good points so far, this will become  n  */
  nd  = 0; /* number deleted */
  for(i=0; i<npoint; i++)
  {
    ip[i] = i;
    xa[i] = x[i];
    ya[i] = y[i];
  }
  if(point != 0)
  {
    ip[0] = point; /* make 'point' be first point, must be included */
    ip[point] = 0;
    tmp = xa[0];
    xa[0] = xa[point];
    xa[point] = tmp;
    tmp = ya[0];
    ya[0] = ya[point];
    ya[point] = tmp;     
  }
  ni = 2; /* trying  ni  points for non singular */
  while(ni<=npoint)
  {
    sing = test_build(ni, xa, ya); /* P and AI global */
    if(sing==0)
    {
      n2i = ni;
      if(debug) printf("test_build added point %d for index %d, nd=%d  \n",
	               ip[ni-1], ni-1, nd);
      if(ni>=36) break;
      ni++;
      continue;
    }
    nd++;
    if(debug)
      printf("test_build deleted point %d for index %d, nd=%d, ni+nd=%d \n",
	     ip[ni-1], ni-1, nd, ni+nd+1);
    if(ni+nd>=npoint) break;
    for(nj=ni-1; nj<npoint-1; nj++)
    {
      ip[nj] = ip[nj+1];
      xa[nj]  = xa[nj+1];
      ya[nj]  = ya[nj+1];
    }
  }
  stat = n2i; /* global */
  n = n2i;

  if(debug)
  {
    printf("\n\nactual points selected\n");
    for(i=0; i<stat; i++) printf("i=%2d, ip[i]=%2d, xa[i]=%f,  ya[i]=%f\n",
                                 i, ip[i], xa[i], ya[i]);
  }
  /* build xpwr, ypwr at point */
  xpwr[0] = 1.0;
  ypwr[0] = 1.0;
  for(k=1; k<8; k++)
  {
    xpwr[k] = xpwr[k-1]*xa[0]; /* 'point' moved to zero */
    ypwr[k] = ypwr[k-1]*ya[0];
  }
  if(n<ordpt[order])
    printf("probable loss of accuracy, order=%d, needs n=%d, has only %d \n",
           order, ordpt[order], n);
} /* end build */


static int test_build(int ni, double xa[], double ya[])
                                                     /* P and AI are global */
{
  int i, j, k, sing;

  for(i=0; i<ni; i++) /* build matrix that will be inverted */
  {
    xpwr[0] = 1.0;
    ypwr[0] = 1.0;
    for(k=1; k<8; k++)
    {
      xpwr[k] = xpwr[k-1]*xa[i];
      ypwr[k] = ypwr[k-1]*ya[i];
    }
    for(j=0; j<ni; j++)
    {
      P[i*ni+j] = xpwr[xpow[j]]*ypwr[ypow[j]];
    }
  }

  sing = inverse(ni, P, AI); /* invert the matrix */

  return sing;
} /* end test_build */


static int inverse(int n, double A[], double AA[])
{
  int *ROW, *COL;
  double *TEMP;
  int HOLD , I_pivot , J_pivot;
  double pivot, abs_pivot;
  int i, j, k;

  ROW = (int *)calloc(n, sizeof(int));
  COL = (int *)calloc(n, sizeof(int));
  TEMP = (double *)calloc(n, sizeof(double));

  /*              set up row and column interchange vectors */
  for (k=0; k<n; k++)
  {
    ROW[k] = k ;
    COL[k] = k ;
  }

  /*               begin main reduction loop */
  for (k=0; k<n; k++)
  {
    /*             find largest element for pivot */
    pivot = A[ROW[k]*n+COL[k]] ;
    I_pivot = k;
    J_pivot = k;
    for(i=k; i<n; i++)
    {
      for(j=k; j<n; j++)
      {
        abs_pivot = abs(pivot) ;
        if( abs(A[ROW[i]*n+COL[j]]) > abs_pivot )
        {
          I_pivot = i ;
          J_pivot = j ;
          pivot = A[ROW[i]*n+COL[j]] ;
        }
      }
    }
    if(abs(pivot) < 1.0E-12)
    {
      free(ROW);
      free(COL);
      free(TEMP);
      return 1;
    }

    HOLD = ROW[k];
    ROW[k]= ROW[I_pivot];
    ROW[I_pivot] = HOLD ;
    HOLD = COL[k];
    COL[k]= COL[J_pivot];
    COL[J_pivot] = HOLD ;

    /*               reduce about pivot */
    A[ROW[k]*n+COL[k]] = 1.0 / pivot ;
    for (j=0; j<n; j++)
    {
      if( j != k )
      {
        A[ROW[k]*n+COL[j]] = A[ROW[k]*n+COL[j]] * A[ROW[k]*n+COL[k]];
      }
    }

    /*               inner reduction loop */
    for (i=0; i<n; i++)
    {
      if( k != i )
      {
        for (j=0; j<n; j++)
        {
          if( k != j )
          {
            A[ROW[i]*n+COL[j]] = A[ROW[i]*n+COL[j]] - A[ROW[i]*n+COL[k]] *
                                 A[ROW[k]*n+COL[j]] ;
          }
        }
        A[ROW[i]*n+COL [k]] = - A[ROW[i]*n+COL[k]] * A[ROW[k]*n+COL[k]] ;
      }
    }
  }
  /*      end of main reduction loop */

  /*      unscramble rows */
  for (j=0; j<n; j++) {
    for (i=0; i<n; i++) {
      TEMP[COL[i]] = A[ROW[i]*n+j];
    }
    for (i=0; i<n; i++) {
      A[i*n+j] = TEMP[i] ;
    }
  }

  /*      unscramble columns */
  for (i=0; i<n; i++)
  {
    for (j=0; j<n; j++)
    {
      TEMP[ROW[j]] = A[i*n+COL[j]] ;
    }
    for (j=0; j<n; j++)
    {
      AA[i*n+j] = TEMP[j]; /* final output */
    }
  }
  free(ROW);
  free(COL);
  free(TEMP);
  return 0;
} /* end inverse */

/* end nuderiv2dg.c  with a few callable functions */