* 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 "config.h"
+
+#include <assert.h>
#include <limits.h>
+#include <math.h>
+
#ifdef HAVE__FINITE
#include <float.h>
#endif
{
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));
}
}
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;
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];
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;
{
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;
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];
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;
return two-tiny;
}
}
+#else
+ /* No IEEE754 information. We need to trust that the OS provides erfc(). */
+ return erfc(x);
+#endif
}
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);
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;
}
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