Streamline branching in Jet copysign and fdim

- Pass the second argument g of copysign(f, g) by const reference
  instead of by value, and return (sa * sb) * f.v directly when g != 0
  to avoid evaluating 0 * g.v.
- Check isgreater(f, g) before isnan in fdim so that normal inputs with
  f > g short-circuit on the first comparison.

Benchmark (Apple M1 Pro):
----------------------------------------------------------------------
Benchmark                                     Before    After   Change
----------------------------------------------------------------------
BM_JetCopysign<3>                            2.27 ns 0.453 ns   -80.0%
BM_JetCopysign<6>                            2.33 ns 0.506 ns   -78.3%
BM_JetCopysign<12>                           4.08 ns  1.23 ns   -69.9%
BM_JetCopysign<30>                           9.76 ns  5.23 ns   -46.4%
BM_JetFdim<3>                                1.98 ns 0.728 ns   -63.2%
BM_JetFdim<6>                                2.07 ns  1.25 ns   -39.6%
BM_JetFdim<12>                               2.78 ns  1.83 ns   -34.2%
BM_JetFdim<30>                               5.28 ns  4.86 ns    -8.0%
----------------------------------------------------------------------

Change-Id: I3b4ceebad76cbb0eddc42be6368872a9e5543b71

include/ceres/jet.h

diff --git a/include/ceres/jet.h b/include/ceres/jet.h
index 312fc49..d2ded9b 100644
--- a/include/ceres/jet.h
+++ b/include/ceres/jet.h
@@ -623,18 +623,19 @@
 //
 // where δ(b) is the Dirac delta function.
 template <typename T, int N>
-EIGEN_STRONG_INLINE Jet<T, N> copysign(const Jet<T, N>& f, const Jet<T, N> g) {
-  // The Dirac delta function  δ(b) is undefined at b=0 (here it's
-  // infinite) and 0 everywhere else.
-  T d = fpclassify(g) == FP_ZERO ? std::numeric_limits<T>::infinity() : T(0);
-  T sa = copysign(T(1), f.a);  // sgn(a)
-  T sb = copysign(T(1), g.a);  // sgn(b)
-  // The second part of the infinitesimal is 2|a|δ(b) which is either infinity
-  // or 0 unless a or any of the values of the b infinitesimal are 0. In the
-  // latter case, the corresponding values become NaNs (multiplying 0 by
-  // infinity gives NaN). We drop the constant factor 2 since it does not change
-  // the result (its values will still be either 0, infinity or NaN).
-  return Jet<T, N>(copysign(f.a, g.a), sa * sb * f.v + abs(f.a) * d * g.v);
+EIGEN_STRONG_INLINE Jet<T, N> copysign(const Jet<T, N>& f,
+                                       const Jet<T, N>& g) {
+  const T sa = copysign(T(1.0), f.a);  // sgn(a)
+  const T sb = copysign(T(1.0), g.a);  // sgn(b)
+  if (fpclassify(g) != FP_ZERO) {
+    return Jet<T, N>(copysign(f.a, g.a), (sa * sb) * f.v);
+  }
+  // The Dirac delta function δ(b) is undefined at b=0 (here it's infinite) and
+  // 0 everywhere else. When a or any of the values of the b infinitesimal are
+  // 0, the corresponding values become NaNs (multiplying 0 by infinity gives
+  // NaN). We drop the constant factor 2 since it does not change the result.
+  const T d = std::numeric_limits<T>::infinity();
+  return Jet<T, N>(copysign(f.a, g.a), (sa * sb) * f.v + (abs(f.a) * d) * g.v);
 }
 
 // log(a + h) ~= log(a) + h / a
@@ -923,10 +924,15 @@
           std::enable_if_t<CompatibleJetOperands_v<Lhs, Rhs>>* = nullptr>
 EIGEN_STRONG_INLINE decltype(auto) fdim(const Lhs& f, const Rhs& g) {
   using J = std::common_type_t<Lhs, Rhs>;
+  // Check isgreater first (which is quiet on NaNs) so that the common f > g
+  // case returns immediately without checking isnan(f) or isnan(g).
+  if (isgreater(f, g)) {
+    return J{f - g};
+  }
   if (isnan(f) || isnan(g)) {
     return std::numeric_limits<J>::quiet_NaN();
   }
-  return isgreater(f, g) ? J{f - g} : J{};
+  return J{};
 }
 
 // erf is defined as an integral that cannot be expressed analytically