diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index 1704c49..96b6dcf 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -573,13 +573,33 @@
 }
 
 // pow -- base is a constant, exponent is a differentiable function.
-// For a > 0 we have: (a)^(p + dp) ~= a^p + a^p log(a) dp
-// For a == 0 and p > 0 we have: (a)^(p + dp) ~= 0
+// We have various special cases, see the comment for pow(Jet, Jet) for
+// analysis:
+//
+// 1. For a > 0 we have: (a)^(p + dp) ~= a^p + a^p log(a) dp
+//
+// 2. For a == 0 and p > 0 we have: (a)^(p + dp) ~= a^p
+//
+// 3. For a < 0 and integer p we have: (a)^(p + dp) ~= a^p
+
 template <typename T, int N> inline
 Jet<T, N> pow(double f, const Jet<T, N>& g) {
   if (f == 0 && g.a > 0) {
+    // Handle case 2.
     return Jet<T, N>(T(0.0));
   }
+  if (f < 0 && g.a == floor(g.a)) {
+    // Handle case 3.
+    Jet<T, N> ret(pow(f, g.a));
+    for (int i = 0; i < N; i++) {
+      if (g.v[i] != T(0.0)) {
+        // Return a NaN when g.v != 0.
+        ret.v[i] = std::numeric_limits<T>::quiet_NaN();
+      }
+    }
+    return ret;
+  }
+  // Handle case 1.
   T const tmp = pow(f, g.a);
   return Jet<T, N>(tmp, log(f) * tmp * g.v);
 }
@@ -587,40 +607,62 @@
 // pow -- both base and exponent are differentiable functions. This has a
 // variety of special cases that require careful handling.
 //
-// * For a > 0: (a + da)^(b + db) ~= a^b + a^(b - 1) * (b*da + a*log(a)*db)
-//   The numerical evaluation of a*log(a) for a > 0 is well behaved, even for
-//   extremely small values (e.g. 1e-99).
+// 1. For a > 0: (a + da)^(b + db) ~= a^b + a^(b - 1) * (b*da + a*log(a)*db)
+//    The numerical evaluation of a*log(a) for a > 0 is well behaved, even for
+//    extremely small values (e.g. 1e-99).
 //
-// * For a == 0 and b > 1: (a + da)^(b + db) ~= 0
-//   This cases is needed because log(0) can not be evaluated in the a > 0
-//   expression. However the function a*log(a) is well behaved around a == 0
-//   and its limit as a-->0 is zero.
+// 2. For a == 0 and b > 1: (a + da)^(b + db) ~= 0
+//    This cases is needed because log(0) can not be evaluated in the a > 0
+//    expression. However the function a*log(a) is well behaved around a == 0
+//    and its limit as a-->0 is zero.
 //
-// * For a == 0 and b == 1: (a + da)^(b + db) ~= 0 + da
+// 3. For a == 0 and b == 1: (a + da)^(b + db) ~= 0 + da
 //
-// * For a == 0 and 0 < b < 1: The value is finite but the derivatives are not.
+// 4. For a == 0 and 0 < b < 1: The value is finite but the derivatives are not.
 //
-// * For a == 0 and b < 0: The value and derivatives of a^b are not finite.
+// 5. For a == 0 and b < 0: The value and derivatives of a^b are not finite.
 //
-// * For a == 0 and b == 0: The C standard incorrectly defines 0^0 to be 1
-//   "because there are applications that can exploit this definition". We
-//   (arbitrarily) decree that derivatives here will be nonfinite, since that
-//   is consistent with the behavior for b < 0 and 0 < b < 1. Practically any
-//   definition could have been justified because mathematical consistency has
-//   been lost at this point.
+// 6. For a == 0 and b == 0: The C standard incorrectly defines 0^0 to be 1
+//    "because there are applications that can exploit this definition". We
+//    (arbitrarily) decree that derivatives here will be nonfinite, since that
+//    is consistent with the behavior for a==0, b < 0 and 0 < b < 1. Practically
+//    any definition could have been justified because mathematical consistency
+//    has been lost at this point.
+//
+// 7. For a < 0, b integer, db == 0: (a + da)^(b + db) ~= a^b + b * a^(b - 1) da
+//    This is equivalent to the case where a is a differentiable function and b
+//    is a constant (to first order).
+//
+// 8. For a < 0, b integer, db != 0: The value is finite but the derivatives are
+//    not, because any change in the value of b moves us away from the point
+//    with a real-valued answer into the region with complex-valued answers.
+//
+// 9. For a < 0, b noninteger: The value and derivatives of a^b are not finite.
 
 template <typename T, int N> inline
 Jet<T, N> pow(const Jet<T, N>& f, const Jet<T, N>& g) {
   if (f.a == 0 && g.a >= 1) {
-    // Handle the special cases when f == 0.
+    // Handle cases 2 and 3.
     if (g.a > 1) {
       return Jet<T, N>(T(0.0));
     }
     return f;
   }
-  // Handle the generic case for f != 0. We also handle f == 0, g < 1 here and
-  // allow the log() function to generate -HUGE_VAL, since this results in a
-  // nonfinite derivative.
+  if (f.a < 0 && g.a == floor(g.a)) {
+    // Handle cases 7 and 8.
+    T const tmp = g.a * pow(f.a, g.a - T(1.0));
+    Jet<T, N> ret(pow(f.a, g.a), tmp * f.v);
+    for (int i = 0; i < N; i++) {
+      if (g.v[i] != T(0.0)) {
+        // Return a NaN when g.v != 0.
+        ret.v[i] = std::numeric_limits<T>::quiet_NaN();
+      }
+    }
+    return ret;
+  }
+  // Handle the remaining cases. For cases 4,5,6,9 we allow the log() function
+  // to generate -HUGE_VAL or NaN, since those cases result in a nonfinite
+  // derivative.
   T const tmp1 = pow(f.a, g.a);
   T const tmp2 = g.a * pow(f.a, g.a - T(1.0));
   T const tmp3 = tmp1 * log(f.a);
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index 589fa56..ec43ef9 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -283,7 +283,7 @@
     }
     for (int i = -10; i < 0; i++) {
       J a = MakeJet(0, 1, 2);
-      J b = MakeJet(i*0.1, 3, 4);       // b = -0.1,-0.9 ... 0
+      J b = MakeJet(i*0.1, 3, 4);       // b = -1,-0.9 ... -0.1
       VL << "a = " << a;
       VL << "b = " << b;
 
@@ -309,6 +309,39 @@
     }
   }
 
+  { // Check that pow(<0, b) is correct for integer b.
+    // This tests special case handling inside pow().
+    J a = MakeJet(-1.5, 3, 4);
+
+    // b integer:
+    for (int i = -10; i <= 10; i++) {
+      J b = MakeJet(i, 0, 5);
+      VL << "a = " << a;
+      VL << "b = " << b;
+
+      J c = pow(a, b);
+      VL << "a^b = " << c;
+      ExpectClose(c.a, pow(-1.5, i), kTolerance);
+      EXPECT_TRUE(IsFinite(c.v[0]));
+      EXPECT_FALSE(IsFinite(c.v[1]));
+      ExpectClose(c.v[0], i * pow(-1.5, i - 1) * 3.0, kTolerance);
+    }
+  }
+
+  { // Check that pow(<0, b) is correct for noninteger b.
+    // This tests special case handling inside pow().
+    J a = MakeJet(-1.5, 3, 4);
+    J b = MakeJet(-2.5, 0, 5);
+    VL << "a = " << a;
+    VL << "b = " << b;
+
+    J c = pow(a, b);
+    VL << "a^b = " << c;
+    EXPECT_FALSE(IsFinite(c.a));
+    EXPECT_FALSE(IsFinite(c.v[0]));
+    EXPECT_FALSE(IsFinite(c.v[1]));
+  }
+
   {
     // Check that pow(0,y) == 0 for y == 2, with the second argument a
     // Jet.  This tests special case handling inside pow().
@@ -322,6 +355,39 @@
     ExpectJetsClose(c, MakeJet(0, 0, 0));
   }
 
+  {
+    // Check that pow(<0,y) is correct for integer y. This tests special case
+    // handling inside pow().
+    double a = -1.5;
+    for (int i = -10; i <= 10; i++) {
+      J b = MakeJet(i, 3, 0);
+      VL << "a = " << a;
+      VL << "b = " << b;
+
+      J c = pow(a, b);
+      VL << "a^b = " << c;
+      ExpectClose(c.a, pow(-1.5, i), kTolerance);
+      EXPECT_FALSE(IsFinite(c.v[0]));
+      EXPECT_TRUE(IsFinite(c.v[1]));
+      ExpectClose(c.v[1], 0, kTolerance);
+    }
+  }
+
+  {
+    // Check that pow(<0,y) is correct for noninteger y. This tests special
+    // case handling inside pow().
+    double a = -1.5;
+    J b = MakeJet(-3.14, 3, 0);
+    VL << "a = " << a;
+    VL << "b = " << b;
+
+    J c = pow(a, b);
+    VL << "a^b = " << c;
+    EXPECT_FALSE(IsFinite(c.a));
+    EXPECT_FALSE(IsFinite(c.v[0]));
+    EXPECT_FALSE(IsFinite(c.v[1]));
+  }
+
   { // Check that 1 + x == x + 1.
     J a = x + 1.0;
     J b = 1.0 + x;