#include <arm_neon.h>
+#include "gromacs/math/utilities.h"
+
namespace gmx
{
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
const int32x4_t exponentBias = vdupq_n_s32(127);
- int32x4_t iExponent;
+ int32x4_t iExponent = vaddq_s32(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = vmaxq_s32(iExponent, vdupq_n_s32(0));
+ }
- iExponent = vshlq_n_s32( vaddq_s32(exponent.simdInternal_, exponentBias), 23);
+ iExponent = vshlq_n_s32( iExponent, 23);
return {
vmulq_f32(value.simdInternal_, vreinterpretq_f32_s32(iExponent))
#include <arm_neon.h>
+#include "gromacs/math/utilities.h"
+
#include "impl_arm_neon_asimd_simd_float.h"
namespace gmx
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble
ldexp(SimdDouble value, SimdDInt32 exponent)
{
- const int64x2_t exponentBias = vdupq_n_s64(1023);
- int64x2_t iExponent;
+ const int32x2_t exponentBias = vdup_n_s32(1023);
+ int32x2_t iExponent = vadd_s32(exponent.simdInternal_, exponentBias);
+ int64x2_t iExponent64;
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = vmax_s32(iExponent, vdup_n_s32(0));
+ }
- iExponent = vmovl_s32(exponent.simdInternal_);
- iExponent = vshlq_n_s64(vaddq_s64(iExponent, exponentBias), 52);
+ iExponent64 = vmovl_s32(iExponent);
+ iExponent64 = vshlq_n_s64(iExponent64, 52);
return {
- vmulq_f64(value.simdInternal_, float64x2_t(iExponent))
+ vmulq_f64(value.simdInternal_, float64x2_t(iExponent64))
};
}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <qpxmath.h>
#endif
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/basedefinitions.h"
#include "impl_ibm_qpx_simd_float.h"
return value;
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble
ldexp(SimdDouble value, SimdDInt32 exponent)
{
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <qpxmath.h>
#endif
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/basedefinitions.h"
namespace gmx
return value;
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat
ldexp(SimdFloat value, SimdFInt32 exponent)
{
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <cstdint>
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/basedefinitions.h"
#include "impl_ibm_vmx_definitions.h"
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
__vector signed int iExponent;
iExponent = vec_add(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = vec_max(iExponent, vec_splat_s32(0));
+ }
+
iExponent = vec_sl(iExponent, vec_add(vec_splat_u32(15), vec_splat_u32(8)));
return {
#include "config.h"
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/basedefinitions.h"
#include "impl_ibm_vsx_definitions.h"
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble
ldexp(SimdDouble value, SimdDInt32 exponent)
{
#endif
iExponent = vec_add(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = vec_max(iExponent, vec_splat_s32(0));
+ }
+
// exponent is now present in pairs of integers; 0011.
// Elements 0/2 already correspond to the upper half of each double,
// so we only need to shift by another 52-32=20 bits.
#include "config.h"
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/basedefinitions.h"
#include "impl_ibm_vsx_definitions.h"
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
const __vector signed int exponentBias = vec_splats(127);
__vector signed int iExponent;
- iExponent = vec_sl( vec_add(exponent.simdInternal_, exponentBias), vec_splats(23U));
+ iExponent = vec_add(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = vec_max(iExponent, vec_splat_s32(0));
+ }
+
+ iExponent = vec_sl( iExponent, vec_splats(23U));
return {
vec_mul(value.simdInternal_, reinterpret_cast<__vector float>(iExponent))
* but if the pointer is not aligned the prefetch might start at the
* lower cache line boundary (meaning fewer bytes are prefetched).
*/
-static inline void
+static inline void gmx_unused
simdPrefetch(void gmx_unused * m)
{
// Do nothing for reference implementation
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <algorithm>
#include <array>
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/fatalerror.h"
#include "impl_reference_definitions.h"
}
/*! \brief Multiply a SIMD double value by the number 2 raised to an exp power.
+ *
+ * \tparam opt By default, this routine will return zero for input arguments
+ * that are so small they cannot be reproduced in the current
+ * precision. If the unsafe math optimization template parameter
+ * setting is used, these tests are skipped, and the result will
+ * be undefined (possible even NaN). This might happen below -127
+ * in single precision or -1023 in double, although some
+ * might use denormal support to extend the range.
*
* \param value Floating-point number to multiply with new exponent
* \param exponent Integer that will not overflow as 2^exponent.
* \return value*2^exponent
*/
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble gmx_simdcall
ldexp(SimdDouble value, SimdDInt32 exponent)
{
for (std::size_t i = 0; i < res.simdInternal_.size(); i++)
{
+ // std::ldexp already takes care of clamping arguments, so we do not
+ // need to do anything in the reference implementation
res.simdInternal_[i] = std::ldexp(value.simdInternal_[i], exponent.simdInternal_[i]);
}
return res;
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <algorithm>
#include <array>
+#include "gromacs/math/utilities.h"
+
#include "impl_reference_definitions.h"
namespace gmx
}
/*! \brief Multiply a SIMD float value by the number 2 raised to an exp power.
+ *
+ * \tparam opt By default, this routine will return zero for input arguments
+ * that are so small they cannot be reproduced in the current
+ * precision. If the unsafe math optimization template parameter
+ * setting is used, these tests are skipped, and the result will
+ * be undefined (possible even NaN). This might happen below -127
+ * in single precision or -1023 in double, although some
+ * might use denormal support to extend the range.
*
* \param value Floating-point number to multiply with new exponent
* \param exponent Integer that will not overflow as 2^exponent.
* \return value*2^exponent
*/
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
for (std::size_t i = 0; i < res.simdInternal_.size(); i++)
{
+ // std::ldexp already takes care of clamping arguments, so we do not
+ // need to do anything in the reference implementation
res.simdInternal_[i] = std::ldexp(value.simdInternal_[i], exponent.simdInternal_[i]);
}
return res;
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <immintrin.h>
+#include "gromacs/math/utilities.h"
#include "gromacs/simd/impl_x86_avx_256/impl_x86_avx_256_simd_double.h"
namespace gmx
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble
ldexp(SimdDouble value, SimdDInt32 exponent)
{
- const __m256i exponentBias = _mm256_set1_epi64x(1023LL);
- __m256i iExponent = _mm256_cvtepi32_epi64(exponent.simdInternal_);
+ const __m128i exponentBias = _mm_set1_epi32(1023);
+ __m128i iExponent = _mm_add_epi32(exponent.simdInternal_, exponentBias);
- iExponent = _mm256_slli_epi64(_mm256_add_epi64(iExponent, exponentBias), 52);
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm_max_epi32(iExponent, _mm_setzero_si128());
+ }
+
+ __m256i iExponent256 = _mm256_slli_epi64(_mm256_cvtepi32_epi64(iExponent), 52);
return {
- _mm256_mul_pd(value.simdInternal_, _mm256_castsi256_pd(iExponent))
+ _mm256_mul_pd(value.simdInternal_, _mm256_castsi256_pd(iExponent256))
};
}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <immintrin.h>
+#include "gromacs/math/utilities.h"
#include "gromacs/simd/impl_x86_avx_256/impl_x86_avx_256_simd_float.h"
namespace gmx
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
const __m256i exponentBias = _mm256_set1_epi32(127);
- __m256i iExponent;
+ __m256i iExponent = _mm256_add_epi32(exponent.simdInternal_, exponentBias);
- iExponent = _mm256_slli_epi32(_mm256_add_epi32(exponent.simdInternal_, exponentBias), 23);
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm256_max_epi32(iExponent, _mm256_setzero_si256());
+ }
+
+ iExponent = _mm256_slli_epi32(iExponent, 23);
return {
_mm256_mul_ps(value.simdInternal_, _mm256_castsi256_ps(iExponent))
};
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <immintrin.h>
+#include "gromacs/math/utilities.h"
+
#include "impl_x86_avx_256_simd_float.h"
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble
ldexp(SimdDouble value, SimdDInt32 exponent)
{
__m256d fExponent;
iExponentLow = _mm_add_epi32(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponentLow = _mm_max_epi32(iExponentLow, _mm_setzero_si128());
+ }
+
iExponentHigh = _mm_shuffle_epi32(iExponentLow, _MM_SHUFFLE(3, 3, 2, 2));
iExponentLow = _mm_shuffle_epi32(iExponentLow, _MM_SHUFFLE(1, 1, 0, 0));
iExponentHigh = _mm_slli_epi64(iExponentHigh, 52);
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <immintrin.h>
+#include "gromacs/math/utilities.h"
+
namespace gmx
{
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
iExponentHigh = _mm256_extractf128_si256(exponent.simdInternal_, 0x1);
iExponentLow = _mm256_castsi256_si128(exponent.simdInternal_);
- iExponentLow = _mm_slli_epi32(_mm_add_epi32(iExponentLow, exponentBias), 23);
- iExponentHigh = _mm_slli_epi32(_mm_add_epi32(iExponentHigh, exponentBias), 23);
+
+ iExponentLow = _mm_add_epi32(iExponentLow, exponentBias);
+ iExponentHigh = _mm_add_epi32(iExponentHigh, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponentLow = _mm_max_epi32(iExponentLow, _mm_setzero_si128());
+ iExponentHigh = _mm_max_epi32(iExponentHigh, _mm_setzero_si128());
+ }
+
+ iExponentLow = _mm_slli_epi32(iExponentLow, 23);
+ iExponentHigh = _mm_slli_epi32(iExponentHigh, 23);
iExponent = _mm256_castsi128_si256(iExponentLow);
iExponent = _mm256_insertf128_si256(iExponent, iExponentHigh, 0x1);
return {
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <immintrin.h>
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/basedefinitions.h"
#include "impl_x86_avx_512_general.h"
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble
ldexp(SimdDouble value, SimdDInt32 exponent)
{
- const __m512i exponentBias = _mm512_set1_epi32(1023);
- __m512i iExponent;
+ const __m256i exponentBias = _mm256_set1_epi32(1023);
+ __m256i iExponent = _mm256_add_epi32(exponent.simdInternal_, exponentBias);
+ __m512i iExponent512;
- iExponent = _mm512_permutexvar_epi32(_mm512_set_epi32(7, 7, 6, 6, 5, 5, 4, 4, 3, 3, 2, 2, 1, 1, 0, 0), _mm512_castsi256_si512(exponent.simdInternal_));
- iExponent = _mm512_mask_slli_epi32(_mm512_setzero_epi32(), avx512Int2Mask(0xAAAA), _mm512_add_epi32(iExponent, exponentBias), 20);
- return _mm512_mul_pd(_mm512_castsi512_pd(iExponent), value.simdInternal_);
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm256_max_epi32(iExponent, _mm256_setzero_si256());
+ }
+
+ iExponent512 = _mm512_permutexvar_epi32(_mm512_set_epi32(7, 7, 6, 6, 5, 5, 4, 4, 3, 3, 2, 2, 1, 1, 0, 0), _mm512_castsi256_si512(iExponent));
+ iExponent512 = _mm512_mask_slli_epi32(_mm512_setzero_epi32(), avx512Int2Mask(0xAAAA), iExponent512, 20);
+ return _mm512_mul_pd(_mm512_castsi512_pd(iExponent512), value.simdInternal_);
}
static inline double gmx_simdcall
#include <immintrin.h>
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/real.h"
#include "impl_x86_avx_512_general.h"
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
const __m512i exponentBias = _mm512_set1_epi32(127);
- __m512i iExponent;
+ __m512i iExponent = _mm512_add_epi32(exponent.simdInternal_, exponentBias);
- iExponent = _mm512_slli_epi32( _mm512_add_epi32(exponent.simdInternal_, exponentBias), 23);
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm512_max_epi32(iExponent, _mm512_setzero_epi32());
+ }
+
+ iExponent = _mm512_slli_epi32(iExponent, 23);
return {
_mm512_mul_ps(value.simdInternal_, _mm512_castsi512_ps(iExponent))
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <immintrin.h>
+#include "gromacs/math/utilities.h"
#include "gromacs/utility/basedefinitions.h"
#include "impl_x86_mic_simd_float.h"
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble
ldexp(SimdDouble value, SimdDInt32 exponent)
{
const __m512i exponentBias = _mm512_set1_epi32(1023);
- __m512i iExponent;
+ __m512i iExponent = _mm512_add_epi32(exponent.simdInternal_, exponentBias);
- iExponent = _mm512_permutevar_epi32(_mm512_set_epi32(7, 7, 6, 6, 5, 5, 4, 4, 3, 3, 2, 2, 1, 1, 0, 0), exponent.simdInternal_);
- iExponent = _mm512_mask_slli_epi32(_mm512_setzero_epi32(), _mm512_int2mask(0xAAAA), _mm512_add_epi32(iExponent, exponentBias), 20);
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm512_max_epi32(iExponent, _mm512_setzero_epi32());
+ }
+
+ iExponent = _mm512_permutevar_epi32(_mm512_set_epi32(7, 7, 6, 6, 5, 5, 4, 4, 3, 3, 2, 2, 1, 1, 0, 0), iExponent);
+ iExponent = _mm512_mask_slli_epi32(_mm512_setzero_epi32(), _mm512_int2mask(0xAAAA), iExponent, 20);
return _mm512_mul_pd(_mm512_castsi512_pd(iExponent), value.simdInternal_);
}
};
}
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
const __m512i exponentBias = _mm512_set1_epi32(127);
- __m512i iExponent;
+ __m512i iExponent = _mm512_add_epi32(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm512_max_epi32(iExponent, _mm512_setzero_epi32());
+ }
- iExponent = _mm512_slli_epi32( _mm512_add_epi32(exponent.simdInternal_, exponentBias), 23);
+ iExponent = _mm512_slli_epi32( iExponent, 23);
return {
_mm512_mul_ps(value.simdInternal_, _mm512_castsi512_ps(iExponent))
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <emmintrin.h>
+#include "gromacs/math/utilities.h"
+
#include "impl_x86_sse2_simd_float.h"
namespace gmx
};
}
+// Override for SSE4.1
+#if GMX_SIMD_X86_SSE2
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdDouble
ldexp(SimdDouble value, SimdDInt32 exponent)
{
const __m128i exponentBias = _mm_set1_epi32(1023);
- __m128i iExponent;
+ __m128i iExponent = _mm_add_epi32(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm_and_si128(iExponent, _mm_cmpgt_epi32(iExponent, _mm_setzero_si128()));
+ }
// After conversion integers will be in slot 0,1. Move them to 0,2 so
// we can do a 64-bit shift and get them to the dp exponents.
- iExponent = _mm_shuffle_epi32(exponent.simdInternal_, _MM_SHUFFLE(3, 1, 2, 0));
- iExponent = _mm_slli_epi64(_mm_add_epi32(iExponent, exponentBias), 52);
+ iExponent = _mm_shuffle_epi32(iExponent, _MM_SHUFFLE(3, 1, 2, 0));
+ iExponent = _mm_slli_epi64(iExponent, 52);
return {
_mm_mul_pd(value.simdInternal_, _mm_castsi128_pd(iExponent))
};
}
+#endif
// Override for AVX-128-FMA and higher
#if GMX_SIMD_X86_SSE2 || GMX_SIMD_X86_SSE4_1
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <emmintrin.h>
+#include "gromacs/math/utilities.h"
+
namespace gmx
{
};
}
+// Override for SSE4.1
+#if GMX_SIMD_X86_SSE2
+template <MathOptimization opt = MathOptimization::Safe>
static inline SimdFloat gmx_simdcall
ldexp(SimdFloat value, SimdFInt32 exponent)
{
const __m128i exponentBias = _mm_set1_epi32(127);
__m128i iExponent;
- iExponent = _mm_slli_epi32( _mm_add_epi32(exponent.simdInternal_, exponentBias), 23);
+ iExponent = _mm_add_epi32(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm_and_si128(iExponent, _mm_cmpgt_epi32(iExponent, _mm_setzero_si128()));
+ }
+
+ iExponent = _mm_slli_epi32( iExponent, 23);
return {
_mm_mul_ps(value.simdInternal_, _mm_castsi128_ps(iExponent))
};
}
+#endif
// Override for AVX-128-FMA and higher
#if GMX_SIMD_X86_SSE2 || GMX_SIMD_X86_SSE4_1
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
};
}
+template <MathOptimization opt = MathOptimization::Safe>
+static inline SimdDouble
+ldexp(SimdDouble value, SimdDInt32 exponent)
+{
+ const __m128i exponentBias = _mm_set1_epi32(1023);
+ __m128i iExponent = _mm_add_epi32(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm_max_epi32(iExponent, _mm_setzero_si128());
+ }
+
+ // After conversion integers will be in slot 0,1. Move them to 0,2 so
+ // we can do a 64-bit shift and get them to the dp exponents.
+ iExponent = _mm_shuffle_epi32(iExponent, _MM_SHUFFLE(3, 1, 2, 0));
+ iExponent = _mm_slli_epi64(iExponent, 52);
+
+ return {
+ _mm_mul_pd(value.simdInternal_, _mm_castsi128_pd(iExponent))
+ };
+}
+
} // namespace gmx
#endif // GMX_SIMD_IMPL_X86_SSE4_1_SIMD_DOUBLE_H
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
};
}
+template <MathOptimization opt = MathOptimization::Safe>
+static inline SimdFloat gmx_simdcall
+ldexp(SimdFloat value, SimdFInt32 exponent)
+{
+ const __m128i exponentBias = _mm_set1_epi32(127);
+ __m128i iExponent;
+
+ iExponent = _mm_add_epi32(exponent.simdInternal_, exponentBias);
+
+ if (opt == MathOptimization::Safe)
+ {
+ // Make sure biased argument is not negative
+ iExponent = _mm_max_epi32(iExponent, _mm_setzero_si128());
+ }
+
+ iExponent = _mm_slli_epi32( iExponent, 23);
+
+ return {
+ _mm_mul_ps(value.simdInternal_, _mm_castsi128_ps(iExponent))
+ };
+}
} // namespace gmx
#include <cmath>
+#include <limits>
+
#include "gromacs/math/utilities.h"
#include "gromacs/simd/simd.h"
#include "gromacs/utility/basedefinitions.h"
static inline SimdFloat gmx_simdcall
exp2(SimdFloat x)
{
- // Lower bound: Clamp args that would lead to an IEEE fp exponent below -127.
- // Note that we use a slightly wider range inside this *implementation* compared
- // to that guaranteed by the API documentation in the comment above (which has
- // to cover all architectures).
- const SimdFloat smallArgLimit(-127.0f);
const SimdFloat CC6(0.0001534581200287996416911311f);
const SimdFloat CC5(0.001339993121934088894618990f);
const SimdFloat CC4(0.009618488957115180159497841f);
// subtraction to get the fractional part loses accuracy, and then we
// can end up with overflows in the polynomial.
//
- // For now, we handle this by clamping smaller arguments to -127. At this
- // point we will already return zero, so we don't need to do anything
- // extra for the exponent.
-
+ // For now, we handle this by forwarding the math optimization setting to
+ // ldexp, where the routine will return zero for very small arguments.
+ //
+ // However, before doing that we need to make sure we do not call cvtR2I
+ // with an argument that is so negative it cannot be converted to an integer.
if (opt == MathOptimization::Safe)
{
- x = max(x, smallArgLimit);
+ x = max(x, SimdFloat(std::numeric_limits<std::int32_t>::lowest()));
}
- fexppart = ldexp(one, cvtR2I(x));
+ fexppart = ldexp<opt>(one, cvtR2I(x));
intpart = round(x);
x = x - intpart;
exp(SimdFloat x)
{
const SimdFloat argscale(1.44269504088896341f);
- // Lower bound: Clamp args that would lead to an IEEE fp exponent below -127.
- // Note that we use a slightly wider range inside this *implementation* compared
- // to that guaranteed by the API documentation in the comment above (which has
- // to cover all architectures).
- const SimdFloat smallArgLimit(-88.0296919311f);
const SimdFloat invargscale0(-0.693145751953125f);
const SimdFloat invargscale1(-1.428606765330187045e-06f);
const SimdFloat CC4(0.00136324646882712841033936f);
// Large negative values are valid arguments to exp2(), so there are two
// things we need to account for:
- // 1. When the IEEE exponents reaches -127, the (biased) exponent field will be
+ // 1. When the exponents reaches -127, the (biased) exponent field will be
// zero and we can no longer multiply with it. There are special IEEE
// formats to handle this range, but for now we have to accept that
// we cannot handle those arguments. If input value becomes even more
// subtraction to get the fractional part loses accuracy, and then we
// can end up with overflows in the polynomial.
//
- // For now, we handle this by clamping smaller arguments to -127*ln(2). At this
- // point we will already return zero, so we don't need to do anything
- // extra for the exponent.
+ // For now, we handle this by forwarding the math optimization setting to
+ // ldexp, where the routine will return zero for very small arguments.
+ //
+ // However, before doing that we need to make sure we do not call cvtR2I
+ // with an argument that is so negative it cannot be converted to an integer
+ // after being multiplied by argscale.
if (opt == MathOptimization::Safe)
{
- x = max(x, smallArgLimit);
+ x = max(x, SimdFloat(std::numeric_limits<std::int32_t>::lowest())/argscale);
}
y = x * argscale;
- fexppart = ldexp(one, cvtR2I(y));
+
+
+ fexppart = ldexp<opt>(one, cvtR2I(y));
intpart = round(y);
// Extended precision arithmetics
static inline SimdDouble gmx_simdcall
exp2(SimdDouble x)
{
- // Lower bound: Clamp args that would lead to an IEEE fp exponent below -1023.
- const SimdDouble smallArgLimit(-1023.0);
const SimdDouble CE11(4.435280790452730022081181e-10);
const SimdDouble CE10(7.074105630863314448024247e-09);
const SimdDouble CE9(1.017819803432096698472621e-07);
// subtraction to get the fractional part loses accuracy, and then we
// can end up with overflows in the polynomial.
//
- // For now, we handle this by clamping smaller arguments to -1023. At this
- // point we will already return zero, so we don't need to do anything
- // extra for the exponent.
-
+ // For now, we handle this by forwarding the math optimization setting to
+ // ldexp, where the routine will return zero for very small arguments.
+ //
+ // However, before doing that we need to make sure we do not call cvtR2I
+ // with an argument that is so negative it cannot be converted to an integer.
if (opt == MathOptimization::Safe)
{
- x = max(x, smallArgLimit);
+ x = max(x, SimdDouble(std::numeric_limits<std::int32_t>::lowest()));
}
- fexppart = ldexp(one, cvtR2I(x));
+ fexppart = ldexp<opt>(one, cvtR2I(x));
intpart = round(x);
x = x - intpart;
exp(SimdDouble x)
{
const SimdDouble argscale(1.44269504088896340735992468100);
- // Lower bound: Clamp args that would lead to an IEEE fp exponent below -1023.
- const SimdDouble smallArgLimit(-709.0895657128);
const SimdDouble invargscale0(-0.69314718055966295651160180568695068359375);
const SimdDouble invargscale1(-2.8235290563031577122588448175013436025525412068e-13);
const SimdDouble CE12(2.078375306791423699350304e-09);
// subtraction to get the fractional part loses accuracy, and then we
// can end up with overflows in the polynomial.
//
- // For now, we handle this by clamping smaller arguments to -1023*ln(2). At this
- // point we will already return zero, so we don't need to do anything
- // extra for the exponent.
+ // For now, we handle this by forwarding the math optimization setting to
+ // ldexp, where the routine will return zero for very small arguments.
+ //
+ // However, before doing that we need to make sure we do not call cvtR2I
+ // with an argument that is so negative it cannot be converted to an integer
+ // after being multiplied by argscale.
if (opt == MathOptimization::Safe)
{
- x = max(x, smallArgLimit);
+ x = max(x, SimdDouble(std::numeric_limits<std::int32_t>::lowest())/argscale);
}
y = x * argscale;
- fexppart = ldexp(one, cvtR2I(y));
+
+ fexppart = ldexp<opt>(one, cvtR2I(y));
intpart = round(y);
// Extended precision arithmetics
static inline SimdDouble gmx_simdcall
exp2SingleAccuracy(SimdDouble x)
{
- // Lower bound: Clamp args that would lead to an IEEE fp exponent below -1023.
- const SimdDouble smallArgLimit(-1023.0);
const SimdDouble CC6(0.0001534581200287996416911311);
const SimdDouble CC5(0.001339993121934088894618990);
const SimdDouble CC4(0.009618488957115180159497841);
// subtraction to get the fractional part loses accuracy, and then we
// can end up with overflows in the polynomial.
//
- // For now, we handle this by clamping smaller arguments to -1023. At this
- // point we will already return zero, so we don't need to do anything
- // extra for the exponent.
-
+ // For now, we handle this by forwarding the math optimization setting to
+ // ldexp, where the routine will return zero for very small arguments.
+ //
+ // However, before doing that we need to make sure we do not call cvtR2I
+ // with an argument that is so negative it cannot be converted to an integer.
if (opt == MathOptimization::Safe)
{
- x = max(x, smallArgLimit);
+ x = max(x, SimdDouble(std::numeric_limits<std::int32_t>::lowest()));
}
ix = cvtR2I(x);
p = fma(p, x, CC2);
p = fma(p, x, CC1);
p = fma(p, x, one);
- x = ldexp(p, ix);
+ x = ldexp<opt>(p, ix);
return x;
}
// Large negative values are valid arguments to exp2(), so there are two
// things we need to account for:
- // 1. When the exponents reaches -127, the (biased) exponent field will be
+ // 1. When the exponents reaches -1023, the (biased) exponent field will be
// zero and we can no longer multiply with it. There are special IEEE
// formats to handle this range, but for now we have to accept that
// we cannot handle those arguments. If input value becomes even more
// subtraction to get the fractional part loses accuracy, and then we
// can end up with overflows in the polynomial.
//
- // For now, we handle this by clamping smaller arguments to -127. At this
- // point we will already return zero, so we don't need to do anything
- // extra for the exponent.
+ // For now, we handle this by forwarding the math optimization setting to
+ // ldexp, where the routine will return zero for very small arguments.
+ //
+ // However, before doing that we need to make sure we do not call cvtR2I
+ // with an argument that is so negative it cannot be converted to an integer
+ // after being multiplied by argscale.
if (opt == MathOptimization::Safe)
{
- x = max(x, smallArgLimit);
+ x = max(x, SimdDouble(std::numeric_limits<std::int32_t>::lowest())/argscale);
}
y = x * argscale;
+
iy = cvtR2I(y);
intpart = round(y); // use same rounding algorithm here
p = fma(p, x, CC0);
p = fma(x*x, p, x);
p = p + one;
- x = ldexp(p, iy);
+ x = ldexp<opt>(p, iy);
return x;
}
SimdReal one = setSimdRealFrom1R(1.0);
GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(pow(2.0, 60.0), pow(2.0, -41.0), pow(2.0, 54.0)),
- ldexp(one, setSimdIntFrom3I(60, -41, 54)));
+ ldexp<MathOptimization::Unsafe>(one, setSimdIntFrom3I(60, -41, 54)));
#if GMX_SIMD_HAVE_DOUBLE && GMX_DOUBLE
GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(pow(2.0, 587.0), pow(2.0, -462.0), pow(2.0, 672.0)),
- ldexp(one, setSimdIntFrom3I(587, -462, 672)));
+ ldexp<MathOptimization::Unsafe>(one, setSimdIntFrom3I(587, -462, 672)));
#endif
/* Rounding mode in conversions must be consistent with simdRound() for SetExponent() to work */
- GMX_EXPECT_SIMD_REAL_EQ(ldexp(one, cvtR2I(round(x0))), ldexp(one, cvtR2I(x0)));
- GMX_EXPECT_SIMD_REAL_EQ(ldexp(one, cvtR2I(round(x1))), ldexp(one, cvtR2I(x1)));
+ GMX_EXPECT_SIMD_REAL_EQ(ldexp<MathOptimization::Unsafe>(one, cvtR2I(round(x0))), ldexp<MathOptimization::Unsafe>(one, cvtR2I(x0)));
+ GMX_EXPECT_SIMD_REAL_EQ(ldexp<MathOptimization::Unsafe>(one, cvtR2I(round(x1))), ldexp<MathOptimization::Unsafe>(one, cvtR2I(x1)));
+
+ // The default safe version must be able to handle very negative arguments too
+ GMX_EXPECT_SIMD_REAL_EQ(setZero(), ldexp(one, setSimdIntFrom3I(-2000, -1000000, -1000000000)));
}
/*
//The test values cannot be computed without FMA for the case that SIMD has no FMA. Even if no explicit FMA were
//used, the compiler could choose to use FMA. This causes up to 6upl error because of the product is up to 6 times
//larger than the final result after the difference.
- setUlpTol(GMX_SIMD_HAVE_FMA ? 0 : 6);
+ setUlpTol(GMX_SIMD_HAVE_FMA ? ulpTol_ : 6);
GMX_EXPECT_SIMD_REAL_NEAR(refcX, cX);
GMX_EXPECT_SIMD_REAL_NEAR(refcY, cY);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff