tdf#113211: fix calculations with big integers

... and munbers with few fractional bits

Change-Id: I86c3e8021e803fed498fae768ded9c9e5337c8bd
Reviewed-on: https://gerrit.libreoffice.org/43477
Tested-by: Jenkins <ci@libreoffice.org>
Reviewed-by: Eike Rathke <erack@redhat.com>

1 files changed, 52 insertions, 1 deletions

diff --git a/sal/rtl/math.cxx b/sal/rtl/math.cxx
index 338f40d1469f..32121b34b2f1 100644
--- a/sal/rtl/math.cxx
+++ b/sal/rtl/math.cxx
@@ -37,6 +37,10 @@
 #include <math.h>
 #include <stdlib.h>
 
+#if !HAVE_GCC_BUILTIN_FFS && !defined _WIN32
+    #include <strings.h>
+#endif
+
 static int const n10Count = 16;
 static double const n10s[2][n10Count] = {
     { 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8,
@@ -169,6 +173,47 @@ bool isRepresentableInteger(double fAbsValue)
     return false;
 }
 
+// Returns 1-based index of least significant bit in a number, or zero if number is zero
+int findFirstSetBit(unsigned n)
+{
+#if HAVE_GCC_BUILTIN_FFS
+    return __builtin_ffs(n);
+#elif defined _WIN32
+    unsigned long pos;
+    unsigned char bNonZero = _BitScanForward(&pos, n);
+    return (bNonZero == 0) ? 0 : pos + 1;
+#else
+    return ffs(n);
+#endif
+}
+
+/** Returns number of binary bits for fractional part of the number
+    Expects a proper non-negative double value, not +-INF, not NAN
+ */
+int getBitsInFracPart(double fAbsValue)
+{
+    assert(rtl::math::isFinite(fAbsValue) && fAbsValue >= 0.0);
+    if (fAbsValue == 0.0)
+        return 0;
+    auto pValParts = reinterpret_cast< const sal_math_Double * >(&fAbsValue);
+    int nExponent = pValParts->inf_parts.exponent - 1023;
+    if (nExponent >= 52)
+        return 0; // All bits in fraction are in integer part of the number
+    int nLeastSignificant = findFirstSetBit(pValParts->inf_parts.fraction_lo);
+    if (nLeastSignificant == 0)
+    {
+        nLeastSignificant = findFirstSetBit(pValParts->inf_parts.fraction_hi);
+        if (nLeastSignificant == 0)
+            nLeastSignificant = 53; // the implied leading 1 is the least significant
+        else
+            nLeastSignificant += 32;
+    }
+    int nFracSignificant = 53 - nLeastSignificant;
+    int nBitsInFracPart = nFracSignificant - nExponent;
+
+    return nBitsInFracPart > 0 ? nBitsInFracPart : 0;
+}
+
 template< typename T >
 inline void doubleToString(typename T::String ** pResult,
                            sal_Int32 * pResultCapacity, sal_Int32 nResultOffset,
@@ -1136,7 +1181,8 @@ double SAL_CALL rtl_math_pow10Exp(double fValue, int nExp) SAL_THROW_EXTERN_C()
 
 double SAL_CALL rtl_math_approxValue( double fValue ) SAL_THROW_EXTERN_C()
 {
-    if (fValue == 0.0 || fValue == HUGE_VAL || !::rtl::math::isFinite( fValue))
+    const double fBigInt = 2199023255552.0; // 2^41 -> only 11 bits left for fractional part, fine as decimal
+    if (fValue == 0.0 || fValue == HUGE_VAL || !::rtl::math::isFinite( fValue) || fValue > fBigInt)
     {
         // We don't handle these conditions.  Bail out.
         return fValue;
@@ -1148,6 +1194,11 @@ double SAL_CALL rtl_math_approxValue( double fValue ) SAL_THROW_EXTERN_C()
     if (bSign)
         fValue = -fValue;
 
+    // If the value is either integer representable as double,
+    // or only has small number of bits in fraction part, then we need not do any approximation
+    if (isRepresentableInteger(fValue) || getBitsInFracPart(fValue) <= 11)
+        return fOrigValue;
+
     int nExp = static_cast< int >(floor(log10(fValue)));
     nExp = 14 - nExp;
     double fExpValue = getN10Exp(nExp);