How do I get the Polynomial Interpolation coefficients using gsl_interp?

Question

So I have the code below. It perfectly calculates all the y-points of the polynomial (and prints them to plot with gnuplot), but how do i get the resulting polynomial (1-x² in this case)?

void twoDegreePoly() {
    int n = 3;
    double x[n],y[n];
    printf ("#m=0,S=16\n");
    for (int i=0; i<n ;i++) {
        x[i] = ((double)2*i)/2 -1;
        y[i] = f(x[i]);
        printf ("%g %g\n", x[i], y[i]);
    }
    printf ("#m=1,S=0\n");

    gsl_interp_accel *acc = gsl_interp_accel_alloc ();
    const gsl_interp_type *t = gsl_interp_polynomial;
    gsl_interp* poly = gsl_interp_alloc(t,n);
    gsl_interp_init (poly, x, y,n);
    for (double xi=x[0]; xi<x[n-1]; xi+= 0.01) {
        double yi = gsl_interp_eval (poly, x, y, xi, acc);
        printf ("%g %g\n", xi, yi);
    }
}

Answer 1

After a quick scan over the documentation, it doesn't seem that such a feature is available in the GSL. This could be caused by two reasons: first, getting polynomial coeffcients is special to this interpolation method doesn't fit well into the general design (which can handle arbitrary functions). Second, citing Numerical Recipes:

Please be certain, however, that the coefficients are what you need. Generally, the coefficients of the interpolating polynomial can be determined much less accurately than its value at a desired abscissa. Therefire, it is not a good idea to determine the coefficients only for use in calculating interpolating values. Values thus calculated will not pass exactly through the tabulated points, for example, ...

The reason for this is that in principle, calculating the coefficients involves solving a linear system with a Vandermonde matrix, which is highly ill-conditioned.

Still, Numerical Recipes gives a routine polcoe by which you can obtain the interpolating polynomial. You can find it in chapter 3.5. in the free second edition.

Answer 2

I have done something similar with the Akima's interpolation. First, define the state as GSL do:

typedef struct
{
  double *b;
  double *c;
  double *d;
  double *_m;
}akima_state_t;

Then, create the interpolant

spline = gsl_spline_alloc (gsl_interp_akima, M_size);
gsl_spline_init (spline, x, y, M_size);

and after that, you can do :

const akima_state_t *state = (const akima_state_t *) ( spline -> interp -> state);
  double _b,_c,_d;
  for (int i = 0; i < M_size; i++)
  {
    _b = state->b[i];
    _c = state->c[i];
    _d = state->d[i];
    std::cout << "(x>"<<x[i]<<")*(x<"<<x[i+1]<<")*("<<y[i]<< "+ (x-"<< x[i]<<")*("<<_b<<"+(x-"<< x[i]<<")*("<<_c<<"+"<<_d<<"*(x-"<<x[i]<<")))) + ";
  }

I do not have tried with a polynomial interpolation, but here the state struct for polynomial, it should be a good starting point.

typedef struct
{
  double *d;
  double *coeff;
  double *work;
}
polynomial_state_t;