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 <dcbw@libreoffice.org>
Tested-by: Jenkins

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).
 

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)

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);

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<int>(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<int>(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:

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;

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<int>(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<int>(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<SCCOL>( nCol + nAddX );
                     nRow = sal::static_int_cast<SCROW>( 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;