# HG changeset patch
# User Marco Caliari <marco.caliari@univr.it>
# Date 1289065386 -3600
# Node ID d17cb8a1271d00dafa044375d5f00d1797ab3e83
# Parent  8f67fe9dd64eab80ac594655455b57ef01abb511
spline.m: Fix algorithm for input with 3 elements (bug #31098)

diff -r 8f67fe9dd64e -r d17cb8a1271d scripts/ChangeLog
--- a/scripts/ChangeLog	Sat Nov 06 09:08:47 2010 -0400
+++ b/scripts/ChangeLog	Sat Nov 06 18:43:06 2010 +0100
@@ -1,3 +1,8 @@
+2010-11-06  Marco Caliari <marco.caliari@univr.it>
+
+	* polynomial/spline.m: Fit a parabola for input with 3
+	elements (bug #31098).
+
 2010-11-04  Rik  <octave@nomad.inbox5.com>
 
 	* plot/__fltk_ginput__.m: Use semicolons to prevent internal
diff -r 8f67fe9dd64e -r d17cb8a1271d scripts/polynomial/spline.m
--- a/scripts/polynomial/spline.m	Sat Nov 06 09:08:47 2010 -0400
+++ b/scripts/polynomial/spline.m	Sat Nov 06 18:43:06 2010 +0100
@@ -144,48 +144,60 @@
 
   else
 
-    g = zeros (n-2, columns (a));
-    g(1,:) = 3 / (h(1) + h(2)) ...
-        * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:)));
-    g(n-2,:) = 3 / (h(n-1) + h(n-2)) ...
-        * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:)));
+    if (n == 3)
 
-    if (n > 4)
+      n = 2;
+      c = (a(1,:) - a(3,:)) / ((x(3) - x(1)) * (x(2) - x(3))) ...
+          + (a(2,:) - a(1,:)) / ((x(2) - x(1)) * (x(2) - x(3)));
+      b = (a(2,:) - a(1,:)) * (x(3) - x(1)) ...
+          / ((x(2) - x(1)) * (x(3) - x(2))) ...
+          + (a(1,:) - a(3,:)) * (x(2) - x(1)) ...
+          / ((x(3) - x(1)) * (x(3) - x(2)));
+      a = a(1,:);
+      d = zeros (size (a));
+      x = [min(x), max(x)];
+
+    else
 
-      g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ...
-          - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx);
+      g = zeros (n-2, columns (a));
+      g(1,:) = 3 / (h(1) + h(2)) ...
+          * (a(3,:) - a(2,:) - h(2) / h(1) * (a(2,:) - a(1,:)));
+      g(n-2,:) = 3 / (h(n-1) + h(n-2)) ...
+          * (h(n-2) / h(n-1) * (a(n,:) - a(n-1,:)) - (a(n-1,:) - a(n-2,:)));
+
+      if (n > 4)
 
-      dg = 2 * (h(1:n-2) .+ h(2:n-1));
-      dg(1) = dg(1) - h(1);
-      dg(n-2) = dg(n-2) - h(n-1);
+        g(2:n - 3,:) = 3 * diff (a(3:n-1,:)) ./ h(3:n-2,idx) ...
+            - 3 * diff (a(2:n-2,:)) ./ h(2:n - 3,idx);
 
-      ldg = udg = h(2:n-2);
-      udg(1) = udg(1) - h(1);
-      ldg(n - 3) = ldg(n-3) - h(n-1);
-      c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]],
+        dg = 2 * (h(1:n-2) .+ h(2:n-1));
+        dg(1) = dg(1) - h(1);
+        dg(n-2) = dg(n-2) - h(n-1);
+
+        ldg = udg = h(2:n-2);
+        udg(1) = udg(1) - h(1);
+        ldg(n - 3) = ldg(n-3) - h(n-1);
+        c(2:n-1,:) = spdiags ([[ldg(:); 0], dg, [0; udg(:)]],
                               [-1, 0, 1], n-2, n-2) \ g;
 
-    elseif (n == 4)
+      elseif (n == 4)
+
+        dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)];
+        ldg = h(2) - h(3);
+        udg = h(2) - h(1);
+        c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]],
+                              [-1, 0, 1], n-2, n-2) \ g;
 
-      dg = [h(1) + 2 * h(2); 2 * h(2) + h(3)];
-      ldg = h(2) - h(3);
-      udg = h(2) - h(1);
-      c(2:n-1,:) = spdiags ([[ldg(:);0], dg, [0; udg(:)]],
-                              [-1, 0, 1], n-2, n-2) \ g;
-      
-    else # n == 3
-            
-      dg = h(1) + 2 * h(2);
-      c(2:n-1,:) = g/dg(1);
+      endif
+
+      c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
+      c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:));
+      b = diff (a) ./ h(1:n-1, idx) ...
+          - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
+      d = diff (c) ./ (3 * h(1:n-1, idx));
 
     endif
 
-    c(1,:) = c(2,:) + h(1) / h(2) * (c(2,:) - c(3,:));
-    c(n,:) = c(n-1,:) + h(n-1) / h(n-2) * (c(n-1,:) - c(n-2,:));
-    b = diff (a) ./ h(1:n-1, idx) ...
-        - h(1:n-1, idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
-    d = diff (c) ./ (3 * h(1:n-1, idx));
-
   endif
 
   d = d(1:n-1,:);