Commits

ssn...@4525493e-7705-40b1-a816-d608a930855b  committed 0424fb1

Add fpcompare.h

  • Participants
  • Parent commits 2357a1d

Comments (0)

Files changed (5)

+2008-09-29  scott snyder  <snyder@bnl.gov>
+
+	* Tagging CxxUtils-00-00-02.
+	* CxxUtils/fpcompare.h: (new)
+	* test/fpcompare_test.cxx: (new)
+	* share/fpcompare_test.ref: (new)
+	* cmt/requirements: Add the test.
+
 2008-09-16  Sebastien Binet  <binet@lblbox>
 
 	* tagging CxxUtils-00-00-01

File CxxUtils/fpcompare.h

+// This file's extension implies that it's C, but it's really -*- C++ -*-.
+// $Id: fpcompare.h,v 1.1 2008-09-29 19:51:59 ssnyder Exp $
+
+/**
+ * @file CxxUtils/fpcompare.h
+ * @author scott snyder
+ * @date Sep 2008
+ * @brief Workaround x86 precision issues for FP inequality comparisons.
+ *
+ * The functions contained here can be used to work around one of the
+ * effects of the brain-damage of the x87 FPU.
+ *
+ * Brief summary: If you're writing a comparison function for sort,
+ * where the comparison depends on computed floating-point values, eg:
+ *
+ *@code
+ *    bool compare (IParticle* a, IParticle* b)
+ *    { return a->pt() > b->bt(); }
+ @endcode
+ *
+ * then you should replace the comparison with a call to one of the
+ * functions in this file:
+ *
+ *@code
+ *    bool compare (IParticle* a, IParticle* b)
+ *    { return CxxUtils::fpcompare::greater (a->pt(), b->bt()); }
+ @endcode
+ *
+ * Longer explanation:
+ *
+ * An expression like this (where pt() returns a double):
+ *
+ *@code
+ *    a->pt() > b->pt()
+ @endcode
+ *
+ * is compiled (on x86) into a sequence like this:
+ *
+ *     call a->pt()
+ *     save result from FPU to a double stack temporary
+ *     call b->pt()
+ *     load the temporary back into the FPU
+ *     do the comparison
+ *
+ * If pt() returns a result with the extra precision bits used
+ * (so that the value changes when rounded to a double), then
+ * it is possible for this comparison to return true for the
+ * case where a==b.  This violates the assumptions that std::sort
+ * makes of the comparison function, and can cause a crash
+ * (possibly even silently wrong results!).
+ *
+ * As a fix, we force both parameters into something that has been declared
+ * @c volatile.  That forces them to be spilled to memory, ensuring
+ * that they are both correctly rounded for the declared data type.
+ * The comparison is then done on these rounded values.
+ *
+ * We condition this on the parameter @c __FLT_EVAL_METHOD__ being 2.
+ * This is defined in the C standard; a value of 2 means that all
+ * FP calculations are done as long double.  For other cases,
+ * we leave out the @c volatile qualifiers; this should result
+ * in the functions being inlined compeltely away.
+ *
+ * In addition to the free functions in the @c CxxUtils::fpcompare
+ * namespace. we define corresponding functionals in the
+ * @c CxxUtils::fpcompare_fn namespace.  Not clear if those will
+ * be useful for anything, though.
+ */
+
+
+#include <cmath>
+#include <functional>
+
+
+#ifndef CXXUTILS_FPCOMPARE_H
+#define CXXUTILS_FPCOMPARE_H
+
+
+// Decide whether we need to use volatile or not.
+#if defined(__FLT_EVAL_METHOD__) && \
+  (__FLT_EVAL_METHOD__ == 2 || __FLT_EVAL_METHOD__ < 0)
+  // __FLT_EVAL_METHOD__ < 0 means unspecified.
+  // Be pessimistic in that case.
+# define CXXUTILS_FPCOMPARE_VOLATILE volatile
+#elif defined(__i386__) && !defined(__SSE2__)
+  // On x86, gcc -msse -mfpmath=sse is observed to _not_ generate
+  // sse fp instructions, but does set __FLT_EVAL_METHOD__ to 0.
+  // -msse2 -mfpmath=sse does seem to work as expected.
+  // Special-case this for now; should follow up with a gcc bug report
+  // if this still happens in current releases.
+# define CXXUTILS_FPCOMPARE_VOLATILE volatile
+#else
+# define CXXUTILS_FPCOMPARE_VOLATILE
+#endif
+
+
+namespace CxxUtils {
+namespace fpcompare {
+
+
+/**
+ * @brief Compare two FP numbers, working around x87 precision issues.
+ * @return @a a > @a b
+ */
+inline
+bool greater (CXXUTILS_FPCOMPARE_VOLATILE double a,
+              CXXUTILS_FPCOMPARE_VOLATILE double b)
+{
+  return a > b;
+}
+
+
+/**
+ * @brief Compare two FP numbers, working around x87 precision issues.
+ * @return @a a < @a b
+ */
+inline
+bool less (CXXUTILS_FPCOMPARE_VOLATILE double a,
+           CXXUTILS_FPCOMPARE_VOLATILE double b)
+{
+  return a < b;
+}
+
+
+/**
+ * @brief Compare two FP numbers, working around x87 precision issues.
+ * @return @a a >= @a b
+ */
+inline
+bool greater_equal (CXXUTILS_FPCOMPARE_VOLATILE double a,
+                    CXXUTILS_FPCOMPARE_VOLATILE double b)
+{
+  return a >= b;
+}
+
+
+/**
+ * @brief Compare two FP numbers, working around x87 precision issues.
+ * @return @a a <= @a b
+ */
+inline
+bool less_equal (CXXUTILS_FPCOMPARE_VOLATILE double a,
+                 CXXUTILS_FPCOMPARE_VOLATILE double b)
+{
+  return a <= b;
+}
+
+
+} // namespace fpcompare
+
+
+namespace fpcompare_fn {
+
+
+/**
+ * @brief Compare two FP numbers, working around x87 precision issues.
+ */
+struct greater
+  : public std::binary_function<double, double, bool>
+{
+  bool
+  operator()(double a, double b) const
+  { return fpcompare::greater (a, b); }
+};
+
+
+/**
+ * @brief Compare two FP numbers, working around x87 precision issues.
+ */
+struct less
+  : public std::binary_function<double, double, bool>
+{
+  bool
+  operator()(double a, double b) const
+  { return fpcompare::less (a, b); }
+};
+
+
+/**
+ * @brief Compare two FP numbers, working around x87 precision issues.
+ */
+struct greater_equal
+  : public std::binary_function<double, double, bool>
+{
+  bool
+  operator()(double a, double b) const
+  { return fpcompare::greater_equal (a, b); }
+};
+
+
+/**
+ * @brief Compare two FP numbers, working around x87 precision issues.
+ */
+struct less_equal
+  : public std::binary_function<double, double, bool>
+{
+  bool
+  operator()(double a, double b) const
+  { return fpcompare::less_equal (a, b); }
+};
+
+
+} // namespace fpcompare_fn
+} // namespace CxxUtils
+
+
+#endif // not CXXUTILS_FPCOMPARE_H

File cmt/requirements

 
 use TestTools	   TestTools-*         AtlasTest
 apply_pattern UnitTest_run unit_test=hashtable
+apply_pattern UnitTest_run unit_test=fpcompare
 
 end_private

File share/fpcompare_test.ref

Empty file added.

File test/fpcompare_test.cxx

+#undef NDEBUG
+/**
+ * @file CxxUtils/test/fpcompare_test.cxx
+ * @author scott snyder <snyder@bnl.gov>
+ * @date Sep, 2008
+ * @brief Regression tests for the fpcompare header.
+ */
+
+
+#include "CxxUtils/fpcompare.h"
+#include <cmath>
+#include <cassert>
+#include <stdio.h>
+
+
+struct tester
+{
+  tester (double x) : m_x (x) {}
+  double pt() const;
+  double m_x;
+};
+
+
+double tester::pt() const
+{
+  return m_x/std::cosh(1.9025873924597605);
+}
+
+
+void test1 (const tester* a, const tester* b, const tester* c)
+{
+  using namespace CxxUtils::fpcompare;
+
+  assert ( ! greater (a->pt(), b->pt()) ); // FAILX
+  assert (   greater (c->pt(), b->pt()) );
+  assert ( ! greater (b->pt(), c->pt()) );
+  
+  assert ( ! less (a->pt(), b->pt()) );
+  assert (   less (a->pt(), c->pt()) );
+  assert ( ! less (c->pt(), a->pt()) );
+
+  assert (   greater_equal (a->pt(), b->pt()) );
+  assert (   greater_equal (c->pt(), b->pt()) );
+  assert ( ! greater_equal (b->pt(), c->pt()) );
+
+  assert (   less_equal (a->pt(), b->pt()) ); // FAILX
+  assert (   less_equal (a->pt(), c->pt()) );
+  assert ( ! less_equal (c->pt(), a->pt()) );
+}
+
+
+void test2 (const tester* a, const tester* b, const tester* c)
+{
+  using namespace CxxUtils::fpcompare_fn;
+
+  greater gt;
+  assert ( ! gt (a->pt(), b->pt()) ); // FAILX
+  assert (   gt (c->pt(), b->pt()) );
+  assert ( ! gt (b->pt(), c->pt()) );
+
+  less lt;
+  assert ( ! lt (a->pt(), b->pt()) );
+  assert (   lt (a->pt(), c->pt()) );
+  assert ( ! lt (c->pt(), a->pt()) );
+
+  greater_equal ge;
+  assert (   ge (a->pt(), b->pt()) );
+  assert (   ge (c->pt(), b->pt()) );
+  assert ( ! ge (b->pt(), c->pt()) );
+
+  less_equal le;
+  assert (   le (a->pt(), b->pt()) ); // FAILX
+  assert (   le (a->pt(), c->pt()) );
+  assert ( ! le (c->pt(), a->pt()) );
+}
+
+
+int main()
+{
+  // chosen so that a->pt() > a->pt is true when compiled
+  // with optimization using x87 instructions.
+  // The assertions marked with FAILX above fail in that case
+  // if the workaround is not operative.
+  tester a (441849.03125);
+  tester c (841849.03125);
+  test1 (&a, &a, &c);
+  test2 (&a, &a, &c);
+  return 0;
+}