Precompute scalar derivatives in Jet elementary functions Compute cbrt(f) using a single cbrt(f.a) call with f.v / (3 * r * r), avoiding both the second cbrt evaluation and f.a * f.a overflow or underflow at extreme magnitudes. Precompute the scalar derivative factor in erf(x) and erfc(x) using a compile-time constant for 2 / sqrt(pi) before multiplying x.v so Eigen evaluates a single scalar-vector product instead of scaling x.v twice. Parenthesize (T(1.0) / log(T(10.0))) / f.a in log10 and (T(1.0) / log(T(2.0))) / f.a in log2 so the reciprocal constant folds at compile time, and use T(0.5) / tmp in sqrt. Standardize on floating-point T(...) constants and left-multiplication by scalar derivatives across elementary functions in jet.h. Benchmark (Apple M1 Pro): ---------------------------------------------------------------------- Benchmark Before After Change ---------------------------------------------------------------------- Cbrt<3>/1000 11.7 us 7.29 us -37.7% Cbrt<10>/1000 15.9 us 11.5 us -27.6% Cbrt<32>/1000 29.9 us 25.7 us -14.0% Erf<3>/1000 23.8 us 21.1 us -11.3% Erf<32>/1000 39.9 us 32.8 us -17.8% Erfc<3>/1000 23.1 us 20.2 us -12.6% Erfc<32>/1000 39.8 us 32.7 us -17.8% Sqrt<15>/1000 4.50 us 4.38 us -2.7% Sqrt<32>/1000 12.5 us 12.4 us -0.8% ---------------------------------------------------------------------- Change-Id: I2a4ceebad76cbb0eddc42be6368872a9e5543b70
diff --git a/include/ceres/constants.h b/include/ceres/constants.h index e4c96da..6c6442d 100644 --- a/include/ceres/constants.h +++ b/include/ceres/constants.h
@@ -34,12 +34,13 @@ // TODO(HSHelson): This header should no longer be necessary once C++20's // <numbers> (e.g. std::numbers::pi_v) becomes usable -// The constant is computed using mpmath by the following shell command: +// The constants are computed using mpmath by the following shell command: // // python3 - <<EOF // from mpmath import mp // mp.dps = 64 // print(mp.pi) +// print(2 / mp.sqrt(mp.pi)) // EOF namespace ceres::constants { // π @@ -47,6 +48,11 @@ inline constexpr T pi_v( 3.141592653589793238462643383279502884197169399375105820974944592L); inline constexpr double pi = pi_v<double>; +// 2 / √π +template <typename T> +inline constexpr T two_over_sqrt_pi_v( + 1.128379167095512573896158903121545171688101258657997713688171443L); +inline constexpr double two_over_sqrt_pi = two_over_sqrt_pi_v<double>; } // namespace ceres::constants #endif // CERES_PUBLIC_CONSTANTS_H_
diff --git a/include/ceres/jet.h b/include/ceres/jet.h index c526148..312fc49 100644 --- a/include/ceres/jet.h +++ b/include/ceres/jet.h
@@ -167,6 +167,7 @@ #include "Eigen/Core" #include "ceres/accurate_norm.h" +#include "ceres/constants.h" #include "ceres/internal/jet_traits.h" #include "ceres/internal/port.h" #include "ceres/jet_fwd.h" @@ -532,7 +533,7 @@ // abs(x + h) ~= abs(x) + sgn(x)h template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> abs(const Jet<T, N>& f) { - return Jet<T, N>(abs(f.a), copysign(T(1), f.a) * f.v); + return Jet<T, N>(abs(f.a), copysign(T(1.0), f.a) * f.v); } // AccurateNorm(x, y, ...) computes the Euclidean norm of its arguments while @@ -640,22 +641,25 @@ template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> log(const Jet<T, N>& f) { const T a_inverse = T(1.0) / f.a; - return Jet<T, N>(log(f.a), f.v * a_inverse); + return Jet<T, N>(log(f.a), a_inverse * f.v); } // log10(a + h) ~= log10(a) + h / (a log(10)) template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> log10(const Jet<T, N>& f) { - // Most compilers will expand log(10) to a constant. - const T a_inverse = T(1.0) / (f.a * log(T(10.0))); - return Jet<T, N>(log10(f.a), f.v * a_inverse); + // Parenthesizing (1 / log(10)) allows the compiler to fold the reciprocal + // constant at compile time and execute a single division at runtime instead + // of a runtime multiplication and division (which IEEE 754 rules forbid + // reassociating). + const T a_inverse = (T(1.0) / log(T(10.0))) / f.a; + return Jet<T, N>(log10(f.a), a_inverse * f.v); } // log1p(a + h) ~= log1p(a) + h / (1 + a) template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> log1p(const Jet<T, N>& f) { const T a_inverse = T(1.0) / (T(1.0) + f.a); - return Jet<T, N>(log1p(f.a), f.v * a_inverse); + return Jet<T, N>(log1p(f.a), a_inverse * f.v); } // exp(a + h) ~= exp(a) + exp(a) h @@ -677,8 +681,10 @@ template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> sqrt(const Jet<T, N>& f) { const T tmp = sqrt(f.a); - const T two_a_inverse = T(1.0) / (T(2.0) * tmp); - return Jet<T, N>(tmp, f.v * two_a_inverse); + // Dividing 0.5 by tmp avoids a runtime multiplication by 2 before the + // division. + const T two_a_inverse = T(0.5) / tmp; + return Jet<T, N>(tmp, two_a_inverse * f.v); } // cos(a + h) ~= cos(a) - sin(a) h @@ -765,23 +771,28 @@ // cbrt(a + h) ~= cbrt(a) + h / (3 a ^ (2/3)) template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> cbrt(const Jet<T, N>& f) { - const T derivative = T(1.0) / (T(3.0) * cbrt(f.a * f.a)); - return Jet<T, N>(cbrt(f.a), f.v * derivative); + // Since a^(2/3) = cbrt(a)^2 for all real a, reusing r = cbrt(f.a) avoids a + // second cbrt evaluation as well as intermediate overflow or underflow of + // f.a * f.a when |f.a| > sqrt(max) or |f.a| < sqrt(min). + const T r = cbrt(f.a); + return Jet<T, N>(r, f.v / (T(3.0) * r * r)); } // exp2(x + h) = 2^(x+h) ~= 2^x + h*2^x*log(2) template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> exp2(const Jet<T, N>& f) { const T tmp = exp2(f.a); - const T derivative = tmp * log(T(2)); - return Jet<T, N>(tmp, f.v * derivative); + const T derivative = tmp * log(T(2.0)); + return Jet<T, N>(tmp, derivative * f.v); } // log2(x + h) ~= log2(x) + h / (x * log(2)) template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> log2(const Jet<T, N>& f) { - const T derivative = T(1.0) / (f.a * log(T(2))); - return Jet<T, N>(log2(f.a), f.v * derivative); + // Parenthesizing (1 / log(2)) folds the reciprocal constant at compile time + // and replaces a runtime multiplication and division with a single division. + const T derivative = (T(1.0) / log(T(2.0))) / f.a; + return Jet<T, N>(log2(f.a), derivative * f.v); } // Like sqrt(x^2 + y^2), @@ -923,23 +934,18 @@ // erf(x + h) = erf(x) + h * 2*exp(-x^2)/sqrt(pi) template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> erf(const Jet<T, N>& x) { - // We evaluate the constant as follows: - // 2 / sqrt(pi) = 1 / sqrt(atan(1.)) - // On POSIX systems it is defined as M_2_SQRTPI, but this is not - // portable and the type may not be T. The above expression - // evaluates to full precision with IEEE arithmetic and, since it's - // constant, the compiler can generate exactly the same code. gcc - // does so even at -O0. - return Jet<T, N>(erf(x.a), x.v * exp(-x.a * x.a) * (T(1) / sqrt(atan(T(1))))); + // Precompute the scalar derivative before multiplying x.v so Eigen evaluates + // a single scalar-vector product instead of scaling x.v twice. + const T derivative = exp(-x.a * x.a) * constants::two_over_sqrt_pi_v<T>; + return Jet<T, N>(erf(x.a), derivative * x.v); } // erfc(x) = 1-erf(x) // erfc(x + h) = erfc(x) + h * (-2*exp(-x^2)/sqrt(pi)) template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> erfc(const Jet<T, N>& x) { - // See in erf() above for the evaluation of the constant in the derivative. - return Jet<T, N>(erfc(x.a), - -x.v * exp(-x.a * x.a) * (T(1) / sqrt(atan(T(1))))); + const T derivative = -exp(-x.a * x.a) * constants::two_over_sqrt_pi_v<T>; + return Jet<T, N>(erfc(x.a), derivative * x.v); } #if defined(CERES_HAS_CPP17_BESSEL_FUNCTIONS) || \ @@ -1224,7 +1230,7 @@ const Jet<T, N>& b, const Jet<T, N>& t) { return Jet<T, N>{lerp(a.a, b.a, t.a), - (T(1) - t.a) * a.v + t.a * b.v + (b.a - a.a) * t.v}; + (T(1.0) - t.a) * a.v + t.a * b.v + (b.a - a.a) * t.v}; } // Computes the midpoint a + (b - a) / 2. @@ -1260,7 +1266,7 @@ // f = a + da // g = b + db - T const tmp = T(1.0) / (f.a * f.a + g.a * g.a); + const T tmp = T(1.0) / (f.a * f.a + g.a * g.a); return Jet<T, N>(atan2(g.a, f.a), tmp * (-g.a * f.v + f.a * g.v)); } @@ -1277,7 +1283,7 @@ // norm(x + h) ~= norm(x) + 2x h template <typename T, int N> EIGEN_STRONG_INLINE Jet<T, N> norm(const Jet<T, N>& f) { - return Jet<T, N>(norm(f.a), T(2) * f.a * f.v); + return Jet<T, N>(norm(f.a), (T(2.0) * f.a) * f.v); } // pow -- base is a differentiable function, exponent is a constant.
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc index a1d2c42..bc982d1 100644 --- a/internal/ceres/jet_test.cc +++ b/internal/ceres/jet_test.cc
@@ -1040,6 +1040,22 @@ } } +TEST(Jet, CbrtExtremeMagnitudes) { + // Verify cbrt avoids f.a * f.a overflow and underflow. + { + J z = cbrt(MakeJet(1e200, 1.0, -2.0)); + const double r = std::cbrt(1e200); + const double d = 3.0 * r * r; + EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(r, 1.0 / d, -2.0 / d))); + } + { + J z = cbrt(MakeJet(-1e-200, 1.0, -2.0)); + const double r = std::cbrt(-1e-200); + const double d = 3.0 * r * r; + EXPECT_THAT(z, IsAlmostEqualTo(MakeJet(r, 1.0 / d, -2.0 / d))); + } +} + TEST(Jet, JetsInEigenMatrices) { J x = MakeJet(2.3, -2.7, 1e-3); J y = MakeJet(1.7, 0.5, 1e+2);