periodic cspline first derivative discontinuity

[David Necas (Yeti) <yeti at physics dot muni dot cz>]

Hello,

I've found a bug (or I think I've found a bug ;-)
in periodic cubic spline implementation -- its
first derivative is not continuous at the edges of
the first segment.

This is because length of next segment (h_ip1
in the patch below) is computed modulo even for
the last segment and so it results as negative
(well, not much clearly explained, perhaps the
patch itself is clearer).

Please note applying this patch actually fixes
only first derivative discontinuity in right edge
of the first segment, some dicontinuity may still
appear at the left edge due to probems in cyclic
tridiagonal linear system solver -- see my next
bugreport.

Yeti

--- gsl-1.1.1.orig/interpolation/cspline.c	Sun Mar 31 19:32:20 2002
+++ gsl-1.1.1/interpolation/cspline.c	Mon Apr  1 14:52:12 2002
@@ -168,13 +168,25 @@
     return GSL_SUCCESS;
   } else {
     
-    state->offdiag[max_index] =  xa[1]-xa[0];
-    
-    for (i = 0; i < sys_size; i++) {
-      const double h_i   = xa[i + 1] - xa[i];
-      const double h_ip1 = xa[(i + 2) % num_points] - xa[i + 1];
+    for (i = 0; i < sys_size-1; i++) {
+      const double h_i       = xa[i + 1] - xa[i];
+      const double h_ip1     = xa[i + 2] - xa[i + 1];
+      const double ydiff_i   = ya[i + 1] - ya[i];
+      const double ydiff_ip1 = ya[i + 2] - ya[i + 1];
+      state->offdiag[i] = h_ip1;
+      state->diag[i] = 2.0 * (h_ip1 + h_i);
+      state->g[i] = 3.0 * (ydiff_ip1 / h_ip1 - ydiff_i / h_i);
+    }
+    /* fixed derivative discontinuity in xa[1], but not in xa[0] --
+     * this one is due to some tridiag solver problem
+     * FIXME
+     */
+    {
+      i = sys_size - 1;
+      const double h_i       = xa[i + 1] - xa[i];
+      const double h_ip1     = xa[1] - xa[0];
       const double ydiff_i   = ya[i + 1] - ya[i];
-      const double ydiff_ip1 = ya[(i + 2) % num_points] - ya[i + 1];
+      const double ydiff_ip1 = ya[1] - ya[0];
       state->offdiag[i] = h_ip1;
       state->diag[i] = 2.0 * (h_ip1 + h_i);
       state->g[i] = 3.0 * (ydiff_ip1 / h_ip1 - ydiff_i / h_i);