[lmi-commits] [lmi] odd/1e24-vs-1e24L b11d029: Extended precision and floating literals of type double

[Greg Chicares]

branch: odd/1e24-vs-1e24L
commit b11d0298502696fe0ff8ed72b05dce13f714e611
Author: Gregory W. Chicares <address@hidden>
Commit: Gregory W. Chicares <address@hidden>

    Extended precision and floating literals of type double
    
    log10(10^24)/3.0 equals 8, but truncating that 8 seemed to make it 7.
    See:
      https://lists.nongnu.org/archive/html/lmi/2018-03/msg00080.html
---
 sandbox_test.cpp | 87 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 87 insertions(+)

diff --git a/sandbox_test.cpp b/sandbox_test.cpp
index 0ed05ae..0ef72a3 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 xyzzy = 999'999'999'999'999.0;
+    double const a = std::log10(xyzzy);
+    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;
 }