Add AccurateRNorm Composing two-argument functions for variadic Jet inputs can produce invalid derivatives because intermediate results do not preserve the full derivative relationship. Introduce overflow and underflow safe Euclidean reciprocal norm with variadic overloads. The Jet overloads compute derivatives from all arguments using the complete chain rule. Change-Id: Id47d4e7a56929548484175cf7a185b351771a67f
diff --git a/docs/source/version_history.rst b/docs/source/version_history.rst index 7217e5e..a987a6b 100644 --- a/docs/source/version_history.rst +++ b/docs/source/version_history.rst
@@ -23,10 +23,10 @@ #. Mixed precision solves + iterative refinement when using ``SUITE_SPARSE`` as the sparse linear algebra library (requires SuiteSparse 7.4.0 or later for single-precision ``CHOLMOD``). (Sameer Agarwal) -#. Added ``AccurateNorm`` for overflow and underflow safe variadic norms - accurate to nearly full working precision. Its Jet derivatives are defined - wherever the norm is nonzero, unlike those of nested two-argument norms. - (Sergiu Deitsch) +#. Added ``AccurateNorm`` and ``AccurateRNorm`` for overflow and underflow safe + variadic norms accurate to nearly full working precision. Their Jet + derivatives are defined wherever the norm is nonzero, unlike those of nested + two-argument norms. (Sergiu Deitsch) #. Simplify instantiation of cost functions and their functors by forwarding functor constructor arguments directly, and add ``std::unique_ptr`` and :enum:`Ownership` constructors across cost and first-order functions.
diff --git a/include/ceres/accurate_norm.h b/include/ceres/accurate_norm.h index 2ace148..37fbf1b 100644 --- a/include/ceres/accurate_norm.h +++ b/include/ceres/accurate_norm.h
@@ -29,26 +29,32 @@ // Author: sergiu.deitsch@gmail.com (Sergiu Deitsch) // // This header implements functions for accurately computing the Euclidean norm -// of two or more arguments while avoiding underflow and overflow. +// of two or more arguments and its reciprocal while avoiding underflow and +// overflow. // -// The functions accumulate the squares of the arguments as an unevaluated sum +// Both functions accumulate the squares of the arguments as an unevaluated sum // of the rounded sum and its rounding error, and correct the square root of -// the rounded sum by a single Newton step that accounts for the rounding -// errors. They share the same implementation for any number of arguments: the -// result is computed without rescaling if the largest magnitude is within a -// range where rescaling cannot change the result, and otherwise after -// rescaling all arguments by a fixed radix power. +// the rounded sum, or its reciprocal, by a single Newton step that accounts for +// the rounding errors. They share the same implementation for any number of +// arguments: the result is computed without rescaling if the largest magnitude +// is within a range where rescaling cannot change the result, and otherwise +// after rescaling all arguments by a fixed radix power. // // Unlike the 2-argument algorithm in [1], the functions do not return the // larger argument if the smaller one is negligible. Neglecting arguments does // not generalize to more arguments since the errors of several neglected -// squares accumulate. +// squares accumulate, and the reciprocal of the larger argument can differ from +// the correctly rounded reciprocal norm. // -// The implementation is derived from the following paper: +// The implementation is derived from the following two papers: // // [1] Borges, C. F. (2021). Algorithm 1014: An Improved Algorithm for // hypot(x,y). ACM Transactions on Mathematical Software, 47(1), 1–12. // https://doi.org/10.1145/3428446 +// +// [2] Borges, C. F. (2021). Fast Compensated Algorithms for the Reciprocal +// Square Root, the Reciprocal Hypotenuse, and Givens Rotations. +// http://arxiv.org/abs/2103.08694 #ifndef CERES_PUBLIC_ACCURATE_NORM_H_ #define CERES_PUBLIC_ACCURATE_NORM_H_ @@ -169,6 +175,16 @@ return (std::numeric_limits<T>::min_exponent - 1) / 2 - std::numeric_limits<T>::digits + 1; } + + // Upper bound of the range for the reciprocal norm of num_arguments + // arguments, replacing UnscaledMaximum(num_arguments). The bound is tighter + // because the reciprocal norm must be at least Tiny() to keep the rounding + // error of its square exact. Larger arguments make the rounding error inexact + // even though the sum of their squares does not overflow. The norm is at most + // the largest magnitude times 𝛽^⌈⌈log₂(num_arguments)⌉/2⌉ ≥ √num_arguments. + static constexpr T ReciprocalUnscaledMaximum(int num_arguments) noexcept { + return T{1} / (Tiny() * PowerOfTwo<T>((CeilLog2(num_arguments) + 1) / 2)); + } }; // Determines the largest magnitude of the arguments. An infinite argument @@ -249,6 +265,42 @@ return fma(tau / h, T(0.5), h); } +// Computes 1 / sqrt(x^2 + y^2 + ...) without checking the arguments. In +// addition to the requirements of UnscaledAccurateNorm, the reciprocal norm +// must be at least Tiny() to keep the rounding error of its square exact. Not +// intended to be invoked by users. +template <typename T, typename... Args> +inline T UnscaledAccurateRNorm(T x, T y, Args... args) { + using std::fma; + using std::sqrt; + + // In the following, σ denotes the sum of squares and σ_e its rounding error. + const auto [sum_of_squares, sum_of_squares_error] = + UnscaledAccurateSquaredNormWithError(x, y, args...); + const T r = T(1) / sum_of_squares; + // Rounding error of the reciprocal, 1 − r⋅(σ + σ_e) + const T residual = + fma(-r, sum_of_squares_error, fma(-r, sum_of_squares, T(1))); + const T rho = sqrt(r); + // Rounding error of the square root, r − ρ² + const T tau = fma(-rho, rho, r); + // To solve 1/ρ² = σ, Newton's correction term g/dg is + // + // ρ⋅(ρ²⋅σ − 1) + // δ_ρ = ──────────── . + // 2 + // + // The update ρ − δ_ρ adds its negation, ρ⋅(1 − ρ²⋅σ) / 2, through the final + // fused multiply-add. Both rounding errors contribute to + // + // 1 − ρ²⋅σ = (1 − r⋅σ) + σ⋅(r − ρ²). + // + // The update combines both terms as in the compensated reciprocal square + // root of [2], here applied to the sum of squares. + const T nu = fma(sum_of_squares, tau, residual) / T(2); + return fma(rho, nu, rho); +} + } // namespace internal // Computes the Euclidean norm of two or more floating-point values of the same @@ -309,6 +361,77 @@ return AccurateNorm(PromotedType(a), PromotedType(b), PromotedType(args)...); } +// Computes the reciprocal of the Euclidean norm of two or more floating-point +// values of the same type while avoiding intermediate underflow and overflow. +// An infinite argument produces zero, even if another argument is NaN. +// Otherwise, a NaN argument produces a NaN with its payload preserved. When all +// arguments are zero, the result is NaN. +template <typename T, typename... Args> +inline auto AccurateRNorm(T a, T b, Args... args) + -> std::enable_if_t<std::is_floating_point_v<T> && + (std::is_same_v<T, Args> && ...), + T> { + using std::fpclassify; + using std::isinf; + using std::isnan; + + const T maximum = internal::MaximumMagnitude(a, b, args...); + + if (isinf(maximum)) { + return 0; + } + + if (isnan(maximum)) { + return maximum; + } + + if (fpclassify(maximum) == FP_ZERO) { + return std::numeric_limits<T>::quiet_NaN(); + } + + using internal::AccurateNormTraits; + using internal::PowerOfTwo; + using internal::UnscaledAccurateRNorm; + + constexpr int num_arguments = 2 + sizeof...(Args); + // Multiplying by a radix power rounds exactly as std::scalbn does. + constexpr int exponent = AccurateNormTraits<T>::ScaleExponent(); + constexpr T down = PowerOfTwo<T>(exponent); + constexpr T up = PowerOfTwo<T>(-exponent); + + // Rescale as in AccurateNorm, but apply the scale to the reciprocal instead + // of its inverse. The reciprocal of the rescaled arguments cannot overflow or + // underflow because the rescaled largest magnitude is bounded. + if (maximum > + AccurateNormTraits<T>::ReciprocalUnscaledMaximum(num_arguments)) { + return UnscaledAccurateRNorm(a * down, b * down, args * down...) * down; + } + + if (maximum < AccurateNormTraits<T>::UnscaledMinimum()) { + // Only a largest magnitude below Tiny() requires the larger power to make + // the rounding errors of the squares of subnormal arguments exact. Scaling + // a larger magnitude by that power could exceed + // ReciprocalUnscaledMaximum(). Scaling by 𝛽^p instead suffices to reach + // UnscaledMinimum(). + const T scale = maximum < AccurateNormTraits<T>::Tiny() + ? up + : PowerOfTwo<T>(std::numeric_limits<T>::digits); + return UnscaledAccurateRNorm(a * scale, b * scale, args * scale...) * scale; + } + + return UnscaledAccurateRNorm(a, b, args...); +} + +// Computes the reciprocal of the Euclidean norm of two or more arithmetic +// values after promoting all arguments to a common floating-point type. +template <typename T, typename U, typename... Args> +inline internal::Promote_t<T, U, Args...> AccurateRNorm(T a, + U b, + Args... args) { + using PromotedType = internal::Promote_t<T, U, Args...>; + return AccurateRNorm(PromotedType(a), PromotedType(b), PromotedType(args)...); +} + } // namespace ceres #endif // CERES_PUBLIC_ACCURATE_NORM_H_
diff --git a/include/ceres/jet.h b/include/ceres/jet.h index fd19f6f..50762e3 100644 --- a/include/ceres/jet.h +++ b/include/ceres/jet.h
@@ -559,6 +559,41 @@ ((x.a / tmp * x.v + y.a / tmp * y.v) + ... + (args.a / tmp * args.v))); } +// AccurateRNorm(x, y, ...) computes the reciprocal of the Euclidean norm of its +// arguments while avoiding underflow and overflow. Therefore, the function can +// be formally defined as +// +// AccurateRNorm(x, y, ...) = 1 / sqrt(x^2 + y^2 + ...) +// +// where +// +// d/dx 1 / sqrt(x^2 + y^2 + ...) = +// -x / (x^2 + y^2 + ...)^3/2 +// d/dy 1 / sqrt(x^2 + y^2 + ...) = +// -y / (x^2 + y^2 + ...)^3/2 +// ... +// +// with the dual representation given by +// +// AccurateRNorm(x + dx, y + dy, ...) ~= +// 1 / sqrt(x^2 + y^2 + ...) - +// (x dx + y dy + ...) / (x^2 + y^2 + ...)^3/2 +// +// The derivatives are undefined when all arguments are zero. +template <typename T, int N, typename... Args> +inline auto AccurateRNorm(const Jet<T, N>& x, + const Jet<T, N>& y, + const Args&... args) + -> std::enable_if_t<(std::is_same_v<Args, Jet<T, N>> && ...), Jet<T, N>> { + const T tmp = AccurateRNorm(x.a, y.a, args.a...); + // Multiplying the argument first avoids overflowing the cube of a large + // reciprocal norm. + return Jet<T, N>( + tmp, + ((-(x.a * tmp * tmp * tmp) * x.v - y.a * tmp * tmp * tmp * y.v) - ... - + (args.a * tmp * tmp * tmp * args.v))); +} + // copysign(a, b) composes a float with the magnitude of a and the sign of b. // Therefore, the function can be formally defined as //
diff --git a/include/ceres/rotation.h b/include/ceres/rotation.h index eede7d2..a5a205a 100644 --- a/include/ceres/rotation.h +++ b/include/ceres/rotation.h
@@ -824,8 +824,7 @@ // 'scale' is 1 / norm(q). const T scale = - T(1) / sqrt(q[Order::kW] * q[Order::kW] + q[Order::kX] * q[Order::kX] + - q[Order::kY] * q[Order::kY] + q[Order::kZ] * q[Order::kZ]); + AccurateRNorm(q[Order::kW], q[Order::kX], q[Order::kY], q[Order::kZ]); // Make unit-norm version of q. const std::array<T, 4> unit = MakeQuaternion<Order>(scale * q[Order::kW],
diff --git a/internal/ceres/accurate_norm_test.cc b/internal/ceres/accurate_norm_test.cc index 8345230..102cd65 100644 --- a/internal/ceres/accurate_norm_test.cc +++ b/internal/ceres/accurate_norm_test.cc
@@ -76,6 +76,7 @@ double zero = 0.0; EXPECT_EQ(ceres::AccurateNorm(integer, floating_point, zero), 5.0); + EXPECT_EQ(ceres::AccurateRNorm(integer, floating_point, zero), 0.2); } TEST(AccurateNorm, PromotesIntegralArgumentsOfTheSameType) { @@ -88,6 +89,17 @@ EXPECT_EQ(ceres::AccurateNorm(3, 4, 0), 5.0); } +TEST(AccurateRNorm, PromotesIntegralArgumentsOfTheSameType) { + static_assert(std::is_same_v<decltype(ceres::AccurateRNorm(3, 4)), double>); + static_assert( + std::is_same_v<decltype(ceres::AccurateRNorm(3, 4, 0)), double>); + static_assert( + std::is_same_v<decltype(ceres::AccurateRNorm(3.0F, 4.0F, 0.0F)), float>); + + EXPECT_EQ(ceres::AccurateRNorm(3, 4), 0.2); + EXPECT_EQ(ceres::AccurateRNorm(3, 4, 0), 0.2); +} + #if GTEST_HAS_TYPED_TEST template <typename T> @@ -152,6 +164,175 @@ ceres::AccurateNorm(0, std::numeric_limits<Scalar>::quiet_NaN()))); } +TYPED_TEST(AccurateNormTest, RNorm) { + using Scalar = TypeParam; + + EXPECT_THAT(ceres::AccurateRNorm(this->kTiny, Scalar{0}), + MaxNumUlp(1 / this->kTiny, 0)); + EXPECT_THAT(ceres::AccurateRNorm(this->kTiny, Scalar{0}, Scalar{0}), + MaxNumUlp(1 / this->kTiny, 0)); + EXPECT_THAT(ceres::AccurateRNorm(Scalar{0}, this->kTiny), + MaxNumUlp(1 / this->kTiny, 0)); + EXPECT_THAT(ceres::AccurateRNorm(Scalar{0}, Scalar{0}, this->kTiny), + MaxNumUlp(1 / this->kTiny, 0)); + + EXPECT_THAT(ceres::AccurateRNorm(this->kTiny, this->kTiny), + MaxNumUlp(1 / (std::sqrt(Scalar{2}) * this->kTiny), 1)); + + const auto tiny3 = std::sqrt(this->kTiny) / Scalar{3}; + EXPECT_THAT( + ceres::AccurateRNorm(tiny3, tiny3, tiny3), + MaxNumUlp(ceres::constants::sqrt_3_v<Scalar> / std::sqrt(this->kTiny), + 1)); + + const auto tiny4 = std::sqrt(this->kTiny) / Scalar{4}; + EXPECT_THAT(ceres::AccurateRNorm(tiny4, tiny4, tiny4, tiny4), + MaxNumUlp(Scalar{2} / std::sqrt(this->kTiny), 1)); + + EXPECT_THAT(ceres::AccurateRNorm(this->kHuge, Scalar{0}), + MaxNumUlp(1 / this->kHuge, 0)); + EXPECT_THAT(ceres::AccurateRNorm(this->kHuge, Scalar{0}, Scalar{0}), + MaxNumUlp(1 / this->kHuge, 0)); + EXPECT_THAT(ceres::AccurateRNorm(Scalar{0}, this->kHuge), + MaxNumUlp(1 / this->kHuge, 0)); + EXPECT_THAT(ceres::AccurateRNorm(Scalar{0}, Scalar{0}, this->kHuge), + MaxNumUlp(1 / this->kHuge, 0)); + + EXPECT_TRUE(std::isnan(ceres::AccurateRNorm(0, 0))); + + const auto large = std::sqrt(this->kHuge / 2); + EXPECT_THAT(ceres::AccurateRNorm(large, large), + MaxNumUlp(1 / std::sqrt(this->kHuge), 1)); + + EXPECT_EQ(ceres::AccurateRNorm(+std::numeric_limits<Scalar>::infinity(), 0), + 0); + EXPECT_EQ(ceres::AccurateRNorm(-std::numeric_limits<Scalar>::infinity(), 0), + 0); + + EXPECT_EQ(ceres::AccurateRNorm(0, +std::numeric_limits<Scalar>::infinity()), + 0); + EXPECT_EQ(ceres::AccurateRNorm(0, -std::numeric_limits<Scalar>::infinity()), + 0); + + const Scalar infinity = std::numeric_limits<Scalar>::infinity(); + const Scalar nan = std::numeric_limits<Scalar>::quiet_NaN(); + EXPECT_EQ(ceres::AccurateRNorm(nan, infinity), Scalar{0}); + EXPECT_EQ(ceres::AccurateRNorm(infinity, nan), Scalar{0}); + EXPECT_EQ(ceres::AccurateRNorm(nan, infinity, Scalar{1}), Scalar{0}); + EXPECT_EQ(ceres::AccurateRNorm(infinity, nan, Scalar{1}), Scalar{0}); + + EXPECT_TRUE(std::isnan( + ceres::AccurateRNorm(std::numeric_limits<Scalar>::quiet_NaN(), 0))); + EXPECT_TRUE(std::isnan( + ceres::AccurateRNorm(0, std::numeric_limits<Scalar>::quiet_NaN()))); +} + +TEST(AccurateNorm, ReciprocalNormSmallArguments) { + constexpr double kArgument = 1e-100; + // Correctly rounded 1/(√2·10^-100) and 1/(√3·10^-100) + constexpr double kExpected = 0x1.9dce3e3d94c83p+331; + constexpr double kExpectedThreeArguments = 0x1.51def28355567p+331; + + EXPECT_EQ(ceres::AccurateRNorm(kArgument, kArgument), kExpected); + EXPECT_EQ(ceres::AccurateRNorm(kArgument, kArgument, kArgument), + kExpectedThreeArguments); +} + +TEST(AccurateNorm, ReciprocalNormIsInvariantToSignsAndOrder) { + // The reciprocal hypotenuse is invariant under signs and argument order. + constexpr double kExpected = 0.2; + + EXPECT_EQ(ceres::AccurateRNorm(3.0, 4.0), kExpected); + EXPECT_EQ(ceres::AccurateRNorm(4.0, 3.0), kExpected); + EXPECT_EQ(ceres::AccurateRNorm(-3.0, 4.0), kExpected); + EXPECT_EQ(ceres::AccurateRNorm(3.0, -4.0), kExpected); +} + +TEST(AccurateNorm, ReciprocalNormIsCorrectlyRoundedForLargeArguments) { + // A reciprocal norm below Tiny() has a square root whose rounding error + // underflows. The arguments must therefore be rescaled even though their + // squares do not overflow. + using DoubleTraits = ceres::internal::AccurateNormTraits<double>; + static_assert(DoubleTraits::ReciprocalUnscaledMaximum(2) == 0x1p+484); + static_assert(DoubleTraits::ReciprocalUnscaledMaximum(3) == 0x1p+484); + static_assert(DoubleTraits::ReciprocalUnscaledMaximum(5) == 0x1p+483); + using FloatTraits = ceres::internal::AccurateNormTraits<float>; + static_assert(FloatTraits::ReciprocalUnscaledMaximum(2) == 0x1p+50f); + + EXPECT_EQ( + ceres::AccurateRNorm(0x1.bf280b7c5d483p+500, 0x1.6cd7ee3e743c3p+510), + 0x1.674152d34023dp-511); + EXPECT_EQ(ceres::AccurateRNorm(0x1.58fa125864086p+509, + 0x1.16fda4679f7bfp+480, + 0x1.8758b55834ccep-36), + 0x1.7bf1c082d7297p-510); + EXPECT_EQ(ceres::AccurateRNorm(0x1.c38472p+59f, 0x1.19b414p+63f), + 0x1.cef76ap-64f); + EXPECT_EQ( + ceres::AccurateRNorm(0x1.4debfp+61f, 0x1.425ed8p+46f, 0x1.edf9c6p+0f), + 0x1.8885f2p-62f); +} + +TEST(AccurateNorm, ReciprocalNormAccountsForNegligibleArgument) { + // The smaller argument does not affect the norm of the arguments. The + // reciprocal of the larger argument, however, is rounded and can differ from + // the correctly rounded reciprocal norm. + EXPECT_EQ(ceres::AccurateRNorm(0x1.5e012ffd8fc0cp+0, 0x1.26a1fb826eb7bp-31), + 0x1.767c89013add5p-1); + EXPECT_EQ(ceres::AccurateRNorm(0x1.63341p+0f, 0x1.346e16p-13f), + 0x1.710168p-1f); + EXPECT_EQ(ceres::AccurateRNorm(0x1.e702212683dccp+700, 0x1.afe84e6a6d27p+672), + 0x1.0d23181d548a3p-701); +} + +TEST(AccurateNorm, ReciprocalNormIsCorrectlyRoundedForSmallArguments) { + // The largest magnitude is between Tiny() and UnscaledMinimum(). Rescaling it + // by the same power as subnormal arguments would exceed + // ReciprocalUnscaledMaximum(). + EXPECT_EQ(ceres::AccurateRNorm( + -0x1.6eaa0ep-28f, -0x1.ed4b4ap-42f, 0x1.8fe444p-122f), + 0x1.6578bp+27f); +} + +TEST(AccurateNorm, ReciprocalNormIsCorrectlyRounded) { + // The expected values are the correctly rounded reciprocal norms. The Newton + // update must account for both the rounding error of the reciprocal sum of + // squares and the rounding error of its square root. + EXPECT_EQ(ceres::AccurateRNorm(0x1.5a32756bbbe09p+0, 0x1.7a1e23a3d7545p+0), + 0x1.ff55f7cf589c9p-2); + EXPECT_EQ(ceres::AccurateRNorm(0x1.4999e2p+0f, 0x1.da58ep+0f), + 0x1.c5d62ap-2f); + EXPECT_EQ( + ceres::AccurateRNorm( + 0x1.bb2e63545507dp+0, 0x1.11f1f5cffc5dfp-1, 0x1.47535dd62c22ap+0), + 0x1.cdbe5e6d38adap-2); + EXPECT_EQ(ceres::AccurateRNorm(0x1.ca3d4p+0f, 0x1.07cbd2p+0f, 0x1.1ac788p-2f), + 0x1.eb697ep-2f); +} + +TEST(AccurateNorm, ReciprocalNormRescalesByRadixPowersWithoutRounding) { + // These operands exercise the radix-power rescaling branches. + constexpr int kLargeExponent = 512; + constexpr int kSmallExponent = -514; + + EXPECT_EQ(ceres::AccurateRNorm(std::scalbn(3.0, kLargeExponent), + std::scalbn(4.0, kLargeExponent)), + std::scalbn(0.2, -kLargeExponent)); + EXPECT_EQ(ceres::AccurateRNorm(std::scalbn(3.0, kSmallExponent), + std::scalbn(4.0, kSmallExponent)), + std::scalbn(0.2, -kSmallExponent)); +} + +TEST(AccurateNorm, VariadicReciprocalNormOfMixedSigns) { + constexpr double kFirst = -1.0; + constexpr double kSecond = 0.5; + constexpr double kThird = 0.25; + // Correctly rounded 1/√1.3125 + constexpr double kExpected = 0x1.bee9056fb9c39p-1; + + EXPECT_EQ(ceres::AccurateRNorm(kFirst, kSecond, kThird), kExpected); +} + TEST(AccurateNorm, VariadicNormAccuracy) { EXPECT_THAT(ceres::AccurateNorm(1.0, 1.0, 1.0), MaxNumUlp(ceres::constants::sqrt_3, 0)); @@ -315,6 +496,57 @@ 0x1.6d979f8e15b5ap-432); } +TEST(AccurateNorm, VariadicRNormIsAccurateAroundUnscaledRange) { + // Inside the range that requires no rescaling + EXPECT_THAT(ceres::AccurateRNorm(0x1p+510, 0x1p+510, 0x1p+510), + MaxNumUlp(0x1.279a74590331cp-511, 0)); + EXPECT_THAT(ceres::AccurateRNorm(0x1p-432, 0x1p-432, 0x1p-432), + MaxNumUlp(0x1.279a74590331cp+431, 0)); + EXPECT_THAT(ceres::AccurateRNorm( + 0x1.3a578a98bcf9dp+509, -0x1.4a15545f86ef4p+510, 0x1.5p-200), + MaxNumUlp(0x1.668504f89c76bp-511, 0)); + + // Just outside the range that requires no rescaling + EXPECT_THAT(ceres::AccurateRNorm(0x1.fffffffffffffp+510, + 0x1.fffffffffffffp+510, + 0x1.fffffffffffffp+510), + MaxNumUlp(0x1.279a74590331dp-512, 0)); + EXPECT_THAT(ceres::AccurateRNorm(0x1.fffffffffffffp-433, + 0x1.fffffffffffffp-433, + 0x1.fffffffffffffp-433), + MaxNumUlp(0x1.279a74590331dp+431, 0)); + EXPECT_THAT(ceres::AccurateRNorm( + 0x1.4a15545f86ef4p-432, -0x1.3a578a98bcf9dp-433, 0x1.5p-500), + MaxNumUlp(0x1.668504f89c76bp+431, 0)); +} + +TEST(AccurateNorm, VariadicRNormIsAccurateForExtremeMagnitudes) { + constexpr double kHuge = 0x1.ep+1023; + + // Subnormal arguments + EXPECT_THAT(ceres::AccurateRNorm(0x0.123456789abcdp-1022, + 0x0.fedcba9876543p-1022, + 0x0.0000000000001p-1022), + MaxNumUlp(0x1.007d48013c608p+1022, 0)); + // Small arguments including a negligible one + EXPECT_THAT(ceres::AccurateRNorm( + 0x1.3a578a98bcf9dp-600, -0x1.4a15545f86ef4p-601, 0x1.5p-650), + MaxNumUlp(0x1.712e4c04dcaa4p+599, 0)); + // Arguments close to the largest finite value + EXPECT_THAT(ceres::AccurateRNorm(0x1.ffffffffffffp+1022, + 0x1.4a15545f86ef4p+1021, + -0x1.3a578a98bcf9dp+1000), + MaxNumUlp(0x0.79d3b767a4e2p-1022, 0)); + // Large arguments including a negligible one + EXPECT_THAT(ceres::AccurateRNorm(0x1.8p+1000, 0x1p-1000, 0x1.4p+999), + MaxNumUlp(0x1.3b13b13b13b14p-1001, 0)); + // The norm exceeds the largest finite value but its reciprocal does not + // underflow to zero. + EXPECT_THAT(ceres::AccurateRNorm( + kHuge, kHuge, kHuge, kHuge, kHuge, kHuge, kHuge, kHuge), + MaxNumUlp(0x0.1822cb17ff2ecp-1022, 0)); +} + TEST(AccurateNorm, RescalesByRadixPowersWithoutRounding) { constexpr int kLargeExponent = 512; constexpr int kSmallExponent = -514; @@ -339,11 +571,16 @@ constexpr double kSeventh = -0x1.5483d7f40eca8p-537; constexpr double kEighth = 0x1.4694b77d1bb38p-643; constexpr double kExpectedNorm = 0x1.5fdef349a2773p+922; + constexpr double kExpectedRNorm = 0x1.748010bb48c7cp-923; EXPECT_THAT( ceres::AccurateNorm( kFirst, kSecond, kThird, kFourth, kFifth, kSixth, kSeventh, kEighth), MaxNumUlp(kExpectedNorm, 0)); + EXPECT_THAT( + ceres::AccurateRNorm( + kFirst, kSecond, kThird, kFourth, kFifth, kSixth, kSeventh, kEighth), + MaxNumUlp(kExpectedRNorm, 0)); } TEST(AccurateNorm, VariadicNormRescalesLargeArguments) { @@ -360,6 +597,36 @@ MaxNumUlp(kExpectedNorm, 0)); } +TEST(AccurateNorm, VariadicRNormRescalesLargeArguments) { + constexpr double kFirst = 0x1.b5bae0f220324p-790; + constexpr double kSecond = -0x1.dbe4984f6662p+922; + constexpr double kThird = -0x1.1d795b706aa6ep-13; + constexpr double kFourth = 0x1.d94bdccd99128p+611; + constexpr double kFifth = -0x1.738b5e5278dcap-365; + constexpr double kSixth = 0x1.49b718b5da944p-191; + constexpr double kSeventh = 0x1.6091a40b7733p+922; + constexpr double kExpected = 0x1.ba9cab0b4f49bp-924; + + EXPECT_THAT(ceres::AccurateRNorm( + kFirst, kSecond, kThird, kFourth, kFifth, kSixth, kSeventh), + MaxNumUlp(kExpected, 0)); +} + +TEST(AccurateNorm, VariadicRNormIsCorrectlyRoundedForMixedMagnitudes) { + constexpr double kFirst = 0x1.435e0c899da16p-163; + constexpr double kSecond = 0x1.61091418e8617p+825; + constexpr double kThird = -0x1.a0645fc0b9c0ap+395; + constexpr double kFourth = 0x1.013f4a26833d2p-597; + constexpr double kFifth = 0x1.6ae50a45142b2p+238; + constexpr double kSixth = -0x1.af001dacbda5dp-87; + constexpr double kSeventh = -0x1.1a68ec954a99ap+163; + constexpr double kExpected = 0x1.73457fef4890ap-826; + + EXPECT_THAT(ceres::AccurateRNorm( + kFirst, kSecond, kThird, kFourth, kFifth, kSixth, kSeventh), + MaxNumUlp(kExpected, 0)); +} + TEST(AccurateNorm, VariadicNormReturnsPositiveInfinity) { constexpr double kInfinity = std::numeric_limits<double>::infinity(); @@ -388,6 +655,10 @@ const double norm = ceres::AccurateNorm(nan, 1.0); const double variadic_norm = ceres::AccurateNorm(nan, 1.0, 2.0); EXPECT_EQ(std::memcmp(&variadic_norm, &norm, sizeof(double)), 0); + + const double rnorm = ceres::AccurateRNorm(nan, 1.0); + const double variadic_rnorm = ceres::AccurateRNorm(nan, 1.0, 2.0); + EXPECT_EQ(std::memcmp(&variadic_rnorm, &rnorm, sizeof(double)), 0); } TEST(AccurateNorm, PreservesNaNPayloadOfAnyArgument) {
diff --git a/internal/ceres/jet_test.cc b/internal/ceres/jet_test.cc index 6e1271c..a1d2c42 100644 --- a/internal/ceres/jet_test.cc +++ b/internal/ceres/jet_test.cc
@@ -1,5 +1,5 @@ // Ceres Solver - A fast non-linear least squares minimizer -// Copyright 2024 Google Inc. All rights reserved. +// Copyright 2026 Google Inc. All rights reserved. // http://ceres-solver.org/ // // Redistribution and use in source and binary forms, with or without @@ -715,6 +715,49 @@ static_assert(!IsAccurateNormCallable<void, J, J, Jet<double, 3>>::value, "AccurateNorm must not accept Jets of different sizes"); +TEST(Jet, AccurateRNorm) { + // The arguments form Pythagorean quadruples to obtain exact norms. + const J a = MakeJet(3.0, 1.0, 0.0); + const J b = MakeJet(4.0, 0.0, 1.0); + const J c = MakeJet(12.0, 2.0, -1.0); + const J d = MakeJet(84.0, -1.0, 3.0); + + // d/dx_i 1 / sqrt(x_1^2 + ... + x_n^2) = -x_i / (x_1^2 + ... + x_n^2)^(3/2) + EXPECT_THAT(AccurateRNorm(a, b), + IsAlmostEqualTo(MakeJet(1.0 / 5.0, -3.0 / 125.0, -4.0 / 125.0))); + EXPECT_THAT( + AccurateRNorm(a, b, c), + IsAlmostEqualTo(MakeJet(1.0 / 13.0, -27.0 / 2197.0, 8.0 / 2197.0))); + EXPECT_THAT( + AccurateRNorm(a, b, c, d), + IsAlmostEqualTo(MakeJet(1.0 / 85.0, 57.0 / 614125.0, -244.0 / 614125.0))); + + // Arguments that are zero do not contribute to the derivative. + const J zero = MakeJet(0.0, 2.0, 3.14); + EXPECT_THAT(AccurateRNorm(a, zero), IsAlmostEqualTo(1.0 / a)); + EXPECT_THAT(AccurateRNorm(zero, b, zero), IsAlmostEqualTo(1.0 / b)); + EXPECT_THAT(AccurateRNorm(zero, zero, zero, d), IsAlmostEqualTo(1.0 / d)); + EXPECT_THAT(AccurateRNorm(-a, zero, zero), IsAlmostEqualTo(1.0 / a)); + + EXPECT_THAT(AccurateRNorm(x, y), IsAlmostEqualTo(1.0 / hypot(x, y))); + EXPECT_THAT(AccurateRNorm(x, y, z), IsAlmostEqualTo(1.0 / hypot(x, y, z))); +} + +template <typename AlwaysVoid, typename... Args> +struct IsAccurateRNormCallable : std::false_type {}; + +template <typename... Args> +struct IsAccurateRNormCallable< + std::void_t<decltype(AccurateRNorm(std::declval<Args>()...))>, + Args...> : std::true_type {}; + +static_assert(IsAccurateRNormCallable<void, J, J>::value); +static_assert(IsAccurateRNormCallable<void, J, J, const J&, J&>::value); +static_assert(!IsAccurateRNormCallable<void, J, J, double>::value, + "AccurateRNorm must not accept mixed Jet and scalar arguments"); +static_assert(!IsAccurateRNormCallable<void, J, J, Jet<double, 3>>::value, + "AccurateRNorm must not accept Jets of different sizes"); + #ifdef CERES_HAS_CPP20 TEST(Jet, Lerp) {