[hist] Add a test for finding bins in TAxis.

The tests are based on two issues brought up in the forum, a few cases
designed by hand to break the old algorithm, github issue #14091, the
tests from PR #14105, and a test sampling 100 random points.

Author: hageboeck

hist/hist/test/test_TH1.cxx

@@ -149,3 +149,132 @@ TEST(TH1, Normalize)
    EXPECT_FLOAT_EQ(v2.Integral("width"), 1.);
    EXPECT_FLOAT_EQ(v2.GetMaximum(), 7.9999990);
 }
+
+TEST(TAxis, BinComputation_FPAccuracy)
+{
+   // Example from 1703c54
+   EXPECT_EQ(TAxis(1, -1., 0.).FindBin(-1e-17), 1);
+
+   {
+      // https://root-forum.cern.ch/t/floating-point-rounding-error-when-filling-the-histogram/35189
+      TAxis axis(128, -0.352, 0.352);
+      constexpr double x = -0.0220;
+      EXPECT_EQ(axis.FindBin(x), 61);
+      EXPECT_EQ(axis.FindFixBin(x), 61);
+      EXPECT_LE(axis.GetBinLowEdge(61), x);
+      EXPECT_GT(axis.GetBinUpEdge(61), x);
+   }
+
+   {
+      // https://github.com/root-project/root/issues/14091
+      constexpr int nBins = 30;
+      constexpr double xMin = 3.0, xMax = 6.0;
+      TAxis ax(nBins, xMin, xMax);
+
+      for (Int_t i = 1; i <= ax.GetNbins(); i++) {
+         EXPECT_EQ(i, ax.FindBin(ax.GetBinLowEdge(i)));
+      }
+   }
+
+   {
+      TAxis axis(2000, -1000., 1000.);
+      EXPECT_EQ(axis.FindFixBin(std::nextafter(-1000., -2000.)), 0);
+      EXPECT_EQ(axis.FindFixBin(-1000.), 1);
+      EXPECT_EQ(axis.FindFixBin(std::nextafter(-1000., 0)), 1);
+
+      EXPECT_EQ(axis.FindFixBin(-500.00000000001), 500);
+      EXPECT_EQ(axis.FindFixBin(-500.), 501);
+      EXPECT_EQ(axis.FindFixBin(-499.9), 501);
+
+      EXPECT_EQ(axis.FindFixBin(-1.E-20), 1000);
+      EXPECT_EQ(axis.FindFixBin(std::nextafter(-0., -1)), 1000);
+      EXPECT_EQ(axis.FindFixBin(-0.), 1001);
+      EXPECT_EQ(axis.FindFixBin(0.), 1001);
+
+      EXPECT_EQ(axis.FindFixBin(499.9), 1500);
+      EXPECT_EQ(axis.FindFixBin(500.), 1501);
+
+      EXPECT_EQ(axis.FindFixBin(1000. - 1.E-13), 2000);
+      EXPECT_EQ(axis.FindFixBin(std::nextafter(1000., 0.)), 2000);
+      EXPECT_EQ(axis.FindFixBin(1000.), 2001);
+
+      EXPECT_EQ(axis.FindBin(std::nextafter(-1000., -2000.)), 0);
+      EXPECT_EQ(axis.FindBin(-1000.), 1);
+      EXPECT_EQ(axis.FindBin(std::nextafter(-1000., 0)), 1);
+
+      EXPECT_EQ(axis.FindBin(-500.00000000001), 500);
+      EXPECT_EQ(axis.FindBin(-500.), 501);
+      EXPECT_EQ(axis.FindBin(-499.9), 501);
+
+      EXPECT_EQ(axis.FindBin(-1.E-20), 1000);
+      EXPECT_EQ(axis.FindBin(std::nextafter(-0., -1)), 1000);
+      EXPECT_EQ(axis.FindBin(-0.), 1001);
+      EXPECT_EQ(axis.FindBin(0.), 1001);
+
+      EXPECT_EQ(axis.FindBin(499.9), 1500);
+      EXPECT_EQ(axis.FindBin(500.), 1501);
+
+      EXPECT_EQ(axis.FindBin(1000. - 1.E-13), 2000);
+      EXPECT_EQ(axis.FindBin(std::nextafter(1000., 0.)), 2000);
+      EXPECT_EQ(axis.FindBin(1000.), 2001);
+   }
+
+   for (const auto &[low, high] : std::initializer_list<std::pair<long double, long double>>{{-10654.1l, 32165.l},
+                                                                                             {-1.0656E23l, -20654.l},
+                                                                                             {1.1234E4l, 4.5678E20l},
+                                                                                             {1.E-60l, 1.E-20l},
+                                                                                             {-1.E-20l, -1.E-60l}}) {
+      constexpr int N = 100;
+      const double width = (high - low) / N;
+      std::mt19937 gen;
+      std::uniform_real_distribution dist{low - width, high + width};
+      TAxis axis(N, low, high);
+      for (unsigned int i = 0; i < 100000; ++i) {
+         long double x = dist(gen);
+         if (i == 0) {
+            x = std::nextafter(double(high), double(low));
+         }
+         const auto bin = axis.FindFixBin(x);
+         EXPECT_EQ(bin, axis.FindBin(x));
+
+         if (x < double(low))
+            EXPECT_EQ(bin, 0);
+         else if (x >= double(high))
+            EXPECT_EQ(bin, N + 1);
+
+         EXPECT_LE(axis.GetBinLowEdge(bin), x);
+         EXPECT_LT(x, axis.GetBinUpEdge(bin));
+      }
+   }
+}
+
+// The tests below are taken from https://github.com/root-project/root/pull/14105
+TEST(TAxis, FindBinExact)
+{
+   // Test the case where bin edges are exactly represented as floating points
+   TAxis ax(88, 1010, 1098);
+   for (int i = 1; i <= ax.GetNbins(); i++) {
+      double x = ax.GetBinLowEdge(i);
+      EXPECT_EQ(i, ax.FindBin(x));
+      EXPECT_EQ(i, ax.FindFixBin(x));
+      x = ax.GetBinUpEdge(i);
+      EXPECT_EQ(i + 1, ax.FindBin(x));
+      EXPECT_EQ(i + 1, ax.FindFixBin(x));
+      x -= x * std::numeric_limits<double>::epsilon();
+      EXPECT_EQ(i, ax.FindBin(x));
+   }
+}
+TEST(TAxis, FindBinApprox)
+{
+   TAxis ax(90, 0., 10.);
+   for (int i = 1; i <= ax.GetNbins(); i++) {
+      double x = ax.GetBinLowEdge(i);
+      EXPECT_EQ(i, ax.FindBin(x));
+      EXPECT_EQ(i, ax.FindFixBin(x));
+      x = ax.GetBinUpEdge(i);
+      EXPECT_EQ(i + 1, ax.FindBin(x));
+      EXPECT_EQ(i + 1, ax.FindFixBin(x));
+      x -= x * std::numeric_limits<double>::epsilon();
+      EXPECT_EQ(i, ax.FindBin(x));
+   }
+}