2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016,2018, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Declares simple math functions
39 * \author Erik Lindahl <erik.lindahl@gmail.com>
41 * \ingroup module_math
43 #ifndef GMX_MATH_FUNCTIONS_H
44 #define GMX_MATH_FUNCTIONS_H
49 #include "gromacs/utility/gmxassert.h"
50 #include "gromacs/utility/real.h"
55 /*! \brief Evaluate log2(n) for integer n statically at compile time.
57 * Use as staticLog2<n>::value, where n must be a positive integer.
58 * Negative n will be reinterpreted as the corresponding unsigned integer,
59 * and you will get a compile-time error if n==0.
60 * The calculation is done by recursively dividing n by 2 (until it is 1),
61 * and incrementing the result by 1 in each step.
63 * \tparam n Value to recursively calculate log2(n) for
65 template<std::uint64_t n
>
68 static const int value
= StaticLog2
<n
/2>::value
+1; //!< Variable value used for recursive static calculation of Log2(int)
71 /*! \brief Specialization of StaticLog2<n> for n==1.
73 * This specialization provides the final value in the recursion; never
74 * call it directly, but use StaticLog2<n>::value.
79 static const int value
= 0; //!< Base value for recursive static calculation of Log2(int)
82 /*! \brief Specialization of StaticLog2<n> for n==0.
84 * This specialization should never actually be used since log2(0) is
85 * negative infinity. However, since Log2() is often used to calculate the number
86 * of bits needed for a number, we might be using the value 0 with a conditional
87 * statement around the logarithm. Depending on the compiler the expansion of
88 * the template can occur before the conditional statement, so to avoid infinite
89 * recursion we need a specialization for the case n==0.
94 static const int value
= -1; //!< Base value for recursive static calculation of Log2(int)
98 /*! \brief Compute floor of logarithm to base 2, 32 bit signed argument
100 * \param x 32-bit signed argument
104 * \note This version of the overloaded function will assert that x is
108 log2I(std::int32_t x
);
110 /*! \brief Compute floor of logarithm to base 2, 64 bit signed argument
112 * \param x 64-bit signed argument
116 * \note This version of the overloaded function will assert that x is
120 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.
132 log2I(std::uint32_t x
);
134 /*! \brief Compute floor of logarithm to base 2, 64 bit unsigned argument
136 * \param x 64-bit unsigned argument
140 * \note This version of the overloaded function uses unsigned arguments to
141 * be able to handle arguments using all 64 bits.
144 log2I(std::uint64_t x
);
146 /*! \brief Find greatest common divisor of two numbers
148 * \param p First number, positive
149 * \param q Second number, positive
151 * \return Greatest common divisor of p and q
154 greatestCommonDivisor(std::int64_t p
, std::int64_t q
);
157 /*! \brief Calculate 1.0/sqrt(x) in single precision
159 * \param x Positive value to calculate inverse square root for
161 * For now this is implemented with std::sqrt(x) since gcc seems to do a
162 * decent job optimizing it. However, we might decide to use instrinsics
163 * or compiler-specific functions in the future.
165 * \return 1.0/sqrt(x)
170 return 1.0f
/std::sqrt(x
);
173 /*! \brief Calculate 1.0/sqrt(x) in double precision, but single range
175 * \param x Positive value to calculate inverse square root for, must be
176 * in the input domain valid for single precision.
178 * For now this is implemented with std::sqrt(x). However, we might
179 * decide to use instrinsics or compiler-specific functions in the future, and
180 * then we want to have the freedom to do the first step in single precision.
182 * \return 1.0/sqrt(x)
187 return 1.0/std::sqrt(x
);
190 /*! \brief Calculate 1.0/sqrt(x) for integer x in double precision.
192 * \param x Positive value to calculate inverse square root for.
194 * \return 1.0/sqrt(x)
199 return invsqrt(static_cast<double>(x
));
202 /*! \brief Calculate inverse cube root of x in single precision
208 * This routine is typically faster than using std::pow().
213 return 1.0f
/std::cbrt(x
);
216 /*! \brief Calculate inverse sixth root of x in double precision
222 * This routine is typically faster than using std::pow().
227 return 1.0/std::cbrt(x
);
230 /*! \brief Calculate inverse sixth root of integer x in double precision
236 * This routine is typically faster than using std::pow().
241 return 1.0/std::cbrt(x
);
244 /*! \brief Calculate sixth root of x in single precision.
246 * \param x Argument, must be greater than or equal to zero.
250 * This routine is typically faster than using std::pow().
255 return std::sqrt(std::cbrt(x
));
258 /*! \brief Calculate sixth root of x in double precision.
260 * \param x Argument, must be greater than or equal to zero.
264 * This routine is typically faster than using std::pow().
269 return std::sqrt(std::cbrt(x
));
272 /*! \brief Calculate sixth root of integer x, return double.
274 * \param x Argument, must be greater than or equal to zero.
278 * This routine is typically faster than using std::pow().
283 return std::sqrt(std::cbrt(x
));
286 /*! \brief Calculate inverse sixth root of x in single precision
288 * \param x Argument, must be greater than zero.
292 * This routine is typically faster than using std::pow().
295 invsixthroot(float x
)
297 return invsqrt(std::cbrt(x
));
300 /*! \brief Calculate inverse sixth root of x in double precision
302 * \param x Argument, must be greater than zero.
306 * This routine is typically faster than using std::pow().
309 invsixthroot(double x
)
311 return invsqrt(std::cbrt(x
));
314 /*! \brief Calculate inverse sixth root of integer x in double precision
316 * \param x Argument, must be greater than zero.
320 * This routine is typically faster than using std::pow().
325 return invsqrt(std::cbrt(x
));
328 /*! \brief calculate x^2
330 * \tparam T Type of argument and return value
335 template <typename T
>
342 /*! \brief calculate x^3
344 * \tparam T Type of argument and return value
349 template <typename T
>
356 /*! \brief calculate x^4
358 * \tparam T Type of argument and return value
363 template <typename T
>
367 return square(square(x
));
370 /*! \brief calculate x^5
372 * \tparam T Type of argument and return value
377 template <typename T
>
384 /*! \brief calculate x^6
386 * \tparam T Type of argument and return value
391 template <typename T
>
395 return square(power3(x
));
398 /*! \brief calculate x^12
400 * \tparam T Type of argument and return value
405 template <typename T
>
409 return square(power6(x
));
412 /*! \brief Maclaurin series for sinh(x)/x.
414 * Used for NH chains and MTTK pressure control.
415 * Here, we compute it to 10th order, which might be an overkill.
416 * 8th is probably enough, but it's not very much more expensive.
418 static inline real
series_sinhx(real x
)
421 return (1 + (x2
/6.0)*(1 + (x2
/20.0)*(1 + (x2
/42.0)*(1 + (x2
/72.0)*(1 + (x2
/110.0))))));
424 /*! \brief Inverse error function, double precision.
426 * \param x Argument, should be in the range -1.0 < x < 1.0
428 * \return The inverse of the error function if the argument is inside the
429 * range, +/- infinity if it is exactly 1.0 or -1.0, and NaN otherwise.
434 /*! \brief Inverse error function, single precision.
436 * \param x Argument, should be in the range -1.0 < x < 1.0
438 * \return The inverse of the error function if the argument is inside the
439 * range, +/- infinity if it is exactly 1.0 or -1.0, and NaN otherwise.
444 /*! \brief Exact integer division, 32bit.
446 * \param a dividend. Function asserts that it is a multiple of divisor
449 * \return quotient of division
451 constexpr int32_t exactDiv(int32_t a
, int32_t b
)
453 return GMX_ASSERT(a
%b
== 0, "exactDiv called with non-divisible arguments"), a
/b
;
456 //! Exact integer division, 64bit.
457 constexpr int64_t exactDiv(int64_t a
, int64_t b
)
459 return GMX_ASSERT(a
%b
== 0, "exactDiv called with non-divisible arguments"), a
/b
;
462 /*! \brief Round float to int
464 * Rounding behavior is round to nearest. Rounding of halfway cases is implemention defined
465 * (either halway to even or halway away from zero).
467 /* Implementation details: It is assumed that FE_TONEAREST is default and not changed by anyone.
468 * Currently the implementation is using rint(f) because 1) on all known HW that is faster than
469 * lround and 2) some compilers (e.g. clang (#22944) and icc) don't optimize (l)lrint(f) well.
470 * GCC(>=4.7) optimizes (l)lrint(f) well but with "-fno-math-errno -funsafe-math-optimizations"
471 * rint(f) is optimized as well. This avoids using intrinsics.
472 * rint(f) followed by float/double to int/int64 conversion produces the same result as directly
473 * rounding to int/int64.
475 static inline int roundToInt(float x
)
477 return static_cast<int>(rintf(x
));
479 //! Round double to int
480 static inline int roundToInt(double x
)
482 return static_cast<int>(rint(x
));
484 //! Round float to int64_t
485 static inline int64_t roundToInt64(float x
)
487 return static_cast<int>(rintf(x
));
489 //! Round double to int64_t
490 static inline int64_t roundToInt64(double x
)
492 return static_cast<int>(rint(x
));
498 #endif // GMX_MATH_FUNCTIONS_H