Manually sort some includes in src/gromacs

[alexxy/gromacs.git] / src / gromacs / math / utilities.c

diff --git a/src/gromacs/math/utilities.c b/src/gromacs/math/utilities.c
index d5459ca6d4bb7bc417fae75c8ca66c90e59ebe56..c2df23ec2ee333d4cadb77154b0205d4c0e56439 100644 (file)
--- a/src/gromacs/math/utilities.c
+++ b/src/gromacs/math/utilities.c
@@ -34,14 +34,16 @@
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
-#include "gromacs/math/utilities.h"
+#include "gmxpre.h"
 
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "utilities.h"
 
-#include <math.h>
+#include <assert.h>
 #include <limits.h>
+#include <math.h>
+
+#include "config.h"
+
 #ifdef HAVE__FINITE
 #include <float.h>
 #endif
@@ -59,11 +61,11 @@ real cuberoot(real x)
 {
     if (x < 0)
     {
-        return (-pow(-x, 1.0/DIM));
+        return (-pow(-x, 1.0/3.0));
     }
     else
     {
-        return (pow(x, 1.0/DIM));
+        return (pow(x, 1.0/3.0));
     }
 }
 
@@ -172,6 +174,7 @@ static const double
 
 double gmx_erfd(double x)
 {
+#ifdef GMX_FLOAT_FORMAT_IEEE754
     gmx_int32_t hx, ix, i;
     double      R, S, P, Q, s, y, z, r;
 
@@ -184,9 +187,7 @@ double gmx_erfd(double x)
 
     conv.d = x;
 
-    /* In release-4-6 and later branches, only the test for
-     * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
-#if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
+#ifdef GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
     hx = conv.i[0];
 #else
     hx = conv.i[1];
@@ -262,9 +263,7 @@ double gmx_erfd(double x)
 
     conv.d = x;
 
-    /* In release-4-6 and later branches, only the test for
-     * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
-#if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
+#ifdef GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
     conv.i[1] = 0;
 #else
     conv.i[0] = 0;
@@ -281,11 +280,16 @@ double gmx_erfd(double x)
     {
         return r/x-one;
     }
+#else
+    /* No IEEE754 information. We need to trust that the OS provides erf(). */
+    return erf(x);
+#endif
 }
 
 
 double gmx_erfcd(double x)
 {
+#ifdef GMX_FLOAT_FORMAT_IEEE754
     gmx_int32_t hx, ix;
     double      R, S, P, Q, s, y, z, r;
 
@@ -298,9 +302,7 @@ double gmx_erfcd(double x)
 
     conv.d = x;
 
-    /* In release-4-6 and later branches, only the test for
-     * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
-#if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
+#ifdef GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
     hx = conv.i[0];
 #else
     hx = conv.i[1];
@@ -378,9 +380,7 @@ double gmx_erfcd(double x)
 
         conv.d = x;
 
-        /* In release-4-6 and later branches, only the test for
-         * GMX_IEEE754_BIG_ENDIAN_WORD_ORDER will be required. */
-#if defined(IEEE754_BIG_ENDIAN_WORD_ORDER) || defined(GMX_IEEE754_BIG_ENDIAN_WORD_ORDER)
+#ifdef GMX_IEEE754_BIG_ENDIAN_WORD_ORDER
         conv.i[1] = 0;
 #else
         conv.i[0] = 0;
@@ -410,6 +410,10 @@ double gmx_erfcd(double x)
             return two-tiny;
         }
     }
+#else
+    /* No IEEE754 information. We need to trust that the OS provides erfc(). */
+    return erfc(x);
+#endif
 }
 
 
@@ -722,9 +726,7 @@ float gmx_erfcf(float x)
 
 gmx_bool gmx_isfinite(real gmx_unused x)
 {
-    gmx_bool returnval = TRUE;
-    /* If no suitable function was found, assume the value is
-     * finite. */
+    gmx_bool returnval;
 
 #ifdef HAVE__FINITE
     returnval = _finite(x);
@@ -732,13 +734,16 @@ gmx_bool gmx_isfinite(real gmx_unused x)
     returnval = isfinite(x);
 #elif defined HAVE__ISFINITE
     returnval = _isfinite(x);
+#else
+    /* If no suitable function was found, assume the value is
+     * finite. */
+    returnval = TRUE;
 #endif
     return returnval;
 }
 
 gmx_bool gmx_isnan(real x)
 {
-    /* cppcheck-suppress duplicateExpression */
     return x != x;
 }
 
@@ -764,12 +769,45 @@ gmx_numzero(double a)
     return gmx_within_tol(a, 0.0, GMX_REAL_MIN/GMX_REAL_EPS);
 }
 
-real
-gmx_log2(real x)
+unsigned int
+gmx_log2i(unsigned int n)
 {
-    const real iclog2 = 1.0/log( 2.0 );
+    assert(n != 0); /* behavior differs for 0 */
+#if defined(__INTEL_COMPILER)
+    return _bit_scan_reverse(n);
+#elif defined(__GNUC__) && UINT_MAX == 4294967295U /*also for clang*/
+    return __builtin_clz(n) ^ 31U;                 /* xor gets optimized out */
+#elif defined(_MSC_VER) && _MSC_VER >= 1400
+    {
+        unsigned long i;
+        _BitScanReverse(&i, n);
+        return i;
+    }
+#elif defined(__xlC__)
+    return 31 - __cntlz4(n);
+#else
+    /* http://graphics.stanford.edu/~seander/bithacks.html#IntegerLogLookup */
+#define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
+    static const char     LogTable256[256] = {
+        -1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
+        LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
+        LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
+    };
+#undef LT
+
+    unsigned int r;     /* r will be lg(n) */
+    unsigned int t, tt; /* temporaries */
 
-    return log( x ) * iclog2;
+    if ((tt = n >> 16) != 0)
+    {
+        r = ((t = tt >> 8) != 0) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
+    }
+    else
+    {
+        r = ((t = n >> 8) != 0) ? 8 + LogTable256[t] : LogTable256[n];
+    }
+    return r;
+#endif
 }
 
 gmx_bool