Add algorithm for difference of squares and move `is_monotonic_increa…

…sing()` to utils (celeritas-project#1082)

* Move is_monotonic_increasing() to utils and add test

* Add algorithm for difference of squares

See  celeritas-project#1080

src/celeritas/em/msc/detail/UrbanMscScatter.hh

@@ -141,7 +141,7 @@ CELER_FUNCTION real_type UrbanMscScatter::calc_displacement(real_type geom_path,
     CELER_EXPECT(true_path >= geom_path);
 
     // true^2 - geo^2
-    real_type rmax2 = (true_path - geom_path) * (true_path + geom_path);
+    real_type rmax2 = diffsq(true_path, geom_path);
 
     // 0.73 is (roughly) the expected value of a distribution of the mean
     // radius given rmax "based on single scattering results"

src/celeritas/grid/ValueGridBuilder.cc

@@ -18,6 +18,7 @@
 #include "corecel/math/SoftEqual.hh"
 
 #include "ValueGridInserter.hh"
+#include "VectorUtils.hh"
 
 namespace celeritas
 {
@@ -65,20 +66,6 @@ bool is_on_grid_point(double value, double lo, double hi, size_type size)
     return soft_mod(value - lo, delta);
 }
 
-bool is_monotonic_increasing(SpanConstDbl grid)
-{
-    CELER_EXPECT(!grid.empty());
-    auto iter = grid.begin();
-    auto prev = *iter++;
-    while (iter != grid.end())
-    {
-        if (*iter <= prev)
-            return false;
-        prev = *iter++;
-    }
-    return true;
-}
-
 //---------------------------------------------------------------------------//
 }  // namespace
 

src/celeritas/grid/VectorUtils.cc

@@ -67,5 +67,23 @@ std::vector<double> logspace(double start, double stop, size_type n)
     return space_impl<Interp::log>(start, stop, n);
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * True if the grid values are monotonically increasing.
+ */
+bool is_monotonic_increasing(Span<double const> grid)
+{
+    CELER_EXPECT(!grid.empty());
+    auto iter = grid.begin();
+    auto prev = *iter++;
+    while (iter != grid.end())
+    {
+        if (*iter <= prev)
+            return false;
+        prev = *iter++;
+    }
+    return true;
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/celeritas/grid/VectorUtils.hh

@@ -11,6 +11,7 @@
 #include <vector>
 
 #include "corecel/Types.hh"
+#include "corecel/cont/Span.hh"
 
 namespace celeritas
 {
@@ -22,5 +23,9 @@ std::vector<double> linspace(double start, double stop, size_type n);
 // Return logarithmically spaced numbers over a specific interval
 std::vector<double> logspace(double start, double stop, size_type n);
 
+//---------------------------------------------------------------------------//
+// True if the grid values are monotonically increasing
+bool is_monotonic_increasing(Span<double const> grid);
+
 //---------------------------------------------------------------------------//
 }  // namespace celeritas

src/corecel/math/Algorithms.hh

@@ -493,6 +493,20 @@ template<class T>
     return T{0} - value;
 }
 
+//---------------------------------------------------------------------------//
+/*!
+ * Calculate the difference of squares \f$ a^2 - b^2 \f$.
+ *
+ * The expression \f$ a^2 - b^2 \f$ exhibits catastrophic cancellation when \f$
+ * a \f$ and \f$ b \f$ are close in magnitude. This can be avoided by
+ * rearranging the formula as \f$ (a - b)(a + b) \f$.
+ */
+template<class T>
+CELER_CONSTEXPR_FUNCTION T diffsq(T a, T b)
+{
+    return (a - b) * (a + b);
+}
+
 //---------------------------------------------------------------------------//
 /*!
  * Double-precision math constant (POSIX derivative).

test/celeritas/grid/VectorUtils.test.cc

@@ -7,6 +7,8 @@
 //---------------------------------------------------------------------------//
 #include "celeritas/grid/VectorUtils.hh"
 
+#include <vector>
+
 #include "celeritas_test.hh"
 
 constexpr double exact_third = 1.0 / 3.0;
@@ -80,6 +82,25 @@ TEST(VectorUtils, logspace)
         EXPECT_DOUBLE_EQ(2 * exact_third, result.back());
     }
 }
+
+//---------------------------------------------------------------------------//
+
+TEST(VectorUtils, monotonic_increasing)
+{
+    {
+        std::vector<double> v{2, 4, 8, 16, 32};
+        EXPECT_TRUE(is_monotonic_increasing(make_span(v)));
+    }
+    {
+        std::vector<double> v{10, 100, 1000, 1000};
+        EXPECT_FALSE(is_monotonic_increasing(make_span(v)));
+    }
+    {
+        std::vector<double> v{1e-16, 0, 1, 2};
+        EXPECT_FALSE(is_monotonic_increasing(make_span(v)));
+    }
+}
+
 //---------------------------------------------------------------------------//
 }  // namespace test
 }  // namespace celeritas

test/corecel/math/Algorithms.test.cc

@@ -328,6 +328,18 @@ TEST(MathTest, negate)
 
 //---------------------------------------------------------------------------//
 
+TEST(MathTest, diffsq)
+{
+    EXPECT_DOUBLE_EQ(9.0, diffsq(5.0, 4.0));
+    EXPECT_DOUBLE_EQ(ipow<2>(std::sin(0.2)), diffsq(1.0, std::cos(0.2)));
+
+    float a{10000.001}, b{10000}, actual{20};
+    EXPECT_FLOAT_EQ(0.46875f, actual - diffsq(a, b));
+    EXPECT_LE(actual - diffsq(a, b), actual - (a * a - b * b));
+}
+
+//---------------------------------------------------------------------------//
+
 TEST(MathTest, sincos)
 {
     {