From 043022622989777b31bea3aa6612c1540416d9dd Mon Sep 17 00:00:00 2001 From: Jakub Trzebiatowski Date: Sun, 5 Apr 2026 21:48:12 +0200 Subject: [PATCH] tdf#129606 move approxDiff to rtl::math and fix ARM64 rounding errors Move the approxDiff heuristic from the sc unnamed namespace to rtl::math and disable floating-point contraction (FMA) for this function. On ARM64 (Apple Silicon), the use of Fused Multiply-Add (FMA) results in higher intermediate precision during floating-point operations. While mathematically more accurate, this prevents legacy rounding heuristics in approxDiff and approxValue from detecting and cleaning rounding artifacts in series generation (e.g., 1.1 - 1.0). Disabling fp-contract forces intermediate rounding to 64-bit IEEE 754 limits. This ensures the noise expected by the approximation logic is present, restoring consistency with x86-64 behavior and fixing issues where values like 6.0 are incorrectly rendered as 6.00000000000001. Change-Id: I54a0f748be6cf817db1f3bdcfd2b355320c3a617 Reviewed-on: https://gerrit.libreoffice.org/c/core/+/203233 Reviewed-by: Dan Williams Tested-by: Jenkins --- include/rtl/math.h | 6 +++++ include/rtl/math.hxx | 5 ++++ sal/qa/rtl/math/test-rtl-math.cxx | 17 ++++++++++++ sal/rtl/math.cxx | 42 +++++++++++++++++++++++++++++ sal/util/sal.map | 5 ++++ sc/source/core/data/table4.cxx | 44 +++---------------------------- 6 files changed, 79 insertions(+), 40 deletions(-) diff --git a/include/rtl/math.h b/include/rtl/math.h index 1938f037c779..149c5cb0f419 100644 --- a/include/rtl/math.h +++ b/include/rtl/math.h @@ -426,6 +426,12 @@ SAL_DLLPUBLIC double SAL_CALL rtl_math_pow10Exp(double fValue, int nExp) SAL_THR */ SAL_DLLPUBLIC double SAL_CALL rtl_math_approxValue(double fValue) SAL_THROW_EXTERN_C(); + +/** Calculate a-b trying to diminish precision errors such as for 0.11-0.12 + not return -0.009999999999999995 but -0.01 instead. + */ +SAL_DLLPUBLIC double SAL_CALL rtl_math_approxDiff(double a, double b) SAL_THROW_EXTERN_C(); + /** Test equality of two values with an accuracy of the magnitude of the given values scaled by 2^-48 (4 bits roundoff stripped). diff --git a/include/rtl/math.hxx b/include/rtl/math.hxx index 604cefcb6b2c..b64e44fad907 100644 --- a/include/rtl/math.hxx +++ b/include/rtl/math.hxx @@ -269,6 +269,11 @@ inline double approxValue(double fValue) return rtl_math_approxValue(fValue); } +inline double approxDiff(double a, double b) +{ + return rtl_math_approxDiff(a, b); +} + /** A wrapper around rtl_math_expm1. */ inline double expm1(double fValue) diff --git a/sal/qa/rtl/math/test-rtl-math.cxx b/sal/qa/rtl/math/test-rtl-math.cxx index 8741623c8648..a537bc00e63d 100644 --- a/sal/qa/rtl/math/test-rtl-math.cxx +++ b/sal/qa/rtl/math/test-rtl-math.cxx @@ -38,6 +38,22 @@ namespace { class Test: public CppUnit::TestFixture { public: + + /** * Tests approxDiff against FMA-induced precision drift on ARM64. + * FMA skips intermediate IEEE 754 rounding, preventing the bit-wise + * divergence needed by the rescaling heuristic to detect and fix + * rounding errors. Accumulating the result 50x makes this drift detectable. + */ + void test_approxDiff() + { + double a = 1.1; + double b = 1.0; + + double diff = rtl::math::approxDiff(a, b); + + CPPUNIT_ASSERT_EQUAL(5.0, diff * 50); + } + void test_stringToDouble_good() { rtl_math_ConversionStatus status; sal_Int32 end; @@ -690,6 +706,7 @@ public: } CPPUNIT_TEST_SUITE(Test); + CPPUNIT_TEST(test_approxDiff); CPPUNIT_TEST(test_stringToDouble_good); CPPUNIT_TEST(test_stringToDouble_bad); CPPUNIT_TEST(test_stringToDouble_exponent_without_digit); diff --git a/sal/rtl/math.cxx b/sal/rtl/math.cxx index 78df29105c25..9198cfe46921 100644 --- a/sal/rtl/math.cxx +++ b/sal/rtl/math.cxx @@ -601,6 +601,48 @@ double SAL_CALL rtl_math_approxValue(double fValue) noexcept return bSign ? -fValue : fValue; } +/* FMA is disabled to force standard IEEE 754 rounding after each multiplication. + On ARM64, FMA's higher precision prevents the bit-level divergence between + 'd' and 'c' that this heuristic relies on. Without this "noise", the logic + fails to detect and fix subtle rounding artifacts, leading to cumulative + errors in Calc series. +*/ +#if defined(__clang__) +#pragma clang fp contract(off) +#elif defined(__GNUC__) +#pragma GCC optimize("fp-contract=off") +#endif +double SAL_CALL rtl_math_approxDiff(double a, double b) noexcept +{ + if (a == b) + return 0.0; + if (a == 0.0) + return -b; + if (b == 0.0) + return a; + const double c = a - b; + const double aa = fabs(a); + const double ab = fabs(b); + if (aa < 1e-16 || aa > 1e+16 || ab < 1e-16 || ab > 1e+16) + // This is going nowhere, live with the result. + return c; + + const double q = aa < ab ? b / a : a / b; + const double d = (a * q - b * q) / q; + if (d == c) + // No differing error, live with the result. + return c; + + // We now have two subtractions with a similar but not equal error. Obtain + // the exponent of the error magnitude and round accordingly. + const double e = fabs(d - c); + const int nExp = static_cast(floor(log10(e))) + 1; + // tdf#129606: Limit precision to the 16th significant digit of the least precise argument. + // Cf. mnMaxGeneralPrecision in sc/source/core/data/column3.cxx. + const int nExpArg = static_cast(floor(log10(std::max(aa, ab)))) - 15; + return rtl::math::round(c, -std::max(nExp, nExpArg)); +} + bool SAL_CALL rtl_math_approxEqual(double a, double b) noexcept { // threshold is last four bits of mantissa: diff --git a/sal/util/sal.map b/sal/util/sal.map index c5c3e4d55641..2c83e0833065 100644 --- a/sal/util/sal.map +++ b/sal/util/sal.map @@ -705,6 +705,11 @@ LIBO_UDK_6.2 { # symbols available in >= LibO 6.2 osl_replaceFile; } LIBO_UDK_5.3; +LIBO_UDK_26.8 { # symbols available in >= LibO 26.8 + global: + rtl_math_approxDiff; +} LIBO_UDK_6.2; + PRIVATE_1.0 { global: osl_detail_ObjectRegistry_storeAddresses; diff --git a/sc/source/core/data/table4.cxx b/sc/source/core/data/table4.cxx index b8bbe4c7bd99..b8c275597be0 100644 --- a/sc/source/core/data/table4.cxx +++ b/sc/source/core/data/table4.cxx @@ -190,42 +190,6 @@ void setSuffixCell( } namespace { -/* TODO: move this to rtl::math::approxDiff() ? Though the name is funny, the - * approx is expected to be more correct than the raw diff. */ -/** Calculate a-b trying to diminish precision errors such as for 0.11-0.12 - not return -0.009999999999999995 but -0.01 instead. - */ -double approxDiff( double a, double b ) -{ - if (a == b) - return 0.0; - if (a == 0.0) - return -b; - if (b == 0.0) - return a; - const double c = a - b; - const double aa = fabs(a); - const double ab = fabs(b); - if (aa < 1e-16 || aa > 1e+16 || ab < 1e-16 || ab > 1e+16) - // This is going nowhere, live with the result. - return c; - - const double q = aa < ab ? b / a : a / b; - const double d = (a * q - b * q) / q; - if (d == c) - // No differing error, live with the result. - return c; - - // We now have two subtractions with a similar but not equal error. Obtain - // the exponent of the error magnitude and round accordingly. - const double e = fabs(d - c); - const int nExp = static_cast(floor(log10(e))) + 1; - // tdf#129606: Limit precision to the 16th significant digit of the least precise argument. - // Cf. mnMaxGeneralPrecision in sc/source/core/data/column3.cxx. - const int nExpArg = static_cast(floor(log10(std::max(aa, ab)))) - 15; - return rtl::math::round(c, -std::max(nExp, nExpArg)); -} - double approxTypedDiff( double a, double b, bool bTime, tools::Duration& rDuration ) { if (bTime) @@ -233,7 +197,7 @@ double approxTypedDiff( double a, double b, bool bTime, tools::Duration& rDurati rDuration = tools::Duration(a - b); return rDuration.GetInDays(); } - return approxDiff( a, b); + return rtl::math::approxDiff(a, b); } } @@ -508,7 +472,7 @@ void ScTable::FillAnalyse( SCCOL nCol1, SCROW nRow1, SCCOL nCol2, SCROW nRow2, HeadNoneTail nFlag2 = lcl_DecompValueString(aStr2, nVal2, &rMinDigits); if (nFlag1 == nFlag2 && aStr == aStr2) { - double nDiff = approxDiff(nVal2, nVal1); + double nDiff = rtl::math::approxDiff(nVal2, nVal1); if (i == 1) rInc = nDiff; else if (!::rtl::math::approxEqual(nDiff, rInc, 13)) @@ -738,7 +702,7 @@ void ScTable::FillAnalyse( SCCOL nCol1, SCROW nRow1, SCCOL nCol2, SCROW nRow2, HeadNoneTail nFlag2 = lcl_DecompValueString(aStr, nVal2, &rMinDigits); if ( nFlag1 == nFlag2 ) { - rInc = approxDiff( nVal2, nVal1); + rInc = rtl::math::approxDiff( nVal2, nVal1); nCol = sal::static_int_cast( nCol + nAddX ); nRow = sal::static_int_cast( nRow + nAddY ); bool bVal = true; @@ -752,7 +716,7 @@ void ScTable::FillAnalyse( SCCOL nCol1, SCROW nRow1, SCCOL nCol2, SCROW nRow2, nFlag2 = lcl_DecompValueString( aStr, nVal2, &rMinDigits ); if ( nFlag1 == nFlag2 ) { - double nDiff = approxDiff( nVal2, nVal1); + double nDiff = rtl::math::approxDiff( nVal2, nVal1); if ( !::rtl::math::approxEqual( nDiff, rInc, 13 ) ) bVal = false; nVal1 = nVal2;