[Bug-gsl] bug in akima interpolation

[ABX]

I tried to implement akima interpolation in splines of POV-Ray. First I compared
how it looks with code from other sources. When I linked my application with
modules from GSL library I recived shape with bug in the most left point:

http://news.povray.org/povray.binaries.images/31227/221922/akima_spline.png

but when I used code from http://www.magic-software.com/Interpolation1D.html
I received almost similiar shape except all control points are "smooth"

http://news.povray.org/povray.binaries.images/31227/222064/akima_spline.png

Since both libraries are builded different way (temporary and helpers variables
are different) I not tried to compare details of implementations. But I suppose
the difference is located here:

GSL:

    m[-2] = 3.0 * m[0] - 2.0 * m[1];
    m[-1] = 2.0 * m[0] - m[1];
    m[size - 1] = 2.0 * m[size - 2] - m[size - 3];
    m[size] = 3.0 * m[size - 2] - 2.0 * m[size - 3];
  
Magic Software:

    afSlope[1] = 2.0f*afSlope[2] - afSlope[3];
    afSlope[0] = 2.0f*afSlope[1] - afSlope[2];
    afSlope[iQuantity+1] = 2.0f*afSlope[iQuantity] - afSlope[iQuantity-1];
    afSlope[iQuantity+2] = 2.0f*afSlope[iQuantity+1] - afSlope[iQuantity];

In my example splines were open, thought first and last point are equal. The
images were created with following set of data:

  t,<x(t),y(t)>

  0,<0,0>
  1,<0,2>
  2,<2,1>
  3,<0,0>
  4,<-1,2>
  5,<-2,1>
  6,<-1,-1>
  7,<0,0>

I do not have original code I tested this cases and even then I could not
deliver it since I linked it with not released code of next POV. But I imagine
to recreate this bug with current gsl code you can take gsl/interpolation/test.c
file and replace two functions with something like:

static int
test_akima (void)
{
  int s;

  double data_t[8] = { 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0 };
  double data_x[8] = { 0.0, 0.0, 2.0, 0.0,-1.0,-2.0,-1.0, 0.0 };
  double data_y[8] = { 0.0, 2.0, 1.0, 0.0, 2.0, 1.0,-1.0, 0.0 };

  xy_table data_table_x = make_xy_table(data_t, data_x, 8);
  xy_table data_table_y = make_xy_table(data_t, data_y, 8);

  s = test_interp (&data_table_x, &data_table_y, gsl_interp_akima);
  gsl_test (s, "akima interpolation");
  return s;
}

static int
test_interp (
  const xy_table * data_table_x,
  const xy_table * data_table_y,
  const gsl_interp_type * T
  )
{
  int status = 0;
  size_t i;

  gsl_interp_accel *a_x = gsl_interp_accel_alloc ();
  gsl_interp_accel *a_y = gsl_interp_accel_alloc ();
  gsl_interp *interp_x = gsl_interp_alloc (T, data_table_x->n);
  gsl_interp *interp_y = gsl_interp_alloc (T, data_table_y->n);
  gsl_interp_init (interp_x, data_table_x->x, data_table_x->y, data_table_x->n);
  gsl_interp_init (interp_y, data_table_y->x, data_table_y->y, data_table_y->n);

  for (DBL i; i<data_table_x->n; i=i+step)  // step has to be defined earlier
    {
      // loop to get values
      double x;
      double y;
      gsl_interp_eval_e (interp_x, data_table_x->x, data_table_x->y, i,a_x,&x);
      gsl_interp_eval_e (interp_y, data_table_y->x, data_table_y->y, i,a_y,&y);
      // output vector f(i)=<x(i),y(i)> in your favorite way
      status++;
    }

  gsl_interp_accel_free (a_x);
  gsl_interp_accel_free (a_y);
  gsl_interp_free (interp_x);
  gsl_interp_free (interp_y);

  return status;
}