X-Git-Url: http://gitweb.michael.orlitzky.com/?p=octave.git;a=blobdiff_plain;f=divided_difference.m;h=261abeded4ccf9eabf468f817d1b9fa9eafe0994;hp=7b8020e95c3b145fe637502802a09113126b329e;hb=5c9754929aaa99f400ee720bb5a4753928898eab;hpb=1a5bceb9de6bd3a511e7edcb6ab08f84dae63948

diff --git a/divided_difference.m b/divided_difference.m
index 7b8020e..261abed 100644
--- a/divided_difference.m
+++ b/divided_difference.m
@@ -1,25 +1,33 @@
 function dd = divided_difference(f, xs)
-  ## Compute divided difference of `f` at points `xs`. The argument `xs`
-  ## is assumed to be a vector containing at least one element. If it
-  ## contains n elements, the (n-1)st divided difference will be
-  ## calculated.
-  order = length(xs);
-  
-  if (order < 1)
-    ## Can't do anything here. Return nothing.
-    dd = NA;
-  elseif (order == 1)
-    ## Our base case.
-    dd = f(xs(1))
-  else
-    ## Order > 1, recurse.
+  % Compute divided difference of `f` at points `xs`. The argument `xs`
+  % is assumed to be a vector containing at least one element. If it
+  % contains n elements, the (n-1)st divided difference will be
+  % calculated.
+  %
+  % INPUTS:
+  %
+  %   * ``f`` - The function whose divided differences we want.
+  %
+  %   * ``xs`` - A vector containing x-coordinates. The length of `xs`
+  %   determines the order of the divided difference.
+  %
+  %
+  % OUTPUTS:
+  %
+  %   * ``dd`` - The divided difference f[xs(1), xs(2),...]
+  %
 
-    ## f[x0,...,x_n-1]
-    f0 = divided_difference(f, xs(1:end-1));
-    ## f[x1,...,x_n]
-    f1 = divided_difference(f, xs(2:end));
+  order = length(xs) - 1;
 
-    # http://mathworld.wolfram.com/DividedDifference.html
-    dd = (f0 - f1)/(xs(1) - xs(end))
+  if (order < 0)
+    % Can't do anything here. Return nothing.
+    dd = NA;
+  elseif (order == 0)
+    % Our base case.
+    dd = f(xs(1));
+  else
+    % Order >= 1.
+    cs = divided_difference_coefficients(xs);
+    dd = dot(cs, f(xs));
   end
 end