Initial commit of compareFloatingPoint utils. (#1443)

* Initial commit of compareFloatingPoint utils.

* Update .gitignore to properly ignore unittest executables.

Makefile

@@ -90,6 +90,7 @@ endif
 ER7_UTILS_OBJS = $(addsuffix /object_$(TRICK_HOST_CPU)/*.o ,$(ER7_UTILS_DIRS))
 
 UTILS_DIRS := \
+	${TRICK_HOME}/trick_source/trick_utils/compareFloatingPoint \
 	${TRICK_HOME}/trick_source/trick_utils/interpolator \
 	${TRICK_HOME}/trick_source/trick_utils/trick_adt \
 	${TRICK_HOME}/trick_source/trick_utils/comm \

include/trick/compareFloatingPoint.hh

@@ -0,0 +1,13 @@
+#ifndef COMPARE_FLOATING_POINT_HH
+#define COMPARE_FLOATING_POINT_HH
+
+/* author: John M. Penn */
+
+namespace Trick {
+
+bool dbl_is_near( double A, double B, double tolerance);
+bool flt_is_near( float A, float B, float tolerance);
+
+}
+
+#endif

trick_source/trick_utils/compareFloatingPoint/Makefile

@@ -0,0 +1,5 @@
+
+include ${TRICK_HOME}/share/trick/makefiles/Makefile.common
+include ${TRICK_HOME}/share/trick/makefiles/Makefile.tricklib
+-include Makefile_deps
+

trick_source/trick_utils/compareFloatingPoint/README.md

@@ -0,0 +1,31 @@
+
+# Compare Floating Point Numbers
+
+The following functions compare floating-point numbers to determine whether they are within a specified tolerance of each other.
+
+These functions are designed to never generate ```FP_SUBNORMAL``` numbers, that could result in a floating point underflow exception, even if the OS doesn't handle floating point underflows by setting their values to zero.
+
+## Header
+
+```#include "trick/compareFloatingPoint.hh" ```
+
+## ```Trick::dbl_is_near```
+```c
+bool Trick::dbl_is_near( double A, double B, double tolerance);
+```
+
+This function compares the values of ```double A``` and ```double B``` to determine whether they are within tolerance of each other. If they are, then the function returns ```true```, otherwise it returns ```false```.
+
+The design of ```Trick::dbl_is_near``` requires that the minimum tolerance be ```DBL_MIN/DBL_EPSILON,``` which is approximately ```1.00208e-292```. That is, any two arguments whose difference is less than or equal to ```1.00208e-292``` are considered to be within tolerance, regardless of the specified tolerance.
+
+Before thinking that doubles should be compared to a tolerance smaller than ```1.00208e-292```, please consider that the ratio of the Planck length to the size of the observable universe is approximately ```1.8e-62```. Also consider that our minimum tolerance is ```5.4e-231``` times smaller than that. So, we think that'll probably be good enough in most cases.
+
+## ```Trick::dbl_is_near```
+
+```c
+bool Trick::flt_is_near( float A, float B, float tolerance);
+```
+
+This function compares the values of ```float A``` and ```float B``` to determine whether they are within tolerance of each other. If they are, then the function returns ```true```, otherwise it returns ```false```.
+
+The minimum tolerance for ```Trick::flt_is_near``` is ```FLT_MIN/FLT_EPSILON```, which is approximately ```9.86076e-32```.

trick_source/trick_utils/compareFloatingPoint/src/compareFloatingPoint.cpp

@@ -0,0 +1,70 @@
+#include <float.h>
+#include <math.h>
+#include "trick/compareFloatingPoint.hh"
+
+bool Trick::dbl_is_near( double A, double B, double tolerance) {
+
+    if (isnan(A) || isnan(B) || isinf(A) || isinf(B)) { return false; }
+
+    // If A or B is a FP_SUBNORMAL that is: less than
+    // DBL_MIN (2.22507e-308) then set it to zero.
+    if ( fpclassify(A) == FP_SUBNORMAL ) { A = 0.0; }
+    if ( fpclassify(B) == FP_SUBNORMAL ) { B = 0.0; }
+
+    // If A and B are identical, then they're close_enough.
+    if (A==B) { return true; }
+
+    // The tolerance must be an FP_NORMAL. Neither of FP_INFINITE and
+    // FP_NAN makes sense. Nor do FP_SUBNORMAL and FP_ZERO, given
+    // the above tests.
+
+    // In order than we not generate an FP_underflow, we must set the minimum
+    // allowable tolerance such that fmin(A,B)+tolerance (or fmax(A,B)-tolerance)
+    // cannot be FP_SUBNORMAL. This is only possible if (tolerance >= DBL_MIN/DBL_EPSILON),
+    // where the gap_size around tolerance >= DBL_MIN, the smallest FP_NORMAL number.
+    // So, if A and B are within 1.00208e-292 of each other, they will always be
+    // considered close_enough.
+    if (( fpclassify(tolerance) != FP_NORMAL ) || (tolerance < (DBL_MIN/DBL_EPSILON) )) {
+       tolerance = (DBL_MIN/DBL_EPSILON);
+    }
+
+    // For A and B to be close enough, the tolerance must be greater than or
+    // equal to the larger of the gaps around A and B.
+    if ( tolerance >= DBL_EPSILON * fmax( fabs(A), fabs(B))) {
+
+        // We want to avoid directly computing the difference between A and B.
+        // We might otherwise generate an FP_SUBNORMAL. For example:
+        // If A = (1.5*DBL_MIN), and B = DBL_MIN, the difference is an
+        // FP_SUBNORMAL value.
+        // When FP_SUBNORMAL values are generated, so are FP_underflows.
+        // By insisting that A and B are FP_NORMAL and that the
+        // gap_size around our tolerance is at least DBL_MIN, then we
+        // can avoid generating FP_SUBNORMALs.
+        if ((fmin(A,B) + tolerance) >= fmax(A,B)) {
+            return true ; // A and B are close enough.
+        } else {
+            return false ;
+        }
+    } else {
+        return false;
+    }
+}
+
+bool Trick::flt_is_near( float A, float B, float tolerance) {
+    if (isnan(A) || isnan(B) || isinf(A) || isinf(B)) { return false; }
+    if ( fpclassify(A) == FP_SUBNORMAL ) { A = 0.0; }
+    if ( fpclassify(B) == FP_SUBNORMAL ) { B = 0.0; }
+    if (A==B) { return true; }
+    if (( fpclassify(tolerance) != FP_NORMAL ) || (tolerance < (FLT_MIN/FLT_EPSILON) )) {
+       tolerance = (FLT_MIN/FLT_EPSILON);
+    }
+    if ( tolerance >= FLT_EPSILON * fmax( fabs(A), fabs(B))) {
+        if ((fmin(A,B) + tolerance) >= fmax(A,B)) {
+            return true ;
+        } else {
+            return false ;
+        }
+    } else {
+        return false;
+    }
+}

trick_source/trick_utils/compareFloatingPoint/test/.gitignore

@@ -0,0 +1,2 @@
+*.o
+*_unittest

trick_source/trick_utils/compareFloatingPoint/test/Makefile

@@ -0,0 +1,47 @@
+
+#SYNOPSIS:
+#
+#   make [all]  - makes everything.
+#   make TARGET - makes the given target.
+#   make clean  - removes all files generated by make.
+
+include ${TRICK_HOME}/share/trick/makefiles/Makefile.common
+
+# Flags passed to the preprocessor.
+TRICK_CXXFLAGS += -I$(GTEST_HOME)/include -I$(TRICK_HOME)/include -g -Wall -Wextra -std=c++11 ${TRICK_SYSTEM_CXXFLAGS}
+
+TRICK_LIBS = ${TRICK_LIB_DIR}/libtrick.a
+TRICK_EXEC_LINK_LIBS += -L${GTEST_HOME}/lib64 -L${GTEST_HOME}/lib -lgtest -lgtest_main -lpthread
+
+# Added for Ubuntu... not required for other systems.
+TRICK_EXEC_LINK_LIBS += -lpthread
+
+# All tests produced by this Makefile.  Remember to add new tests you
+# created to the list.
+TESTS = dbl_is_near_unittest flt_is_near_unittest
+
+OTHER_OBJECTS =
+
+# House-keeping build targets.
+
+all : $(TESTS)
+
+test: $(TESTS)
+	./dbl_is_near_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/dbl_is_near.xml
+	./flt_is_near_unittest --gtest_output=xml:${TRICK_HOME}/trick_test/flt_is_near.xml
+
+clean :
+	rm -f $(TESTS) *.o
+	rm -rf io_src xml
+
+dbl_is_near_unittest.o : dbl_is_near_unittest.cc
+	$(TRICK_CXX) $(TRICK_CXXFLAGS) -c $<
+
+dbl_is_near_unittest : dbl_is_near_unittest.o
+	$(TRICK_CXX) $(TRICK_CXXFLAGS) -o $@ $^ $(OTHER_OBJECTS) -L${TRICK_HOME}/lib_${TRICK_HOST_CPU} $(TRICK_LIBS) $(TRICK_EXEC_LINK_LIBS)
+
+flt_is_near_unittest.o : flt_is_near_unittest.cc
+	$(TRICK_CXX) $(TRICK_CXXFLAGS) -c $<
+
+flt_is_near_unittest : flt_is_near_unittest.o
+	$(TRICK_CXX) $(TRICK_CXXFLAGS) -o $@ $^ $(OTHER_OBJECTS) -L${TRICK_HOME}/lib_${TRICK_HOST_CPU} $(TRICK_LIBS) $(TRICK_EXEC_LINK_LIBS)

trick_source/trick_utils/compareFloatingPoint/test/dbl_is_near_unittest.cc

@@ -0,0 +1,160 @@
+
+#include <gtest/gtest.h>
+#include <iostream>
+#include "trick/compareFloatingPoint.hh"
+#include <math.h>
+#include <float.h>
+
+TEST(dbl_is_near_unittest, Simple_1) {
+   bool result;
+   double A = 1.0;
+   double B = 1.1;
+   double tolerance = 0.2;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(dbl_is_near_unittest, Simple_2) {
+   bool result;
+   double A = 1234.567891;
+   double B = 1234.567882;
+   double tolerance = 0.00001;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(dbl_is_near_unittest, Simple_3) {
+   bool result;
+   double A = -1.562154;
+   double B =  0.435837;
+   double tolerance = 2.0;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(dbl_is_near_unittest, Simple_4) {
+   bool result;
+   double A = -1.562154;
+   double B =  0.435837;
+   double tolerance = 1.8;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(dbl_is_near_unittest, Simple_5) {
+   bool result;
+   double A = -1.562154;
+   double B = -0.435837;
+   double tolerance = 1.2;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(dbl_is_near_unittest, A_is_FP_NAN) {
+   bool result;
+   double A = NAN;
+   double B = 0.0;
+   double tolerance = DBL_MAX;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+TEST(dbl_is_near_unittest, B_is_FP_NAN) {
+   bool result;
+   double A = 0.0;
+   double B = NAN;
+   double tolerance = DBL_MAX;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(dbl_is_near_unittest, A_is_FP_INFINITE) {
+   bool result;
+   double A = HUGE_VAL;
+   double B = DBL_MAX;
+   double tolerance = 2 * DBL_EPSILON * DBL_MAX;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(dbl_is_near_unittest, B_is_FP_INFINITE) {
+   bool result;
+   double A = DBL_MAX;
+   double B = HUGE_VAL;
+   double tolerance = 2 * DBL_EPSILON * DBL_MAX;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(dbl_is_near_unittest, A_and_B_are_identical) {
+   // Tolerance is irrelavant because A and B are identical.
+   bool result;
+   double A = DBL_MIN;
+   double B = DBL_MIN;
+   double tolerance = 0.0;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(dbl_is_near_unittest, Within_MinimumTolerance) {
+  // The specified tolerance is < DBL_MIN/DBL_EPSILON, and so defaults to the
+  // minimum. Since the difference between A and B is within the tolerance,
+  // they are "near".
+   bool result;
+   double A = DBL_MIN;
+   double B = 1.5 * DBL_MIN;
+   double tolerance = 0.0;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(dbl_is_near_unittest, Exactly_Minimum_Tolerance) {
+   // The specified tolerance is < DBL_MIN/DBL_EPSILON, and so defaults to the
+   // minimum (DBL_MIN/DBL_EPSILON). Since the difference between A and B is
+   // exactly equal to the tolerance, they are "near".
+   bool result;
+   double A = 0.0;
+   double B = DBL_MIN/DBL_EPSILON;
+   double tolerance = 0.0;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(dbl_is_near_unittest, Greater_than_Minimum_Tolerance) {
+  // The specified tolerance is < DBL_MIN/DBL_EPSILON, and so defaults to the
+  // minimum (DBL_MIN/DBL_EPSILON). Since the difference between A and B is
+  // slightly greater than the tolerance, they are not "near".
+   bool result;
+   double A = 0.0;
+   double B = DBL_MIN/DBL_EPSILON + DBL_MIN;
+   double tolerance = 0.0;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(dbl_is_near_unittest, Tolerance_greater_than_minimum) {
+   // This test is like the previous, but specifies the tolerance to be slightly
+   // larger than the minimum, so that A and B are near.
+   bool result;
+   double A = 0.0;
+   double B = DBL_MIN/DBL_EPSILON + DBL_MIN;
+   double tolerance = DBL_MIN/DBL_EPSILON + DBL_MIN;
+   result = Trick::dbl_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(dbl_is_near_unittest, Tolerance_is_small_enough) {
+   // Is the tolerance small enough?
+   double min_tolerance = DBL_MIN/DBL_EPSILON;
+   const double planck_length = 1.6e-35; // meters
+   const double speed_of_light = 2.99e8; // meters/second
+   const double seconds_per_year = 3.1e7; // seconds/year
+   const double light_years_per_known_universe = 9.3e10; // lightyears/known_universe
+   double size_of_known_universe = speed_of_light * seconds_per_year * light_years_per_known_universe;
+   double universe_min_to_max_ratio = planck_length / size_of_known_universe;
+   std::cout << "=========================================" << std::endl;
+   std::cout << "minimum tolerance = " << min_tolerance << std::endl;
+   std::cout << "planck_length / size_of_known_universe = " << universe_min_to_max_ratio << std::endl;
+   std::cout << "=========================================" << std::endl;
+   bool result = (min_tolerance < universe_min_to_max_ratio);
+   EXPECT_TRUE(result);
+}

trick_source/trick_utils/compareFloatingPoint/test/flt_is_near_unittest.cc

@@ -0,0 +1,151 @@
+
+#include <gtest/gtest.h>
+#include <iostream>
+#include "trick/compareFloatingPoint.hh"
+#include <math.h>
+#include <float.h>
+
+TEST(flt_is_near_unittest, Simple_1) {
+   bool result;
+   float A = 1.0;
+   float B = 1.1;
+   float tolerance = 0.2;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(flt_is_near_unittest, Simple_2) {
+   bool result;
+   float A = 1234.567891;
+   float B = 1234.567882;
+   float tolerance = 0.00001;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(flt_is_near_unittest, Simple_3) {
+   bool result;
+   float A = -1.562154;
+   float B =  0.435837;
+   float tolerance = 2.0;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(flt_is_near_unittest, Simple_4) {
+   bool result;
+   float A = -1.562154;
+   float B =  0.435837;
+   float tolerance = 1.8;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(flt_is_near_unittest, Simple_5) {
+   bool result;
+   float A = -1.562154;
+   float B = -0.435837;
+   float tolerance = 1.2;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(flt_is_near_unittest, A_is_FP_NAN) {
+   bool result;
+   float A = NAN;
+   float B = 0.0;
+   float tolerance = FLT_MAX;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+TEST(flt_is_near_unittest, B_is_FP_NAN) {
+   bool result;
+   float A = 0.0;
+   float B = NAN;
+   float tolerance = FLT_MAX;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(flt_is_near_unittest, A_is_FP_INFINITE) {
+   bool result;
+   float A = HUGE_VAL;
+   float B = FLT_MAX;
+   float tolerance = 2 * FLT_EPSILON * FLT_MAX;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(flt_is_near_unittest, B_is_FP_INFINITE) {
+   bool result;
+   float A = FLT_MAX;
+   float B = HUGE_VAL;
+   float tolerance = 2 * FLT_EPSILON * FLT_MAX;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(flt_is_near_unittest, A_and_B_are_identical) {
+   // Tolerance is irrelavant because A and B are identical.
+   bool result;
+   float A = FLT_MIN;
+   float B = FLT_MIN;
+   float tolerance = 0.0;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(flt_is_near_unittest, Within_MinimumTolerance) {
+  // The specified tolerance is < FLT_MIN/FLT_EPSILON, and so defaults to the
+  // minimum. Since the difference between A and B is within the tolerance,
+  // they are "near".
+   bool result;
+   float A = FLT_MIN;
+   float B = 1.5 * FLT_MIN;
+   float tolerance = 0.0;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(flt_is_near_unittest, Exactly_Minimum_Tolerance) {
+   // The specified tolerance is < FLT_MIN/FLT_EPSILON, and so defaults to the
+   // minimum (FLT_MIN/FLT_EPSILON). Since the difference between A and B is
+   // exactly equal to the tolerance, they are "near".
+   bool result;
+   float A = 0.0;
+   float B = FLT_MIN/FLT_EPSILON;
+   float tolerance = 0.0;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(flt_is_near_unittest, Greater_than_Minimum_Tolerance) {
+  // The specified tolerance is < FLT_MIN/FLT_EPSILON, and so defaults to the
+  // minimum (FLT_MIN/FLT_EPSILON). Since the difference between A and B is
+  // slightly greater than the tolerance, they are not "near".
+   bool result;
+   float A = 0.0;
+   float B = FLT_MIN/FLT_EPSILON + FLT_MIN;
+   float tolerance = 0.0;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_FALSE(result);
+}
+
+TEST(flt_is_near_unittest, Tolerance_greater_than_minimum) {
+   // This test is like the previous, but specifies the tolerance to be slightly
+   // larger than the minimum, so that A and B are near.
+   bool result;
+   float A = 0.0;
+   float B = FLT_MIN/FLT_EPSILON + FLT_MIN;
+   float tolerance = FLT_MIN/FLT_EPSILON + FLT_MIN;
+   result = Trick::flt_is_near(A, B, tolerance);
+   EXPECT_TRUE(result);
+}
+
+TEST(flt_is_near_unittest, PrintNumbers) {
+// This isn't really a test. It's purpose is to print interesting values.
+   bool result;
+   float min_tolerance = FLT_MIN/FLT_EPSILON;
+   std::cout << "Minimum tolerance = " << min_tolerance << std::endl;
+   EXPECT_TRUE(true);
+}