[lmi] Converting cents to dollars [Was: Contradictory performance measurements]

[Greg Chicares]

On 4/8/21 1:10 PM, Vadim Zeitlin wrote:
[...]
> ---------------------------------- >8 --------------------------------------
> diff --git a/currency.hpp b/currency.hpp
> index ce04d75a2..a73a2b627 100644
> --- a/currency.hpp
> +++ b/currency.hpp
> @@ -40,6 +40,7 @@ class currency
> 
>      static constexpr int    cents_digits     = 2;
>      static constexpr double cents_per_dollar = 100.0;
> +    static constexpr double cents_per_dollar_inv = 0.01;
> 
>    public:
>      using data_type = double;
> @@ -58,8 +59,7 @@ class currency
> 
>      data_type cents() const {return m_;}
>      // CURRENCY !! add a unit test for possible underflow
> -    // CURRENCY !! is multiplication by reciprocal faster or more accurate?
> -    double d() const {return m_ / cents_per_dollar;}
> +    double d() const {return m_ * cents_per_dollar_inv;}
> 
>    private:
>      explicit currency(data_type z, raw_cents) : m_ {z} {}
> ---------------------------------- >8 --------------------------------------

As I awoke this morning, I thought we might perhaps combine this
idea with something nearly similar in 'round_to.hpp':

/// Division by an exact integer value would afford a stronger
/// guarantee of accuracy, particularly when the value to be rounded
/// is already rounded--e.g., if 3.00 is to be rounded to hundredths,
/// then
///   (100.0 * 3.00) / 100
/// is preferable to
///   (100.0 * 3.00) * 0.01
/// especially if the rounding style is anything but to-nearest.
///
/// However, reciprocal multiplication is faster than division--see:
///   https://lists.nongnu.org/archive/html/lmi/2021-04/msg00010.html
/// et seqq., which demonstrates an average three-percent speedup,
/// with no observed difference in a comprehensive system test as
/// long as the multiplications are performed in extended precision.

IOW, instead of

> +    static constexpr double cents_per_dollar_inv = 0.01;
...
> -    double d() const {return m_ / cents_per_dollar;}
> +    double d() const {return m_ * cents_per_dollar_inv;}

above, use a 'long double' reciprocal.

It does cause regressions, so of course I'm not disposed to commit it;
but it does raise several questions, enumerated below, for discussion.

(1) I'm testing this in our i686+x87 production environment. What are
your thoughts about migrating (later) to x86_64 and using the x87
hardware only for 'long double'? I emphatically do not mean
  -mfpmath=sse+387  // NOT
because:

  https://gcc.gnu.org/onlinedocs/gcc/x86-Options.html
| Use this option with care, as it is still experimental

Instead, I mean using type 'long double' along with
  -mfpmath=sse // implicit for x86_64
which happens to use the x87 hardware for 'long double':
i.e., exactly what our current x86_64 builds do.

Now consider this very experimental patch:

--8<----8<----8<----8<----8<----8<----8<----8<----8<----8<--
diff --git a/currency.hpp b/currency.hpp
index ce04d75a2..db6d020df 100644
--- a/currency.hpp
+++ b/currency.hpp
@@ -29,8 +29,16 @@
 #include <stdexcept>                    // runtime_error
 #include <vector>
 
+#include "materially_equal.hpp"
+
+#include <cfloat>
+#include <iomanip>
+#include <iostream>
+
 class raw_cents {}; // Tag class.
 
+static constexpr long double dollars_per_cent = {0.01L};
+
 class currency
 {
     friend class currency_test;
@@ -38,8 +46,13 @@ class currency
     template<typename> friend class round_to; // explicit ctor
     friend class round_to_test;               // currency::cents_digits
 
-    static constexpr int    cents_digits     = 2;
-    static constexpr double cents_per_dollar = 100.0;
+    static constexpr      int    cents_digits     = 2;
+    static constexpr      double cents_per_dollar = 100.0;
+//  static constexpr long double dollars_per_cent = 1.0L / cents_per_dollar;
+//  static constexpr long double dollars_per_cent = 1.0L / 100.0L;
+//  static long double const     dollars_per_cent = 1.0L / cents_per_dollar;
+//  static long double const     dollars_per_cent = 1.0L / 100.0L;
+//  static long double const     dollars_per_cent = {0.01L};
 
   public:
     using data_type = double;
@@ -58,8 +71,27 @@ class currency
 
     data_type cents() const {return m_;}
     // CURRENCY !! add a unit test for possible underflow
-    // CURRENCY !! is multiplication by reciprocal faster or more accurate?
-    double d() const {return m_ / cents_per_dollar;}
+//  // CURRENCY !! is multiplication by reciprocal faster or more accurate?
+//  double d() const {return m_ / cents_per_dollar;}
+//  double d() const {return static_cast<double>(m_ * dollars_per_cent);}
+    double d() const
+        {
+        double volatile x = m_ / cents_per_dollar;
+        double volatile y = static_cast<double>(m_ * dollars_per_cent);
+//      if(x != y)
+        if(!materially_equal(x, y, 1.0E-18L))
+            {
+            std::cout
+                << m_ << " m_\n"
+                << std::setprecision(DECIMAL_DIG)
+                << static_cast<long double>(x) << " x\n"
+                << static_cast<long double>(y) << " y\n"
+                << static_cast<long double>(x-y) << " x-y\n"
+                << std::endl
+                ;
+            }
+        return x;
+        }
 
   private:
     explicit currency(data_type z, raw_cents) : m_ {z} {}
--8<----8<----8<----8<----8<----8<----8<----8<----8<----8<--

(2) In the first part of that patch, I couldn't find a way to write
a 'long double' constexpr inside the class. I imagine this is a
gcc-8.3 or C++17 limitation that will probably be lifted in the near
future; if not, writing that constexpr at namespace scope is a good
enough expedient.

(3) Over the course of this experiment, that patch evolved into
something different from the original intention. In effect, it's
equivalent to just this:

-    double d() const {return              m_ / cents_per_dollar;}
+    double d() const {double volatile z = m_ / cents_per_dollar; return z;}

IOW, compared to HEAD, it just forces rounding to 'double'
instead of using an 80-bit x87 register for the return value.
Thus, it does what I suppose I imagined HEAD to do.

It does engender regressions:

*** System test failed ***
  1504 system-test files compared
  1117 system-test files match
  387 system-test files differ
  0 system-test files missing
...system test completed.

and it decreases accuracy, so I hesitate to commit it, but it
is illuminating. And the "if(!materially_equal" code is never
executed in that system test, which I take to mean that (fast)
reciprocal multiplication is equivalent (for this rather broad
system test) to (slow) division: that is, the regressions are
caused only by a forced 80- to 64-bit register spill.

So should we just do whatever testing is needed to accept
something like this patch? No. This is only a symptom. We
should instead spend time on the underlying problem: i.e.,
that lmi was written in terms of fractional dollars, and then
only partly migrated to integral cents. We don't want d() to
be faster--we want to render it unnecessary.