[lmi-commits] [lmi] master 1718c057 1/2: For immediate reversion: long double idea

[Greg Chicares]

branch: master
commit 1718c057cb25acf319b8665a9a8ca5fc4383b4c2
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>

    For immediate reversion: long double idea
    
    Instead of testing lmi's fdlibm adaptation (which uses double precision)
    against the RTL's double-precision expm1() and log1p(), why not test it
    against the RTL's extended-precision expm1() and log1p()?
    
    At best, this could be useful, assuming that the extended-precision
    implementation is more accurate, and assuming that casting the values
    it returns to double (once and only once, as the final step) satisfies
    the technical definition of "correctly rounded", which C99 does not
    guarantee.
    
    At worst, if those assumptions aren't both true, it's not useful.
    But it is an intriguing idea nonetheless.
    
    x86_64, sse, gcc, posix
    
    test_expm1_log1p: 1000000 trials
      181250 errors [worst 0.99998232458484737072979 ulp] in expm1()
      151400 errors [worst 0.99999911531551366472570 ulp] in log1p()
      705317 errors [worst 1.56317871631794114151148 ulp] in monthly int
      729473 errors [worst 1.49460122354889279883139 ulp] in daily int
    
      lmi::expm1()     7.535e-04 s mean;        750 us least of 100 runs
      std::expm1()     7.806e-04 s mean;        752 us least of 100 runs
      lmi::log1p()     7.717e-04 s mean;        770 us least of 100 runs
      std::log1p()     7.999e-04 s mean;        797 us least of 100 runs
    
    As previously demonstrated, lmi's fdlibm adaptation reproduces glibc's
    results exactly for all 1,000,000 trials, with nearly identical speed.
    Here, comparing to glibc's long double implementation instead, the
    speed is still nearly identical: does that suggest that glibc uses the
    same code for both precisions, since extended precision generally takes
    longer? And the numerical results differ, but only in 15-20% of the
    trials, and never by more than one ulp; could it be that the apparent
    errors really just show that static_cast<double>(long double) does not
    return a "correctly rounded" result?
    
    Results for MinGW-w64 look very much like results before this revision.
    The code is easily modified to compare double- and extended-precision
    RTL (ignoring lmi's fdlibm); that experiment shows discrepancies like
    those above for pc-linux-gnu, but zero discrepancies for MinGW-w64.
    Does that suggest that they coded something inaccurate like
      long double expm1(long double z) {return (double)expm1((double)z);}
    ? No:
      
https://github.com/mingw-w64/mingw-w64/blob/master/mingw-w64-crt/math/x86/expm1.def.h
    They use the same extended-precision code for all floating-point types,
    and for finite non-NaNs it's
          x /= __FLT_LOGE2;
          __asm__ __volatile__ ("f2xm1" : "=t" (x) : "0" (x));
          return x;
    That's almost twice as slow as fdlibm, and log1p() is about five times
    as slow. It was a brilliant idea in 1980 to implement almost correctly-
    rounded math functions in silicon; unfortunately, they use slow silicon.
    
    And even a brilliant tool can do harm in the wrong hands: the original
    release used the wrong arithmetic operation, requiring this correction:
    -     x *= __FLT_LOGE2;
    +     x /= __FLT_LOGE2;
    as shown here:
      https://sourceforge.net/p/mingw-w64/mailman/message/27298327/
      
https://mingw-w64-public.narkive.com/CxmokaLr/expm1-is-giving-incorrect-results
    so one might doubt the rigor of their testing.
    
    Anyway, this seemed too interesting to put on a branch where it'll get
    lost, though it'll be reverted immediately because it's not useful to
    keep it in master.
---
 math_functions_test.cpp | 18 ++++++++++++------
 1 file changed, 12 insertions(+), 6 deletions(-)

diff --git a/math_functions_test.cpp b/math_functions_test.cpp
index 354c60d0..e6254541 100644
--- a/math_functions_test.cpp
+++ b/math_functions_test.cpp
@@ -595,6 +595,9 @@ void test_expm1_log1p()
         return std::fabs(t - u) / std::min(std::fabs(t), std::fabs(u));
         };
 
+    // Keeps lines short.
+    auto dbl = [](long double z) {return static_cast<double>(z);};
+
     // Test fdlibm vs. C RTL for many parameters.
     int    err_count0 {0};
     int    err_count1 {0};
@@ -610,17 +613,20 @@ void test_expm1_log1p()
     for(int i = 1 - how_many; i < how_many; ++i)
         {
         // interest rate
-        double const irate = i / (1.0 * how_many);
+             double const  irate = i / (1.0  * how_many);
+        long double const lirate = irate; // identical value
+        // this would not be identical:
+        //                       = i / (1.0L * how_many);
         // fdlibm
         double const a0 = lmi::expm1(irate);
         double const a1 = lmi::log1p(irate);
         double const a2 = lmi::expm1(lmi::log1p(irate) / 12);
         double const a3 = lmi::expm1(lmi::log1p(irate) / 365);
-        // RTL
-        double const b0 = std::expm1(irate);
-        double const b1 = std::log1p(irate);
-        double const b2 = std::expm1(std::log1p(irate) / 12);
-        double const b3 = std::expm1(std::log1p(irate) / 365);
+        // RTL, long double precision, stored as double
+        double const b0 = dbl(std::expm1(lirate));
+        double const b1 = dbl(std::log1p(lirate));
+        double const b2 = dbl(std::expm1(std::log1p(lirate) / 12));
+        double const b3 = dbl(std::expm1(std::log1p(lirate) / 365));
         // relative error
         double const e0 = rel_err(a0, b0);
         double const e1 = rel_err(a1, b1);