a gsl_interp_poly type candidate submited [code maturity is low]

["Dan, Ho-Jin" <hjdan at sys713 dot kaist dot ac dot kr>]

Hi there,
I couldn't seek the global polynomial interpolation  for distinct n 
points similar
to POLINT in the slatec library. The global polynomial interpolation 
generally
is not good for interpolating data since its error behavior is not so 
good. But
it is frequently used as a basic building block in other routines. (I 
need also
another for complex data set). I have to submit the low mature code in 
order to
confirm  my work.

1. I am a novice in developing a free software, so I confuse how to 
submit codes,
    patches and manuals(I can handle a document of the latex format).
2. in the vim editor, how to set the gnu-recommended c-indenting.
3. If a submitted code doesn't seem to be useful, please suggest another 
proxy
4. Is the name of "polynomial" type good or too general?
5. One of Higher-level interfaces will generate coefs. of a polynomial for
    gsl_poly_eval(),  gsl_poly_solve() etc.
6. complex counterparts will be submited.

note. GNU GPL lisence will be attached.



#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_interp.h>

/* Note 1: gsl_interp_polynomial type added 
 * 	is polynomial name(analogous to slatec) right? 
 * 	another candidate is a divdif (divided difference) 
 * Note 2: This code is under development
 * Note 3: ref: K. E. Atkinson, An introduction to Numerical Analysis 2nd ed., 
 * 	John Wiley & Sons, Inc., Singapore, pp. 138--147
 * Dan, Ho-Jin ( mailto:hjdan@sys713.kaist.ac.kr )
 * Nov. 26, 2001
 */

/* for ctags ref. 
 * gsl_interp_alloc
 * gsl_spline_alloc
 * gsl_interp
 * cspline_state_t
 * gsl_interp_cspline
 */

typedef struct {
  double* d;
} polynomial_state_t;

void* polynomial_alloc(size_t size)
{
  polynomial_state_t* state = (polynomial_state_t*)malloc(
      sizeof(polynomial_state_t));
  if (state == 0) {
    GSL_ERROR_NULL("failed to allocate space for polynomial_state_t", 
	GSL_ENOMEM);
  }
  state->d = (double*)malloc(sizeof(double)*size);
  if (state->d == 0) {
    free(state);
    GSL_ERROR_NULL("failed to allocate space for polynomial_state_t.d", 
	GSL_ENOMEM);
  }
  return state;
}
    
int polynomial_init(void *vstate, 
    const double xa[], const double ya[], size_t size)
{
  polynomial_state_t* state = (polynomial_state_t*) vstate;
  int i, j;
  /* newton's divide difference */
  state->d[0] = ya[0];
  for (j=size-1; j>=1; j--) {
    state->d[j] = (ya[j] - ya[j-1])/(xa[j] - xa[j-1]);
  }
  for (i=2; i<size; i++) {
    for (j = size-1; j>=i; j--) {
      state->d[j] = (state->d[j] - state->d[j-1])/(xa[j] - xa[j-i]);
    }
  }
  /*
  for (i=0; i<size; i++) 
    printf("[%d]%e\n", i, state->d[i]);
    */
  return GSL_SUCCESS;
}

int polynomial_eval(const void* vstate, 
    const double xa[], const double ya[], size_t size, double x, 
    gsl_interp_accel *acc, double * y)
{
  polynomial_state_t* state = (polynomial_state_t*) vstate;
  int i;
  *y = state->d[size - 1];
  for (i=size-2; i>=0; i--) {
    *y = state->d[i] + (x - xa[i])*(*y);
  }
  return GSL_SUCCESS;
}

void polynomial_free(void* vstate)
{
  polynomial_state_t* state;
  free(state->d);
  free(state);
}

static const gsl_interp_type polynomial_type = 
{
  "polynomial", 
  2,
  &polynomial_alloc, /* alloc, not applicable */
  &polynomial_init,
  &polynomial_eval,
  NULL, /* currenly unavailable */
  NULL, /* currenly unavailable */
  NULL, /* currenly unavailable */
  &polynomial_free, /* free, not applicable */
};

const gsl_interp_type * gsl_interp_polynomial = &polynomial_type;

int main()
{
  int i;
  double xi, yi;
  double x[5] = {2., 2.1, 2.2, 2.3, 2.4};
  double y[5] = {1.414214, 1.449138, 1.483240, 1.516575, 1.549193};

  printf("#m=0,S=2\n");
  for (i=0; i<5; i++) 
  {
    printf("%g %g\n", x[i], y[i]);
  }
  printf("#m=1,S=0\n");
  {
    gsl_interp* divdif = gsl_interp_alloc(gsl_interp_polynomial, 5);
    gsl_interp_init(divdif, x, y, 5);
    for (xi = x[0]; xi < x[4]; xi+=0.01)
    {
      yi = gsl_interp_eval(divdif, x, y, xi, 0);
      printf("%g %g\n", xi, yi);
    }
    gsl_interp_free(divdif);
  }
  return 0;
}