Resolves: tdf#138360 better accuracy in rtl_math_round()

Decimal negative exponents (powers) are imprecise
1e-1   0.10000000000000001
1e-2   0.01
1e-3   0.001
1e-4   0.0001
1e-5   1.0000000000000001e-05
1e-6   9.9999999999999995e-07
1e-7   9.9999999999999995e-08
1e-8   1e-08
1e-9   1.0000000000000001e-09
1e-10  1e-10
1e-11  9.9999999999999994e-12
1e-12  9.9999999999999998e-13
1e-13  1e-13
1e-14  1e-14
1e-15  1.0000000000000001e-15
1e-16  9.9999999999999998e-17
1e-17  1.0000000000000001e-17
1e-18  1.0000000000000001e-18
1e-19  9.9999999999999998e-20
1e-20  9.9999999999999995e-21

so use the positive exponents instead and swap multiplication and
division when scaling and scaling back the value, which multiplies
the rounded value with a precise integer instead of dividing it by
an imprecise fraction.

For a large (absolute) value check if it is roundable to the
desired decimals at all with the binade's distance precision gap
and if not then diminish the decimals parameter. This prevents
possible inaccuracies due to overly scaling the value.

Change-Id: I41a113078031a552cf98d72f8cb2b10bdc88dea4
Reviewed-on: https://gerrit.libreoffice.org/c/core/+/106830
Reviewed-by: Eike Rathke <erack@redhat.com>
Tested-by: Jenkins

1 files changed, 36 insertions, 3 deletions

diff --git a/sal/rtl/math.cxx b/sal/rtl/math.cxx
index 6ed4906270e0..1d6cb88327f9 100644
--- a/sal/rtl/math.cxx
+++ b/sal/rtl/math.cxx
@@ -1144,16 +1144,44 @@ double SAL_CALL rtl_math_round(double fValue, int nDecPlaces,
     if (bSign)
         fValue = -fValue;
 
+    // Rounding to decimals between integer distance precision (gaps) does not
+    // make sense, do not even try to multiply/divide and introduce inaccuracy.
+    if (nDecPlaces >= 0 && fValue >= (static_cast<sal_Int64>(1) << 52))
+        return bSign ? -fValue : fValue;
+
     double fFac = 0;
     if (nDecPlaces != 0)
     {
+        if (nDecPlaces > 1 && fValue > 4294967296.0)
+        {
+            // 4294967296 is 2^32 with room for at least 20 decimals, checking
+            // smaller values is not necessary. Lower the limit if more than 20
+            // decimals were to be allowed.
+
+            // Determine how many decimals are representable in the precision.
+            // Anything greater 2^52 and 0.0 was already ruled out above.
+            // Theoretically 0.5, 0.25, 0.125, 0.0625, 0.03125, ...
+            const double fDec = 52 - log2(fValue) + 1;
+            if (fDec < nDecPlaces)
+                nDecPlaces = static_cast<sal_Int32>(fDec);
+        }
+
+        /* TODO: this was without the inverse factor and determining max
+         * possible decimals, it could now be adjusted to be more lenient. */
         // max 20 decimals, we don't have unlimited precision
         // #38810# and no overflow on fValue*=fFac
         if (nDecPlaces < -20 || 20 < nDecPlaces || fValue > (DBL_MAX / 1e20))
             return bSign ? -fValue : fValue;
 
-        fFac = getN10Exp(nDecPlaces);
-        fValue *= fFac;
+        // Avoid 1e-5 (1.0000000000000001e-05) and such inaccurate fractional
+        // factors that later when dividing back spoil things. For negative
+        // decimals divide first with the inverse, then multiply the rounded
+        // value back.
+        fFac = getN10Exp(abs(nDecPlaces));
+        if (nDecPlaces < 0)
+            fValue /= fFac;
+        else
+            fValue *= fFac;
     }
 
     switch ( eMode )
@@ -1245,7 +1273,12 @@ double SAL_CALL rtl_math_round(double fValue, int nDecPlaces,
     }
 
     if (nDecPlaces != 0)
-        fValue /= fFac;
+    {
+        if (nDecPlaces < 0)
+            fValue *= fFac;
+        else
+            fValue /= fFac;
+    }
 
     return bSign ? -fValue : fValue;
 }