* 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));
}
}
* ====================================================
*/
-#if ( (defined SIZEOF_INT && SIZEOF_INT == 4) || (SIZEOF_INT_MAX == 2147483647) )
-typedef int erf_int32_t;
-typedef unsigned int erf_u_int32_t;
-#elif (LONG_MAX == 2147483647L)
-typedef long erf_int32_t;
-typedef unsigned long erf_u_int32_t;
-#elif (SHRT_MAX == 2147483647)
-typedef short erf_int32_t;
-typedef unsigned short erf_u_int32_t;
-#else
-# error ERROR: No 32 bit wide integer type found!
-#endif
-
-
static const double
tiny = 1e-300,
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
double gmx_erfd(double x)
{
-
- erf_int32_t hx, ix, i;
+#ifdef GMX_FLOAT_FORMAT_IEEE754
+ gmx_int32_t hx, ix, i;
double R, S, P, Q, s, y, z, r;
union
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];
if (ix >= 0x7ff00000)
{
/* erf(nan)=nan */
- i = ((erf_u_int32_t)hx>>31)<<1;
+ i = ((gmx_uint32_t)hx>>31)<<1;
return (double)(1-i)+one/x; /* erf(+-inf)=+-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)
{
- erf_int32_t hx, ix;
+#ifdef GMX_FLOAT_FORMAT_IEEE754
+ gmx_int32_t hx, ix;
double R, S, P, Q, s, y, z, r;
union
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];
{
/* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
- return (double)(((erf_u_int32_t)hx>>31)<<1)+one/x;
+ return (double)(((gmx_uint32_t)hx>>31)<<1)+one/x;
}
if (ix < 0x3feb0000)
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
}
typedef union
{
float value;
- erf_u_int32_t word;
+ gmx_uint32_t word;
} ieee_float_shape_type;
#define GET_FLOAT_WORD(i, d) \
float gmx_erff(float x)
{
- erf_int32_t hx, ix, i;
+ gmx_int32_t hx, ix, i;
float R, S, P, Q, s, y, z, r;
union
if (ix >= 0x7f800000)
{
/* erf(nan)=nan */
- i = ((erf_u_int32_t)hx>>31)<<1;
+ i = ((gmx_uint32_t)hx>>31)<<1;
return (float)(1-i)+onef/x; /* erf(+-inf)=+-1 */
}
float gmx_erfcf(float x)
{
- erf_int32_t hx, ix;
+ gmx_int32_t hx, ix;
float R, S, P, Q, s, y, z, r;
union
{
/* erfc(nan)=nan */
/* erfc(+-inf)=0,2 */
- return (float)(((erf_u_int32_t)hx>>31)<<1)+onef/x;
+ return (float)(((gmx_uint32_t)hx>>31)<<1)+onef/x;
}
if (ix < 0x3f580000)
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
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff