Fix bug where pow(JetA,JetB) returned wrong result for JetA==0.

Change-Id: Ife7f6f3aca06322fa414ef747008d2b4dc468a57

include/ceres/jet.h

diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index b476d80..87197ad 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -573,22 +573,57 @@
 }
 
 // pow -- base is a constant, exponent is a differentiable function.
-// (a)^(p+dp) ~= a^p + a^p log(a) dp
+// 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
 template <typename T, int N> inline
 Jet<T, N> pow(double f, const Jet<T, N>& g) {
+  if (f == 0 && g.a > 0) {
+    return Jet<T, N>(T(0.0));
+  }
   T const tmp = pow(f, g.a);
   return Jet<T, N>(tmp, log(f) * tmp * g.v);
 }
 
+// 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).
+//
+// * 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
+//
+// * 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.
+//
+// * 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.
 
-// pow -- both base and exponent are differentiable functions.
-// (a+da)^(b+db) ~= a^b + b * a^(b-1) da + a^b log(a) * db
 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.
+    if (g.a > 1) {
+      return Jet<T, N>(T(0.0));
+    }
+    return Jet<T, N>(T(0.0), f.v);
+  }
+  // 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.
   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);
-
   return Jet<T, N>(tmp1, tmp2 * f.v + tmp3 * g.v);
 }
 

internal/ceres/jet_test.cc

diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index 4598099..394dbb8 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -244,6 +244,82 @@
     ExpectJetsClose(v, u);
   }
 
+  { // Check that pow(0,y) == 1 for y>1, with both arguments Jets.
+    // This tests special case handling inside pow().
+    J a = MakeJet(0, 1, 2);
+    J b = MakeJet(2, 3, 4);
+    VL << "a = " << a;
+    VL << "b = " << b;
+
+    J c = pow(a, b);
+    VL << "a^b = " << c;
+    ExpectJetsClose(c, MakeJet(0, 0, 0));
+  }
+
+  { // Check that pow(0,y) == 1 for y==1, with both arguments Jets.
+    // This tests special case handling inside pow().
+    J a = MakeJet(0, 1, 2);
+    J b = MakeJet(1, 3, 4);
+    VL << "a = " << a;
+    VL << "b = " << b;
+
+    J c = pow(a, b);
+    VL << "a^b = " << c;
+    ExpectJetsClose(c, MakeJet(0, 1, 2));
+  }
+
+  { // Check that pow(0,<1) is not finite, with both arguments Jets.
+    for (int i = 1; i < 10; i++) {
+      J a = MakeJet(0, 1, 2);
+      J b = MakeJet(i*0.1, 3, 4);       // b=0.1 ... 0.9
+      VL << "a = " << a;
+      VL << "b = " << b;
+
+      J c = pow(a, b);
+      VL << "a^b = " << c;
+      EXPECT_EQ(c.a, 0.0);
+      EXPECT_FALSE(IsFinite(c.v[0]));
+      EXPECT_FALSE(IsFinite(c.v[1]));
+    }
+    for (int i = -10; i < 0; i++) {
+      J a = MakeJet(0, 1, 2);
+      J b = MakeJet(i*0.1, 3, 4);       // b=-1,-0.9 ... 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]));
+    }
+    {
+      // The special case of 0^0=1 defined by the C standard.
+      J a = MakeJet(0, 1, 2);
+      J b = MakeJet(0, 3, 4);
+      VL << "a = " << a;
+      VL << "b = " << b;
+
+      J c = pow(a, b);
+      VL << "a^b = " << c;
+      EXPECT_EQ(c.a, 1.0);
+      EXPECT_FALSE(IsFinite(c.v[0]));
+      EXPECT_FALSE(IsFinite(c.v[1]));
+    }
+  }
+
+  { // Check that pow(0,y) == 1 for y==2, with the second argument a Jets.
+    // This tests special case handling inside pow().
+    double a = 0;
+    J b = MakeJet(2, 3, 4);
+    VL << "a = " << a;
+    VL << "b = " << b;
+
+    J c = pow(a, b);
+    VL << "a^b = " << c;
+    ExpectJetsClose(c, MakeJet(0, 0, 0));
+  }
+
   { // Check that 1 + x == x + 1.
     J a = x + 1.0;
     J b = 1.0 + x;