diff --git a/include/ceres/rotation.h b/include/ceres/rotation.h
index e014e04..397787d 100644
--- a/include/ceres/rotation.h
+++ b/include/ceres/rotation.h
@@ -395,7 +395,7 @@
     const MatrixAdapter<T, row_stride, col_stride>& R) {
   static const T kOne = T(1.0);
   const T theta2 = DotProduct(angle_axis, angle_axis);
-  if (theta2 > 0.0) {
+  if (theta2 > T(std::numeric_limits<double>::epsilon())) {
     // We want to be careful to only evaluate the square root if the
     // norm of the angle_axis vector is greater than zero. Otherwise
     // we get a division by zero.
@@ -419,13 +419,13 @@
   } else {
     // At zero, we switch to using the first order Taylor expansion.
     R(0, 0) =  kOne;
-    R(1, 0) = -angle_axis[2];
-    R(2, 0) =  angle_axis[1];
-    R(0, 1) =  angle_axis[2];
+    R(1, 0) =  angle_axis[2];
+    R(2, 0) = -angle_axis[1];
+    R(0, 1) = -angle_axis[2];
     R(1, 1) =  kOne;
-    R(2, 1) = -angle_axis[0];
-    R(0, 2) = -angle_axis[1];
-    R(1, 2) =  angle_axis[0];
+    R(2, 1) =  angle_axis[0];
+    R(0, 2) =  angle_axis[1];
+    R(1, 2) = -angle_axis[0];
     R(2, 2) = kOne;
   }
 }
@@ -581,7 +581,7 @@
 template<typename T> inline
 void AngleAxisRotatePoint(const T angle_axis[3], const T pt[3], T result[3]) {
   const T theta2 = DotProduct(angle_axis, angle_axis);
-  if (theta2 > 0.0) {
+  if (theta2 > T(std::numeric_limits<double>::epsilon())) {
     // Away from zero, use the rodriguez formula
     //
     //   result = pt costheta +
@@ -625,9 +625,9 @@
     // and actually performing multiplication with the point pt, gives us
     // R * pt = pt + w x pt.
     //
-    // Switching to the Taylor expansion at zero helps avoid all sorts
-    // of numerical nastiness.
-
+    // Switching to the Taylor expansion at zero provides meaningful
+    // derivatives when evaluated using Jets.
+    //
     // Explicitly inlined evaluation of the cross product for
     // performance reasons.
     const T w_cross_pt[3] = { angle_axis[1] * pt[2] - angle_axis[2] * pt[1],
diff --git a/internal/ceres/rotation_test.cc b/internal/ceres/rotation_test.cc
index 8de1bbd..cc4cc92 100644
--- a/internal/ceres/rotation_test.cc
+++ b/internal/ceres/rotation_test.cc
@@ -548,6 +548,40 @@
   }
 }
 
+// Takes a bunch of random axis/angle values near zero, converts them
+// to rotation matrices, and back again.
+TEST(Rotation, AngleAxisToRotationMatrixAndBackNearZero) {
+  srand(5);
+  for (int i = 0; i < kNumTrials; i++) {
+    double axis_angle[3];
+    // Make an axis by choosing three random numbers in [-1, 1) and
+    // normalizing.
+    double norm = 0;
+    for (int i = 0; i < 3; i++) {
+      axis_angle[i] = RandDouble() * 2 - 1;
+      norm += axis_angle[i] * axis_angle[i];
+    }
+    norm = sqrt(norm);
+
+    // Tiny theta.
+    double theta = 1e-16 * (kPi * 2 * RandDouble() - kPi);
+    for (int i = 0; i < 3; i++) {
+      axis_angle[i] = axis_angle[i] * theta / norm;
+    }
+
+    double matrix[9];
+    double round_trip[3];
+    AngleAxisToRotationMatrix(axis_angle, matrix);
+    ASSERT_THAT(matrix, IsOrthonormal());
+    RotationMatrixToAngleAxis(matrix, round_trip);
+
+    for (int i = 0; i < 3; ++i) {
+      EXPECT_NEAR(round_trip[i], axis_angle[i], kTolerance);
+    }
+  }
+}
+
+
 // Transposes a 3x3 matrix.
 static void Transpose3x3(double m[9]) {
   std::swap(m[1], m[3]);