2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016,2018,2019,2020 by the GROMACS development team.
5 * Copyright (c) 2021, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 * Declares simple math functions
40 * \author Erik Lindahl <erik.lindahl@gmail.com>
42 * \ingroup module_math
44 #ifndef GMX_MATH_FUNCTIONS_H
45 #define GMX_MATH_FUNCTIONS_H
50 #include "gromacs/utility/gmxassert.h"
51 #include "gromacs/utility/real.h"
56 /*! \brief Evaluate log2(n) for integer n statically at compile time.
58 * Use as staticLog2<n>::value, where n must be a positive integer.
59 * Negative n will be reinterpreted as the corresponding unsigned integer,
60 * and you will get a compile-time error if n==0.
61 * The calculation is done by recursively dividing n by 2 (until it is 1),
62 * and incrementing the result by 1 in each step.
64 * \tparam n Value to recursively calculate log2(n) for
66 template<std::uint64_t n>
69 static const int value = StaticLog2<n / 2>::value
70 + 1; //!< Variable value used for recursive static calculation of Log2(int)
73 /*! \brief Specialization of StaticLog2<n> for n==1.
75 * This specialization provides the final value in the recursion; never
76 * call it directly, but use StaticLog2<n>::value.
81 static const int value = 0; //!< Base value for recursive static calculation of Log2(int)
84 /*! \brief Specialization of StaticLog2<n> for n==0.
86 * This specialization should never actually be used since log2(0) is
87 * negative infinity. However, since Log2() is often used to calculate the number
88 * of bits needed for a number, we might be using the value 0 with a conditional
89 * statement around the logarithm. Depending on the compiler the expansion of
90 * the template can occur before the conditional statement, so to avoid infinite
91 * recursion we need a specialization for the case n==0.
96 static const int value = -1; //!< Base value for recursive static calculation of Log2(int)
100 /*! \brief Compute floor of logarithm to base 2, 32 bit signed argument
102 * \param x 32-bit signed argument
106 * \note This version of the overloaded function will assert that x is
109 unsigned int log2I(std::int32_t x);
111 /*! \brief Compute floor of logarithm to base 2, 64 bit signed argument
113 * \param x 64-bit signed argument
117 * \note This version of the overloaded function will assert that x is
120 unsigned int log2I(std::int64_t x);
122 /*! \brief Compute floor of logarithm to base 2, 32 bit unsigned argument
124 * \param x 32-bit unsigned argument
128 * \note This version of the overloaded function uses unsigned arguments to
129 * be able to handle arguments using all 32 bits.
131 unsigned int log2I(std::uint32_t x);
133 /*! \brief Compute floor of logarithm to base 2, 64 bit unsigned argument
135 * \param x 64-bit unsigned argument
139 * \note This version of the overloaded function uses unsigned arguments to
140 * be able to handle arguments using all 64 bits.
142 unsigned int log2I(std::uint64_t x);
144 /*! \brief Find greatest common divisor of two numbers
146 * \param p First number, positive
147 * \param q Second number, positive
149 * \return Greatest common divisor of p and q
151 std::int64_t greatestCommonDivisor(std::int64_t p, std::int64_t q);
154 /*! \brief Calculate 1.0/sqrt(x) in single precision
156 * \param x Positive value to calculate inverse square root for
158 * For now this is implemented with std::sqrt(x) since gcc seems to do a
159 * decent job optimizing it. However, we might decide to use instrinsics
160 * or compiler-specific functions in the future.
162 * \return 1.0/sqrt(x)
164 static inline float invsqrt(float x)
166 return 1.0F / std::sqrt(x);
169 /*! \brief Calculate 1.0/sqrt(x) in double precision, but single range
171 * \param x Positive value to calculate inverse square root for, must be
172 * in the input domain valid for single precision.
174 * For now this is implemented with std::sqrt(x). However, we might
175 * decide to use instrinsics or compiler-specific functions in the future, and
176 * then we want to have the freedom to do the first step in single precision.
178 * \return 1.0/sqrt(x)
180 static inline double invsqrt(double x)
182 return 1.0 / std::sqrt(x);
185 /*! \brief Calculate 1.0/sqrt(x) for integer x in double precision.
187 * \param x Positive value to calculate inverse square root for.
189 * \return 1.0/sqrt(x)
191 static inline double invsqrt(int x)
193 return invsqrt(static_cast<double>(x));
196 /*! \brief Calculate inverse cube root of x in single precision
202 * This routine is typically faster than using std::pow().
204 static inline float invcbrt(float x)
206 return 1.0F / std::cbrt(x);
209 /*! \brief Calculate inverse sixth root of x in double precision
215 * This routine is typically faster than using std::pow().
217 static inline double invcbrt(double x)
219 return 1.0 / std::cbrt(x);
222 /*! \brief Calculate inverse sixth root of integer x in double precision
228 * This routine is typically faster than using std::pow().
230 static inline double invcbrt(int x)
232 return 1.0 / std::cbrt(x);
235 /*! \brief Calculate sixth root of x in single precision.
237 * \param x Argument, must be greater than or equal to zero.
241 * This routine is typically faster than using std::pow().
243 static inline float sixthroot(float x)
245 return std::sqrt(std::cbrt(x));
248 /*! \brief Calculate sixth root of x in double precision.
250 * \param x Argument, must be greater than or equal to zero.
254 * This routine is typically faster than using std::pow().
256 static inline double sixthroot(double x)
258 return std::sqrt(std::cbrt(x));
261 /*! \brief Calculate sixth root of integer x, return double.
263 * \param x Argument, must be greater than or equal to zero.
267 * This routine is typically faster than using std::pow().
269 static inline double sixthroot(int x)
271 return std::sqrt(std::cbrt(x));
274 /*! \brief Calculate inverse sixth root of x in single precision
276 * \param x Argument, must be greater than zero.
280 * This routine is typically faster than using std::pow().
282 static inline float invsixthroot(float x)
284 return invsqrt(std::cbrt(x));
287 /*! \brief Calculate inverse sixth root of x in double precision
289 * \param x Argument, must be greater than zero.
293 * This routine is typically faster than using std::pow().
295 static inline double invsixthroot(double x)
297 return invsqrt(std::cbrt(x));
300 /*! \brief Calculate inverse sixth root of integer x in double precision
302 * \param x Argument, must be greater than zero.
306 * This routine is typically faster than using std::pow().
308 static inline double invsixthroot(int x)
310 return invsqrt(std::cbrt(x));
313 /*! \brief calculate x^2
315 * \tparam T Type of argument and return value
326 /*! \brief calculate x^3
328 * \tparam T Type of argument and return value
336 return x * square(x);
339 /*! \brief calculate x^4
341 * \tparam T Type of argument and return value
349 return square(square(x));
352 /*! \brief calculate x^5
354 * \tparam T Type of argument and return value
362 return x * power4(x);
365 /*! \brief calculate x^6
367 * \tparam T Type of argument and return value
375 return square(power3(x));
378 /*! \brief calculate x^12
380 * \tparam T Type of argument and return value
388 return square(power6(x));
391 /*! \brief Maclaurin series for sinh(x)/x.
393 * Used for NH chains and MTTK pressure control.
394 * Here, we compute it to 10th order, which might be an overkill.
395 * 8th is probably enough, but it's not very much more expensive.
397 static inline real series_sinhx(real x)
404 * (1 + (x2 / 42.0_real) * (1 + (x2 / 72.0_real) * (1 + (x2 / 110.0_real))))));
407 /*! \brief Inverse error function, double precision.
409 * \param x Argument, should be in the range -1.0 < x < 1.0
411 * \return The inverse of the error function if the argument is inside the
412 * range, +/- infinity if it is exactly 1.0 or -1.0, and NaN otherwise.
414 double erfinv(double x);
416 /*! \brief Inverse error function, single precision.
418 * \param x Argument, should be in the range -1.0 < x < 1.0
420 * \return The inverse of the error function if the argument is inside the
421 * range, +/- infinity if it is exactly 1.0 or -1.0, and NaN otherwise.
423 float erfinv(float x);
425 /*! \brief Exact integer division, 32bit.
427 * \param a dividend. Function asserts that it is a multiple of divisor
430 * \return quotient of division
432 constexpr int32_t exactDiv(int32_t a, int32_t b)
434 return GMX_ASSERT(a % b == 0, "exactDiv called with non-divisible arguments"), a / b;
437 //! Exact integer division, 64bit.
438 constexpr int64_t exactDiv(int64_t a, int64_t b)
440 return GMX_ASSERT(a % b == 0, "exactDiv called with non-divisible arguments"), a / b;
443 /*! \brief Round float to int
445 * Rounding behavior is round to nearest. Rounding of halfway cases is implemention defined
446 * (either halway to even or halway away from zero).
448 /* Implementation details: It is assumed that FE_TONEAREST is default and not changed by anyone.
449 * Currently the implementation is using rint(f) because 1) on all known HW that is faster than
450 * lround and 2) some compilers (e.g. clang (#22944) and icc) don't optimize (l)lrint(f) well.
451 * GCC(>=4.7) optimizes (l)lrint(f) well but with "-fno-math-errno -funsafe-math-optimizations"
452 * rint(f) is optimized as well. This avoids using intrinsics.
453 * rint(f) followed by float/double to int/int64 conversion produces the same result as directly
454 * rounding to int/int64.
456 static inline int roundToInt(float x)
458 return static_cast<int>(rintf(x));
460 //! Round double to int
461 static inline int roundToInt(double x)
463 return static_cast<int>(rint(x));
465 //! Round float to int64_t
466 static inline int64_t roundToInt64(float x)
468 return static_cast<int>(rintf(x));
470 //! Round double to int64_t
471 static inline int64_t roundToInt64(double x)
473 return static_cast<int>(rint(x));
476 //! \brief Check whether \p v is an integer power of 2.
477 template<typename T, typename = std::enable_if_t<std::is_integral<T>::value>>
478 #if defined(__NVCC__) && !defined(__CUDACC_RELAXED_CONSTEXPR__)
479 /* In CUDA 11, a constexpr function cannot be called from a function with incompatible execution
480 * space, unless --expt-relaxed-constexpr flag is set */
483 static inline constexpr bool
484 isPowerOfTwo(const T v)
486 return (v > 0) && ((v & (v - 1)) == 0);
492 #endif // GMX_MATH_FUNCTIONS_H