diff --git a/include/ceres/rotation.h b/include/ceres/rotation.h
index 367a7d0..dc04a16 100644
--- a/include/ceres/rotation.h
+++ b/include/ceres/rotation.h
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2023 Google Inc. All rights reserved.
+// Copyright 2025 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -318,27 +318,28 @@
   const T& a0 = angle_axis[0];
   const T& a1 = angle_axis[1];
   const T& a2 = angle_axis[2];
-  const T theta = hypot(a0, a1, a2);
+
+  T k;
 
   // For points not at the origin, the full conversion is numerically stable.
-  if (fpclassify(theta) != FP_ZERO) {
+  if (fpclassify(a0) != FP_ZERO || fpclassify(a1) != FP_ZERO ||
+      fpclassify(a2) != FP_ZERO) {
+    const T theta = hypot(a0, a1, a2);
     const T half_theta = theta * T(0.5);
-    const T k = sin(half_theta) / theta;
+    k = sin(half_theta) / theta;
     quaternion[0] = cos(half_theta);
-    quaternion[1] = a0 * k;
-    quaternion[2] = a1 * k;
-    quaternion[3] = a2 * k;
   } else {
     // At the origin, sqrt() will produce NaN in the derivative since
     // the argument is zero.  By approximating with a Taylor series,
     // and truncating at one term, the value and first derivatives will be
     // computed correctly when Jets are used.
-    const T k(0.5);
+    k = T(0.5);
     quaternion[0] = T(1.0);
-    quaternion[1] = a0 * k;
-    quaternion[2] = a1 * k;
-    quaternion[3] = a2 * k;
   }
+
+  quaternion[1] = a0 * k;
+  quaternion[2] = a1 * k;
+  quaternion[3] = a2 * k;
 }
 
 template <typename T>
@@ -348,11 +349,14 @@
   const T& q1 = quaternion[1];
   const T& q2 = quaternion[2];
   const T& q3 = quaternion[3];
-  const T sin_theta = hypot(q1, q2, q3);
+
+  T k;
 
   // For quaternions representing non-zero rotation, the conversion
   // is numerically stable.
-  if (fpclassify(sin_theta) != FP_ZERO) {
+  if (fpclassify(q1) != FP_ZERO || fpclassify(q2) != FP_ZERO ||
+      fpclassify(q3) != FP_ZERO) {
+    const T sin_theta = hypot(q1, q2, q3);
     const T& cos_theta = quaternion[0];
 
     // If cos_theta is negative, theta is greater than pi/2, which
@@ -368,23 +372,20 @@
     //   theta - pi = atan(sin(theta - pi), cos(theta - pi))
     //              = atan(-sin(theta), -cos(theta))
     //
-    const T two_theta =
-        T(2.0) * ((cos_theta < T(0.0)) ? atan2(-sin_theta, -cos_theta)
-                                       : atan2(sin_theta, cos_theta));
-    const T k = two_theta / sin_theta;
-    angle_axis[0] = q1 * k;
-    angle_axis[1] = q2 * k;
-    angle_axis[2] = q3 * k;
+    const T sign = copysign(T(1), cos_theta);
+    const T two_theta = T(2.0) * atan2(sign * sin_theta, sign * cos_theta);
+    k = two_theta / sin_theta;
   } else {
     // For zero rotation, sqrt() will produce NaN in the derivative since
     // the argument is zero.  By approximating with a Taylor series,
     // and truncating at one term, the value and first derivatives will be
     // computed correctly when Jets are used.
-    const T k(2.0);
-    angle_axis[0] = q1 * k;
-    angle_axis[1] = q2 * k;
-    angle_axis[2] = q3 * k;
+    k = T(2.0);
   }
+
+  angle_axis[0] = q1 * k;
+  angle_axis[1] = q2 * k;
+  angle_axis[2] = q3 * k;
 }
 
 template <typename T>
diff --git a/internal/ceres/autodiff_manifold_test.cc b/internal/ceres/autodiff_manifold_test.cc
index e691da6..f4b6156 100644
--- a/internal/ceres/autodiff_manifold_test.cc
+++ b/internal/ceres/autodiff_manifold_test.cc
@@ -142,12 +142,10 @@
 struct QuaternionFunctor {
   template <typename T>
   bool Plus(const T* x, const T* delta, T* x_plus_delta) const {
-    const T squared_norm_delta =
-        delta[0] * delta[0] + delta[1] * delta[1] + delta[2] * delta[2];
-
     T q_delta[4];
-    if (squared_norm_delta > T(0.0)) {
-      T norm_delta = sqrt(squared_norm_delta);
+    if (fpclassify(delta[0]) != FP_ZERO || fpclassify(delta[1]) != FP_ZERO ||
+        fpclassify(delta[2]) != FP_ZERO) {
+      T norm_delta = hypot(delta[0], delta[1], delta[2]);
       const T sin_delta_by_delta = sin(norm_delta) / norm_delta;
       q_delta[0] = cos(norm_delta);
       q_delta[1] = sin_delta_by_delta * delta[0];
@@ -173,10 +171,11 @@
     T minus_x[4] = {x[0], -x[1], -x[2], -x[3]};
     T ambient_y_minus_x[4];
     QuaternionProduct(y, minus_x, ambient_y_minus_x);
-    T u_norm = sqrt(ambient_y_minus_x[1] * ambient_y_minus_x[1] +
-                    ambient_y_minus_x[2] * ambient_y_minus_x[2] +
-                    ambient_y_minus_x[3] * ambient_y_minus_x[3]);
-    if (u_norm > 0.0) {
+    if (fpclassify(ambient_y_minus_x[1]) != FP_ZERO ||
+        fpclassify(ambient_y_minus_x[2]) != FP_ZERO ||
+        fpclassify(ambient_y_minus_x[3]) != FP_ZERO) {
+      T u_norm = hypot(
+          ambient_y_minus_x[1], ambient_y_minus_x[2], ambient_y_minus_x[3]);
       T theta = atan2(u_norm, ambient_y_minus_x[0]);
       y_minus_x[0] = theta * ambient_y_minus_x[1] / u_norm;
       y_minus_x[1] = theta * ambient_y_minus_x[2] / u_norm;
diff --git a/internal/ceres/is_close.cc b/internal/ceres/is_close.cc
index 575918b..c652b68 100644
--- a/internal/ceres/is_close.cc
+++ b/internal/ceres/is_close.cc
@@ -1,5 +1,5 @@
 // Ceres Solver - A fast non-linear least squares minimizer
-// Copyright 2023 Google Inc. All rights reserved.
+// Copyright 2025 Google Inc. All rights reserved.
 // http://ceres-solver.org/
 //
 // Redistribution and use in source and binary forms, with or without
@@ -41,18 +41,19 @@
              double* absolute_error) {
   double local_absolute_error;
   double local_relative_error;
-  if (!absolute_error) {
+  if (absolute_error == nullptr) {
     absolute_error = &local_absolute_error;
   }
-  if (!relative_error) {
+  if (relative_error == nullptr) {
     relative_error = &local_relative_error;
   }
   *absolute_error = std::fabs(x - y);
-  *relative_error = *absolute_error / std::max(std::fabs(x), std::fabs(y));
-  if (x == 0 || y == 0) {
+  if (std::fpclassify(x) == FP_ZERO || std::fpclassify(y) == FP_ZERO) {
     // If x or y is exactly zero, then relative difference doesn't have any
     // meaning. Take the absolute difference instead.
     *relative_error = *absolute_error;
+  } else {
+    *relative_error = *absolute_error / std::max(std::fabs(x), std::fabs(y));
   }
   return *relative_error < std::fabs(relative_precision);
 }
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc
index a9de40f..67e5580 100644
--- a/internal/ceres/jet_test.cc
+++ b/internal/ceres/jet_test.cc
@@ -103,13 +103,15 @@
 
   Fenv env;  // Do not leak floating-point exceptions to the caller
   double absolute_difference = std::abs(x - y);
-  double relative_difference =
-      absolute_difference / std::max(std::abs(x), std::abs(y));
+  double relative_difference;
 
   if (std::fpclassify(x) == FP_ZERO || std::fpclassify(y) == FP_ZERO) {
     // If x or y is exactly zero, then relative difference doesn't have any
     // meaning. Take the absolute difference instead.
     relative_difference = absolute_difference;
+  } else {
+    relative_difference =
+        absolute_difference / std::max(std::abs(x), std::abs(y));
   }
   return std::islessequal(relative_difference, max_abs_relative_difference);
 }
diff --git a/internal/ceres/polynomial.cc b/internal/ceres/polynomial.cc
index db3a404..a0596c6 100644
--- a/internal/ceres/polynomial.cc
+++ b/internal/ceres/polynomial.cc
@@ -70,26 +70,30 @@
     scaling_has_changed = false;
 
     for (int i = 0; i < degree; ++i) {
-      const double row_norm = companion_matrix_offdiagonal.row(i).lpNorm<1>();
       const double col_norm = companion_matrix_offdiagonal.col(i).lpNorm<1>();
 
-      // Decompose row_norm/col_norm into mantissa * 2^exponent,
-      // where 0.5 <= mantissa < 1. Discard mantissa (return value
-      // of frexp), as only the exponent is needed.
-      int exponent = 0;
-      std::frexp(row_norm / col_norm, &exponent);
-      exponent /= 2;
+      // Avoid division by zero
+      if (std::fpclassify(col_norm) != FP_ZERO) {
+        const double row_norm = companion_matrix_offdiagonal.row(i).lpNorm<1>();
+        // Decompose row_norm/col_norm into mantissa * 2^exponent,
+        // where 0.5 <= mantissa < 1. Discard mantissa (return value
+        // of frexp), as only the exponent is needed.
+        int exponent = 0;
+        std::frexp(row_norm / col_norm, &exponent);
+        exponent /= 2;
 
-      if (exponent != 0) {
-        const double scaled_col_norm = std::ldexp(col_norm, exponent);
-        const double scaled_row_norm = std::ldexp(row_norm, -exponent);
-        if (scaled_col_norm + scaled_row_norm < gamma * (col_norm + row_norm)) {
-          // Accept the new scaling. (Multiplication by powers of 2 should not
-          // introduce rounding errors (ignoring non-normalized numbers and
-          // over- or underflow))
-          scaling_has_changed = true;
-          companion_matrix_offdiagonal.row(i) *= std::ldexp(1.0, -exponent);
-          companion_matrix_offdiagonal.col(i) *= std::ldexp(1.0, exponent);
+        if (exponent != 0) {
+          const double scaled_col_norm = std::ldexp(col_norm, exponent);
+          const double scaled_row_norm = std::ldexp(row_norm, -exponent);
+          if (scaled_col_norm + scaled_row_norm <
+              gamma * (col_norm + row_norm)) {
+            // Accept the new scaling. (Multiplication by powers of 2 should not
+            // introduce rounding errors (ignoring non-normalized numbers and
+            // over- or underflow))
+            scaling_has_changed = true;
+            companion_matrix_offdiagonal.row(i) *= std::ldexp(1.0, -exponent);
+            companion_matrix_offdiagonal.col(i) *= std::ldexp(1.0, exponent);
+          }
         }
       }
     }
diff --git a/internal/ceres/test_util.cc b/internal/ceres/test_util.cc
index fc73764..9ae661e 100644
--- a/internal/ceres/test_util.cc
+++ b/internal/ceres/test_util.cc
@@ -64,11 +64,14 @@
   }
 
   double absolute_difference = fabs(x - y);
-  double relative_difference = absolute_difference / std::max(fabs(x), fabs(y));
-  if (x == 0 || y == 0) {
+  double relative_difference;
+  if (std::fpclassify(x) == FP_ZERO || std::fpclassify(y) == FP_ZERO) {
     // If x or y is exactly zero, then relative difference doesn't have any
     // meaning. Take the absolute difference instead.
     relative_difference = absolute_difference;
+  } else {
+    relative_difference =
+        absolute_difference / std::max(std::fabs(x), std::fabs(y));
   }
   if (relative_difference > max_abs_relative_difference) {
     VLOG(1) << absl::StrFormat("x=%17g y=%17g abs=%17g rel=%17g",