[lmi] 7 == log10(10^24)/3.0

[Greg Chicares]

TL;DR:
     1'000'000'000'000'000'000'000'000.0
  != 1'000'000'000'000'000'000'000'000.0L
(but there's an actual git-style question below).

I was unable to create a simple test case for the oddity
in commit 1c1bafa402cb625a465ce6d87ab72185bcf5eca1
(which, as it turned out, was due to a vain hope that
because the floating value 1000.0 is an exact integer,
1.0/(1.0/1000.0) would also be exact)...

    Further document a MinGW-w64 gcc-7.2.0 anomaly
    
    floor() and trunc() yield different results for a particular floating-
    point argument slightly different from positive three. To reproduce:

...so it was interesting to come across a similar symptom
while developing unit tests for related code (which doesn't
have the same 1.0/(1.0/1000.0) problem), for which I _was_
able to create a simple test case--and which turned out to
be attributable to the different vain hope TL;DR'd above.

It cost me some time to figure this out, so, not wanting to
throw the work away, I thought I'd post it here. Vadim, would
it be in poor taste to commit something like this to the
savannah repository as a forlorn branch, to be abandoned
immediately after committing it?

---------8<--------8<--------8<--------8<--------8<--------8<--------8<-------
diff --git a/sandbox_test.cpp b/sandbox_test.cpp
index 0ed05aed..801d29fe 100644
--- a/sandbox_test.cpp
+++ b/sandbox_test.cpp
@@ -23,8 +23,95 @@
 
 #include "test_tools.hpp"
 
+#include <cmath>                        // floor(), log10()
+
+#include <iomanip>
+#include <ios>
+#include <sstream>
+#include <string>
+
+template<typename T>
+std::string floating_rep(T t)
+{
+    std::ostringstream oss;
+    oss << std::hex << std::setfill('0');
+    int realsize = sizeof(T);
+#if defined __GNUC__
+    if(12 == realsize) realsize = 10;
+#endif // defined __GNUC__
+    unsigned char const* u = reinterpret_cast<unsigned char const*>(&t);
+    for(int j = realsize - 1; 0 <= j; --j)
+        oss << std::setw(2) << static_cast<int>(u[j]);
+    return oss.str();
+}
+
 int test_main(int, char*[])
 {
+    std::cout
+        << "Unsurprisingly, log10(10^24)/3.0 equals 8,"
+        << "\nbut truncating that 8 in any of several ways makes it 7:"
+        << "\n" << std::endl
+        ;
+    double const ten24 = 1'000'000'000'000'000'000'000'000.0;
+    double const a = std::log10(ten24);
+    double const b = a / 3.0;
+    double const c = std::floor(b);
+    double const d = std::trunc(b);
+    double const e = static_cast<int>(b);
+    std::cout << "a: " << a << " ; b: " << b << " ; c: " << c << std::endl;
+    std::cout << "a: " << a << ' ' << floating_rep(a) << std::endl;
+    std::cout << "b: " << b << ' ' << floating_rep(b) << std::endl;
+    std::cout << "c: " << c << ' ' << floating_rep(c) << std::endl;
+    std::cout << "d: " << d << ' ' << floating_rep(d) << std::endl;
+    std::cout << "e: " << e << ' ' << floating_rep(e) << std::endl;
+
+    {
+    std::cout
+        << "\nThat seemed to be an excess-precision problem;"
+        << "\nyet it still occurs with type 'long double'"
+        << "\n(but the hex representation of that 8 looks weird):"
+        << "\n" << std::endl
+        ;
+    long double const t = 1'000'000'000'000'000'000'000'000.0;
+    long double const x =           (std::log10(t) / 3.0);
+    long double const y = std::trunc(std::log10(t) / 3.0);
+    std::cout << "x: " << x << ' ' << floating_rep(x) << std::endl;
+    std::cout << "y: " << y << ' ' << floating_rep(y) << std::endl;
+    }
+    {
+    std::cout
+        << "\nIt still occurs with even with 'long double' C functions"
+        << "\n(notably using log10l() instead of log10()):"
+        << "\n" << std::endl
+        ;
+    long double const t = 1'000'000'000'000'000'000'000'000.0;
+    long double const u = 1'000'000'000'000'000'000'000'000.0L;
+    long double const x =            (   log10l(t) / 3.0);
+    long double const y =      truncl(   log10l(t) / 3.0);
+    std::cout << "t: " << t << ' ' << floating_rep(t) << std::endl;
+    std::cout << "u: " << u << ' ' << floating_rep(u) << std::endl;
+    std::cout << "x: " << x << ' ' << floating_rep(x) << std::endl;
+    std::cout << "y: " << y << ' ' << floating_rep(y) << std::endl;
+    }
+    {
+    std::cout
+        << "\nThe explanation is that 10^24 is inexact: everything"
+        << "\nseems to work as hoped if all arithmetic is performed"
+        << "\nin long double precision, with long double library"
+        << "\nfunctions, on a long double value initialized from"
+        << "\na long double floating literal:"
+        << "\n" << std::endl
+        ;
+    long double const t = 1'000'000'000'000'000'000'000'000.0;
+    long double const u = 1'000'000'000'000'000'000'000'000.0L;
+    long double const x =            (   log10l(u) / 3.0);
+    long double const y =      truncl(   log10l(u) / 3.0);
+    std::cout << "t: " << t << ' ' << floating_rep(t) << std::endl;
+    std::cout << "u: " << u << ' ' << floating_rep(u) << std::endl;
+    std::cout << "x: " << x << ' ' << floating_rep(x) << std::endl;
+    std::cout << "y: " << y << ' ' << floating_rep(y) << std::endl;
+    }
+
     return 0;
 }
 
--------->8-------->8-------->8-------->8-------->8-------->8-------->8-------