// lsfit.c  for 1, 2, 3, 4, 5, 6, 7 dimensions
//              6, 6, 6, 6, 6, 6, 4 th power
//          Compute a multiple dimensional polynomial
//          that is a least square fit of some data.
//          Evaluate the polynomial at desired points.
// use:   fit_init(n,m);  n is dimension, number of independent variables
//                        m is max power of polynomial
//                       
//        fit2_bpoly(nx, ny, xg, yg, u) build coefficients
//                   nx number of x coordinates i=0,nx-1
//                   ny number of y coordinates j=0,ny-1
//                   xg array of x coordinate values
//                   yg array of y coordinate values
//                   u[i*ny+j] value at point i,j
// 3D example
//    fit_init(3,6); // 3D to 6th power u = poly(x,y,z)
//    last = false;
//    for(i=0; i<n; i++){ // must provide enough independent data
//      if(i==n-1) last = true;
//      fit3_data(x, y, z, u, last); } // now have fit
//    ufit = eval3_poly(x, y, z);      // get intermediate vals 
//    write_fit3(fname.c);

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "lsfit.h"

#undef  abs
#define abs(x) (((x)<0.0)?(-(x)):(x))
#undef  max
#define max(x,y) (((x)>(y))?((x)):(y))
#undef  min
#define min(x,y) (((x)<(y))?((x)):(y))
#include "lsfit.h"

static int debug = 0;
static int nnn;   /* number of coefficients */
static int n, m;  /* max power, number of variables */
static double *A; /* matrix of simultaneous equations*/
static double *C; /* polynomial coefficients */
static double *Y; // RHS of equations
static int *Apwr; // powers of each variable in a term
static int *Bpwr;
static int *Cpwr;
static int *Dpwr;
static int *Epwr;
static int *Fpwr;
static int *Gpwr;
static double *Av; // values of powers
static double *Bv;
static double *Cv;
static double *Dv;
static double *Ev;
static double *Fv;
static double *Gv;

static void simeq(int nn, double X[]); /* A and Y global */
static void gen_1d_powers();
static void gen_2d_powers();
static void gen_3d_powers();
static void gen_4d_powers();
static void gen_5d_powers();
static void gen_6d_powers();
static void gen_7d_powers();
static void check_2d(int nx, int ny, double x[], double y[], double u[], 
                     double *rms_err, double *avg_err, double *max_err);

void fit_init(int an, int am)  /* dimension, power, only one fit at a time */
{
  int i, j, k;

  n = an;
  m = am;
  nnn = 0; /* number of terms to be computed*/
  if(n>7)
  {
    printf("this lsfit.c can handle only handle to dimension 7 , not %d\n",n);
    exit(1);
  }
  if(m>6 || (n==7 && m>4))
  {
    printf("this lsfit.c can handle only handle 6 power, 4 at 7D , not %d\n",
           m);
    exit(1);
  }

  if(n==1)
  {
    gen_1d_powers();
  }
  else if(n==2)
  {
    gen_2d_powers();
  }
  else if(n==3)
  {
    gen_3d_powers();
  }
  else if(n==4)
  {
    gen_4d_powers();
  }
  else if(n==5)
  {
    gen_5d_powers();
  }
  else if(n==6)
  {
    gen_6d_powers();
  }
  else if(n==7)
  {
    gen_7d_powers();
  }
  else
  {
   printf("this lsfit.c can handle only 1,2,3,4,5,6,7 dimensions, not %d\n",n);
   exit(1);
  }
  if(debug) printf("nnn= %d polynomial terms \n", nnn);
  if(debug) fflush(stdin);
  A = (double *)malloc(nnn*nnn*sizeof(double));
  C = (double *)malloc(nnn*sizeof(double));
  Y = (double *)malloc(nnn*sizeof(double));
  for(i=0; i<nnn; i++)
  {
    for(j=0; j<nnn; j++) A[i*nnn+j] = 0.0;
    C[i] = 0.0;
    Y[i] = 0.0;
  }
  Av = (double *)malloc((m+1)*sizeof(double));
  Bv = (double *)malloc((m+1)*sizeof(double));
  Cv = (double *)malloc((m+1)*sizeof(double));
  Dv = (double *)malloc((m+1)*sizeof(double));
  Ev = (double *)malloc((m+1)*sizeof(double));
  Fv = (double *)malloc((m+1)*sizeof(double));
  Gv = (double *)malloc((m+1)*sizeof(double));
  if(debug) printf("nnn= %d space allocated \n", nnn);
} /* end fit_init */

void fit_free()  /* free allocated ram */
{
  free(A);
  free(C);
  free(Y);
  free(Av);
  free(Bv);
  free(Cv);
  free(Dv);
  free(Ev);
  free(Fv);
  free(Gv);
  free(Apwr);
  if(n>1) free(Bpwr);
  if(n>2) free(Cpwr);
  if(n>3) free(Dpwr);
  if(n>4) free(Epwr);
  if(n>5) free(Fpwr);
  if(n>6) free(Gpwr);
} /* end fit_free */ 

static void gen_1d_powers()
{
  int i;
  int pwrsix[9]={0,1,2,3,4,5,6,7,8}; // start of each order
  int pwrs[8]={0, 1, 2, 3, 4, 5, 6, 7};

  nnn = pwrsix[m+1];
  Apwr = (int *)malloc(nnn*sizeof(int));

  if(debug) printf("terms used to find fit \n");
  for(i=0; i<=m; i++) // power requested
  {
    Apwr[i]=pwrs[i];
    if(debug) printf("%d   a^%d \n", i, Apwr[i]);
  }
  if(debug) printf(" \n");
} /* end gen_1d_powers */

static void gen_2d_powers()
{
  int i, ii;
  int pwrsix[8]={0,1,3,6,10,15,21,28}; // start of each order
  int pwrs[28][2]={{0,0},
		{1,0},{0,1},
		{2,0},{1,1},{0,2},
		{3,0},{2,1},{1,2},{0,3},
		{4,0},{3,1},{2,2},{1,3},{0,4},
		{5,0},{4,1},{3,2},{2,3},{1,4},{0,5},
		{6,0},{5,1},{4,2},{3,3},{2,4},{1,5},{0,6}};
  nnn = pwrsix[m+1];
  Apwr = (int *)malloc(nnn*sizeof(int));
  Bpwr = (int *)malloc(nnn*sizeof(int));

  if(debug) printf("terms used to find fit \n");
  for(i=0; i<=m; i++) // power requested
  {
    for(ii=pwrsix[i]; ii<pwrsix[i+1]; ii++)
    {
      Apwr[ii]=pwrs[ii][0];
      Bpwr[ii]=pwrs[ii][1];
      if(debug) printf("%d   a^%d * b^%d\n", ii, Apwr[ii], Bpwr[ii]);
    }
    if(debug) printf(" \n");
  }
} /* end gen_2d_powers */

static void gen_3d_powers()
{
  int i, ii;
  int pwrsix[8]={0,1,4,10,20,35,56,84}; // start of each order                 
  int pwrs[84][3]={{0,0,0},{1,0,0},{0,1,0},{0,0,1},
		{2,0,0},{1,1,0},{1,0,1},{0,2,0},
		{0,1,1},{0,0,2},{3,0,0},{2,1,0},
		{2,0,1},{1,2,0},{1,1,1},{1,0,2},
		{0,3,0},{0,2,1},{0,1,2},{0,0,3},
		{4,0,0},{3,1,0},{3,0,1},{2,2,0},
		{2,1,1},{2,0,2},{1,3,0},{1,2,1},
		{1,1,2},{1,0,3},{0,4,0},{0,3,1},
		{0,2,2},{0,1,3},{0,0,4},{5,0,0},
		{4,1,0},{3,2,0},{2,3,0},{1,4,0},
		{0,5,0},{4,0,1},{3,1,1},{2,2,1},
		{1,3,1},{0,4,1},{3,0,2},{2,1,2},
		{1,2,2},{0,3,2},{2,0,3},{1,1,3},
		{0,2,3},{1,0,4},{0,1,4},{0,0,5},
		{6,0,0},{5,1,0},{4,2,0},{3,3,0},
		{2,4,0},{1,5,0},{0,6,0},{5,0,1},
		{4,1,1},{3,2,1},{2,3,1},{1,4,1},
		{0,5,1},{4,0,2},{3,1,2},{2,2,2},
		{1,3,2},{0,4,2},{3,0,3},{2,1,3},
		{1,2,3},{0,3,3},{2,0,4},{1,1,4},
		{0,2,4},{1,0,5},{0,1,5},{0,0,6} };
  nnn = pwrsix[m+1];
  Apwr = (int *)malloc(nnn*sizeof(int));
  Bpwr = (int *)malloc(nnn*sizeof(int));
  Cpwr = (int *)malloc(nnn*sizeof(int));

  if(debug) printf("terms used to find fit\n");
  for(i=0; i<=m; i++) // power requested                                  
  {
    for(ii=pwrsix[i]; ii<pwrsix[i+1]; ii++)
    {
      Apwr[ii]=pwrs[ii][0];
      Bpwr[ii]=pwrs[ii][1];
      Cpwr[ii]=pwrs[ii][2];
      if(debug) printf("%d   a^%d * b^%d * c^%d \n",
                       ii, Apwr[ii], Bpwr[ii], Cpwr[ii]);
    }
    if(debug) printf(" \n");
  }
} /* end gen_3d_powers */

static void gen_4d_powers()
{
  int i, ii;
  int pwrsix[8]={0,1,5,15,35,70,126,210}; // start of each order               
  int pwrs[210][4]={{0,0,0,0},
		{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1},
		{2,0,0,0},{1,1,0,0},{0,2,0,0},{1,0,1,0},{0,1,1,0},
		{0,0,2,0},{1,0,0,1},{0,1,0,1},{0,0,1,1},{0,0,0,2},
		{3,0,0,0},{2,1,0,0},{1,2,0,0},{0,3,0,0},{2,0,1,0},
		{1,1,1,0},{0,2,1,0},{1,0,2,0},{0,1,2,0},{0,0,3,0},
		{2,0,0,1},{1,1,0,1},{0,2,0,1},{1,0,1,1},{0,1,1,1},
		{0,0,2,1},{1,0,0,2},{0,1,0,2},{0,0,1,2},{0,0,0,3},
		{4,0,0,0},{3,1,0,0},{2,2,0,0},{1,3,0,0},{0,4,0,0},
		{3,0,1,0},{2,1,1,0},{1,2,1,0},{0,3,1,0},{2,0,2,0},
		{1,1,2,0},{0,2,2,0},{1,0,3,0},{0,1,3,0},{0,0,4,0},
		{3,0,0,1},{2,1,0,1},{1,2,0,1},{0,3,0,1},{2,0,1,1},
		{1,1,1,1},{0,2,1,1},{1,0,2,1},{0,1,2,1},{0,0,3,1},
		{2,0,0,2},{1,1,0,2},{0,2,0,2},{1,0,1,2},{0,1,1,2},
		{0,0,2,2},{1,0,0,3},{0,1,0,3},{0,0,1,3},{0,0,0,4},
		{5,0,0,0},{4,1,0,0},{3,2,0,0},{2,3,0,0},{1,4,0,0},
		{0,5,0,0},{4,0,1,0},{3,1,1,0},{2,2,1,0},{1,3,1,0},
		{0,4,1,0},{3,0,2,0},{2,1,2,0},{1,2,2,0},{0,3,2,0},
		{2,0,3,0},{1,1,3,0},{0,2,3,0},{1,0,4,0},{0,1,4,0},
		{0,0,5,0},{4,0,0,1},{3,1,0,1},{2,2,0,1},{1,3,0,1},
		{0,4,0,1},{3,0,1,1},{2,1,1,1},{1,2,1,1},{0,3,1,1},
		{2,0,2,1},{1,1,2,1},{0,2,2,1},{1,0,3,1},{0,1,3,1},
		{0,0,4,1},{3,0,0,2},{2,1,0,2},{1,2,0,2},{0,3,0,2},
		{2,0,1,2},{1,1,1,2},{0,2,1,2},{1,0,2,2},{0,1,2,2},
		{0,0,3,2},{2,0,0,3},{1,1,0,3},{0,2,0,3},{1,0,1,3},
		{0,1,1,3},{0,0,2,3},{1,0,0,4},{0,1,0,4},{0,0,1,4},
		{0,0,0,5},
		{6,0,0,0},{5,1,0,0},{4,2,0,0},{3,3,0,0},{2,4,0,0},
		{1,5,0,0},{0,6,0,0},{5,0,1,0},{4,1,1,0},{3,2,1,0},
		{2,3,1,0},{1,4,1,0},{0,5,1,0},{4,0,2,0},{3,1,2,0},
		{2,2,2,0},{1,3,2,0},{0,4,2,0},{3,0,3,0},{2,1,3,0},
		{1,2,3,0},{0,3,3,0},{2,0,4,0},{1,1,4,0},{0,2,4,0},
		{1,0,5,0},{0,1,5,0},{0,0,6,0},{5,0,0,1},{4,1,0,1},
		{3,2,0,1},{2,3,0,1},{1,4,0,1},{0,5,0,1},{4,0,1,1},
		{3,1,1,1},{2,2,1,1},{1,3,1,1},{0,4,1,1},{3,0,2,1},
		{2,1,2,1},{1,2,2,1},{0,3,2,1},{2,0,3,1},{1,1,3,1},
		{0,2,3,1},{1,0,4,1},{0,1,4,1},{0,0,5,1},{4,0,0,2},
		{3,1,0,2},{2,2,0,2},{1,3,0,2},{0,4,0,2},{3,0,1,2},
		{2,1,1,2},{1,2,1,2},{0,3,1,2},{2,0,2,2},{1,1,2,2},
		{0,2,2,2},{1,0,3,2},{0,1,3,2},{0,0,4,2},{3,0,0,3},
		{2,1,0,3},{1,2,0,3},{0,3,0,3},{2,0,1,3},{1,1,1,3},
		{0,2,1,3},{1,0,2,3},{0,1,2,3},{0,0,3,3},{2,0,0,4},
		{1,1,0,4},{0,2,0,4},{1,0,1,4},{0,1,1,4},{0,0,2,4},
		{1,0,0,5},{0,1,0,5},{0,0,1,5},{0,0,0,6}};
  nnn = pwrsix[m+1];
  Apwr = (int *)malloc(nnn*sizeof(int));
  Bpwr = (int *)malloc(nnn*sizeof(int));
  Cpwr = (int *)malloc(nnn*sizeof(int));
  Dpwr = (int *)malloc(nnn*sizeof(int));

  if(debug) printf("terms used to find fit \n");
  for(i=0; i<=m; i++) // power requested                                  
  {
    for(ii=pwrsix[i]; ii<pwrsix[i+1]; ii++)
    {
      Apwr[ii]=pwrs[ii][0];
      Bpwr[ii]=pwrs[ii][1];
      Cpwr[ii]=pwrs[ii][2];
      Dpwr[ii]=pwrs[ii][3];
      if(debug) printf("%d   a^%d * b^%d * c^%d * d^%d\n",
                       ii, Apwr[ii],Bpwr[ii],Cpwr[ii],Dpwr[ii]);
    }
    if(debug) printf(" \n");
  }
} /* end gen_4d_powers */

static void gen_5d_powers()
{
  int i, ii;
  int pwrsix[8]={0,1,6,21,56,126,252,462}; // start of each order
  int pwrs[462][5]={{ 0, 0, 0, 0, 0},
                    { 1, 0, 0, 0, 0},{ 0, 1, 0, 0, 0},{ 0, 0, 1, 0, 0},
                    { 0, 0, 0, 1, 0},{ 0, 0, 0, 0, 1},
                    { 2, 0, 0, 0, 0},{ 1, 1, 0, 0, 0},{ 0, 2, 0, 0, 0},
                    { 1, 0, 1, 0, 0},{ 0, 1, 1, 0, 0},{ 0, 0, 2, 0, 0},
                    { 1, 0, 0, 1, 0},{ 0, 1, 0, 1, 0},{ 0, 0, 1, 1, 0},
                    { 0, 0, 0, 2, 0},{ 1, 0, 0, 0, 1},{ 0, 1, 0, 0, 1},
                    { 0, 0, 1, 0, 1},{ 0, 0, 0, 1, 1},{ 0, 0, 0, 0, 2},
                    { 3, 0, 0, 0, 0},{ 2, 1, 0, 0, 0},{ 1, 2, 0, 0, 0},
                    { 0, 3, 0, 0, 0},{ 2, 0, 1, 0, 0},{ 1, 1, 1, 0, 0},
                    { 0, 2, 1, 0, 0},{ 1, 0, 2, 0, 0},{ 0, 1, 2, 0, 0},
                    { 0, 0, 3, 0, 0},{ 2, 0, 0, 1, 0},{ 1, 1, 0, 1, 0},
                    { 0, 2, 0, 1, 0},{ 1, 0, 1, 1, 0},{ 0, 1, 1, 1, 0},
                    { 0, 0, 2, 1, 0},{ 1, 0, 0, 2, 0},{ 0, 1, 0, 2, 0},
                    { 0, 0, 1, 2, 0},{ 0, 0, 0, 3, 0},{ 2, 0, 0, 0, 1},
                    { 1, 1, 0, 0, 1},{ 0, 2, 0, 0, 1},{ 1, 0, 1, 0, 1},
                    { 0, 1, 1, 0, 1},{ 0, 0, 2, 0, 1},{ 1, 0, 0, 1, 1},
                    { 0, 1, 0, 1, 1},{ 0, 0, 1, 1, 1},{ 0, 0, 0, 2, 1},
                    { 1, 0, 0, 0, 2},{ 0, 1, 0, 0, 2},{ 0, 0, 1, 0, 2},
                    { 0, 0, 0, 1, 2},{ 0, 0, 0, 0, 3},
                    { 4, 0, 0, 0, 0},{ 3, 1, 0, 0, 0},{ 2, 2, 0, 0, 0},
                    { 1, 3, 0, 0, 0},{ 0, 4, 0, 0, 0},{ 3, 0, 1, 0, 0},
                    { 2, 1, 1, 0, 0},{ 1, 2, 1, 0, 0},{ 0, 3, 1, 0, 0},
                    { 2, 0, 2, 0, 0},{ 1, 1, 2, 0, 0},{ 0, 2, 2, 0, 0},
                    { 1, 0, 3, 0, 0},{ 0, 1, 3, 0, 0},{ 0, 0, 4, 0, 0},
                    { 3, 0, 0, 1, 0},{ 2, 1, 0, 1, 0},{ 1, 2, 0, 1, 0},
                    { 0, 3, 0, 1, 0},{ 2, 0, 1, 1, 0},{ 1, 1, 1, 1, 0},
                    { 0, 2, 1, 1, 0},{ 1, 0, 2, 1, 0},{ 0, 1, 2, 1, 0},
                    { 0, 0, 3, 1, 0},{ 2, 0, 0, 2, 0},{ 1, 1, 0, 2, 0},
                    { 0, 2, 0, 2, 0},{ 1, 0, 1, 2, 0},{ 0, 1, 1, 2, 0},
                    { 0, 0, 2, 2, 0},{ 1, 0, 0, 3, 0},{ 0, 1, 0, 3, 0},
                    { 0, 0, 1, 3, 0},{ 0, 0, 0, 4, 0},{ 3, 0, 0, 0, 1},
                    { 2, 1, 0, 0, 1},{ 1, 2, 0, 0, 1},{ 0, 3, 0, 0, 1},
                    { 2, 0, 1, 0, 1},{ 1, 1, 1, 0, 1},{ 0, 2, 1, 0, 1},
                    { 1, 0, 2, 0, 1},{ 0, 1, 2, 0, 1},{ 0, 0, 3, 0, 1},
                    { 2, 0, 0, 1, 1},{ 1, 1, 0, 1, 1},{ 0, 2, 0, 1, 1},
                    { 1, 0, 1, 1, 1},{ 0, 1, 1, 1, 1},{ 0, 0, 2, 1, 1},
                    { 1, 0, 0, 2, 1},{ 0, 1, 0, 2, 1},{ 0, 0, 1, 2, 1},
                    { 0, 0, 0, 3, 1},{ 2, 0, 0, 0, 2},{ 1, 1, 0, 0, 2},
                    { 0, 2, 0, 0, 2},{ 1, 0, 1, 0, 2},{ 0, 1, 1, 0, 2},
                    { 0, 0, 2, 0, 2},{ 1, 0, 0, 1, 2},{ 0, 1, 0, 1, 2},
                    { 0, 0, 1, 1, 2},{ 0, 0, 0, 2, 2},{ 1, 0, 0, 0, 3},
                    { 0, 1, 0, 0, 3},{ 0, 0, 1, 0, 3},{ 0, 0, 0, 1, 3},
                    { 0, 0, 0, 0, 4},
                    { 5, 0, 0, 0, 0},{ 4, 1, 0, 0, 0},{ 3, 2, 0, 0, 0},
                    { 2, 3, 0, 0, 0},{ 1, 4, 0, 0, 0},{ 0, 5, 0, 0, 0},
                    { 4, 0, 1, 0, 0},{ 3, 1, 1, 0, 0},{ 2, 2, 1, 0, 0},
                    { 1, 3, 1, 0, 0},{ 0, 4, 1, 0, 0},{ 3, 0, 2, 0, 0},
                    { 2, 1, 2, 0, 0},{ 1, 2, 2, 0, 0},{ 0, 3, 2, 0, 0},
                    { 2, 0, 3, 0, 0},{ 1, 1, 3, 0, 0},{ 0, 2, 3, 0, 0},
                    { 1, 0, 4, 0, 0},{ 0, 1, 4, 0, 0},{ 0, 0, 5, 0, 0},
                    { 4, 0, 0, 1, 0},{ 3, 1, 0, 1, 0},{ 2, 2, 0, 1, 0},
                    { 1, 3, 0, 1, 0},{ 0, 4, 0, 1, 0},{ 3, 0, 1, 1, 0},
                    { 2, 1, 1, 1, 0},{ 1, 2, 1, 1, 0},{ 0, 3, 1, 1, 0},
                    { 2, 0, 2, 1, 0},{ 1, 1, 2, 1, 0},{ 0, 2, 2, 1, 0},
                    { 1, 0, 3, 1, 0},{ 0, 1, 3, 1, 0},{ 0, 0, 4, 1, 0},
                    { 3, 0, 0, 2, 0},{ 2, 1, 0, 2, 0},{ 1, 2, 0, 2, 0},
                    { 0, 3, 0, 2, 0},{ 2, 0, 1, 2, 0},{ 1, 1, 1, 2, 0},
                    { 0, 2, 1, 2, 0},{ 1, 0, 2, 2, 0},{ 0, 1, 2, 2, 0},
                    { 0, 0, 3, 2, 0},{ 2, 0, 0, 3, 0},{ 1, 1, 0, 3, 0},
                    { 0, 2, 0, 3, 0},{ 1, 0, 1, 3, 0},{ 0, 1, 1, 3, 0},
                    { 0, 0, 2, 3, 0},{ 1, 0, 0, 4, 0},{ 0, 1, 0, 4, 0},
                    { 0, 0, 1, 4, 0},{ 0, 0, 0, 5, 0},{ 4, 0, 0, 0, 1},
                    { 3, 1, 0, 0, 1},{ 2, 2, 0, 0, 1},{ 1, 3, 0, 0, 1},
                    { 0, 4, 0, 0, 1},{ 3, 0, 1, 0, 1},{ 2, 1, 1, 0, 1},
                    { 1, 2, 1, 0, 1},{ 0, 3, 1, 0, 1},{ 2, 0, 2, 0, 1},
                    { 1, 1, 2, 0, 1},{ 0, 2, 2, 0, 1},{ 1, 0, 3, 0, 1},
                    { 0, 1, 3, 0, 1},{ 0, 0, 4, 0, 1},{ 3, 0, 0, 1, 1},
                    { 2, 1, 0, 1, 1},{ 1, 2, 0, 1, 1},{ 0, 3, 0, 1, 1},
                    { 2, 0, 1, 1, 1},{ 1, 1, 1, 1, 1},{ 0, 2, 1, 1, 1},
                    { 1, 0, 2, 1, 1},{ 0, 1, 2, 1, 1},{ 0, 0, 3, 1, 1},
                    { 2, 0, 0, 2, 1},{ 1, 1, 0, 2, 1},{ 0, 2, 0, 2, 1},
                    { 1, 0, 1, 2, 1},{ 0, 1, 1, 2, 1},{ 0, 0, 2, 2, 1},
                    { 1, 0, 0, 3, 1},{ 0, 1, 0, 3, 1},{ 0, 0, 1, 3, 1},
                    { 0, 0, 0, 4, 1},{ 3, 0, 0, 0, 2},{ 2, 1, 0, 0, 2},
                    { 1, 2, 0, 0, 2},{ 0, 3, 0, 0, 2},{ 2, 0, 1, 0, 2},
                    { 1, 1, 1, 0, 2},{ 0, 2, 1, 0, 2},{ 1, 0, 2, 0, 2},
                    { 0, 1, 2, 0, 2},{ 0, 0, 3, 0, 2},{ 2, 0, 0, 1, 2},
                    { 1, 1, 0, 1, 2},{ 0, 2, 0, 1, 2},{ 1, 0, 1, 1, 2},
                    { 0, 1, 1, 1, 2},{ 0, 0, 2, 1, 2},{ 1, 0, 0, 2, 2},
                    { 0, 1, 0, 2, 2},{ 0, 0, 1, 2, 2},{ 0, 0, 0, 3, 2},
                    { 2, 0, 0, 0, 3},{ 1, 1, 0, 0, 3},{ 0, 2, 0, 0, 3},
                    { 1, 0, 1, 0, 3},{ 0, 1, 1, 0, 3},{ 0, 0, 2, 0, 3},
                    { 1, 0, 0, 1, 3},{ 0, 1, 0, 1, 3},{ 0, 0, 1, 1, 3},
                    { 0, 0, 0, 2, 3},{ 1, 0, 0, 0, 4},{ 0, 1, 0, 0, 4},
                    { 0, 0, 1, 0, 4},{ 0, 0, 0, 1, 4},{ 0, 0, 0, 0, 5},
                    { 6, 0, 0, 0, 0},{ 5, 1, 0, 0, 0},{ 4, 2, 0, 0, 0},
                    { 3, 3, 0, 0, 0},{ 2, 4, 0, 0, 0},{ 1, 5, 0, 0, 0},
                    { 0, 6, 0, 0, 0},{ 5, 0, 1, 0, 0},{ 4, 1, 1, 0, 0},
                    { 3, 2, 1, 0, 0},{ 2, 3, 1, 0, 0},{ 1, 4, 1, 0, 0},
                    { 0, 5, 1, 0, 0},{ 4, 0, 2, 0, 0},{ 3, 1, 2, 0, 0},
                    { 2, 2, 2, 0, 0},{ 1, 3, 2, 0, 0},{ 0, 4, 2, 0, 0},
                    { 3, 0, 3, 0, 0},{ 2, 1, 3, 0, 0},{ 1, 2, 3, 0, 0},
                    { 0, 3, 3, 0, 0},{ 2, 0, 4, 0, 0},{ 1, 1, 4, 0, 0},
                    { 0, 2, 4, 0, 0},{ 1, 0, 5, 0, 0},{ 0, 1, 5, 0, 0},
                    { 0, 0, 6, 0, 0},{ 5, 0, 0, 1, 0},{ 4, 1, 0, 1, 0},
                    { 3, 2, 0, 1, 0},{ 2, 3, 0, 1, 0},{ 1, 4, 0, 1, 0},
                    { 0, 5, 0, 1, 0},{ 4, 0, 1, 1, 0},{ 3, 1, 1, 1, 0},
                    { 2, 2, 1, 1, 0},{ 1, 3, 1, 1, 0},{ 0, 4, 1, 1, 0},
                    { 3, 0, 2, 1, 0},{ 2, 1, 2, 1, 0},{ 1, 2, 2, 1, 0},
                    { 0, 3, 2, 1, 0},{ 2, 0, 3, 1, 0},{ 1, 1, 3, 1, 0},
                    { 0, 2, 3, 1, 0},{ 1, 0, 4, 1, 0},{ 0, 1, 4, 1, 0},
                    { 0, 0, 5, 1, 0},{ 4, 0, 0, 2, 0},{ 3, 1, 0, 2, 0},
                    { 2, 2, 0, 2, 0},{ 1, 3, 0, 2, 0},{ 0, 4, 0, 2, 0},
                    { 3, 0, 1, 2, 0},{ 2, 1, 1, 2, 0},{ 1, 2, 1, 2, 0},
                    { 0, 3, 1, 2, 0},{ 2, 0, 2, 2, 0},{ 1, 1, 2, 2, 0},
                    { 0, 2, 2, 2, 0},{ 1, 0, 3, 2, 0},{ 0, 1, 3, 2, 0},
                    { 0, 0, 4, 2, 0},{ 3, 0, 0, 3, 0},{ 2, 1, 0, 3, 0},
                    { 1, 2, 0, 3, 0},{ 0, 3, 0, 3, 0},{ 2, 0, 1, 3, 0},
                    { 1, 1, 1, 3, 0},{ 0, 2, 1, 3, 0},{ 1, 0, 2, 3, 0},
                    { 0, 1, 2, 3, 0},{ 0, 0, 3, 3, 0},{ 2, 0, 0, 4, 0},
                    { 1, 1, 0, 4, 0},{ 0, 2, 0, 4, 0},{ 1, 0, 1, 4, 0},
                    { 0, 1, 1, 4, 0},{ 0, 0, 2, 4, 0},{ 1, 0, 0, 5, 0},
                    { 0, 1, 0, 5, 0},{ 0, 0, 1, 5, 0},{ 0, 0, 0, 6, 0},
                    { 5, 0, 0, 0, 1},{ 4, 1, 0, 0, 1},{ 3, 2, 0, 0, 1},
                    { 2, 3, 0, 0, 1},{ 1, 4, 0, 0, 1},{ 0, 5, 0, 0, 1},
                    { 4, 0, 1, 0, 1},{ 3, 1, 1, 0, 1},{ 2, 2, 1, 0, 1},
                    { 1, 3, 1, 0, 1},{ 0, 4, 1, 0, 1},{ 3, 0, 2, 0, 1},
                    { 2, 1, 2, 0, 1},{ 1, 2, 2, 0, 1},{ 0, 3, 2, 0, 1},
                    { 2, 0, 3, 0, 1},{ 1, 1, 3, 0, 1},{ 0, 2, 3, 0, 1},
                    { 1, 0, 4, 0, 1},{ 0, 1, 4, 0, 1},{ 0, 0, 5, 0, 1},
                    { 4, 0, 0, 1, 1},{ 3, 1, 0, 1, 1},{ 2, 2, 0, 1, 1},
                    { 1, 3, 0, 1, 1},{ 0, 4, 0, 1, 1},{ 3, 0, 1, 1, 1},
                    { 2, 1, 1, 1, 1},{ 1, 2, 1, 1, 1},{ 0, 3, 1, 1, 1},
                    { 2, 0, 2, 1, 1},{ 1, 1, 2, 1, 1},{ 0, 2, 2, 1, 1},
                    { 1, 0, 3, 1, 1},{ 0, 1, 3, 1, 1},{ 0, 0, 4, 1, 1},
                    { 3, 0, 0, 2, 1},{ 2, 1, 0, 2, 1},{ 1, 2, 0, 2, 1},
                    { 0, 3, 0, 2, 1},{ 2, 0, 1, 2, 1},{ 1, 1, 1, 2, 1},
                    { 0, 2, 1, 2, 1},{ 1, 0, 2, 2, 1},{ 0, 1, 2, 2, 1},
                    { 0, 0, 3, 2, 1},{ 2, 0, 0, 3, 1},{ 1, 1, 0, 3, 1},
                    { 0, 2, 0, 3, 1},{ 1, 0, 1, 3, 1},{ 0, 1, 1, 3, 1},
                    { 0, 0, 2, 3, 1},{ 1, 0, 0, 4, 1},{ 0, 1, 0, 4, 1},
                    { 0, 0, 1, 4, 1},{ 0, 0, 0, 5, 1},{ 4, 0, 0, 0, 2},
                    { 3, 1, 0, 0, 2},{ 2, 2, 0, 0, 2},{ 1, 3, 0, 0, 2},
                    { 0, 4, 0, 0, 2},{ 3, 0, 1, 0, 2},{ 2, 1, 1, 0, 2},
                    { 1, 2, 1, 0, 2},{ 0, 3, 1, 0, 2},{ 2, 0, 2, 0, 2},
                    { 1, 1, 2, 0, 2},{ 0, 2, 2, 0, 2},{ 1, 0, 3, 0, 2},
                    { 0, 1, 3, 0, 2},{ 0, 0, 4, 0, 2},{ 3, 0, 0, 1, 2},
                    { 2, 1, 0, 1, 2},{ 1, 2, 0, 1, 2},{ 0, 3, 0, 1, 2},
                    { 2, 0, 1, 1, 2},{ 1, 1, 1, 1, 2},{ 0, 2, 1, 1, 2},
                    { 1, 0, 2, 1, 2},{ 0, 1, 2, 1, 2},{ 0, 0, 3, 1, 2},
                    { 2, 0, 0, 2, 2},{ 1, 1, 0, 2, 2},{ 0, 2, 0, 2, 2},
                    { 1, 0, 1, 2, 2},{ 0, 1, 1, 2, 2},{ 0, 0, 2, 2, 2},
                    { 1, 0, 0, 3, 2},{ 0, 1, 0, 3, 2},{ 0, 0, 1, 3, 2},
                    { 0, 0, 0, 4, 2},{ 3, 0, 0, 0, 3},{ 2, 1, 0, 0, 3},
                    { 1, 2, 0, 0, 3},{ 0, 3, 0, 0, 3},{ 2, 0, 1, 0, 3},
                    { 1, 1, 1, 0, 3},{ 0, 2, 1, 0, 3},{ 1, 0, 2, 0, 3},
                    { 0, 1, 2, 0, 3},{ 0, 0, 3, 0, 3},{ 2, 0, 0, 1, 3},
                    { 1, 1, 0, 1, 3},{ 0, 2, 0, 1, 3},{ 1, 0, 1, 1, 3},
                    { 0, 1, 1, 1, 3},{ 0, 0, 2, 1, 3},{ 1, 0, 0, 2, 3},
                    { 0, 1, 0, 2, 3},{ 0, 0, 1, 2, 3},{ 0, 0, 0, 3, 3},
                    { 2, 0, 0, 0, 4},{ 1, 1, 0, 0, 4},{ 0, 2, 0, 0, 4},
                    { 1, 0, 1, 0, 4},{ 0, 1, 1, 0, 4},{ 0, 0, 2, 0, 4},
                    { 1, 0, 0, 1, 4},{ 0, 1, 0, 1, 4},{ 0, 0, 1, 1, 4},
                    { 0, 0, 0, 2, 4},{ 1, 0, 0, 0, 5},{ 0, 1, 0, 0, 5},
                    { 0, 0, 1, 0, 5},{ 0, 0, 0, 1, 5},{ 0, 0, 0, 0, 6}};

  nnn = pwrsix[m+1];
  Apwr = (int *)malloc(nnn*sizeof(int));
  Bpwr = (int *)malloc(nnn*sizeof(int));
  Cpwr = (int *)malloc(nnn*sizeof(int));
  Dpwr = (int *)malloc(nnn*sizeof(int));
  Epwr = (int *)malloc(nnn*sizeof(int));

  if(debug) printf("terms used to find fit \n");
  for(i=0; i<=m; i++) // power requested                                  
  {
    for(ii=pwrsix[i]; ii<pwrsix[i+1]; ii++)
    {
      Apwr[ii]=pwrs[ii][0];
      Bpwr[ii]=pwrs[ii][1];
      Cpwr[ii]=pwrs[ii][2];
      Dpwr[ii]=pwrs[ii][3];
      Epwr[ii]=pwrs[ii][4];
      if(debug) printf("%d   a^%d * b^%d * c^%d * d^%d * e^%d\n",
                       ii, Apwr[ii],Bpwr[ii],Cpwr[ii],Dpwr[ii],Epwr[ii]);
    }
    if(debug) printf(" \n");
  }
} /* end gen_5d_powers */

static void gen_6d_powers()
{
  int i, ii;
  int pwrsix[8]={0,1,7,28,84,210,462,924}; // start of each order
  int pwrs[924][6]={
{0,0,0,0,0,0},
{1,0,0,0,0,0},{0,1,0,0,0,0},{0,0,1,0,0,0},{0,0,0,1,0,0},
{0,0,0,0,1,0},{0,0,0,0,0,1},
{2,0,0,0,0,0},{1,1,0,0,0,0},{0,2,0,0,0,0},{1,0,1,0,0,0},
{0,1,1,0,0,0},{0,0,2,0,0,0},{1,0,0,1,0,0},{0,1,0,1,0,0},
{0,0,1,1,0,0},{0,0,0,2,0,0},{1,0,0,0,1,0},{0,1,0,0,1,0},
{0,0,1,0,1,0},{0,0,0,1,1,0},{0,0,0,0,2,0},{1,0,0,0,0,1},
{0,1,0,0,0,1},{0,0,1,0,0,1},{0,0,0,1,0,1},{0,0,0,0,1,1},
{0,0,0,0,0,2},
{3,0,0,0,0,0},{2,1,0,0,0,0},{1,2,0,0,0,0},{0,3,0,0,0,0},
{2,0,1,0,0,0},{1,1,1,0,0,0},{0,2,1,0,0,0},{1,0,2,0,0,0},
{0,1,2,0,0,0},{0,0,3,0,0,0},{2,0,0,1,0,0},{1,1,0,1,0,0},
{0,2,0,1,0,0},{1,0,1,1,0,0},{0,1,1,1,0,0},{0,0,2,1,0,0},
{1,0,0,2,0,0},{0,1,0,2,0,0},{0,0,1,2,0,0},{0,0,0,3,0,0},
{2,0,0,0,1,0},{1,1,0,0,1,0},{0,2,0,0,1,0},{1,0,1,0,1,0},
{0,1,1,0,1,0},{0,0,2,0,1,0},{1,0,0,1,1,0},{0,1,0,1,1,0},
{0,0,1,1,1,0},{0,0,0,2,1,0},{1,0,0,0,2,0},{0,1,0,0,2,0},
{0,0,1,0,2,0},{0,0,0,1,2,0},{0,0,0,0,3,0},{2,0,0,0,0,1},
{1,1,0,0,0,1},{0,2,0,0,0,1},{1,0,1,0,0,1},{0,1,1,0,0,1},
{0,0,2,0,0,1},{1,0,0,1,0,1},{0,1,0,1,0,1},{0,0,1,1,0,1},
{0,0,0,2,0,1},{1,0,0,0,1,1},{0,1,0,0,1,1},{0,0,1,0,1,1},
{0,0,0,1,1,1},{0,0,0,0,2,1},{1,0,0,0,0,2},{0,1,0,0,0,2},
{0,0,1,0,0,2},{0,0,0,1,0,2},{0,0,0,0,1,2},{0,0,0,0,0,3},
{4,0,0,0,0,0},{3,1,0,0,0,0},{2,2,0,0,0,0},{1,3,0,0,0,0},
{0,4,0,0,0,0},{3,0,1,0,0,0},{2,1,1,0,0,0},{1,2,1,0,0,0},
{0,3,1,0,0,0},{2,0,2,0,0,0},{1,1,2,0,0,0},{0,2,2,0,0,0},
{1,0,3,0,0,0},{0,1,3,0,0,0},{0,0,4,0,0,0},{3,0,0,1,0,0},
{2,1,0,1,0,0},{1,2,0,1,0,0},{0,3,0,1,0,0},{2,0,1,1,0,0},
{1,1,1,1,0,0},{0,2,1,1,0,0},{1,0,2,1,0,0},{0,1,2,1,0,0},
{0,0,3,1,0,0},{2,0,0,2,0,0},{1,1,0,2,0,0},{0,2,0,2,0,0},
{1,0,1,2,0,0},{0,1,1,2,0,0},{0,0,2,2,0,0},{1,0,0,3,0,0},
{0,1,0,3,0,0},{0,0,1,3,0,0},{0,0,0,4,0,0},{3,0,0,0,1,0},
{2,1,0,0,1,0},{1,2,0,0,1,0},{0,3,0,0,1,0},{2,0,1,0,1,0},
{1,1,1,0,1,0},{0,2,1,0,1,0},{1,0,2,0,1,0},{0,1,2,0,1,0},
{0,0,3,0,1,0},{2,0,0,1,1,0},{1,1,0,1,1,0},{0,2,0,1,1,0},
{1,0,1,1,1,0},{0,1,1,1,1,0},{0,0,2,1,1,0},{1,0,0,2,1,0},
{0,1,0,2,1,0},{0,0,1,2,1,0},{0,0,0,3,1,0},{2,0,0,0,2,0},
{1,1,0,0,2,0},{0,2,0,0,2,0},{1,0,1,0,2,0},{0,1,1,0,2,0},
{0,0,2,0,2,0},{1,0,0,1,2,0},{0,1,0,1,2,0},{0,0,1,1,2,0},
{0,0,0,2,2,0},{1,0,0,0,3,0},{0,1,0,0,3,0},{0,0,1,0,3,0},
{0,0,0,1,3,0},{0,0,0,0,4,0},{3,0,0,0,0,1},{2,1,0,0,0,1},
{1,2,0,0,0,1},{0,3,0,0,0,1},{2,0,1,0,0,1},{1,1,1,0,0,1},
{0,2,1,0,0,1},{1,0,2,0,0,1},{0,1,2,0,0,1},{0,0,3,0,0,1},
{2,0,0,1,0,1},{1,1,0,1,0,1},{0,2,0,1,0,1},{1,0,1,1,0,1},
{0,1,1,1,0,1},{0,0,2,1,0,1},{1,0,0,2,0,1},{0,1,0,2,0,1},
{0,0,1,2,0,1},{0,0,0,3,0,1},{2,0,0,0,1,1},{1,1,0,0,1,1},
{0,2,0,0,1,1},{1,0,1,0,1,1},{0,1,1,0,1,1},{0,0,2,0,1,1},
{1,0,0,1,1,1},{0,1,0,1,1,1},{0,0,1,1,1,1},{0,0,0,2,1,1},
{1,0,0,0,2,1},{0,1,0,0,2,1},{0,0,1,0,2,1},{0,0,0,1,2,1},
{0,0,0,0,3,1},{2,0,0,0,0,2},{1,1,0,0,0,2},{0,2,0,0,0,2},
{1,0,1,0,0,2},{0,1,1,0,0,2},{0,0,2,0,0,2},{1,0,0,1,0,2},
{0,1,0,1,0,2},{0,0,1,1,0,2},{0,0,0,2,0,2},{1,0,0,0,1,2},
{0,1,0,0,1,2},{0,0,1,0,1,2},{0,0,0,1,1,2},{0,0,0,0,2,2},
{1,0,0,0,0,3},{0,1,0,0,0,3},{0,0,1,0,0,3},{0,0,0,1,0,3},
{0,0,0,0,1,3},{0,0,0,0,0,4},
{5,0,0,0,0,0},{4,1,0,0,0,0},{3,2,0,0,0,0},{2,3,0,0,0,0},
{1,4,0,0,0,0},{0,5,0,0,0,0},{4,0,1,0,0,0},{3,1,1,0,0,0},
{2,2,1,0,0,0},{1,3,1,0,0,0},{0,4,1,0,0,0},{3,0,2,0,0,0},
{2,1,2,0,0,0},{1,2,2,0,0,0},{0,3,2,0,0,0},{2,0,3,0,0,0},
{1,1,3,0,0,0},{0,2,3,0,0,0},{1,0,4,0,0,0},{0,1,4,0,0,0},
{0,0,5,0,0,0},{4,0,0,1,0,0},{3,1,0,1,0,0},{2,2,0,1,0,0},
{1,3,0,1,0,0},{0,4,0,1,0,0},{3,0,1,1,0,0},{2,1,1,1,0,0},
{1,2,1,1,0,0},{0,3,1,1,0,0},{2,0,2,1,0,0},{1,1,2,1,0,0},
{0,2,2,1,0,0},{1,0,3,1,0,0},{0,1,3,1,0,0},{0,0,4,1,0,0},
{3,0,0,2,0,0},{2,1,0,2,0,0},{1,2,0,2,0,0},{0,3,0,2,0,0},
{2,0,1,2,0,0},{1,1,1,2,0,0},{0,2,1,2,0,0},{1,0,2,2,0,0},
{0,1,2,2,0,0},{0,0,3,2,0,0},{2,0,0,3,0,0},{1,1,0,3,0,0},
{0,2,0,3,0,0},{1,0,1,3,0,0},{0,1,1,3,0,0},{0,0,2,3,0,0},
{1,0,0,4,0,0},{0,1,0,4,0,0},{0,0,1,4,0,0},{0,0,0,5,0,0},
{4,0,0,0,1,0},{3,1,0,0,1,0},{2,2,0,0,1,0},{1,3,0,0,1,0},
{0,4,0,0,1,0},{3,0,1,0,1,0},{2,1,1,0,1,0},{1,2,1,0,1,0},
{0,3,1,0,1,0},{2,0,2,0,1,0},{1,1,2,0,1,0},{0,2,2,0,1,0},
{1,0,3,0,1,0},{0,1,3,0,1,0},{0,0,4,0,1,0},{3,0,0,1,1,0},
{2,1,0,1,1,0},{1,2,0,1,1,0},{0,3,0,1,1,0},{2,0,1,1,1,0},
{1,1,1,1,1,0},{0,2,1,1,1,0},{1,0,2,1,1,0},{0,1,2,1,1,0},
{0,0,3,1,1,0},{2,0,0,2,1,0},{1,1,0,2,1,0},{0,2,0,2,1,0},
{1,0,1,2,1,0},{0,1,1,2,1,0},{0,0,2,2,1,0},{1,0,0,3,1,0},
{0,1,0,3,1,0},{0,0,1,3,1,0},{0,0,0,4,1,0},{3,0,0,0,2,0},
{2,1,0,0,2,0},{1,2,0,0,2,0},{0,3,0,0,2,0},{2,0,1,0,2,0},
{1,1,1,0,2,0},{0,2,1,0,2,0},{1,0,2,0,2,0},{0,1,2,0,2,0},
{0,0,3,0,2,0},{2,0,0,1,2,0},{1,1,0,1,2,0},{0,2,0,1,2,0},
{1,0,1,1,2,0},{0,1,1,1,2,0},{0,0,2,1,2,0},{1,0,0,2,2,0},
{0,1,0,2,2,0},{0,0,1,2,2,0},{0,0,0,3,2,0},{2,0,0,0,3,0},
{1,1,0,0,3,0},{0,2,0,0,3,0},{1,0,1,0,3,0},{0,1,1,0,3,0},
{0,0,2,0,3,0},{1,0,0,1,3,0},{0,1,0,1,3,0},{0,0,1,1,3,0},
{0,0,0,2,3,0},{1,0,0,0,4,0},{0,1,0,0,4,0},{0,0,1,0,4,0},
{0,0,0,1,4,0},{0,0,0,0,5,0},{4,0,0,0,0,1},{3,1,0,0,0,1},
{2,2,0,0,0,1},{1,3,0,0,0,1},{0,4,0,0,0,1},{3,0,1,0,0,1},
{2,1,1,0,0,1},{1,2,1,0,0,1},{0,3,1,0,0,1},{2,0,2,0,0,1},
{1,1,2,0,0,1},{0,2,2,0,0,1},{1,0,3,0,0,1},{0,1,3,0,0,1},
{0,0,4,0,0,1},{3,0,0,1,0,1},{2,1,0,1,0,1},{1,2,0,1,0,1},
{0,3,0,1,0,1},{2,0,1,1,0,1},{1,1,1,1,0,1},{0,2,1,1,0,1},
{1,0,2,1,0,1},{0,1,2,1,0,1},{0,0,3,1,0,1},{2,0,0,2,0,1},
{1,1,0,2,0,1},{0,2,0,2,0,1},{1,0,1,2,0,1},{0,1,1,2,0,1},
{0,0,2,2,0,1},{1,0,0,3,0,1},{0,1,0,3,0,1},{0,0,1,3,0,1},
{0,0,0,4,0,1},{3,0,0,0,1,1},{2,1,0,0,1,1},{1,2,0,0,1,1},
{0,3,0,0,1,1},{2,0,1,0,1,1},{1,1,1,0,1,1},{0,2,1,0,1,1},
{1,0,2,0,1,1},{0,1,2,0,1,1},{0,0,3,0,1,1},{2,0,0,1,1,1},
{1,1,0,1,1,1},{0,2,0,1,1,1},{1,0,1,1,1,1},{0,1,1,1,1,1},
{0,0,2,1,1,1},{1,0,0,2,1,1},{0,1,0,2,1,1},{0,0,1,2,1,1},
{0,0,0,3,1,1},{2,0,0,0,2,1},{1,1,0,0,2,1},{0,2,0,0,2,1},
{1,0,1,0,2,1},{0,1,1,0,2,1},{0,0,2,0,2,1},{1,0,0,1,2,1},
{0,1,0,1,2,1},{0,0,1,1,2,1},{0,0,0,2,2,1},{1,0,0,0,3,1},
{0,1,0,0,3,1},{0,0,1,0,3,1},{0,0,0,1,3,1},{0,0,0,0,4,1},
{3,0,0,0,0,2},{2,1,0,0,0,2},{1,2,0,0,0,2},{0,3,0,0,0,2},
{2,0,1,0,0,2},{1,1,1,0,0,2},{0,2,1,0,0,2},{1,0,2,0,0,2},
{0,1,2,0,0,2},{0,0,3,0,0,2},{2,0,0,1,0,2},{1,1,0,1,0,2},
{0,2,0,1,0,2},{1,0,1,1,0,2},{0,1,1,1,0,2},{0,0,2,1,0,2},
{1,0,0,2,0,2},{0,1,0,2,0,2},{0,0,1,2,0,2},{0,0,0,3,0,2},
{2,0,0,0,1,2},{1,1,0,0,1,2},{0,2,0,0,1,2},{1,0,1,0,1,2},
{0,1,1,0,1,2},{0,0,2,0,1,2},{1,0,0,1,1,2},{0,1,0,1,1,2},
{0,0,1,1,1,2},{0,0,0,2,1,2},{1,0,0,0,2,2},{0,1,0,0,2,2},
{0,0,1,0,2,2},{0,0,0,1,2,2},{0,0,0,0,3,2},{2,0,0,0,0,3},
{1,1,0,0,0,3},{0,2,0,0,0,3},{1,0,1,0,0,3},{0,1,1,0,0,3},
{0,0,2,0,0,3},{1,0,0,1,0,3},{0,1,0,1,0,3},{0,0,1,1,0,3},
{0,0,0,2,0,3},{1,0,0,0,1,3},{0,1,0,0,1,3},{0,0,1,0,1,3},
{0,0,0,1,1,3},{0,0,0,0,2,3},{1,0,0,0,0,4},{0,1,0,0,0,4},
{0,0,1,0,0,4},{0,0,0,1,0,4},{0,0,0,0,1,4},{0,0,0,0,0,5},
{6,0,0,0,0,0},{5,1,0,0,0,0},{4,2,0,0,0,0},{3,3,0,0,0,0},
{2,4,0,0,0,0},{1,5,0,0,0,0},{0,6,0,0,0,0},{5,0,1,0,0,0},
{4,1,1,0,0,0},{3,2,1,0,0,0},{2,3,1,0,0,0},{1,4,1,0,0,0},
{0,5,1,0,0,0},{4,0,2,0,0,0},{3,1,2,0,0,0},{2,2,2,0,0,0},
{1,3,2,0,0,0},{0,4,2,0,0,0},{3,0,3,0,0,0},{2,1,3,0,0,0},
{1,2,3,0,0,0},{0,3,3,0,0,0},{2,0,4,0,0,0},{1,1,4,0,0,0},
{0,2,4,0,0,0},{1,0,5,0,0,0},{0,1,5,0,0,0},{0,0,6,0,0,0},
{5,0,0,1,0,0},{4,1,0,1,0,0},{3,2,0,1,0,0},{2,3,0,1,0,0},
{1,4,0,1,0,0},{0,5,0,1,0,0},{4,0,1,1,0,0},{3,1,1,1,0,0},
{2,2,1,1,0,0},{1,3,1,1,0,0},{0,4,1,1,0,0},{3,0,2,1,0,0},
{2,1,2,1,0,0},{1,2,2,1,0,0},{0,3,2,1,0,0},{2,0,3,1,0,0},
{1,1,3,1,0,0},{0,2,3,1,0,0},{1,0,4,1,0,0},{0,1,4,1,0,0},
{0,0,5,1,0,0},{4,0,0,2,0,0},{3,1,0,2,0,0},{2,2,0,2,0,0},
{1,3,0,2,0,0},{0,4,0,2,0,0},{3,0,1,2,0,0},{2,1,1,2,0,0},
{1,2,1,2,0,0},{0,3,1,2,0,0},{2,0,2,2,0,0},{1,1,2,2,0,0},
{0,2,2,2,0,0},{1,0,3,2,0,0},{0,1,3,2,0,0},{0,0,4,2,0,0},
{3,0,0,3,0,0},{2,1,0,3,0,0},{1,2,0,3,0,0},{0,3,0,3,0,0},
{2,0,1,3,0,0},{1,1,1,3,0,0},{0,2,1,3,0,0},{1,0,2,3,0,0},
{0,1,2,3,0,0},{0,0,3,3,0,0},{2,0,0,4,0,0},{1,1,0,4,0,0},
{0,2,0,4,0,0},{1,0,1,4,0,0},{0,1,1,4,0,0},{0,0,2,4,0,0},
{1,0,0,5,0,0},{0,1,0,5,0,0},{0,0,1,5,0,0},{0,0,0,6,0,0},
{5,0,0,0,1,0},{4,1,0,0,1,0},{3,2,0,0,1,0},{2,3,0,0,1,0},
{1,4,0,0,1,0},{0,5,0,0,1,0},{4,0,1,0,1,0},{3,1,1,0,1,0},
{2,2,1,0,1,0},{1,3,1,0,1,0},{0,4,1,0,1,0},{3,0,2,0,1,0},
{2,1,2,0,1,0},{1,2,2,0,1,0},{0,3,2,0,1,0},{2,0,3,0,1,0},
{1,1,3,0,1,0},{0,2,3,0,1,0},{1,0,4,0,1,0},{0,1,4,0,1,0},
{0,0,5,0,1,0},{4,0,0,1,1,0},{3,1,0,1,1,0},{2,2,0,1,1,0},
{1,3,0,1,1,0},{0,4,0,1,1,0},{3,0,1,1,1,0},{2,1,1,1,1,0},
{1,2,1,1,1,0},{0,3,1,1,1,0},{2,0,2,1,1,0},{1,1,2,1,1,0},
{0,2,2,1,1,0},{1,0,3,1,1,0},{0,1,3,1,1,0},{0,0,4,1,1,0},
{3,0,0,2,1,0},{2,1,0,2,1,0},{1,2,0,2,1,0},{0,3,0,2,1,0},
{2,0,1,2,1,0},{1,1,1,2,1,0},{0,2,1,2,1,0},{1,0,2,2,1,0},
{0,1,2,2,1,0},{0,0,3,2,1,0},{2,0,0,3,1,0},{1,1,0,3,1,0},
{0,2,0,3,1,0},{1,0,1,3,1,0},{0,1,1,3,1,0},{0,0,2,3,1,0},
{1,0,0,4,1,0},{0,1,0,4,1,0},{0,0,1,4,1,0},{0,0,0,5,1,0},
{4,0,0,0,2,0},{3,1,0,0,2,0},{2,2,0,0,2,0},{1,3,0,0,2,0},
{0,4,0,0,2,0},{3,0,1,0,2,0},{2,1,1,0,2,0},{1,2,1,0,2,0},
{0,3,1,0,2,0},{2,0,2,0,2,0},{1,1,2,0,2,0},{0,2,2,0,2,0},
{1,0,3,0,2,0},{0,1,3,0,2,0},{0,0,4,0,2,0},{3,0,0,1,2,0},
{2,1,0,1,2,0},{1,2,0,1,2,0},{0,3,0,1,2,0},{2,0,1,1,2,0},
{1,1,1,1,2,0},{0,2,1,1,2,0},{1,0,2,1,2,0},{0,1,2,1,2,0},
{0,0,3,1,2,0},{2,0,0,2,2,0},{1,1,0,2,2,0},{0,2,0,2,2,0},
{1,0,1,2,2,0},{0,1,1,2,2,0},{0,0,2,2,2,0},{1,0,0,3,2,0},
{0,1,0,3,2,0},{0,0,1,3,2,0},{0,0,0,4,2,0},{3,0,0,0,3,0},
{2,1,0,0,3,0},{1,2,0,0,3,0},{0,3,0,0,3,0},{2,0,1,0,3,0},
{1,1,1,0,3,0},{0,2,1,0,3,0},{1,0,2,0,3,0},{0,1,2,0,3,0},
{0,0,3,0,3,0},{2,0,0,1,3,0},{1,1,0,1,3,0},{0,2,0,1,3,0},
{1,0,1,1,3,0},{0,1,1,1,3,0},{0,0,2,1,3,0},{1,0,0,2,3,0},
{0,1,0,2,3,0},{0,0,1,2,3,0},{0,0,0,3,3,0},{2,0,0,0,4,0},
{1,1,0,0,4,0},{0,2,0,0,4,0},{1,0,1,0,4,0},{0,1,1,0,4,0},
{0,0,2,0,4,0},{1,0,0,1,4,0},{0,1,0,1,4,0},{0,0,1,1,4,0},
{0,0,0,2,4,0},{1,0,0,0,5,0},{0,1,0,0,5,0},{0,0,1,0,5,0},
{0,0,0,1,5,0},{0,0,0,0,6,0},{5,0,0,0,0,1},{4,1,0,0,0,1},
{3,2,0,0,0,1},{2,3,0,0,0,1},{1,4,0,0,0,1},{0,5,0,0,0,1},
{4,0,1,0,0,1},{3,1,1,0,0,1},{2,2,1,0,0,1},{1,3,1,0,0,1},
{0,4,1,0,0,1},{3,0,2,0,0,1},{2,1,2,0,0,1},{1,2,2,0,0,1},
{0,3,2,0,0,1},{2,0,3,0,0,1},{1,1,3,0,0,1},{0,2,3,0,0,1},
{1,0,4,0,0,1},{0,1,4,0,0,1},{0,0,5,0,0,1},{4,0,0,1,0,1},
{3,1,0,1,0,1},{2,2,0,1,0,1},{1,3,0,1,0,1},{0,4,0,1,0,1},
{3,0,1,1,0,1},{2,1,1,1,0,1},{1,2,1,1,0,1},{0,3,1,1,0,1},
{2,0,2,1,0,1},{1,1,2,1,0,1},{0,2,2,1,0,1},{1,0,3,1,0,1},
{0,1,3,1,0,1},{0,0,4,1,0,1},{3,0,0,2,0,1},{2,1,0,2,0,1},
{1,2,0,2,0,1},{0,3,0,2,0,1},{2,0,1,2,0,1},{1,1,1,2,0,1},
{0,2,1,2,0,1},{1,0,2,2,0,1},{0,1,2,2,0,1},{0,0,3,2,0,1},
{2,0,0,3,0,1},{1,1,0,3,0,1},{0,2,0,3,0,1},{1,0,1,3,0,1},
{0,1,1,3,0,1},{0,0,2,3,0,1},{1,0,0,4,0,1},{0,1,0,4,0,1},
{0,0,1,4,0,1},{0,0,0,5,0,1},{4,0,0,0,1,1},{3,1,0,0,1,1},
{2,2,0,0,1,1},{1,3,0,0,1,1},{0,4,0,0,1,1},{3,0,1,0,1,1},
{2,1,1,0,1,1},{1,2,1,0,1,1},{0,3,1,0,1,1},{2,0,2,0,1,1},
{1,1,2,0,1,1},{0,2,2,0,1,1},{1,0,3,0,1,1},{0,1,3,0,1,1},
{0,0,4,0,1,1},{3,0,0,1,1,1},{2,1,0,1,1,1},{1,2,0,1,1,1},
{0,3,0,1,1,1},{2,0,1,1,1,1},{1,1,1,1,1,1},{0,2,1,1,1,1},
{1,0,2,1,1,1},{0,1,2,1,1,1},{0,0,3,1,1,1},{2,0,0,2,1,1},
{1,1,0,2,1,1},{0,2,0,2,1,1},{1,0,1,2,1,1},{0,1,1,2,1,1},
{0,0,2,2,1,1},{1,0,0,3,1,1},{0,1,0,3,1,1},{0,0,1,3,1,1},
{0,0,0,4,1,1},{3,0,0,0,2,1},{2,1,0,0,2,1},{1,2,0,0,2,1},
{0,3,0,0,2,1},{2,0,1,0,2,1},{1,1,1,0,2,1},{0,2,1,0,2,1},
{1,0,2,0,2,1},{0,1,2,0,2,1},{0,0,3,0,2,1},{2,0,0,1,2,1},
{1,1,0,1,2,1},{0,2,0,1,2,1},{1,0,1,1,2,1},{0,1,1,1,2,1},
{0,0,2,1,2,1},{1,0,0,2,2,1},{0,1,0,2,2,1},{0,0,1,2,2,1},
{0,0,0,3,2,1},{2,0,0,0,3,1},{1,1,0,0,3,1},{0,2,0,0,3,1},
{1,0,1,0,3,1},{0,1,1,0,3,1},{0,0,2,0,3,1},{1,0,0,1,3,1},
{0,1,0,1,3,1},{0,0,1,1,3,1},{0,0,0,2,3,1},{1,0,0,0,4,1},
{0,1,0,0,4,1},{0,0,1,0,4,1},{0,0,0,1,4,1},{0,0,0,0,5,1},
{4,0,0,0,0,2},{3,1,0,0,0,2},{2,2,0,0,0,2},{1,3,0,0,0,2},
{0,4,0,0,0,2},{3,0,1,0,0,2},{2,1,1,0,0,2},{1,2,1,0,0,2},
{0,3,1,0,0,2},{2,0,2,0,0,2},{1,1,2,0,0,2},{0,2,2,0,0,2},
{1,0,3,0,0,2},{0,1,3,0,0,2},{0,0,4,0,0,2},{3,0,0,1,0,2},
{2,1,0,1,0,2},{1,2,0,1,0,2},{0,3,0,1,0,2},{2,0,1,1,0,2},
{1,1,1,1,0,2},{0,2,1,1,0,2},{1,0,2,1,0,2},{0,1,2,1,0,2},
{0,0,3,1,0,2},{2,0,0,2,0,2},{1,1,0,2,0,2},{0,2,0,2,0,2},
{1,0,1,2,0,2},{0,1,1,2,0,2},{0,0,2,2,0,2},{1,0,0,3,0,2},
{0,1,0,3,0,2},{0,0,1,3,0,2},{0,0,0,4,0,2},{3,0,0,0,1,2},
{2,1,0,0,1,2},{1,2,0,0,1,2},{0,3,0,0,1,2},{2,0,1,0,1,2},
{1,1,1,0,1,2},{0,2,1,0,1,2},{1,0,2,0,1,2},{0,1,2,0,1,2},
{0,0,3,0,1,2},{2,0,0,1,1,2},{1,1,0,1,1,2},{0,2,0,1,1,2},
{1,0,1,1,1,2},{0,1,1,1,1,2},{0,0,2,1,1,2},{1,0,0,2,1,2},
{0,1,0,2,1,2},{0,0,1,2,1,2},{0,0,0,3,1,2},{2,0,0,0,2,2},
{1,1,0,0,2,2},{0,2,0,0,2,2},{1,0,1,0,2,2},{0,1,1,0,2,2},
{0,0,2,0,2,2},{1,0,0,1,2,2},{0,1,0,1,2,2},{0,0,1,1,2,2},
{0,0,0,2,2,2},{1,0,0,0,3,2},{0,1,0,0,3,2},{0,0,1,0,3,2},
{0,0,0,1,3,2},{0,0,0,0,4,2},{3,0,0,0,0,3},{2,1,0,0,0,3},
{1,2,0,0,0,3},{0,3,0,0,0,3},{2,0,1,0,0,3},{1,1,1,0,0,3},
{0,2,1,0,0,3},{1,0,2,0,0,3},{0,1,2,0,0,3},{0,0,3,0,0,3},
{2,0,0,1,0,3},{1,1,0,1,0,3},{0,2,0,1,0,3},{1,0,1,1,0,3},
{0,1,1,1,0,3},{0,0,2,1,0,3},{1,0,0,2,0,3},{0,1,0,2,0,3},
{0,0,1,2,0,3},{0,0,0,3,0,3},{2,0,0,0,1,3},{1,1,0,0,1,3},
{0,2,0,0,1,3},{1,0,1,0,1,3},{0,1,1,0,1,3},{0,0,2,0,1,3},
{1,0,0,1,1,3},{0,1,0,1,1,3},{0,0,1,1,1,3},{0,0,0,2,1,3},
{1,0,0,0,2,3},{0,1,0,0,2,3},{0,0,1,0,2,3},{0,0,0,1,2,3},
{0,0,0,0,3,3},{2,0,0,0,0,4},{1,1,0,0,0,4},{0,2,0,0,0,4},
{1,0,1,0,0,4},{0,1,1,0,0,4},{0,0,2,0,0,4},{1,0,0,1,0,4},
{0,1,0,1,0,4},{0,0,1,1,0,4},{0,0,0,2,0,4},{1,0,0,0,1,4},
{0,1,0,0,1,4},{0,0,1,0,1,4},{0,0,0,1,1,4},{0,0,0,0,2,4},
{1,0,0,0,0,5},{0,1,0,0,0,5},{0,0,1,0,0,5},{0,0,0,1,0,5},
{0,0,0,0,1,5},{0,0,0,0,0,6}};

  nnn = pwrsix[m+1];
  Apwr = (int *)malloc(nnn*sizeof(int));
  Bpwr = (int *)malloc(nnn*sizeof(int));
  Cpwr = (int *)malloc(nnn*sizeof(int));
  Dpwr = (int *)malloc(nnn*sizeof(int));
  Epwr = (int *)malloc(nnn*sizeof(int));
  Fpwr = (int *)malloc(nnn*sizeof(int));

  if(debug) printf("terms used to find fit \n");
  for(i=0; i<=m; i++) // power requested                                  
  {
    for(ii=pwrsix[i]; ii<pwrsix[i+1]; ii++)
    {
      Apwr[ii]=pwrs[ii][0];
      Bpwr[ii]=pwrs[ii][1];
      Cpwr[ii]=pwrs[ii][2];
      Dpwr[ii]=pwrs[ii][3];
      Epwr[ii]=pwrs[ii][4];
      Fpwr[ii]=pwrs[ii][5];
      if(debug) printf("%d   a^%d * b^%d * c^%d * d^%d * e^%d * f^%d\n",
                   ii, Apwr[ii],Bpwr[ii],Cpwr[ii],Dpwr[ii],Epwr[ii],Fpwr[ii]);
    }
    if(debug) printf(" \n");
  }
} /* end gen_6d_powers */


static void gen_7d_powers()
{
  int i, ii;
  int pwrsix[6]={0,1,8,36,120,330}; // start of each order
  int pwrs[330][7]={
{0,0,0,0,0,0,0},

{1,0,0,0,0,0,0},
{0,1,0,0,0,0,0},
{0,0,1,0,0,0,0},
{0,0,0,1,0,0,0},
{0,0,0,0,1,0,0},
{0,0,0,0,0,1,0},
{0,0,0,0,0,0,1},

{2,0,0,0,0,0,0},
{1,1,0,0,0,0,0},
{0,2,0,0,0,0,0},
{1,0,1,0,0,0,0},
{0,1,1,0,0,0,0},
{0,0,2,0,0,0,0},
{1,0,0,1,0,0,0},
{0,1,0,1,0,0,0},
{0,0,1,1,0,0,0},
{0,0,0,2,0,0,0},
{1,0,0,0,1,0,0},
{0,1,0,0,1,0,0},
{0,0,1,0,1,0,0},
{0,0,0,1,1,0,0},
{0,0,0,0,2,0,0},
{1,0,0,0,0,1,0},
{0,1,0,0,0,1,0},
{0,0,1,0,0,1,0},
{0,0,0,1,0,1,0},
{0,0,0,0,1,1,0},
{0,0,0,0,0,2,0},
{1,0,0,0,0,0,1},
{0,1,0,0,0,0,1},
{0,0,1,0,0,0,1},
{0,0,0,1,0,0,1},
{0,0,0,0,1,0,1},
{0,0,0,0,0,1,1},
{0,0,0,0,0,0,2},

{3,0,0,0,0,0,0},
{2,1,0,0,0,0,0},
{1,2,0,0,0,0,0},
{0,3,0,0,0,0,0},
{2,0,1,0,0,0,0},
{1,1,1,0,0,0,0},
{0,2,1,0,0,0,0},
{1,0,2,0,0,0,0},
{0,1,2,0,0,0,0},
{0,0,3,0,0,0,0},
{2,0,0,1,0,0,0},
{1,1,0,1,0,0,0},
{0,2,0,1,0,0,0},
{1,0,1,1,0,0,0},
{0,1,1,1,0,0,0},
{0,0,2,1,0,0,0},
{1,0,0,2,0,0,0},
{0,1,0,2,0,0,0},
{0,0,1,2,0,0,0},
{0,0,0,3,0,0,0},
{2,0,0,0,1,0,0},
{1,1,0,0,1,0,0},
{0,2,0,0,1,0,0},
{1,0,1,0,1,0,0},
{0,1,1,0,1,0,0},
{0,0,2,0,1,0,0},
{1,0,0,1,1,0,0},
{0,1,0,1,1,0,0},
{0,0,1,1,1,0,0},
{0,0,0,2,1,0,0},
{1,0,0,0,2,0,0},
{0,1,0,0,2,0,0},
{0,0,1,0,2,0,0},
{0,0,0,1,2,0,0},
{0,0,0,0,3,0,0},
{2,0,0,0,0,1,0},
{1,1,0,0,0,1,0},
{0,2,0,0,0,1,0},
{1,0,1,0,0,1,0},
{0,1,1,0,0,1,0},
{0,0,2,0,0,1,0},
{1,0,0,1,0,1,0},
{0,1,0,1,0,1,0},
{0,0,1,1,0,1,0},
{0,0,0,2,0,1,0},
{1,0,0,0,1,1,0},
{0,1,0,0,1,1,0},
{0,0,1,0,1,1,0},
{0,0,0,1,1,1,0},
{0,0,0,0,2,1,0},
{1,0,0,0,0,2,0},
{0,1,0,0,0,2,0},
{0,0,1,0,0,2,0},
{0,0,0,1,0,2,0},
{0,0,0,0,1,2,0},
{0,0,0,0,0,3,0},
{2,0,0,0,0,0,1},
{1,1,0,0,0,0,1},
{0,2,0,0,0,0,1},
{1,0,1,0,0,0,1},
{0,1,1,0,0,0,1},
{0,0,2,0,0,0,1},
{1,0,0,1,0,0,1},
{0,1,0,1,0,0,1},
{0,0,1,1,0,0,1},
{0,0,0,2,0,0,1},
{1,0,0,0,1,0,1},
{0,1,0,0,1,0,1},
{0,0,1,0,1,0,1},
{0,0,0,1,1,0,1},
{0,0,0,0,2,0,1},
{1,0,0,0,0,1,1},
{0,1,0,0,0,1,1},
{0,0,1,0,0,1,1},
{0,0,0,1,0,1,1},
{0,0,0,0,1,1,1},
{0,0,0,0,0,2,1},
{1,0,0,0,0,0,2},
{0,1,0,0,0,0,2},
{0,0,1,0,0,0,2},
{0,0,0,1,0,0,2},
{0,0,0,0,1,0,2},
{0,0,0,0,0,1,2},
{0,0,0,0,0,0,3},

{4,0,0,0,0,0,0},
{3,1,0,0,0,0,0},
{2,2,0,0,0,0,0},
{1,3,0,0,0,0,0},
{0,4,0,0,0,0,0},
{3,0,1,0,0,0,0},
{2,1,1,0,0,0,0},
{1,2,1,0,0,0,0},
{0,3,1,0,0,0,0},
{2,0,2,0,0,0,0},
{1,1,2,0,0,0,0},
{0,2,2,0,0,0,0},
{1,0,3,0,0,0,0},
{0,1,3,0,0,0,0},
{0,0,4,0,0,0,0},
{3,0,0,1,0,0,0},
{2,1,0,1,0,0,0},
{1,2,0,1,0,0,0},
{0,3,0,1,0,0,0},
{2,0,1,1,0,0,0},
{1,1,1,1,0,0,0},
{0,2,1,1,0,0,0},
{1,0,2,1,0,0,0},
{0,1,2,1,0,0,0},
{0,0,3,1,0,0,0},
{2,0,0,2,0,0,0},
{1,1,0,2,0,0,0},
{0,2,0,2,0,0,0},
{1,0,1,2,0,0,0},
{0,1,1,2,0,0,0},
{0,0,2,2,0,0,0},
{1,0,0,3,0,0,0},
{0,1,0,3,0,0,0},
{0,0,1,3,0,0,0},
{0,0,0,4,0,0,0},
{3,0,0,0,1,0,0},
{2,1,0,0,1,0,0},
{1,2,0,0,1,0,0},
{0,3,0,0,1,0,0},
{2,0,1,0,1,0,0},
{1,1,1,0,1,0,0},
{0,2,1,0,1,0,0},
{1,0,2,0,1,0,0},
{0,1,2,0,1,0,0},
{0,0,3,0,1,0,0},
{2,0,0,1,1,0,0},
{1,1,0,1,1,0,0},
{0,2,0,1,1,0,0},
{1,0,1,1,1,0,0},
{0,1,1,1,1,0,0},
{0,0,2,1,1,0,0},
{1,0,0,2,1,0,0},
{0,1,0,2,1,0,0},
{0,0,1,2,1,0,0},
{0,0,0,3,1,0,0},
{2,0,0,0,2,0,0},
{1,1,0,0,2,0,0},
{0,2,0,0,2,0,0},
{1,0,1,0,2,0,0},
{0,1,1,0,2,0,0},
{0,0,2,0,2,0,0},
{1,0,0,1,2,0,0},
{0,1,0,1,2,0,0},
{0,0,1,1,2,0,0},
{0,0,0,2,2,0,0},
{1,0,0,0,3,0,0},
{0,1,0,0,3,0,0},
{0,0,1,0,3,0,0},
{0,0,0,1,3,0,0},
{0,0,0,0,4,0,0},
{3,0,0,0,0,1,0},
{2,1,0,0,0,1,0},
{1,2,0,0,0,1,0},
{0,3,0,0,0,1,0},
{2,0,1,0,0,1,0},
{1,1,1,0,0,1,0},
{0,2,1,0,0,1,0},
{1,0,2,0,0,1,0},
{0,1,2,0,0,1,0},
{0,0,3,0,0,1,0},
{2,0,0,1,0,1,0},
{1,1,0,1,0,1,0},
{0,2,0,1,0,1,0},
{1,0,1,1,0,1,0},
{0,1,1,1,0,1,0},
{0,0,2,1,0,1,0},
{1,0,0,2,0,1,0},
{0,1,0,2,0,1,0},
{0,0,1,2,0,1,0},
{0,0,0,3,0,1,0},
{2,0,0,0,1,1,0},
{1,1,0,0,1,1,0},
{0,2,0,0,1,1,0},
{1,0,1,0,1,1,0},
{0,1,1,0,1,1,0},
{0,0,2,0,1,1,0},
{1,0,0,1,1,1,0},
{0,1,0,1,1,1,0},
{0,0,1,1,1,1,0},
{0,0,0,2,1,1,0},
{1,0,0,0,2,1,0},
{0,1,0,0,2,1,0},
{0,0,1,0,2,1,0},
{0,0,0,1,2,1,0},
{0,0,0,0,3,1,0},
{2,0,0,0,0,2,0},
{1,1,0,0,0,2,0},
{0,2,0,0,0,2,0},
{1,0,1,0,0,2,0},
{0,1,1,0,0,2,0},
{0,0,2,0,0,2,0},
{1,0,0,1,0,2,0},
{0,1,0,1,0,2,0},
{0,0,1,1,0,2,0},
{0,0,0,2,0,2,0},
{1,0,0,0,1,2,0},
{0,1,0,0,1,2,0},
{0,0,1,0,1,2,0},
{0,0,0,1,1,2,0},
{0,0,0,0,2,2,0},
{1,0,0,0,0,3,0},
{0,1,0,0,0,3,0},
{0,0,1,0,0,3,0},
{0,0,0,1,0,3,0},
{0,0,0,0,1,3,0},
{0,0,0,0,0,4,0},
{3,0,0,0,0,0,1},
{2,1,0,0,0,0,1},
{1,2,0,0,0,0,1},
{0,3,0,0,0,0,1},
{2,0,1,0,0,0,1},
{1,1,1,0,0,0,1},
{0,2,1,0,0,0,1},
{1,0,2,0,0,0,1},
{0,1,2,0,0,0,1},
{0,0,3,0,0,0,1},
{2,0,0,1,0,0,1},
{1,1,0,1,0,0,1},
{0,2,0,1,0,0,1},
{1,0,1,1,0,0,1},
{0,1,1,1,0,0,1},
{0,0,2,1,0,0,1},
{1,0,0,2,0,0,1},
{0,1,0,2,0,0,1},
{0,0,1,2,0,0,1},
{0,0,0,3,0,0,1},
{2,0,0,0,1,0,1},
{1,1,0,0,1,0,1},
{0,2,0,0,1,0,1},
{1,0,1,0,1,0,1},
{0,1,1,0,1,0,1},
{0,0,2,0,1,0,1},
{1,0,0,1,1,0,1},
{0,1,0,1,1,0,1},
{0,0,1,1,1,0,1},
{0,0,0,2,1,0,1},
{1,0,0,0,2,0,1},
{0,1,0,0,2,0,1},
{0,0,1,0,2,0,1},
{0,0,0,1,2,0,1},
{0,0,0,0,3,0,1},
{2,0,0,0,0,1,1},
{1,1,0,0,0,1,1},
{0,2,0,0,0,1,1},
{1,0,1,0,0,1,1},
{0,1,1,0,0,1,1},
{0,0,2,0,0,1,1},
{1,0,0,1,0,1,1},
{0,1,0,1,0,1,1},
{0,0,1,1,0,1,1},
{0,0,0,2,0,1,1},
{1,0,0,0,1,1,1},
{0,1,0,0,1,1,1},
{0,0,1,0,1,1,1},
{0,0,0,1,1,1,1},
{0,0,0,0,2,1,1},
{1,0,0,0,0,2,1},
{0,1,0,0,0,2,1},
{0,0,1,0,0,2,1},
{0,0,0,1,0,2,1},
{0,0,0,0,1,2,1},
{0,0,0,0,0,3,1},
{2,0,0,0,0,0,2},
{1,1,0,0,0,0,2},
{0,2,0,0,0,0,2},
{1,0,1,0,0,0,2},
{0,1,1,0,0,0,2},
{0,0,2,0,0,0,2},
{1,0,0,1,0,0,2},
{0,1,0,1,0,0,2},
{0,0,1,1,0,0,2},
{0,0,0,2,0,0,2},
{1,0,0,0,1,0,2},
{0,1,0,0,1,0,2},
{0,0,1,0,1,0,2},
{0,0,0,1,1,0,2},
{0,0,0,0,2,0,2},
{1,0,0,0,0,1,2},
{0,1,0,0,0,1,2},
{0,0,1,0,0,1,2},
{0,0,0,1,0,1,2},
{0,0,0,0,1,1,2},
{0,0,0,0,0,2,2},
{1,0,0,0,0,0,3},
{0,1,0,0,0,0,3},
{0,0,1,0,0,0,3},
{0,0,0,1,0,0,3},
{0,0,0,0,1,0,3},
{0,0,0,0,0,1,3},
{0,0,0,0,0,0,4}};

  nnn = pwrsix[m+1];
  Apwr = (int *)malloc(nnn*sizeof(int));
  Bpwr = (int *)malloc(nnn*sizeof(int));
  Cpwr = (int *)malloc(nnn*sizeof(int));
  Dpwr = (int *)malloc(nnn*sizeof(int));
  Epwr = (int *)malloc(nnn*sizeof(int));
  Fpwr = (int *)malloc(nnn*sizeof(int));
  Gpwr = (int *)malloc(nnn*sizeof(int));

  if(debug) printf("terms used to find fit \n");
  for(i=0; i<=m; i++) // power requested                                  
  {
    for(ii=pwrsix[i]; ii<pwrsix[i+1]; ii++)
    {
      Apwr[ii]=pwrs[ii][0];
      Bpwr[ii]=pwrs[ii][1];
      Cpwr[ii]=pwrs[ii][2];
      Dpwr[ii]=pwrs[ii][3];
      Epwr[ii]=pwrs[ii][4];
      Fpwr[ii]=pwrs[ii][5];
      Gpwr[ii]=pwrs[ii][6];
      if(debug) printf("%d   a^%d * b^%d * c^%d * d^%d * e^%d * f^%d * g^%d\n",
		       ii, Apwr[ii],Bpwr[ii],Cpwr[ii],Dpwr[ii],
                       Epwr[ii],Fpwr[ii],Gpwr[ii]);
    }
    if(debug) printf(" \n");
  }
} /* end gen_7d_powers */


void fit1_data(double x, double u, int last)
{
  double term_i, term_j;
  int i, k, ia, ja;

  Av[0] = 1.0;
  for(k=1; k<=m; k++)
  {
    Av[k] = Av[k-1]*x;
    if(last) printf("Av[%d]=%f \n", k, Av[k]);
  }
  for(ia=0; ia<nnn; ia++)
  {
    term_i = Av[Apwr[ia]];
    for(ja=0; ja<nnn; ja++)
    {
      term_j = Av[Apwr[ja]];
      A[ia*nnn+ja] = A[ia*nnn+ja] + term_i * term_j;
    } // ja
    Y[ia] = Y[ia] + u * term_i;
  } // ia
  if(last)
  {
    printf("A[first]=%f, Y[first]=%f, A[last,last]=%f, Y[last]=%f \n",
           A[0], Y[0], A[(nnn-1)*nnn+(nnn-1)], Y[nnn-1]);
    simeq(nnn, C);
    if(debug)
      for(i=0; i<nnn; i++)
      {
	printf("C[%d]=%f x^%d \n",i,C[i],Apwr[i]);
      }
  }
} /* end fit1_data */

void fit2_data(double x, double y, double u, int last)
{
  double term_i, term_j;
  int i, k, ia, ja;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  for(k=1; k<=m; k++)
  {
    Av[k] = Av[k-1]*x;
    Bv[k] = Bv[k-1]*y;
    if(last) printf("Av[%d]=%f, Bv[%d]=%f \n", k, Av[k], k, Bv[k]);
  }
  for(ia=0; ia<nnn; ia++)
  {
    term_i = Av[Apwr[ia]] * Bv[Bpwr[ia]];
    for(ja=0; ja<nnn; ja++)
    {
      term_j = Av[Apwr[ja]] * Bv[Bpwr[ja]];
      A[ia*nnn+ja] = A[ia*nnn+ja] + term_i * term_j;
    } // ja
    Y[ia] = Y[ia] + u * term_i;
  } // ia
  if(last)
  {
    printf("A[first]=%f, Y[first]=%f, A[last,last]=%f, Y[last]=%f \n",
           A[0], Y[0], A[(nnn-1)*nnn+(nnn-1)], Y[nnn-1]);
    simeq(nnn, C);
    if(debug)
      for(i=0; i<nnn; i++)
      {
	printf("C[%d]=%f x^%d y^%d\n",i,C[i],Apwr[i],Bpwr[i]);
      }
  }
} /* end fit2_data */

void fit3_data(double x, double y, double z, double u, int last)
{
  double term_i, term_j;
  int i, k, ia, ja;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  for(k=1; k<=m; k++)
  {
    Av[k] = Av[k-1]*x;
    Bv[k] = Bv[k-1]*y;
    Cv[k] = Cv[k-1]*z;
  }
  for(ia=0; ia<nnn; ia++)
  {
    term_i = Av[Apwr[ia]] * Bv[Bpwr[ia]] * Cv[Cpwr[ia]];
    for(ja=0; ja<nnn; ja++)
    {
      term_j = Av[Apwr[ja]] * Bv[Bpwr[ja]] * Cv[Cpwr[ja]];
      A[ia*nnn+ja] = A[ia*nnn+ja] + term_i * term_j;
    } // ja                                                                   
    Y[ia] = Y[ia] + u * term_i;
  } // ia                                                                     
  if(last)
  {
    simeq(nnn, C);
    if(debug)
      for(i=0; i<nnn; i++)
      {
	printf("C[%d]=%f x^%d y^%d z^%d\n",i, C[i], Apwr[i], Bpwr[i], Cpwr[i]);
      }
  }
} /* end fit3_data */

void fit4_data(double x, double y, double z, double t, double u, int last)
{
  double term_i, term_j;
  int i, k, ia, ja;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  Dv[0] = 1.0;
  for(k=1; k<=m; k++)
  {
    Av[k] = Av[k-1]*x;
    Bv[k] = Bv[k-1]*y;
    Cv[k] = Cv[k-1]*z;
    Dv[k] = Dv[k-1]*t;
  }
  for(ia=0; ia<nnn; ia++)
  {
    term_i = Av[Apwr[ia]] * Bv[Bpwr[ia]] * Cv[Cpwr[ia]] * Dv[Dpwr[ia]];
    for(ja=0; ja<nnn; ja++)
    {
      term_j = Av[Apwr[ja]] * Bv[Bpwr[ja]] * Cv[Cpwr[ja]] * Dv[Dpwr[ja]];
      A[ia*nnn+ja] = A[ia*nnn+ja] + term_i * term_j;
    } // ja                                                                   
    Y[ia] = Y[ia] + u * term_i;
  } // ia                                                                     
  if(last)
  {
    simeq(nnn, C);
    if(debug)
      for(i=0; i<nnn; i++)
      {
	printf("C[%d]=%f  x^%d * y^%d * z^%d * t^%d\n",
               i,C[i],Apwr[i],Bpwr[i],Cpwr[i],Dpwr[i]);
      }
  }
} /* end fit4_data */

void fit5_data(double x, double y, double z, double t, double u, 
               double w, int last)
{
  double term_i, term_j;
  int i, k, ia, ja;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  Dv[0] = 1.0;
  Ev[0] = 1.0;
  for(k=1; k<=m; k++)
  {
    Av[k] = Av[k-1]*x;
    Bv[k] = Bv[k-1]*y;
    Cv[k] = Cv[k-1]*z;
    Dv[k] = Dv[k-1]*t;
    Ev[k] = Ev[k-1]*u;
  }
  for(ia=0; ia<nnn; ia++)
  {
    term_i = Av[Apwr[ia]] * Bv[Bpwr[ia]] * Cv[Cpwr[ia]] * Dv[Dpwr[ia]] *
             Ev[Epwr[ia]];
    for(ja=0; ja<nnn; ja++)
    {
      term_j = Av[Apwr[ja]] * Bv[Bpwr[ja]] * Cv[Cpwr[ja]] * Dv[Dpwr[ja]] *
               Ev[Epwr[ja]];
      A[ia*nnn+ja] = A[ia*nnn+ja] + term_i * term_j;
    } // ja                                                                   
    Y[ia] = Y[ia] + w * term_i;
  } // ia                                                                     
  if(last)
  {
    simeq(nnn, C);
    if(debug)
      for(i=0; i<nnn; i++)
      {
	printf("C[%d]=%f  x^%d * y^%d * z^%d * t^%d * u^%d \n",
               i,C[i],Apwr[i],Bpwr[i],Cpwr[i],Dpwr[i],Epwr[i]);
      }
  }
} /* end fit5_data */


void fit6_data(double x, double y, double z, double t, double u, 
               double v, double w, int last)
{
  double term_i, term_j;
  int i, k, ia, ja;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  Dv[0] = 1.0;
  Ev[0] = 1.0;
  Fv[0] = 1.0;
  for(k=1; k<=m; k++)
  {
    Av[k] = Av[k-1]*x;
    Bv[k] = Bv[k-1]*y;
    Cv[k] = Cv[k-1]*z;
    Dv[k] = Dv[k-1]*t;
    Ev[k] = Ev[k-1]*u;
    Fv[k] = Fv[k-1]*v;
  }
  for(ia=0; ia<nnn; ia++)
  {
    term_i = Av[Apwr[ia]] * Bv[Bpwr[ia]] * Cv[Cpwr[ia]] * Dv[Dpwr[ia]] *
             Ev[Epwr[ia]] * Fv[Fpwr[ia]];
    for(ja=0; ja<nnn; ja++)
    {
      term_j = Av[Apwr[ja]] * Bv[Bpwr[ja]] * Cv[Cpwr[ja]] * Dv[Dpwr[ja]] *
               Ev[Epwr[ja]] * Fv[Fpwr[ja]];
      A[ia*nnn+ja] = A[ia*nnn+ja] + term_i * term_j;
    } // ja                                                                   
    Y[ia] = Y[ia] + w * term_i;
  } // ia                                                                     
  if(last)
  {
    simeq(nnn, C);
    if(debug)
      for(i=0; i<nnn; i++)
      {
	printf("C[%d]=%f  x^%d * y^%d * z^%d * t^%d * u^%d * v^%d \n",
               i,C[i],Apwr[i],Bpwr[i],Cpwr[i],Dpwr[i],Epwr[i],Fpwr[i]);
      }
  }
} /* end fit6_data */

void fit7_data(double x, double y, double z, double t, double u, 
               double v, double w, double val, int last)
{
  double term_i, term_j;
  int i, k, ia, ja;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  Dv[0] = 1.0;
  Ev[0] = 1.0;
  Fv[0] = 1.0;
  Gv[0] = 1.0;
  for(k=1; k<=m; k++)
  {
    Av[k] = Av[k-1]*x;
    Bv[k] = Bv[k-1]*y;
    Cv[k] = Cv[k-1]*z;
    Dv[k] = Dv[k-1]*t;
    Ev[k] = Ev[k-1]*u;
    Fv[k] = Fv[k-1]*v;
    Gv[k] = Gv[k-1]*w;
  }
  for(ia=0; ia<nnn; ia++)
  {
    term_i = Av[Apwr[ia]] * Bv[Bpwr[ia]] * Cv[Cpwr[ia]] * Dv[Dpwr[ia]] *
             Ev[Epwr[ia]] * Fv[Fpwr[ia]] * Gv[Gpwr[ia]];
    for(ja=0; ja<nnn; ja++)
    {
      term_j = Av[Apwr[ja]] * Bv[Bpwr[ja]] * Cv[Cpwr[ja]] * Dv[Dpwr[ja]] *
               Ev[Epwr[ja]] * Fv[Fpwr[ja]] * Gv[Gpwr[ja]];
      A[ia*nnn+ja] = A[ia*nnn+ja] + term_i * term_j;
    } // ja                                                                   
    Y[ia] = Y[ia] + val * term_i;
  } // ia                                                                     
  if(last)
  {
    simeq(nnn, C);
    if(debug)
      for(i=0; i<nnn; i++)
      {
	printf("C[%d]=%f  x^%d * y^%d * z^%d * t^%d * u^%d * v^%d w^%d\n",
               i,C[i],Apwr[i],Bpwr[i],Cpwr[i],Dpwr[i],Epwr[i],Fpwr[i], Gpwr[i]);
      }
  }
} /* end fit7_data */


static void fit_2d(int nx, int ny, double x[], double y[], double u[])
{
  int i, j, ii, jj, k;
  double Av[15]; /* at least n */
  double Bv[15]; /* powers of variables */
  double u1, a1, b1;
  double term_i, term_j;
  
  if(debug) printf("fit_2d n=%d, m=%d, nnn=%d \n", n, m, nnn);
  gen_2d_powers();
  if(debug) printf("after gen_2d_powers n=%d, m=%d, nnn=%d \n", n, m, nnn);
  
  for(i=0; i<nnn; i++)
  {
    for(j=0; j<nnn; j++)
    {
      A[i*nnn+j] = 0.0;
    }
    Y[i] = 0.0;
  }
  if(debug) printf("A initialized\n");

  for(i=0; i<nx; i++)
  {
    for(j=0; j<ny; j++)
    {
      a1 = x[i];
      b1 = y[j];
      u1 = u[i*ny+j];
      Av[0] = 1.0;
      Bv[0] = 1.0;
      if(debug) printf("Av[0]=%g, Bv[0]=%g \n",
                        Av[0], Bv[0]);
      for(k=1; k<=n; k++)
      {
        Av[k] = Av[k-1]*a1;
        Bv[k] = Bv[k-1]*b1;
        if(debug) printf("Av[%d]=%g, Bv[%d]=%g \n",
                          k, Av[k], k, Bv[k]);
      } /* end k */
      for(ii=0; ii<nnn; ii++)
      {
        term_i = Av[Apwr[ii]] * Bv[Bpwr[ii]];
        printf("ii=%d, term_i=%g \n",ii,term_i);
        for(jj=0; jj<nnn; jj++)
        {
          term_j = Av[Apwr[jj]] * Bv[Bpwr[jj]];
          printf("jj=%d, term_j=%g \n",jj,term_j);
          A[ii*nnn+jj] = A[ii*nnn+jj] + term_i * term_j;
        } /* end jj */
        Y[ii] = Y[ii] + u1 * term_i;
      } /* end ii */
    } /* end y[j] data loop */
  } /* end x[i] data loop */
  if(debug)
  {
    for(i=0; i<nnn; i++)
    {
      for(j=0; j<nnn; j++)
      {
        printf("A[%d][%d]=%g \n", i, j, A[i*nnn+j]);
      }
      printf("Y[%d]=%g \n", i, Y[i]);
    }
  }
  if(debug) printf("solve for C \n");
  simeq(nnn, C);
} /* end fit_2d */

void fit2_bpoly(int nx, int ny, double x[], double y[], double u[])
{
  double rms_err, avg_err, max_err;

  fit_2d(nx, ny, x, y, u); /* C coefficients computed */
  check_2d(nx, ny, x, y, u, &rms_err, &avg_err, &max_err);
} /* end fit2_bpoly */

double eval2_poly(double x, double y)
{ 
  double u;
  int i;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  for(i=1; i<=m; i++)
  {
    Av[i] = Av[i-1]*x;
    Bv[i] = Bv[i-1]*y;
  }
  u = 0.0;
  for(i=0; i<nnn; i++)
  {
    u = u + Av[Apwr[i]] * Bv[Bpwr[i]] * C[i];
  }
  return u;
} /* end eval2_poly */

double eval1_poly(double x)
{ 
  double u;
  int i;

  Av[0] = 1.0;
  for(i=1; i<=m; i++)
  {
    Av[i] = Av[i-1]*x;
  }
  u = 0.0;
  for(i=0; i<nnn; i++)
  {
    u = u + Av[Apwr[i]] * C[i];
  }
  return u;
} /* end eval1_poly */

double eval3_poly(double x, double y, double z)
{ 
  double u;
  int i;
  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  for(i=1; i<=m; i++)
  {
    Av[i] = Av[i-1]*x;
    Bv[i] = Bv[i-1]*y;
    Cv[i] = Cv[i-1]*z;
  }
  u = 0.0;
  for(i=0; i<nnn; i++)
  {
    u = u + C[i] * Av[Apwr[i]] * Bv[Bpwr[i]] * Cv[Cpwr[i]];
  }
  return u;
} /* end eval3_poly */

double eval4_poly(double x, double y, double z, double t)
{ 
  double u;
  int i;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  Dv[0] = 1.0;
  for(i=1; i<=m; i++)
  {
    Av[i] = Av[i-1]*x;
    Bv[i] = Bv[i-1]*y;
    Cv[i] = Cv[i-1]*z;
    Dv[i] = Dv[i-1]*t;
  }
  u = 0.0;
  for(i=0; i<nnn; i++)
  {
    u = u + C[i] * Av[Apwr[i]] * Bv[Bpwr[i]] * Cv[Cpwr[i]] * Dv[Dpwr[i]];
  }
  return u;
} /* end eval4_poly */

double eval5_poly(double x, double y, double z, double t, double u)
{ 
  double w;
  int i;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  Dv[0] = 1.0;
  Ev[0] = 1.0;
  for(i=1; i<=m; i++)
  {
    Av[i] = Av[i-1]*x;
    Bv[i] = Bv[i-1]*y;
    Cv[i] = Cv[i-1]*z;
    Dv[i] = Dv[i-1]*t;
    Ev[i] = Ev[i-1]*u;
  }
  w = 0.0;
  for(i=0; i<nnn; i++)
  {
    w = w + C[i] * Av[Apwr[i]] * Bv[Bpwr[i]] * Cv[Cpwr[i]] * Dv[Dpwr[i]] *
                   Ev[Epwr[i]];
  }
  return w;
} /* end eval5_poly */

double eval6_poly(double x, double y, double z, double t, double u, double v)
{ 
  double w;
  int i;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  Dv[0] = 1.0;
  Ev[0] = 1.0;
  Fv[0] = 1.0;
  for(i=1; i<=m; i++)
  {
    Av[i] = Av[i-1]*x;
    Bv[i] = Bv[i-1]*y;
    Cv[i] = Cv[i-1]*z;
    Dv[i] = Dv[i-1]*t;
    Ev[i] = Ev[i-1]*u;
    Fv[i] = Fv[i-1]*v;
  }
  w = 0.0;
  for(i=0; i<nnn; i++)
  {
    w = w + C[i] * Av[Apwr[i]] * Bv[Bpwr[i]] * Cv[Cpwr[i]] * Dv[Dpwr[i]] *
                   Ev[Epwr[i]] * Fv[Fpwr[i]];
  }
  return w;
} /* end eval6_poly */


double eval7_poly(double x, double y, double z, double t, double u, double v, double w)
{ 
  double val;
  int i;

  Av[0] = 1.0;
  Bv[0] = 1.0;
  Cv[0] = 1.0;
  Dv[0] = 1.0;
  Ev[0] = 1.0;
  Fv[0] = 1.0;
  Gv[0] = 1.0;
  for(i=1; i<=m; i++)
  {
    Av[i] = Av[i-1]*x;
    Bv[i] = Bv[i-1]*y;
    Cv[i] = Cv[i-1]*z;
    Dv[i] = Dv[i-1]*t;
    Ev[i] = Ev[i-1]*u;
    Fv[i] = Fv[i-1]*v;
    Gv[i] = Gv[i-1]*w;
  }
  val = 0.0;
  for(i=0; i<nnn; i++)
  {
    val = val + C[i] * Av[Apwr[i]] * Bv[Bpwr[i]] * Cv[Cpwr[i]] * Dv[Dpwr[i]] *
                   Ev[Epwr[i]] * Fv[Fpwr[i]] * Gv[Gpwr[i]];
  }
  return val;
} /* end eval7_poly */

static void check_2d(int nx, int ny, double x[], double y[], double u[], 
                     double *rms_err, double *avg_err, double *max_err)
{
  double a1, b1, u1, ua, diff;
  double sumsq = 0.0;
  double sum = 0.0;
  double maxe = 0.0;
  double amin, amax, bmin, bmax, umin, umax;
  int i, j, k;
  
  for(i=0; i<nx; i++)
  {
    for(j=0; j<ny; j++)
    {
      if(i==0 && j==0)
      {
        u1 = eval2_poly(x[i], y[j]);
        amin=x[i]; amax=x[i];
        bmin=y[j]; bmax=y[j];
        umin=u[i*ny+j]; umax=u[i*ny+j];
      }
      if(x[i]>amax) amax=x[i];
      if(a1<amin) amin=x[i];
      if(b1>bmax) bmax=y[j];
      if(b1<bmin) bmin=y[j];
      if(u1>umax) umax=u1;
      if(u1<umin) umin=u1;

      ua = u[i*ny+j];
      diff = abs(u1-ua);
      if(debug) printf("sample=%d, u_fit=%g, u_data=%g, diff=%g \n",
                     i*nnn+j, u1, ua, diff);
      if(diff>maxe) maxe=diff;
      sum = sum + diff;
      sumsq = sumsq + diff*diff;
    } /* end j */
  } /* end i */
  printf("check_2d k=%d, amin=%g, amax=%g, bmin=%g, bmax=%g \n",
          k, amin, amax, bmin, bmax);
  printf("               umin=%g, umax=%g \n",
             umin, umax);
  *max_err = maxe;
  *avg_err = sum/(double)(nnn);
  *rms_err = sqrt(sumsq/(double)(nnn));
} /* end check_2d */

static void simeq(int nn, double X[]) /* A and Y global */
{

/*      PURPOSE : SOLVE THE LINEAR SYSTEM OF EQUATIONS WITH REAL     */
/*                COEFFICIENTS   [A] * |X| = |Y|                     */
/*                                                                   */
/*      INPUT  : THE NUMBER OF EQUATIONS  n                          */
/*               THE REAL MATRIX  A   should be A[i][j] but A[i*n+j] */
/*               THE REAL VECTOR  Y                                  */
/*      OUTPUT : THE REAL VECTOR  X                                  */
/*                                                                   */
/*      METHOD : GAUSS-JORDAN ELIMINATION USING MAXIMUM ELEMENT      */
/*               FOR PIVOT.                                          */
/*                                                                   */
/*      USAGE  :     simeq(n,X); A and Y global here                 */
/*                                                                   */
/*                                                                   */
/*    WRITTEN BY : JON SQUIRE , 28 MAY 1983                          */
/*    ORIGONAL DEC 1959 for IBM 650, TRANSLATED TO OTHER LANGUAGES   */
/*    e.g. FORTRAN converted to Ada converted to C                   */

    double *B;           /* [nn][nn+1]  WORKING MATRIX */
    int *ROW;            /* ROW INTERCHANGE INDICIES */
    int HOLD , I_PIVOT;  /* PIVOT INDICIES */
    double PIVOT;        /* PIVOT ELEMENT VALUE */
    double ABS_PIVOT;
    int i,j,k,m;

    B = (double *)calloc(nn*(nn+1), sizeof(double));
    ROW = (int *)calloc(nn, sizeof(int));
    m = nn+1;

    /* BUILD WORKING DATA STRUCTURE */
    for(i=0; i<nn; i++){
      for(j=0; j<nn; j++){
        B[i*m+j] = A[i*nn+j];
      }
      B[i*m+nn] = Y[i];
    }
    /* SET UP ROW  INTERCHANGE VECTORS */
    for(k=0; k<nn; k++){
      ROW[k] = k;
    }

    /* BEGIN MAIN REDUCTION LOOP */
    for(k=0; k<nn; k++){

      /* FIND LARGEST ELEMENT FOR PIVOT */
      PIVOT = B[ROW[k]*m+k];
      ABS_PIVOT = abs(PIVOT);
      I_PIVOT = k;
      for(i=k; i<nn; i++){
        if( abs(B[ROW[i]*m+k]) > ABS_PIVOT){
          I_PIVOT = i;
          PIVOT = B[ROW[i]*m+k];
          ABS_PIVOT = abs ( PIVOT );
        }
      }

      /* HAVE PIVOT, INTERCHANGE ROW POINTERS */
      HOLD = ROW[k];
      ROW[k] = ROW[I_PIVOT];
      ROW[I_PIVOT] = HOLD;

      /* CHECK FOR NEAR SINGULAR */
      if( ABS_PIVOT < 1.0E-12 ){
        for(j=k+1; j<nn+1; j++){
          B[ROW[k]*m+j] = 0.0;
        }
        printf("redundant row (singular) %d \n", ROW[k]);
      } /* singular, delete row */
      else{

        /* REDUCE ABOUT PIVOT */
        for(j=k+1; j<nn+1; j++){
          B[ROW[k]*m+j] = B[ROW[k]*m+j] / B[ROW[k]*m+k];
        }

        /* INNER REDUCTION LOOP */
        for(i=0; i<nn; i++){
          if( i != k){
            for(j=k+1; j<nn+1; j++){
              B[ROW[i]*m+j] = B[ROW[i]*m+j] - B[ROW[i]*m+k] * B[ROW[k]*m+j];
            }
          }
        }
      }
      /* FINISHED INNER REDUCTION */
    }

    /* END OF MAIN REDUCTION LOOP */
    /* BUILD  X  FOR RETURN, UNSCRAMBLING ROWS */
    for(i=0; i<nn; i++){
      X[i] = B[ROW[i]*m+nn];
    }
    free(B);
    free(ROW);
} /* end simeq */

/* end lsfit.c */

int get_fitnnn()
{
  return nnn; // user allocates CC[nnn] coefficients, pwr[n*nnn]
} /* end get_fitnnn */

void get_fit(double CC[], int pwr[])
{
  int i, j;
  for(i=0; i<nnn; i++)
  {
    CC[i] = C[i];                 // coefficient of term
    pwr[i*n+0] = Apwr[i];         // power of first variable
    if(n>1) pwr[i*n+1] = Bpwr[i]; // power of second variable
    if(n>2) pwr[i*n+2] = Cpwr[i];
    if(n>3) pwr[i*n+3] = Dpwr[i];
    if(n>4) pwr[i*n+4] = Epwr[i];
    if(n>5) pwr[i*n+5] = Fpwr[i];
    if(n>6) pwr[i*n+6] = Gpwr[i];
  }
} /* end get_fit */

void write_fit7(char* fname) // .c will ne added
{
  FILE * outp;
  int i, j, nn; // non zero coefficients
  char dest[100];

  strcpy(dest, fname);
  outp = fopen(strcat(dest,".c"), "w");
  if(outp==NULL) 
  {
    printf("file %s could not be opened for write\n", fname);
    return;
  }
  fprintf(outp, "double %s(double a, double b, double c,\n", fname);
  fprintf(outp, "       double d, double e, double f, double g) \n");
  fprintf(outp, "{ \n");
  fprintf(outp, "  int i, j; \n");
  fprintf(outp, "  int m = %d; \n", m);
  fprintf(outp, "  double t, val=0.0; \n");
  fprintf(outp, "  double vpwr[%d]; \n", nnn*n); // variable powers computed
  fprintf(outp, "  double C[%d] = { \n", nnn); // coefficients
  fprintf(outp, "    ");
  for(i=0; i<nnn-1; i++) // write C coefficients
  {
    fprintf(outp, " %19.16e,", C[i]);
    if(i%3==0)
    {
      fprintf(outp, " \n");
      fprintf(outp, "    ");
    }
  }
  fprintf(outp, " %19.15e}; \n", C[nnn-1]);
  fprintf(outp, " \n");
  fprintf(outp, "  int pwrs[%d] = { \n", nnn*n);
  for(i=0; i<nnn-1; i++) // write powers of variables
  {
    fprintf(outp, "    %d, %d, %d, %d, %d, %d, %d, \n", Apwr[i], Bpwr[i],
	    Cpwr[i], Dpwr[i], Epwr[i], Fpwr[i], Gpwr[i]);
  }
  i = nnn-1;
  fprintf(outp, "    %d, %d, %d, %d, %d, %d, %d}; \n", Apwr[i], Bpwr[i],
	    Cpwr[i], Dpwr[i], Epwr[i], Fpwr[i], Gpwr[i]);

  fprintf(outp, " \n");
  
  fprintf(outp, "  for(j=0; j<%d; j++) \n", n);
  fprintf(outp, "    vpwr[j] = 1.0; \n");
  fprintf(outp, "  for(i=1; i<=%d; i++) \n", m);
  fprintf(outp, "  { \n");
  fprintf(outp, "    vpwr[i*%d+0] = vpwr[(i-1)*%d+0] * a; \n", n, n);
  fprintf(outp, "    vpwr[i*%d+1] = vpwr[(i-1)*%d+1] * b; \n", n, n);
  fprintf(outp, "    vpwr[i*%d+2] = vpwr[(i-1)*%d+2] * c; \n", n, n);
  fprintf(outp, "    vpwr[i*%d+3] = vpwr[(i-1)*%d+3] * d; \n", n, n);
  fprintf(outp, "    vpwr[i*%d+4] = vpwr[(i-1)*%d+4] * e; \n", n, n);
  fprintf(outp, "    vpwr[i*%d+5] = vpwr[(i-1)*%d+5] * f; \n", n, n);
  fprintf(outp, "    vpwr[i*%d+6] = vpwr[(i-1)*%d+6] * g; \n", n, n);
  fprintf(outp, "  } \n");
  fprintf(outp, "  for(i=0; i<%d; i++) \n", nnn);
  fprintf(outp, "  { \n");
  fprintf(outp, "    t = C[i]; \n");
  fprintf(outp, "    for(j=0; j<%d; j++) \n", n);
  fprintf(outp, "    { \n");
  fprintf(outp, "      t = t * vpwr[pwrs[i*%d+j]*%d+j]; \n", n, n);
  fprintf(outp, "    } \n");
  fprintf(outp, "    val = val + t; \n");
  fprintf(outp, "  } \n");
  fprintf(outp, "  return val; \n");
  fprintf(outp, "} \n");
  fclose(outp);
} /* end write_fit7 */