From c2dbf1f8827b39b5f5916675f7d319ec1e728e96 Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Fri, 14 Sep 2012 23:50:50 -0400
Subject: [PATCH] Add the divided_difference_coefficients() function and
 reimplement divided_difference() in terms of that.

---
 divided_difference.m              | 14 ++++----------
 divided_difference_coefficients.m | 26 ++++++++++++++++++++++++++
 2 files changed, 30 insertions(+), 10 deletions(-)
 create mode 100644 divided_difference_coefficients.m

diff --git a/divided_difference.m b/divided_difference.m
index 2b8f8df..a9039b1 100644
--- a/divided_difference.m
+++ b/divided_difference.m
@@ -15,7 +15,7 @@ function dd = divided_difference(f, xs)
   ## OUTPUTS:
   ##
   ##   * ``dd`` - The divided difference f[xs(1), xs(2),...]
-  ##
+  ## 
   order = length(xs) - 1;
   
   if (order < 0)
@@ -25,14 +25,8 @@ function dd = divided_difference(f, xs)
     ## Our base case.
     dd = f(xs(1))
   else
-    ## Order >= 1, recurse.
-
-    ## f[x0,...,x_n-1]
-    f0 = divided_difference(f, xs(1:end-1));
-    ## f[x1,...,x_n]
-    f1 = divided_difference(f, xs(2:end));
-
-    # http://mathworld.wolfram.com/DividedDifference.html
-    dd = (f0 - f1)/(xs(1) - xs(end));
+    ## Order >= 1.
+    cs = divided_difference_coefficients(xs)
+    dd = dot(cs, f(xs));
   end
 end
diff --git a/divided_difference_coefficients.m b/divided_difference_coefficients.m
new file mode 100644
index 0000000..53d5799
--- /dev/null
+++ b/divided_difference_coefficients.m
@@ -0,0 +1,26 @@
+function coefficients = divided_difference_coefficients(xs)
+  ## Compute divided difference coefficients of `f` at points `xs`.
+  ##
+  ## INPUTS:
+  ##
+  ##   * ``xs`` - A vector containing x-coordinates.
+  ##
+  ## OUTPUTS:
+  ##
+  ##   * ``coefficients`` - The vector of coefficients such that
+  ##   dot(coefficients, f(xs)) == dd. Used to solve linear systems.
+  ##
+
+  coefficients = [];
+    
+  for xj = xs
+    this_coeff = 1;
+    for xi = xs
+      if (xi != xj)
+	## Append (xj - xi) to the vector of coefficients.
+	this_coeff = this_coeff * (1 / (xj - xi));
+      end
+    end
+    coefficients(end+1) = this_coeff;
+  end
+end
-- 
2.43.2