// Ceres Solver - A fast non-linear least squares minimizer
// Copyright 2026 Google Inc. All rights reserved.
// http://ceres-solver.org/
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice,
//   this list of conditions and the following disclaimer.
// * Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
// * Neither the name of Google Inc. nor the names of its contributors may be
//   used to endorse or promote products derived from this software without
//   specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
// ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
// LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
// CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
// SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
// INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
// CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
// ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
// POSSIBILITY OF SUCH DAMAGE.
//
// Author: sergiu.deitsch@gmail.com (Sergiu Deitsch)

#include "ceres/accurate_norm.h"

#include <cfloat>
#include <cmath>
#include <cstring>
#include <limits>
#include <type_traits>

#include "gmock/gmock.h"
#include "gtest/gtest.h"

// The constant is computed using mpmath by the following shell command:
//
// python3 - <<EOF
// from mpmath import mp
// mp.dps = 64
// print(mp.sqrt(3))
// EOF
// √3
template <typename T>
inline constexpr T kSqrt3(
    1.732050807568877293527446341505872366942805253810380628055806979L);

// Matches if the argument is at most n floating-point values away from the
// expected value.
MATCHER_P2(MaxNumUlp, expected, n, "") {
  using Scalar = std::decay_t<decltype(arg)>;
  const Scalar target = static_cast<Scalar>(expected);
  Scalar value = arg;

  for (int distance = 0; value != target && distance < n; ++distance) {
    value = std::nextafter(value, target);
  }

  *result_listener << "actual " << arg << " is not within " << n
                   << " ULP of expected " << target;
  return value == target;
}

TEST(AccurateNorm, Promote) {
  static_assert(std::is_same_v<ceres::internal::Promote_t<int>, double>,
                "Promotion of an int must be a double");
  static_assert(std::is_same_v<ceres::internal::Promote_t<int, int>, double>,
                "Promotion of multuple ints must be a double");
  static_assert(std::is_same_v<ceres::internal::Promote_t<unsigned>, double>,
                "Promotion of an unsigned int must be double");
  static_assert(std::is_same_v<ceres::internal::Promote_t<long>, double>,
                "Promotion of a long must be double");
  static_assert(
      std::is_same_v<ceres::internal::Promote_t<int, long, float>, double>,
      "Promotion of arithmetic types must be double");
}

TEST(AccurateNorm, PromotesLvalueArguments) {
  int integer = 3;
  double floating_point = 4.0;
  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) {
  static_assert(std::is_same_v<decltype(ceres::AccurateNorm(3, 4)), double>);
  static_assert(std::is_same_v<decltype(ceres::AccurateNorm(3, 4, 0)), double>);
  static_assert(
      std::is_same_v<decltype(ceres::AccurateNorm(3.0F, 4.0F, 0.0F)), float>);

  EXPECT_EQ(ceres::AccurateNorm(3, 4), 5.0);
  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>
class AccurateNormTest : public testing::Test {
 public:
  static constexpr auto kTiny = std::numeric_limits<T>::min();
  static constexpr auto kHuge = std::numeric_limits<T>::max();
};

using Types = testing::Types<float, double, long double>;

TYPED_TEST_SUITE(AccurateNormTest, Types);

TEST(AccurateNorm, ScaleIsRadixExponent) {
  using Traits = ceres::internal::AccurateNormTraits<double>;

  // ulp(√F_min) = ulp(2^-511) = 2^-563 for double
  static_assert(std::is_same_v<decltype(Traits::ScaleExponent()), int>);
  static_assert(Traits::ScaleExponent() == -563);
}

TYPED_TEST(AccurateNormTest, Norm) {
  using Scalar = TypeParam;

  EXPECT_THAT(ceres::AccurateNorm(this->kTiny, Scalar{0}),
              MaxNumUlp(this->kTiny, 0));
  EXPECT_THAT(ceres::AccurateNorm(this->kTiny, Scalar{0}, Scalar{0}),
              MaxNumUlp(this->kTiny, 0));
  EXPECT_THAT(ceres::AccurateNorm(Scalar{0}, this->kTiny),
              MaxNumUlp(this->kTiny, 0));
  EXPECT_THAT(ceres::AccurateNorm(Scalar{0}, Scalar{0}, this->kTiny),
              MaxNumUlp(this->kTiny, 0));

  EXPECT_THAT(ceres::AccurateNorm(this->kHuge, Scalar{0}),
              MaxNumUlp(this->kHuge, 0));
  EXPECT_THAT(ceres::AccurateNorm(this->kHuge, Scalar{0}, Scalar{0}),
              MaxNumUlp(this->kHuge, 0));
  EXPECT_THAT(ceres::AccurateNorm(Scalar{0}, this->kHuge),
              MaxNumUlp(this->kHuge, 0));
  EXPECT_THAT(ceres::AccurateNorm(Scalar{0}, Scalar{0}, this->kHuge),
              MaxNumUlp(this->kHuge, 0));

  EXPECT_THAT(ceres::AccurateNorm(this->kTiny, this->kTiny),
              MaxNumUlp(this->kTiny * std::sqrt(Scalar{2}), 1));

  EXPECT_THAT(ceres::AccurateNorm(Scalar{0}, Scalar{0}),
              MaxNumUlp(Scalar{0}, 0));

  EXPECT_TRUE(std::isinf(
      ceres::AccurateNorm(+std::numeric_limits<Scalar>::infinity(), 0)));
  EXPECT_TRUE(std::isinf(
      ceres::AccurateNorm(-std::numeric_limits<Scalar>::infinity(), 0)));

  EXPECT_TRUE(std::isinf(
      ceres::AccurateNorm(0, +std::numeric_limits<Scalar>::infinity())));
  EXPECT_TRUE(std::isinf(
      ceres::AccurateNorm(0, -std::numeric_limits<Scalar>::infinity())));

  EXPECT_TRUE(std::isnan(
      ceres::AccurateNorm(std::numeric_limits<Scalar>::quiet_NaN(), 0)));
  EXPECT_TRUE(std::isnan(
      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(kSqrt3<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(kSqrt3<double>, 0));

  // Combination exposing a difference of at least two ULPs in inaccurate
  // implementations found by random search.
  constexpr double kFirst = 0.0;
  constexpr double kSecond = -0x1.c4a46e8d5e9f3p-940;
  constexpr double kThird = 0x1.2870a0a1f3fa3p-943;
  constexpr double kExpected = 0x1.c628110110bf1p-940;
  EXPECT_THAT(ceres::AccurateNorm(kFirst, kSecond, kThird),
              MaxNumUlp(kExpected, 0));
}

TEST(AccurateNorm, HandlesTableMakerDilemma) {
  // Values from Borges, Algorithm 1014, Section 6.
  constexpr double kFirst = 0x1.a308e1455f447p+0;
  constexpr double kSecond = 0x1.9d931a83ef879p+0;
  constexpr double kExpected = 0x1.2660d009d54f9p+1;

  EXPECT_THAT(ceres::AccurateNorm(kFirst, kSecond), MaxNumUlp(kExpected, 1));
}

TEST(AccurateNorm, ReturnsLargerArgumentWhenSmallerIsNegligible) {
  constexpr double kLarger = 1.0;
  // √(ε/2) is the largest ratio of the arguments for which the norm equals the
  // larger argument.
  constexpr double kSmaller = 0x1.6a09e667f3bcdp-27;

  EXPECT_EQ(ceres::AccurateNorm(kLarger, kSmaller), kLarger);
}

TYPED_TEST(AccurateNormTest, AccountsForSmallerArgumentAboveCutoff) {
  using Scalar = TypeParam;
  // The squared ratio of the arguments equals the machine epsilon which raises
  // the norm by three quarters of an ULP above the larger argument.
  constexpr Scalar kLarger{1.5};
  const Scalar kSmaller =
      kLarger * std::sqrt(std::numeric_limits<Scalar>::epsilon());
  const Scalar kExpected =
      std::nextafter(kLarger, std::numeric_limits<Scalar>::infinity());

  EXPECT_EQ(ceres::AccurateNorm(kLarger, kSmaller), kExpected);
}

TEST(AccurateNorm, RecoversSquaringErrorsOfSmallArguments) {
  // The arguments are not rescaled but their squares are small enough for the
  // rounding errors of the squares to underflow unless the scaling threshold
  // accounts for them.
  EXPECT_EQ(ceres::AccurateNorm(0x1.23342ep-63f, 0x1.271008p-63f),
            0x1.9e8fdcp-63f);
  EXPECT_EQ(ceres::AccurateNorm(0x1.4a15545f86ef4p-510, 0x1.3a578a98bcf9dp-510),
            0x1.c7d03b992e6a8p-510);
#if LDBL_MANT_DIG == 64 && LDBL_MIN_EXP == -16381
  EXPECT_EQ(ceres::AccurateNorm(0x8.b912469c50fab7ap-8194L,
                                0x8.c54f2ded5ba8268p-8194L),
            0xc.5eb4864ce876f86p-8194L);
#endif
}

TEST(AccurateNorm, IsAtLeastAsAccurateAsHypot) {
  // The expected values are correctly rounded. Implementations of std::hypot,
  // e.g., the one shipped with libstdc++ 16, are off by one ULP for these
  // arguments.
  constexpr double kFirst = 0x1.424bf2ed916bfp+0;
  constexpr double kSecond = 0x1.1435107c5d458p+0;
  constexpr double kExpected = 0x1.a8758df39043fp+0;

  EXPECT_EQ(ceres::AccurateNorm(kFirst, kSecond), kExpected);
  EXPECT_LE(std::fabs(ceres::AccurateNorm(kFirst, kSecond) - kExpected),
            std::fabs(std::hypot(kFirst, kSecond) - kExpected));

  constexpr double kThreeFirst = 0x1.c11f6531eb66ep+0;
  constexpr double kThreeSecond = 0x1.f30567547a34cp+0;
  constexpr double kThreeThird = 0x1.1e0edcc120696p+0;
  constexpr double kThreeExpected = 0x1.6ce2663da03a7p+1;

  EXPECT_EQ(ceres::AccurateNorm(kThreeFirst, kThreeSecond, kThreeThird),
            kThreeExpected);
  EXPECT_LE(
      std::fabs(ceres::AccurateNorm(kThreeFirst, kThreeSecond, kThreeThird) -
                kThreeExpected),
      std::fabs(std::hypot(kThreeFirst, kThreeSecond, kThreeThird) -
                kThreeExpected));
}

TYPED_TEST(AccurateNormTest, PowerOfTwoIsExact) {
  using Scalar = TypeParam;
  using std::scalbn;
  constexpr int kMinExponent = std::numeric_limits<Scalar>::min_exponent -
                               std::numeric_limits<Scalar>::digits;
  constexpr int kMaxExponent = std::numeric_limits<Scalar>::max_exponent - 1;

  for (int exponent = kMinExponent; exponent <= kMaxExponent; ++exponent) {
    EXPECT_EQ(ceres::internal::PowerOfTwo<Scalar>(exponent),
              scalbn(Scalar{1}, exponent))
        << "exponent " << exponent;
  }
}

TEST(AccurateNorm, TraitsAreConstantExpressions) {
  // Compilers that do not fold std::scalbn would otherwise compute the
  // thresholds at runtime on every invocation.
  using DoubleTraits = ceres::internal::AccurateNormTraits<double>;
  static_assert(DoubleTraits::Tiny() == 0x1p-485);
  static_assert(DoubleTraits::UnscaledMinimum() == 0x1p-432);
  static_assert(DoubleTraits::UnscaledMaximum(2) == 0x1p+511);
  static_assert(DoubleTraits::UnscaledMaximum(3) == 0x1p+510);

  using FloatTraits = ceres::internal::AccurateNormTraits<float>;
  static_assert(FloatTraits::Tiny() == 0x1p-51f);
  static_assert(FloatTraits::UnscaledMinimum() == 0x1p-27f);
  static_assert(FloatTraits::UnscaledMaximum(3) == 0x1p+62f);
}

TEST(AccurateNorm, VariadicNormIsAccurateForExtremeMagnitudes) {
  // Subnormal arguments
  EXPECT_EQ(ceres::AccurateNorm(0x0.123456789abcdp-1022,
                                0x0.fedcba9876543p-1022,
                                0x0.0000000000001p-1022),
            0x0.ff82f53036b9cp-1022);
  // Small arguments including a negligible one
  EXPECT_EQ(ceres::AccurateNorm(
                0x1.3a578a98bcf9dp-600, -0x1.4a15545f86ef4p-601, 0x1.5p-650),
            0x1.6308e012498e5p-600);
  // Arguments close to the largest finite value
  EXPECT_EQ(ceres::AccurateNorm(0x1.ffffffffffffp+1022,
                                0x1.4a15545f86ef4p+1021,
                                -0x1.3a578a98bcf9dp+1000),
            0x1.0cf8b69a0aff1p+1023);
  // Large arguments including a negligible one
  EXPECT_EQ(ceres::AccurateNorm(0x1.8p+1000, 0x1p-1000, 0x1.4p+999),
            0x1.ap+1000);
  // The norm exceeds the largest finite value
  EXPECT_EQ(ceres::AccurateNorm(0x1.ep+1023, 0x1.ep+1023, 0x1.ep+1023),
            std::numeric_limits<double>::infinity());
}

TEST(AccurateNorm, VariadicNormIsAccurateAroundUnscaledRange) {
  constexpr double kSqrt3Largest = 0x1.bb67ae8584caap+510;
  constexpr double kSqrt3Smallest = 0x1.bb67ae8584caap-432;

  // Inside the range that requires no rescaling
  EXPECT_EQ(ceres::AccurateNorm(0x1p+510, 0x1p+510, 0x1p+510), kSqrt3Largest);
  EXPECT_EQ(ceres::AccurateNorm(0x1p-432, 0x1p-432, 0x1p-432), kSqrt3Smallest);
  EXPECT_EQ(ceres::AccurateNorm(
                0x1.3a578a98bcf9dp+509, -0x1.4a15545f86ef4p+510, 0x1.5p-200),
            0x1.6d979f8e15b5ap+510);

  // Just outside the range that requires no rescaling
  EXPECT_EQ(ceres::AccurateNorm(0x1.fffffffffffffp+510,
                                0x1.fffffffffffffp+510,
                                0x1.fffffffffffffp+510),
            0x1.bb67ae8584caap+511);
  EXPECT_EQ(ceres::AccurateNorm(0x1.fffffffffffffp-433,
                                0x1.fffffffffffffp-433,
                                0x1.fffffffffffffp-433),
            kSqrt3Smallest);
  EXPECT_EQ(ceres::AccurateNorm(
                0x1.4a15545f86ef4p-432, -0x1.3a578a98bcf9dp-433, 0x1.5p-500),
            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;

  EXPECT_EQ(ceres::AccurateNorm(std::scalbn(3.0, kLargeExponent),
                                std::scalbn(4.0, kLargeExponent)),
            std::scalbn(5.0, kLargeExponent));
  EXPECT_EQ(ceres::AccurateNorm(std::scalbn(3.0, kSmallExponent),
                                std::scalbn(4.0, kSmallExponent)),
            std::scalbn(5.0, kSmallExponent));
}

TEST(AccurateNorm, VariadicNormHandlesWideDynamicRange) {
  // Combination found by random search that exposes a difference of more than
  // two ULPs when normalized values are not stored.
  constexpr double kFirst = -0x1.5fdef349a2773p+922;
  constexpr double kSecond = -0x1.72ec46b66e1d9p-114;
  constexpr double kThird = -0x1.eeb9ef28337eep-462;
  constexpr double kFourth = -0x1.137462de2cf44p+205;
  constexpr double kFifth = -0x1.3a6cb8edf6622p+264;
  constexpr double kSixth = 0x1.0f5b36fe2970ap-326;
  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) {
  constexpr double kFirst = -0x1.5cf602c1b383ep-360;
  constexpr double kSecond = 0x1.e7855aa96a0c7p+850;
  constexpr double kThird = -0x1.e3bf678284e4p+876;
  constexpr double kFourth = 0x1.7d18d55d5723fp+173;
  constexpr double kFifth = -0x1.1a6c3b4a7f1f8p+437;
  constexpr double kSixth = -0x1.14b8d80b68457p-561;
  constexpr double kExpectedNorm = 0x1.e3bf678284e41p+876;

  EXPECT_THAT(
      ceres::AccurateNorm(kFirst, kSecond, kThird, kFourth, kFifth, kSixth),
      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();

  EXPECT_EQ(ceres::AccurateNorm(-kInfinity, 1.0, 2.0), kInfinity);
}

TEST(AccurateNorm, NonfiniteArgumentHandling) {
  constexpr double kInfinity = std::numeric_limits<double>::infinity();
  constexpr double kNaN = std::numeric_limits<double>::quiet_NaN();

  EXPECT_EQ(ceres::AccurateNorm(kNaN, kInfinity), kInfinity);
  EXPECT_EQ(ceres::AccurateNorm(kInfinity, kNaN), kInfinity);
  EXPECT_EQ(ceres::AccurateNorm(1.0, kInfinity, kNaN), kInfinity);
  EXPECT_EQ(ceres::AccurateNorm(kNaN, 1.0, kInfinity), kInfinity);
  EXPECT_EQ(ceres::AccurateNorm(kNaN, kInfinity, 1.0), kInfinity);
  EXPECT_EQ(ceres::AccurateNorm(kInfinity, kNaN, 1.0), kInfinity);

  EXPECT_TRUE(std::isnan(ceres::AccurateNorm(kNaN, kNaN, 1.0)));
  EXPECT_TRUE(std::isnan(ceres::AccurateNorm(1.0, kNaN, kNaN)));
  EXPECT_TRUE(std::isnan(ceres::AccurateNorm(kNaN, 1.0, kNaN)));
}

TEST(AccurateNorm, PreservesNaNPayloadAcrossArity) {
  const double nan = std::copysign(std::nan("12345"), -1.0);

  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) {
  const double nan = std::copysign(std::nan("12345"), -1.0);
  const double expected = std::fabs(nan);

  const double first = ceres::AccurateNorm(nan, 1.0, 2.0);
  const double last = ceres::AccurateNorm(1.0, 2.0, nan);
  EXPECT_EQ(std::memcmp(&first, &expected, sizeof(double)), 0);
  EXPECT_EQ(std::memcmp(&last, &expected, sizeof(double)), 0);
}

#endif