[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[lmi-commits] [lmi] master 3ebab874 8/9: Test fdlibm against C RTL
From: |
Greg Chicares |
Subject: |
[lmi-commits] [lmi] master 3ebab874 8/9: Test fdlibm against C RTL |
Date: |
Sat, 21 May 2022 20:13:52 -0400 (EDT) |
branch: master
commit 3ebab8745288a786a377de83bdc1aa742223791d
Author: Gregory W. Chicares <gchicares@sbcglobal.net>
Commit: Gregory W. Chicares <gchicares@sbcglobal.net>
Test fdlibm against C RTL
glibc's expm1() and log1p() are adapted from fdlibm, but they don't
always give exactly the same answer. This could be due to differences
in gcc version or options, or to some change glibc made. Summary:
x86_64, sse, gcc, posix
test_expm1_log1p: 1000000 trials
1129 errors [worst 0.99847825656009148165282 ulp] in expm1()
2155 errors [worst 0.99987376777519854087473 ulp] in log1p()
2031 errors [worst 1.99543305292237360681895 ulp] in monthly int
1818 errors [worst 1.94732362992901619769270 ulp] in daily int
The MinGW-w64 implementation disagrees with fdlibm more often:
x86_64, sse, gcc, msw
test_expm1_log1p: 1000000 trials
627369 errors [worst 1.60640547797468680180089 ulp] in expm1()
151389 errors [worst 0.99999911531551366472570 ulp] in log1p()
705552 errors [worst 2.46837813563044816689285 ulp] in monthly int
711194 errors [worst 2.14128061677744607749219 ulp] in daily int
but the worst discrepancy seems to be...one and a half ulps for
expm1(), which suggests that the casual ulp calculation is imperfect.
Of course, the monthly and daily interest calculations involve calls
to both expm1() and log1p(), so their results are more discrepant.
---
math_functions_test.cpp | 98 ++++++++++++++++++++++++++++++++++++++++++++++++-
1 file changed, 97 insertions(+), 1 deletion(-)
diff --git a/math_functions_test.cpp b/math_functions_test.cpp
index 986852ee..bf1e8128 100644
--- a/math_functions_test.cpp
+++ b/math_functions_test.cpp
@@ -31,7 +31,8 @@
#include "timer.hpp"
#include <algorithm> // min()
-#include <cmath> // isnan(), pow()
+#include <cfloat> // DBL_EPSILON
+#include <cmath> // fabs(), isnan(), pow()
#include <iomanip>
#include <limits>
#include <type_traits>
@@ -528,6 +529,101 @@ void test_expm1_log1p()
LMI_TEST_EQUAL(0.009950330853168083, y);
// digits 12345678901234567
LMI_TEST_EQUAL(0.0032737397821988637, z);
+
+ // For sampled (integer/1000000.0) arguments in the open range
+ // ]-0.043348, +0.042151[
+ // lmi's fdlibm implementation of expm1() and log1p() matches glibc's
+ // except for this single example:
+ double const g0 {25610 / 1000000.0};
+ double const g1 = lmi::expm1(g0); // 0.02594075354662067622868
+ double const g2 = std::expm1(g0); // 0.02594075354662067275924
+// LMI_TEST_EQUAL(g1, g2); fails
+ LMI_TEST(materially_equal(g1, g2));
+ // and in that example glibc is correct:
+ // https://www.wolframalpha.com/input?i=exp(25610/1000000)-1
+ // 0.0259407535466206736037231992174016233440736931692437771090797988...
+ // https://babbage.cs.qc.cuny.edu/IEEE-754.old/Decimal.html
+ // 0.025940753546620673 = 3F9A903680771FAF correctly rounded
+ // 0.025940753546620676 = 3F9A903680771FB0 fdlibm
+ // 0.025940753546620673 = 3F9A903680771FAF glibc
+
+ // Absolute value of relative error.
+ auto rel_err = [](double t, double u)
+ {
+ return std::fabs(t - u) / std::min(std::fabs(t), std::fabs(u));
+ };
+
+ // Test fdlibm vs. C RTL for many parameters.
+ int err_count0 {0};
+ int err_count1 {0};
+ int err_count2 {0};
+ int err_count3 {0};
+ double err_max0 {0.0};
+ double err_max1 {0.0};
+ double err_max2 {0.0};
+ double err_max3 {0.0};
+ std::cout.setf(std::ios_base::fixed, std::ios_base::floatfield);
+ std::cout.precision(23);
+ int const how_many = 1000000;
+ for(int i = 1 - how_many; i < how_many; ++i)
+ {
+ // interest rate
+ double const irate = i / (1.0 * 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);
+ // relative error
+ double const e0 = rel_err(a0, b0);
+ double const e1 = rel_err(a1, b1);
+ double const e2 = rel_err(a2, b2);
+ double const e3 = rel_err(a3, b3);
+ // comparison
+ if(a0 != b0 || a1 != b1 || a2 != b2 || a3 != b3)
+ {
+ err_count0 += a0 != b0;
+ err_count1 += a1 != b1;
+ err_count2 += a2 != b2;
+ err_count3 += a3 != b3;
+ err_max0 = std::max(err_max0, e0);
+ err_max1 = std::max(err_max1, e1);
+ err_max2 = std::max(err_max2, e2);
+ err_max3 = std::max(err_max3, e3);
+// Enable this to show each one:
+#if 0
+ std::cout << "unequal" << std::endl;
+ std::cout
+ << i << ' ' << irate << ' ' << how_many << '\n'
+ << a0 << ' ' << a1 << ' ' << a2 << ' ' << a3 << '\n'
+ << b0 << ' ' << b1 << ' ' << b2 << ' ' << b3 << '\n'
+ << e0 << ' ' << e1 << ' ' << e2 << ' ' << e3 << '\n'
+ << std::endl
+ ;
+#endif // 0
+ }
+ }
+ std::cout
+ << __func__ << ": " << how_many << " trials\n"
+ << " " << err_count0 << " errors"
+ << " [worst " << err_max0 / DBL_EPSILON << " ulp]"
+ << " in expm1()\n"
+ << " " << err_count1 << " errors"
+ << " [worst " << err_max1 / DBL_EPSILON << " ulp]"
+ << " in log1p()\n"
+ << " " << err_count2 << " errors"
+ << " [worst " << err_max2 / DBL_EPSILON << " ulp]"
+ << " in monthly int\n"
+ << " " << err_count3 << " errors"
+ << " [worst " << err_max3 / DBL_EPSILON << " ulp]"
+ << " in daily int\n"
+ << std::endl
+ ;
}
/// This function isn't a unit test per se. Its purpose is to show
- [lmi-commits] [lmi] master updated (926d9f91 -> 82468e66), Greg Chicares, 2022/05/21
- [lmi-commits] [lmi] master b52923fa 2/9: Make further use of values known to be correctly rounded, Greg Chicares, 2022/05/21
- [lmi-commits] [lmi] master fdf41a42 4/9: Reorganize unit tests, Greg Chicares, 2022/05/21
- [lmi-commits] [lmi] master 82468e66 9/9: Show architectural context at start of main test function, Greg Chicares, 2022/05/21
- [lmi-commits] [lmi] master e59c42d3 7/9: Rearrange unit-test functions in invocation order, Greg Chicares, 2022/05/21
- [lmi-commits] [lmi] master d393cd47 3/9: Move a block to make an incipient refactoring simpler, Greg Chicares, 2022/05/21
- [lmi-commits] [lmi] master ff833df3 1/9: Verify against independent calculations, Greg Chicares, 2022/05/21
- [lmi-commits] [lmi] master 0c667ac5 5/9: Move an output line, Greg Chicares, 2022/05/21
- [lmi-commits] [lmi] master 3ebab874 8/9: Test fdlibm against C RTL,
Greg Chicares <=
- [lmi-commits] [lmi] master 4e300fbc 6/9: Rearrange a unit test's output, Greg Chicares, 2022/05/21