tdf#113211: fix calculations with big integers

... and munbers with few fractional bits

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

config_host.mk.in

@@ -238,6 +238,7 @@ export HAMCREST_JAR=@HAMCREST_JAR@
 export HAVE_GCC_AVX=@HAVE_GCC_AVX@
 export HAVE_GCC_STACK_PROTECTOR_STRONG=@HAVE_GCC_STACK_PROTECTOR_STRONG@
 export HAVE_GCC_BUILTIN_ATOMIC=@HAVE_GCC_BUILTIN_ATOMIC@
+export HAVE_GCC_BUILTIN_FFS=@HAVE_GCC_BUILTIN_FFS@
 export HAVE_GCC_FINLINE_LIMIT=@HAVE_GCC_FINLINE_LIMIT@
 export HAVE_GCC_FNO_DEFAULT_INLINE=@HAVE_GCC_FNO_DEFAULT_INLINE@
 export HAVE_GCC_FNO_ENFORCE_EH_SPECS=@HAVE_GCC_FNO_ENFORCE_EH_SPECS@

config_host/config_global.h.in

@@ -14,6 +14,7 @@ Any change in this header will cause a rebuild of almost everything.
 
 #define HAVE_CXX14_CONSTEXPR 0
 #define HAVE_GCC_BUILTIN_ATOMIC 0
+#define HAVE_GCC_BUILTIN_FFS 0
 /* _Pragma */
 #define HAVE_GCC_PRAGMA_OPERATOR 0
 #define HAVE_GCC_DEPRECATED_MESSAGE 0

configure.ac

@@ -5795,6 +5795,15 @@ if test "$GCC" = "yes" -o "$COM_IS_CLANG" = TRUE; then
         AC_MSG_RESULT([no])
     fi
 
+    AC_MSG_CHECKING([whether $CC supports __builtin_ffs])
+    AC_LINK_IFELSE([AC_LANG_PROGRAM([[]], [[ return __builtin_ffs(1); ]])],[HAVE_GCC_BUILTIN_FFS=TRUE],[])
+    if test "$HAVE_GCC_BUILTIN_FFS" = "TRUE"; then
+        AC_MSG_RESULT([yes])
+        AC_DEFINE(HAVE_GCC_BUILTIN_FFS)
+    else
+        AC_MSG_RESULT([no])
+    fi
+
     AC_MSG_CHECKING([whether $CC supports __attribute__((deprecated(message)))])
     save_CFLAGS=$CFLAGS
     CFLAGS="$CFLAGS -Werror"
@@ -5943,6 +5952,7 @@ fi
 AC_SUBST(HAVE_GCC_AVX)
 AC_SUBST(HAVE_GCC_STACK_PROTECTOR_STRONG)
 AC_SUBST(HAVE_GCC_BUILTIN_ATOMIC)
+AC_SUBST(HAVE_GCC_BUILTIN_FFS)
 
 dnl ===================================================================
 dnl Identify the C++ library

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