diff --git a/internal/ceres/CMakeLists.txt b/internal/ceres/CMakeLists.txt
index adad6dc..0960b95 100644
--- a/internal/ceres/CMakeLists.txt
+++ b/internal/ceres/CMakeLists.txt
@@ -66,7 +66,7 @@
     normal_prior.cc
     parameter_block_ordering.cc
     partitioned_matrix_view.cc
-    polynomial_solver.cc
+    polynomial.cc
     problem.cc
     problem_impl.cc
     program.cc
@@ -238,7 +238,7 @@
   CERES_TEST(parameter_block)
   CERES_TEST(parameter_block_ordering)
   CERES_TEST(partitioned_matrix_view)
-  CERES_TEST(polynomial_solver)
+  CERES_TEST(polynomial)
   CERES_TEST(problem)
   CERES_TEST(residual_block)
   CERES_TEST(residual_block_utils)
diff --git a/internal/ceres/dogleg_strategy.cc b/internal/ceres/dogleg_strategy.cc
index 668fa54..da861fe 100644
--- a/internal/ceres/dogleg_strategy.cc
+++ b/internal/ceres/dogleg_strategy.cc
@@ -35,7 +35,7 @@
 #include "ceres/array_utils.h"
 #include "ceres/internal/eigen.h"
 #include "ceres/linear_solver.h"
-#include "ceres/polynomial_solver.h"
+#include "ceres/polynomial.h"
 #include "ceres/sparse_matrix.h"
 #include "ceres/trust_region_strategy.h"
 #include "ceres/types.h"
diff --git a/internal/ceres/polynomial_solver.cc b/internal/ceres/polynomial.cc
similarity index 62%
rename from internal/ceres/polynomial_solver.cc
rename to internal/ceres/polynomial.cc
index 0ece7bc..3b88471 100644
--- a/internal/ceres/polynomial_solver.cc
+++ b/internal/ceres/polynomial.cc
@@ -27,11 +27,14 @@
 // POSSIBILITY OF SUCH DAMAGE.
 //
 // Author: moll.markus@arcor.de (Markus Moll)
+//         sameeragarwal@google.com (Sameer Agarwal)
 
-#include "ceres/polynomial_solver.h"
+#include "ceres/polynomial.h"
 
 #include <cmath>
 #include <cstddef>
+#include <vector>
+
 #include "Eigen/Dense"
 #include "ceres/internal/port.h"
 #include "glog/logging.h"
@@ -179,5 +182,132 @@
   return true;
 }
 
+Vector DifferentiatePolynomial(const Vector& polynomial) {
+  const int degree = polynomial.rows() - 1;
+  CHECK_GE(degree, 0);
+  Vector derivative(degree);
+  for (int i = 0; i < degree; ++i) {
+    derivative(i) = (degree - i) * polynomial(i);
+  }
+
+  return derivative;
+}
+
+void MinimizePolynomial(const Vector& polynomial,
+                        const double x_min,
+                        const double x_max,
+                        double* optimal_x,
+                        double* optimal_value) {
+  // Find the minimum of the polynomial at the two ends.
+  //
+  // We start by inspecting the middle of the interval. Technically
+  // this is not needed, but we do this to make this code as close to
+  // the minFunc package as possible.
+  *optimal_x = (x_min + x_max) / 2.0;
+  *optimal_value = EvaluatePolynomial(polynomial, *optimal_x);
+
+  const double x_min_value = EvaluatePolynomial(polynomial, x_min);
+  if (x_min_value < *optimal_value) {
+    *optimal_value = x_min_value;
+    *optimal_x = x_min;
+  }
+
+  const double x_max_value = EvaluatePolynomial(polynomial, x_max);
+  if (x_max_value < *optimal_value) {
+    *optimal_value = x_max_value;
+    *optimal_x = x_max;
+  }
+
+  // If the polynomial is linear or constant, we are done.
+  if (polynomial.rows() <= 2) {
+    return;
+  }
+
+  const Vector derivative = DifferentiatePolynomial(polynomial);
+  Vector roots_real;
+  if (!FindPolynomialRoots(derivative, &roots_real, NULL)) {
+    LOG(WARNING) << "Unable to find the critical points of "
+                 << "the interpolating polynomial.";
+    return;
+  }
+
+  // This is a bit of an overkill, as some of the roots may actually
+  // have a complex part, but its simpler to just check these values.
+  for (int i = 0; i < roots_real.rows(); ++i) {
+    const double root = roots_real(i);
+    if ((root < x_min) || (root > x_max)) {
+      continue;
+    }
+
+    const double value = EvaluatePolynomial(polynomial, root);
+    if (value < *optimal_value) {
+      *optimal_value = value;
+      *optimal_x = root;
+    }
+  }
+}
+
+Vector FindInterpolatingPolynomial(const vector<FunctionSample>& samples) {
+  const int num_samples = samples.size();
+  int num_constraints = 0;
+  for (int i = 0; i < num_samples; ++i) {
+    if (samples[i].value_is_valid) {
+      ++num_constraints;
+    }
+    if (samples[i].gradient_is_valid) {
+      ++num_constraints;
+    }
+  }
+
+  const int degree = num_constraints - 1;
+  Matrix lhs = Matrix::Zero(num_constraints, num_constraints);
+  Vector rhs = Vector::Zero(num_constraints);
+
+  int row = 0;
+  for (int i = 0; i < num_samples; ++i) {
+    const FunctionSample& sample = samples[i];
+    if (sample.value_is_valid) {
+      LOG(INFO) << "value constraint";
+      for (int j = 0; j <= degree; ++j) {
+        lhs(row, j) = pow(sample.x, degree - j);
+      }
+      rhs(row) = sample.value;
+      ++row;
+    }
+
+    if (sample.gradient_is_valid) {
+      for (int j = 0; j < degree; ++j) {
+        LOG(INFO) << "gradient constraint";
+        lhs(row, j) = (degree - j) * pow(sample.x, degree - j - 1);
+      }
+      rhs(row) = sample.gradient;
+      ++row;
+    }
+  }
+
+  return lhs.fullPivLu().solve(rhs);
+}
+
+void MinimizeInterpolatingPolynomial(const vector<FunctionSample>& samples,
+                                     double x_min,
+                                     double x_max,
+                                     double* optimal_x,
+                                     double* optimal_value) {
+  const Vector polynomial = FindInterpolatingPolynomial(samples);
+  MinimizePolynomial(polynomial, x_min, x_max, optimal_x, optimal_value);
+  for (int i = 0; i < samples.size(); ++i) {
+    const FunctionSample& sample = samples[i];
+    if ((sample.x < x_min) || (sample.x > x_max)) {
+      continue;
+    }
+
+    const double value = EvaluatePolynomial(polynomial, sample.x);
+    if (value < *optimal_value) {
+      *optimal_x = sample.x;
+      *optimal_value = value;
+    }
+  }
+}
+
 }  // namespace internal
 }  // namespace ceres
diff --git a/internal/ceres/polynomial.h b/internal/ceres/polynomial.h
new file mode 100644
index 0000000..42ffdcb
--- /dev/null
+++ b/internal/ceres/polynomial.h
@@ -0,0 +1,134 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: moll.markus@arcor.de (Markus Moll)
+//         sameeragarwal@google.com (Sameer Agarwal)
+
+#ifndef CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
+#define CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
+
+#include <vector>
+#include "ceres/internal/eigen.h"
+#include "ceres/internal/port.h"
+
+namespace ceres {
+namespace internal {
+
+// All polynomials are assumed to be the form
+//
+//   sum_{i=0}^N polynomial(i) x^{N-i}.
+//
+// and are given by a vector of coefficients of size N + 1.
+
+// Evaluate the polynomial at x using the Horner scheme.
+inline double EvaluatePolynomial(const Vector& polynomial, double x) {
+  double v = 0.0;
+  for (int i = 0; i < polynomial.size(); ++i) {
+    v = v * x + polynomial(i);
+  }
+  return v;
+}
+
+// Use the companion matrix eigenvalues to determine the roots of the
+// polynomial.
+//
+// This function returns true on success, false otherwise.
+// Failure indicates that the polynomial is invalid (of size 0) or
+// that the eigenvalues of the companion matrix could not be computed.
+// On failure, a more detailed message will be written to LOG(ERROR).
+// If real is not NULL, the real parts of the roots will be returned in it.
+// Likewise, if imaginary is not NULL, imaginary parts will be returned in it.
+bool FindPolynomialRoots(const Vector& polynomial,
+                         Vector* real,
+                         Vector* imaginary);
+
+// Return the derivative of the given polynomial. It is assumed that
+// the input polynomial is at least of degree zero.
+Vector DifferentiatePolynomial(const Vector& polynomial);
+
+// Find the minimum value of the polynomial in the interval [x_min,
+// x_max]. The minimum is obtained by computing all the roots of the
+// derivative of the input polynomial. All real roots within the
+// interval [x_min, x_max] are considered as well as the end points
+// x_min and x_max. Since polynomials are differentiable functions,
+// this ensures that the true minimum is found.
+void MinimizePolynomial(const Vector& polynomial,
+                        double x_min,
+                        double x_max,
+                        double* optimal_x,
+                        double* optimal_value);
+
+// Structure for storing sample values of a function.
+//
+// Clients can use this struct to communicate the value of the
+// function and or its gradient at a given point x.
+struct FunctionSample {
+  FunctionSample()
+      : x(0.0),
+        value(0.0),
+        value_is_valid(false),
+        gradient(0.0),
+        gradient_is_valid(false) {
+  }
+
+  double x;
+  double value;      // value = f(x)
+  bool value_is_valid;
+  double gradient;   // gradient = f'(x)
+  bool gradient_is_valid;
+};
+
+// Given a set of function value and/or gradient samples, find a
+// polynomial whose value and gradients are exactly equal to the ones
+// in samples.
+//
+// Generally speaking,
+//
+// degree = # values + # gradients - 1
+//
+// Of course its possible to sample a polynomial any number of times,
+// in which case, generally speaking the spurious higher order
+// coefficients will be zero.
+Vector FindInterpolatingPolynomial(const vector<FunctionSample>& samples);
+
+// Interpolate the function described by samples with a polynomial,
+// and minimize it on the interval [x_min, x_max]. Depending on the
+// input samples, it is possible that the interpolation or the root
+// finding algorithms may fail due to numerical difficulties. But the
+// function is guaranteed to return its best guess of an answer, by
+// considering the samples and the end points as possible solutions.
+void MinimizeInterpolatingPolynomial(const vector<FunctionSample>& samples,
+                                     double x_min,
+                                     double x_max,
+                                     double* optimal_x,
+                                     double* optimal_value);
+
+}  // namespace internal
+}  // namespace ceres
+
+#endif  // CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
diff --git a/internal/ceres/polynomial_solver.h b/internal/ceres/polynomial_solver.h
deleted file mode 100644
index 1cf07dd..0000000
--- a/internal/ceres/polynomial_solver.h
+++ /dev/null
@@ -1,65 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2012 Google Inc. All rights reserved.
-// http://code.google.com/p/ceres-solver/
-//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are met:
-//
-// * Redistributions of source code must retain the above copyright notice,
-//   this list of conditions and the following disclaimer.
-// * Redistributions in binary form must reproduce the above copyright notice,
-//   this list of conditions and the following disclaimer in the documentation
-//   and/or other materials provided with the distribution.
-// * Neither the name of Google Inc. nor the names of its contributors may be
-//   used to endorse or promote products derived from this software without
-//   specific prior written permission.
-//
-// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
-// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-// POSSIBILITY OF SUCH DAMAGE.
-//
-// Author: moll.markus@arcor.de (Markus Moll)
-
-#ifndef CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
-#define CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
-
-#include "ceres/internal/eigen.h"
-
-namespace ceres {
-namespace internal {
-
-// Use the companion matrix eigenvalues to determine the roots of the polynomial
-//
-//   sum_{i=0}^N polynomial(i) x^{N-i}.
-//
-// This function returns true on success, false otherwise.
-// Failure indicates that the polynomial is invalid (of size 0) or
-// that the eigenvalues of the companion matrix could not be computed.
-// On failure, a more detailed message will be written to LOG(ERROR).
-// If real is not NULL, the real parts of the roots will be returned in it.
-// Likewise, if imaginary is not NULL, imaginary parts will be returned in it.
-bool FindPolynomialRoots(const Vector& polynomial,
-                         Vector* real,
-                         Vector* imaginary);
-
-// Evaluate the polynomial at x using the Horner scheme.
-inline double EvaluatePolynomial(const Vector& polynomial, double x) {
-  double v = 0.0;
-  for (int i = 0; i < polynomial.size(); ++i) {
-    v = v * x + polynomial(i);
-  }
-  return v;
-}
-
-}  // namespace internal
-}  // namespace ceres
-
-#endif  // CERES_INTERNAL_POLYNOMIAL_SOLVER_H_
diff --git a/internal/ceres/polynomial_solver_test.cc b/internal/ceres/polynomial_solver_test.cc
deleted file mode 100644
index ae347fb..0000000
--- a/internal/ceres/polynomial_solver_test.cc
+++ /dev/null
@@ -1,223 +0,0 @@
-// Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2012 Google Inc. All rights reserved.
-// http://code.google.com/p/ceres-solver/
-//
-// Redistribution and use in source and binary forms, with or without
-// modification, are permitted provided that the following conditions are met:
-//
-// * Redistributions of source code must retain the above copyright notice,
-//   this list of conditions and the following disclaimer.
-// * Redistributions in binary form must reproduce the above copyright notice,
-//   this list of conditions and the following disclaimer in the documentation
-//   and/or other materials provided with the distribution.
-// * Neither the name of Google Inc. nor the names of its contributors may be
-//   used to endorse or promote products derived from this software without
-//   specific prior written permission.
-//
-// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
-// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-// POSSIBILITY OF SUCH DAMAGE.
-//
-// Author: moll.markus@arcor.de (Markus Moll)
-
-#include "ceres/polynomial_solver.h"
-
-#include <limits>
-#include <cmath>
-#include <cstddef>
-#include <algorithm>
-#include "gtest/gtest.h"
-#include "ceres/test_util.h"
-
-namespace ceres {
-namespace internal {
-namespace {
-
-// For IEEE-754 doubles, machine precision is about 2e-16.
-const double kEpsilon = 1e-13;
-const double kEpsilonLoose = 1e-9;
-
-// Return the constant polynomial p(x) = 1.23.
-Vector ConstantPolynomial(double value) {
-  Vector poly(1);
-  poly(0) = value;
-  return poly;
-}
-
-// Return the polynomial p(x) = poly(x) * (x - root).
-Vector AddRealRoot(const Vector& poly, double root) {
-  Vector poly2(poly.size() + 1);
-  poly2.setZero();
-  poly2.head(poly.size()) += poly;
-  poly2.tail(poly.size()) -= root * poly;
-  return poly2;
-}
-
-// Return the polynomial
-// p(x) = poly(x) * (x - real - imag*i) * (x - real + imag*i).
-Vector AddComplexRootPair(const Vector& poly, double real, double imag) {
-  Vector poly2(poly.size() + 2);
-  poly2.setZero();
-  // Multiply poly by x^2 - 2real + abs(real,imag)^2
-  poly2.head(poly.size()) += poly;
-  poly2.segment(1, poly.size()) -= 2 * real * poly;
-  poly2.tail(poly.size()) += (real*real + imag*imag) * poly;
-  return poly2;
-}
-
-// Sort the entries in a vector.
-// Needed because the roots are not returned in sorted order.
-Vector SortVector(const Vector& in) {
-  Vector out(in);
-  std::sort(out.data(), out.data() + out.size());
-  return out;
-}
-
-// Run a test with the polynomial defined by the N real roots in roots_real.
-// If use_real is false, NULL is passed as the real argument to
-// FindPolynomialRoots. If use_imaginary is false, NULL is passed as the
-// imaginary argument to FindPolynomialRoots.
-template<int N>
-void RunPolynomialTestRealRoots(const double (&real_roots)[N],
-                                bool use_real,
-                                bool use_imaginary,
-                                double epsilon) {
-  Vector real;
-  Vector imaginary;
-  Vector poly = ConstantPolynomial(1.23);
-  for (int i = 0; i < N; ++i) {
-    poly = AddRealRoot(poly, real_roots[i]);
-  }
-  Vector* const real_ptr = use_real ? &real : NULL;
-  Vector* const imaginary_ptr = use_imaginary ? &imaginary : NULL;
-  bool success = FindPolynomialRoots(poly, real_ptr, imaginary_ptr);
-
-  EXPECT_EQ(success, true);
-  if (use_real) {
-    EXPECT_EQ(real.size(), N);
-    real = SortVector(real);
-    ExpectArraysClose(N, real.data(), real_roots, epsilon);
-  }
-  if (use_imaginary) {
-    EXPECT_EQ(imaginary.size(), N);
-    const Vector zeros = Vector::Zero(N);
-    ExpectArraysClose(N, imaginary.data(), zeros.data(), epsilon);
-  }
-}
-}  // namespace
-
-TEST(PolynomialSolver, InvalidPolynomialOfZeroLengthIsRejected) {
-  // Vector poly(0) is an ambiguous constructor call, so
-  // use the constructor with explicit column count.
-  Vector poly(0, 1);
-  Vector real;
-  Vector imag;
-  bool success = FindPolynomialRoots(poly, &real, &imag);
-
-  EXPECT_EQ(success, false);
-}
-
-TEST(PolynomialSolver, ConstantPolynomialReturnsNoRoots) {
-  Vector poly = ConstantPolynomial(1.23);
-  Vector real;
-  Vector imag;
-  bool success = FindPolynomialRoots(poly, &real, &imag);
-
-  EXPECT_EQ(success, true);
-  EXPECT_EQ(real.size(), 0);
-  EXPECT_EQ(imag.size(), 0);
-}
-
-TEST(PolynomialSolver, LinearPolynomialWithPositiveRootWorks) {
-  const double roots[1] = { 42.42 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
-}
-
-TEST(PolynomialSolver, LinearPolynomialWithNegativeRootWorks) {
-  const double roots[1] = { -42.42 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
-}
-
-TEST(PolynomialSolver, QuadraticPolynomialWithPositiveRootsWorks) {
-  const double roots[2] = { 1.0, 42.42 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
-}
-
-TEST(PolynomialSolver, QuadraticPolynomialWithOneNegativeRootWorks) {
-  const double roots[2] = { -42.42, 1.0 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
-}
-
-TEST(PolynomialSolver, QuadraticPolynomialWithTwoNegativeRootsWorks) {
-  const double roots[2] = { -42.42, -1.0 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
-}
-
-TEST(PolynomialSolver, QuadraticPolynomialWithCloseRootsWorks) {
-  const double roots[2] = { 42.42, 42.43 };
-  RunPolynomialTestRealRoots(roots, true, false, kEpsilonLoose);
-}
-
-TEST(PolynomialSolver, QuadraticPolynomialWithComplexRootsWorks) {
-  Vector real;
-  Vector imag;
-
-  Vector poly = ConstantPolynomial(1.23);
-  poly = AddComplexRootPair(poly, 42.42, 4.2);
-  bool success = FindPolynomialRoots(poly, &real, &imag);
-
-  EXPECT_EQ(success, true);
-  EXPECT_EQ(real.size(), 2);
-  EXPECT_EQ(imag.size(), 2);
-  ExpectClose(real(0), 42.42, kEpsilon);
-  ExpectClose(real(1), 42.42, kEpsilon);
-  ExpectClose(std::abs(imag(0)), 4.2, kEpsilon);
-  ExpectClose(std::abs(imag(1)), 4.2, kEpsilon);
-  ExpectClose(std::abs(imag(0) + imag(1)), 0.0, kEpsilon);
-}
-
-TEST(PolynomialSolver, QuarticPolynomialWorks) {
-  const double roots[4] = { 1.23e-4, 1.23e-1, 1.23e+2, 1.23e+5 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
-}
-
-TEST(PolynomialSolver, QuarticPolynomialWithTwoClustersOfCloseRootsWorks) {
-  const double roots[4] = { 1.23e-1, 2.46e-1, 1.23e+5, 2.46e+5 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilonLoose);
-}
-
-TEST(PolynomialSolver, QuarticPolynomialWithTwoZeroRootsWorks) {
-  const double roots[4] = { -42.42, 0.0, 0.0, 42.42 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilonLoose);
-}
-
-TEST(PolynomialSolver, QuarticMonomialWorks) {
-  const double roots[4] = { 0.0, 0.0, 0.0, 0.0 };
-  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
-}
-
-TEST(PolynomialSolver, NullPointerAsImaginaryPartWorks) {
-  const double roots[4] = { 1.23e-4, 1.23e-1, 1.23e+2, 1.23e+5 };
-  RunPolynomialTestRealRoots(roots, true, false, kEpsilon);
-}
-
-TEST(PolynomialSolver, NullPointerAsRealPartWorks) {
-  const double roots[4] = { 1.23e-4, 1.23e-1, 1.23e+2, 1.23e+5 };
-  RunPolynomialTestRealRoots(roots, false, true, kEpsilon);
-}
-
-TEST(PolynomialSolver, BothOutputArgumentsNullWorks) {
-  const double roots[4] = { 1.23e-4, 1.23e-1, 1.23e+2, 1.23e+5 };
-  RunPolynomialTestRealRoots(roots, false, false, kEpsilon);
-}
-
-}  // namespace internal
-}  // namespace ceres
diff --git a/internal/ceres/polynomial_test.cc b/internal/ceres/polynomial_test.cc
new file mode 100644
index 0000000..f56352f
--- /dev/null
+++ b/internal/ceres/polynomial_test.cc
@@ -0,0 +1,512 @@
+// Ceres Solver - A fast non-linear least squares minimizer
+// Copyright 2012 Google Inc. All rights reserved.
+// http://code.google.com/p/ceres-solver/
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are met:
+//
+// * Redistributions of source code must retain the above copyright notice,
+//   this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above copyright notice,
+//   this list of conditions and the following disclaimer in the documentation
+//   and/or other materials provided with the distribution.
+// * Neither the name of Google Inc. nor the names of its contributors may be
+//   used to endorse or promote products derived from this software without
+//   specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+// POSSIBILITY OF SUCH DAMAGE.
+//
+// Author: moll.markus@arcor.de (Markus Moll)
+//         sameeragarwal@google.com (Sameer Agarwal)
+
+#include "ceres/polynomial.h"
+
+#include <limits>
+#include <cmath>
+#include <cstddef>
+#include <algorithm>
+#include "gtest/gtest.h"
+#include "ceres/test_util.h"
+
+namespace ceres {
+namespace internal {
+namespace {
+
+// For IEEE-754 doubles, machine precision is about 2e-16.
+const double kEpsilon = 1e-13;
+const double kEpsilonLoose = 1e-9;
+
+// Return the constant polynomial p(x) = 1.23.
+Vector ConstantPolynomial(double value) {
+  Vector poly(1);
+  poly(0) = value;
+  return poly;
+}
+
+// Return the polynomial p(x) = poly(x) * (x - root).
+Vector AddRealRoot(const Vector& poly, double root) {
+  Vector poly2(poly.size() + 1);
+  poly2.setZero();
+  poly2.head(poly.size()) += poly;
+  poly2.tail(poly.size()) -= root * poly;
+  return poly2;
+}
+
+// Return the polynomial
+// p(x) = poly(x) * (x - real - imag*i) * (x - real + imag*i).
+Vector AddComplexRootPair(const Vector& poly, double real, double imag) {
+  Vector poly2(poly.size() + 2);
+  poly2.setZero();
+  // Multiply poly by x^2 - 2real + abs(real,imag)^2
+  poly2.head(poly.size()) += poly;
+  poly2.segment(1, poly.size()) -= 2 * real * poly;
+  poly2.tail(poly.size()) += (real*real + imag*imag) * poly;
+  return poly2;
+}
+
+// Sort the entries in a vector.
+// Needed because the roots are not returned in sorted order.
+Vector SortVector(const Vector& in) {
+  Vector out(in);
+  std::sort(out.data(), out.data() + out.size());
+  return out;
+}
+
+// Run a test with the polynomial defined by the N real roots in roots_real.
+// If use_real is false, NULL is passed as the real argument to
+// FindPolynomialRoots. If use_imaginary is false, NULL is passed as the
+// imaginary argument to FindPolynomialRoots.
+template<int N>
+void RunPolynomialTestRealRoots(const double (&real_roots)[N],
+                                bool use_real,
+                                bool use_imaginary,
+                                double epsilon) {
+  Vector real;
+  Vector imaginary;
+  Vector poly = ConstantPolynomial(1.23);
+  for (int i = 0; i < N; ++i) {
+    poly = AddRealRoot(poly, real_roots[i]);
+  }
+  Vector* const real_ptr = use_real ? &real : NULL;
+  Vector* const imaginary_ptr = use_imaginary ? &imaginary : NULL;
+  bool success = FindPolynomialRoots(poly, real_ptr, imaginary_ptr);
+
+  EXPECT_EQ(success, true);
+  if (use_real) {
+    EXPECT_EQ(real.size(), N);
+    real = SortVector(real);
+    ExpectArraysClose(N, real.data(), real_roots, epsilon);
+  }
+  if (use_imaginary) {
+    EXPECT_EQ(imaginary.size(), N);
+    const Vector zeros = Vector::Zero(N);
+    ExpectArraysClose(N, imaginary.data(), zeros.data(), epsilon);
+  }
+}
+}  // namespace
+
+TEST(Polynomial, InvalidPolynomialOfZeroLengthIsRejected) {
+  // Vector poly(0) is an ambiguous constructor call, so
+  // use the constructor with explicit column count.
+  Vector poly(0, 1);
+  Vector real;
+  Vector imag;
+  bool success = FindPolynomialRoots(poly, &real, &imag);
+
+  EXPECT_EQ(success, false);
+}
+
+TEST(Polynomial, ConstantPolynomialReturnsNoRoots) {
+  Vector poly = ConstantPolynomial(1.23);
+  Vector real;
+  Vector imag;
+  bool success = FindPolynomialRoots(poly, &real, &imag);
+
+  EXPECT_EQ(success, true);
+  EXPECT_EQ(real.size(), 0);
+  EXPECT_EQ(imag.size(), 0);
+}
+
+TEST(Polynomial, LinearPolynomialWithPositiveRootWorks) {
+  const double roots[1] = { 42.42 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
+}
+
+TEST(Polynomial, LinearPolynomialWithNegativeRootWorks) {
+  const double roots[1] = { -42.42 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
+}
+
+TEST(Polynomial, QuadraticPolynomialWithPositiveRootsWorks) {
+  const double roots[2] = { 1.0, 42.42 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
+}
+
+TEST(Polynomial, QuadraticPolynomialWithOneNegativeRootWorks) {
+  const double roots[2] = { -42.42, 1.0 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
+}
+
+TEST(Polynomial, QuadraticPolynomialWithTwoNegativeRootsWorks) {
+  const double roots[2] = { -42.42, -1.0 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
+}
+
+TEST(Polynomial, QuadraticPolynomialWithCloseRootsWorks) {
+  const double roots[2] = { 42.42, 42.43 };
+  RunPolynomialTestRealRoots(roots, true, false, kEpsilonLoose);
+}
+
+TEST(Polynomial, QuadraticPolynomialWithComplexRootsWorks) {
+  Vector real;
+  Vector imag;
+
+  Vector poly = ConstantPolynomial(1.23);
+  poly = AddComplexRootPair(poly, 42.42, 4.2);
+  bool success = FindPolynomialRoots(poly, &real, &imag);
+
+  EXPECT_EQ(success, true);
+  EXPECT_EQ(real.size(), 2);
+  EXPECT_EQ(imag.size(), 2);
+  ExpectClose(real(0), 42.42, kEpsilon);
+  ExpectClose(real(1), 42.42, kEpsilon);
+  ExpectClose(std::abs(imag(0)), 4.2, kEpsilon);
+  ExpectClose(std::abs(imag(1)), 4.2, kEpsilon);
+  ExpectClose(std::abs(imag(0) + imag(1)), 0.0, kEpsilon);
+}
+
+TEST(Polynomial, QuarticPolynomialWorks) {
+  const double roots[4] = { 1.23e-4, 1.23e-1, 1.23e+2, 1.23e+5 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
+}
+
+TEST(Polynomial, QuarticPolynomialWithTwoClustersOfCloseRootsWorks) {
+  const double roots[4] = { 1.23e-1, 2.46e-1, 1.23e+5, 2.46e+5 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilonLoose);
+}
+
+TEST(Polynomial, QuarticPolynomialWithTwoZeroRootsWorks) {
+  const double roots[4] = { -42.42, 0.0, 0.0, 42.42 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilonLoose);
+}
+
+TEST(Polynomial, QuarticMonomialWorks) {
+  const double roots[4] = { 0.0, 0.0, 0.0, 0.0 };
+  RunPolynomialTestRealRoots(roots, true, true, kEpsilon);
+}
+
+TEST(Polynomial, NullPointerAsImaginaryPartWorks) {
+  const double roots[4] = { 1.23e-4, 1.23e-1, 1.23e+2, 1.23e+5 };
+  RunPolynomialTestRealRoots(roots, true, false, kEpsilon);
+}
+
+TEST(Polynomial, NullPointerAsRealPartWorks) {
+  const double roots[4] = { 1.23e-4, 1.23e-1, 1.23e+2, 1.23e+5 };
+  RunPolynomialTestRealRoots(roots, false, true, kEpsilon);
+}
+
+TEST(Polynomial, BothOutputArgumentsNullWorks) {
+  const double roots[4] = { 1.23e-4, 1.23e-1, 1.23e+2, 1.23e+5 };
+  RunPolynomialTestRealRoots(roots, false, false, kEpsilon);
+}
+
+TEST(Polynomial, DifferentiateConstantPolynomial) {
+  // p(x) = 1;
+  Vector polynomial(1);
+  polynomial(0) = 1.0;
+  const Vector derivative = DifferentiatePolynomial(polynomial);
+  EXPECT_EQ(derivative.rows(), 0);
+}
+
+TEST(Polynomial, DifferentiateQuadraticPolynomial) {
+  // p(x) = x^2 + 2x + 3;
+  Vector polynomial(3);
+  polynomial(0) = 1.0;
+  polynomial(1) = 2.0;
+  polynomial(2) = 3.0;
+
+  const Vector derivative = DifferentiatePolynomial(polynomial);
+  EXPECT_EQ(derivative.rows(), 2);
+  EXPECT_EQ(derivative(0), 2.0);
+  EXPECT_EQ(derivative(1), 2.0);
+}
+
+TEST(Polynomial, MinimizeConstantPolynomial) {
+  // p(x) = 1;
+  Vector polynomial(1);
+  polynomial(0) = 1.0;
+
+  double optimal_x = 0.0;
+  double optimal_value = 0.0;
+  double min_x = 0.0;
+  double max_x = 1.0;
+  MinimizePolynomial(polynomial, min_x, max_x, &optimal_x, &optimal_value);
+
+  EXPECT_EQ(optimal_value, 1.0);
+  EXPECT_LE(optimal_x, max_x);
+  EXPECT_GE(optimal_x, min_x);
+}
+
+TEST(Polynomial, MinimizeLinearPolynomial) {
+  // p(x) = x - 2
+  Vector polynomial(2);
+
+  polynomial(0) = 1.0;
+  polynomial(1) = 2.0;
+
+  double optimal_x = 0.0;
+  double optimal_value = 0.0;
+  double min_x = 0.0;
+  double max_x = 1.0;
+  MinimizePolynomial(polynomial, min_x, max_x, &optimal_x, &optimal_value);
+
+  EXPECT_EQ(optimal_x, 0.0);
+  EXPECT_EQ(optimal_value, 2.0);
+}
+
+
+TEST(Polynomial, MinimizeQuadraticPolynomial) {
+  // p(x) = x^2 - 3 x + 2
+  // min_x = 3/2
+  // min_value = -1/4;
+  Vector polynomial(3);
+  polynomial(0) = 1.0;
+  polynomial(1) = -3.0;
+  polynomial(2) = 2.0;
+
+  double optimal_x = 0.0;
+  double optimal_value = 0.0;
+  double min_x = -2.0;
+  double max_x = 2.0;
+  MinimizePolynomial(polynomial, min_x, max_x, &optimal_x, &optimal_value);
+  EXPECT_EQ(optimal_x, 3.0/2.0);
+  EXPECT_EQ(optimal_value, -1.0/4.0);
+
+  min_x = -2.0;
+  max_x = 1.0;
+  MinimizePolynomial(polynomial, min_x, max_x, &optimal_x, &optimal_value);
+  EXPECT_EQ(optimal_x, 1.0);
+  EXPECT_EQ(optimal_value, 0.0);
+
+  min_x = 2.0;
+  max_x = 3.0;
+  MinimizePolynomial(polynomial, min_x, max_x, &optimal_x, &optimal_value);
+  EXPECT_EQ(optimal_x, 2.0);
+  EXPECT_EQ(optimal_value, 0.0);
+}
+
+TEST(Polymomial, ConstantInterpolatingPolynomial) {
+  // p(x) = 1.0
+  Vector true_polynomial(1);
+  true_polynomial << 1.0;
+
+  vector<FunctionSample> samples;
+  FunctionSample sample;
+  sample.x = 1.0;
+  sample.value = 1.0;
+  sample.value_is_valid = true;
+  samples.push_back(sample);
+
+  const Vector polynomial = FindInterpolatingPolynomial(samples);
+  EXPECT_NEAR((true_polynomial - polynomial).norm(), 0.0, 1e-15);
+}
+
+TEST(Polynomial, LinearInterpolatingPolynomial) {
+  // p(x) = 2x - 1
+  Vector true_polynomial(2);
+  true_polynomial << 2.0, -1.0;
+
+  vector<FunctionSample> samples;
+  FunctionSample sample;
+  sample.x = 1.0;
+  sample.value = 1.0;
+  sample.value_is_valid = true;
+  sample.gradient = 2.0;
+  sample.gradient_is_valid = true;
+  samples.push_back(sample);
+
+  const Vector polynomial = FindInterpolatingPolynomial(samples);
+  EXPECT_NEAR((true_polynomial - polynomial).norm(), 0.0, 1e-15);
+}
+
+TEST(Polynomial, QuadraticInterpolatingPolynomial) {
+  // p(x) = 2x^2 + 3x + 2
+   Vector true_polynomial(3);
+   true_polynomial << 2.0, 3.0, 2.0;
+
+  vector<FunctionSample> samples;
+  {
+    FunctionSample sample;
+    sample.x = 1.0;
+    sample.value = 7.0;
+    sample.value_is_valid = true;
+    sample.gradient = 7.0;
+    sample.gradient_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  {
+    FunctionSample sample;
+    sample.x = -3.0;
+    sample.value = 11.0;
+    sample.value_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  Vector polynomial = FindInterpolatingPolynomial(samples);
+  EXPECT_NEAR((true_polynomial - polynomial).norm(), 0.0, 1e-15);
+}
+
+TEST(Polynomial, DeficientCubicInterpolatingPolynomial) {
+  // p(x) = 2x^2 + 3x + 2
+  Vector true_polynomial(4);
+  true_polynomial << 0.0, 2.0, 3.0, 2.0;
+
+  vector<FunctionSample> samples;
+  {
+    FunctionSample sample;
+    sample.x = 1.0;
+    sample.value = 7.0;
+    sample.value_is_valid = true;
+    sample.gradient = 7.0;
+    sample.gradient_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  {
+    FunctionSample sample;
+    sample.x = -3.0;
+    sample.value = 11.0;
+    sample.value_is_valid = true;
+    sample.gradient = -9;
+    sample.gradient_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  const Vector polynomial = FindInterpolatingPolynomial(samples);
+  EXPECT_NEAR((true_polynomial - polynomial).norm(), 0.0, 1e-14);
+}
+
+
+TEST(Polynomial, CubicInterpolatingPolynomialFromValues) {
+  // p(x) = x^3 + 2x^2 + 3x + 2
+ Vector true_polynomial(4);
+ true_polynomial << 1.0, 2.0, 3.0, 2.0;
+
+ vector<FunctionSample> samples;
+  {
+    FunctionSample sample;
+    sample.x = 1.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  {
+    FunctionSample sample;
+    sample.x = -3.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  {
+    FunctionSample sample;
+    sample.x = 2.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  {
+    FunctionSample sample;
+    sample.x = 0.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  const Vector polynomial = FindInterpolatingPolynomial(samples);
+  EXPECT_NEAR((true_polynomial - polynomial).norm(), 0.0, 1e-14);
+}
+
+TEST(Polynomial, CubicInterpolatingPolynomialFromValuesAndOneGradient) {
+  // p(x) = x^3 + 2x^2 + 3x + 2
+ Vector true_polynomial(4);
+ true_polynomial << 1.0, 2.0, 3.0, 2.0;
+ Vector true_gradient_polynomial = DifferentiatePolynomial(true_polynomial);
+
+ vector<FunctionSample> samples;
+  {
+    FunctionSample sample;
+    sample.x = 1.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  {
+    FunctionSample sample;
+    sample.x = -3.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  {
+    FunctionSample sample;
+    sample.x = 2.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    sample.gradient = EvaluatePolynomial(true_gradient_polynomial, sample.x);
+    sample.gradient_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  const Vector polynomial = FindInterpolatingPolynomial(samples);
+  EXPECT_NEAR((true_polynomial - polynomial).norm(), 0.0, 1e-14);
+}
+
+TEST(Polynomial, CubicInterpolatingPolynomialFromValuesAndGradients) {
+  // p(x) = x^3 + 2x^2 + 3x + 2
+ Vector true_polynomial(4);
+ true_polynomial << 1.0, 2.0, 3.0, 2.0;
+ Vector true_gradient_polynomial = DifferentiatePolynomial(true_polynomial);
+
+ vector<FunctionSample> samples;
+  {
+    FunctionSample sample;
+    sample.x = -3.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    sample.gradient = EvaluatePolynomial(true_gradient_polynomial, sample.x);
+    sample.gradient_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  {
+    FunctionSample sample;
+    sample.x = 2.0;
+    sample.value = EvaluatePolynomial(true_polynomial, sample.x);
+    sample.value_is_valid = true;
+    sample.gradient = EvaluatePolynomial(true_gradient_polynomial, sample.x);
+    sample.gradient_is_valid = true;
+    samples.push_back(sample);
+  }
+
+  const Vector polynomial = FindInterpolatingPolynomial(samples);
+  EXPECT_NEAR((true_polynomial - polynomial).norm(), 0.0, 1e-14);
+}
+
+}  // namespace internal
+}  // namespace ceres