Fix bug where pow(JetA,JetB) returned wrong result for JetA==0.
Change-Id: Ife7f6f3aca06322fa414ef747008d2b4dc468a57
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);
}
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;
