api/legacy/include/gromacs/math/functions.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   9  *
  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.
  14  *
  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.
  19  *
  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.
  24  *
  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.
  32  *
  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.
  35  */
  36 /*! \file
  37  * \brief
  38  * Declares simple math functions
  39  *
  40  * \author Erik Lindahl <erik.lindahl@gmail.com>
  41  * \inpublicapi
  42  * \ingroup module_math
  43  */
  44 #ifndef GMX_MATH_FUNCTIONS_H
  45 #define GMX_MATH_FUNCTIONS_H
  46
  47 #include <cmath>
  48 #include <cstdint>
  49
  50 #include "gromacs/utility/gmxassert.h"
  51 #include "gromacs/utility/real.h"
  52
  53 namespace gmx
  54 {
  55
  56 /*! \brief Evaluate log2(n) for integer n statically at compile time.
  57  *
  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.
  63  *
  64  * \tparam n Value to recursively calculate log2(n) for
  65  */
  66 template<std::uint64_t n>
  67 struct StaticLog2
  68 {
  69     static const int value = StaticLog2<n / 2>::value
  70                              + 1; //!< Variable value used for recursive static calculation of Log2(int)
  71 };
  72
  73 /*! \brief Specialization of StaticLog2<n> for n==1.
  74  *
  75  *  This specialization provides the final value in the recursion; never
  76  *  call it directly, but use StaticLog2<n>::value.
  77  */
  78 template<>
  79 struct StaticLog2<1>
  80 {
  81     static const int value = 0; //!< Base value for recursive static calculation of Log2(int)
  82 };
  83
  84 /*! \brief Specialization of StaticLog2<n> for n==0.
  85  *
  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.
  92  */
  93 template<>
  94 struct StaticLog2<0>
  95 {
  96     static const int value = -1; //!< Base value for recursive static calculation of Log2(int)
  97 };
  98
  99
 100 /*! \brief Compute floor of logarithm to base 2, 32 bit signed argument
 101  *
 102  *  \param x 32-bit signed argument
 103  *
 104  *  \return log2(x)
 105  *
 106  *  \note This version of the overloaded function will assert that x is
 107  *        not negative.
 108  */
 109 unsigned int log2I(std::int32_t x);
 110
 111 /*! \brief Compute floor of logarithm to base 2, 64 bit signed argument
 112  *
 113  *  \param x 64-bit signed argument
 114  *
 115  *  \return log2(x)
 116  *
 117  *  \note This version of the overloaded function will assert that x is
 118  *        not negative.
 119  */
 120 unsigned int log2I(std::int64_t x);
 121
 122 /*! \brief Compute floor of logarithm to base 2, 32 bit unsigned argument
 123  *
 124  *  \param x 32-bit unsigned argument
 125  *
 126  *  \return log2(x)
 127  *
 128  *  \note This version of the overloaded function uses unsigned arguments to
 129  *        be able to handle arguments using all 32 bits.
 130  */
 131 unsigned int log2I(std::uint32_t x);
 132
 133 /*! \brief Compute floor of logarithm to base 2, 64 bit unsigned argument
 134  *
 135  *  \param x 64-bit unsigned argument
 136  *
 137  *  \return log2(x)
 138  *
 139  *  \note This version of the overloaded function uses unsigned arguments to
 140  *        be able to handle arguments using all 64 bits.
 141  */
 142 unsigned int log2I(std::uint64_t x);
 143
 144 /*! \brief Find greatest common divisor of two numbers
 145  *
 146  *  \param p First number, positive
 147  *  \param q Second number, positive
 148  *
 149  * \return Greatest common divisor of p and q
 150  */
 151 std::int64_t greatestCommonDivisor(std::int64_t p, std::int64_t q);
 152
 153
 154 /*! \brief Calculate 1.0/sqrt(x) in single precision
 155  *
 156  * \param  x  Positive value to calculate inverse square root for
 157  *
 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.
 161  *
 162  * \return 1.0/sqrt(x)
 163  */
 164 static inline float invsqrt(float x)
 165 {
 166     return 1.0F / std::sqrt(x);
 167 }
 168
 169 /*! \brief Calculate 1.0/sqrt(x) in double precision, but single range
 170  *
 171  * \param  x  Positive value to calculate inverse square root for, must be
 172  *            in the input domain valid for single precision.
 173  *
 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.
 177  *
 178  * \return 1.0/sqrt(x)
 179  */
 180 static inline double invsqrt(double x)
 181 {
 182     return 1.0 / std::sqrt(x);
 183 }
 184
 185 /*! \brief Calculate 1.0/sqrt(x) for integer x in double precision.
 186  *
 187  * \param  x  Positive value to calculate inverse square root for.
 188  *
 189  * \return 1.0/sqrt(x)
 190  */
 191 static inline double invsqrt(int x)
 192 {
 193     return invsqrt(static_cast<double>(x));
 194 }
 195
 196 /*! \brief Calculate inverse cube root of x in single precision
 197  *
 198  *  \param  x  Argument
 199  *
 200  *  \return x^(-1/3)
 201  *
 202  *  This routine is typically faster than using std::pow().
 203  */
 204 static inline float invcbrt(float x)
 205 {
 206     return 1.0F / std::cbrt(x);
 207 }
 208
 209 /*! \brief Calculate inverse sixth root of x in double precision
 210  *
 211  *  \param  x  Argument
 212  *
 213  *  \return x^(-1/3)
 214  *
 215  *  This routine is typically faster than using std::pow().
 216  */
 217 static inline double invcbrt(double x)
 218 {
 219     return 1.0 / std::cbrt(x);
 220 }
 221
 222 /*! \brief Calculate inverse sixth root of integer x in double precision
 223  *
 224  *  \param  x  Argument
 225  *
 226  *  \return x^(-1/3)
 227  *
 228  *  This routine is typically faster than using std::pow().
 229  */
 230 static inline double invcbrt(int x)
 231 {
 232     return 1.0 / std::cbrt(x);
 233 }
 234
 235 /*! \brief Calculate sixth root of x in single precision.
 236  *
 237  *  \param  x  Argument, must be greater than or equal to zero.
 238  *
 239  *  \return x^(1/6)
 240  *
 241  *  This routine is typically faster than using std::pow().
 242  */
 243 static inline float sixthroot(float x)
 244 {
 245     return std::sqrt(std::cbrt(x));
 246 }
 247
 248 /*! \brief Calculate sixth root of x in double precision.
 249  *
 250  *  \param  x  Argument, must be greater than or equal to zero.
 251  *
 252  *  \return x^(1/6)
 253  *
 254  *  This routine is typically faster than using std::pow().
 255  */
 256 static inline double sixthroot(double x)
 257 {
 258     return std::sqrt(std::cbrt(x));
 259 }
 260
 261 /*! \brief Calculate sixth root of integer x, return double.
 262  *
 263  *  \param  x  Argument, must be greater than or equal to zero.
 264  *
 265  *  \return x^(1/6)
 266  *
 267  *  This routine is typically faster than using std::pow().
 268  */
 269 static inline double sixthroot(int x)
 270 {
 271     return std::sqrt(std::cbrt(x));
 272 }
 273
 274 /*! \brief Calculate inverse sixth root of x in single precision
 275  *
 276  *  \param  x  Argument, must be greater than zero.
 277  *
 278  *  \return x^(-1/6)
 279  *
 280  *  This routine is typically faster than using std::pow().
 281  */
 282 static inline float invsixthroot(float x)
 283 {
 284     return invsqrt(std::cbrt(x));
 285 }
 286
 287 /*! \brief Calculate inverse sixth root of x in double precision
 288  *
 289  *  \param  x  Argument, must be greater than zero.
 290  *
 291  *  \return x^(-1/6)
 292  *
 293  *  This routine is typically faster than using std::pow().
 294  */
 295 static inline double invsixthroot(double x)
 296 {
 297     return invsqrt(std::cbrt(x));
 298 }
 299
 300 /*! \brief Calculate inverse sixth root of integer x in double precision
 301  *
 302  *  \param  x  Argument, must be greater than zero.
 303  *
 304  *  \return x^(-1/6)
 305  *
 306  *  This routine is typically faster than using std::pow().
 307  */
 308 static inline double invsixthroot(int x)
 309 {
 310     return invsqrt(std::cbrt(x));
 311 }
 312
 313 /*! \brief calculate x^2
 314  *
 315  *  \tparam T  Type of argument and return value
 316  *  \param  x  argument
 317  *
 318  *  \return x^2
 319  */
 320 template<typename T>
 321 T square(T x)
 322 {
 323     return x * x;
 324 }
 325
 326 /*! \brief calculate x^3
 327  *
 328  *  \tparam T  Type of argument and return value
 329  *  \param  x  argument
 330  *
 331  *  \return x^3
 332  */
 333 template<typename T>
 334 T power3(T x)
 335 {
 336     return x * square(x);
 337 }
 338
 339 /*! \brief calculate x^4
 340  *
 341  *  \tparam T  Type of argument and return value
 342  *  \param  x  argument
 343  *
 344  *  \return x^4
 345  */
 346 template<typename T>
 347 T power4(T x)
 348 {
 349     return square(square(x));
 350 }
 351
 352 /*! \brief calculate x^5
 353  *
 354  *  \tparam T  Type of argument and return value
 355  *  \param  x  argument
 356  *
 357  *  \return x^5
 358  */
 359 template<typename T>
 360 T power5(T x)
 361 {
 362     return x * power4(x);
 363 }
 364
 365 /*! \brief calculate x^6
 366  *
 367  *  \tparam T  Type of argument and return value
 368  *  \param  x  argument
 369  *
 370  *  \return x^6
 371  */
 372 template<typename T>
 373 T power6(T x)
 374 {
 375     return square(power3(x));
 376 }
 377
 378 /*! \brief calculate x^12
 379  *
 380  *  \tparam T  Type of argument and return value
 381  *  \param  x  argument
 382  *
 383  *  \return x^12
 384  */
 385 template<typename T>
 386 T power12(T x)
 387 {
 388     return square(power6(x));
 389 }
 390
 391 /*! \brief Maclaurin series for sinh(x)/x.
 392  *
 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.
 396  */
 397 static inline real series_sinhx(real x)
 398 {
 399     real x2 = x * x;
 400     return (1
 401             + (x2 / 6.0_real)
 402                       * (1
 403                          + (x2 / 20.0_real)
 404                                    * (1 + (x2 / 42.0_real) * (1 + (x2 / 72.0_real) * (1 + (x2 / 110.0_real))))));
 405 }
 406
 407 /*! \brief Inverse error function, double precision.
 408  *
 409  *  \param x Argument, should be in the range -1.0 < x < 1.0
 410  *
 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.
 413  */
 414 double erfinv(double x);
 415
 416 /*! \brief Inverse error function, single precision.
 417  *
 418  *  \param x Argument, should be in the range -1.0 < x < 1.0
 419  *
 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.
 422  */
 423 float erfinv(float x);
 424
 425 /*! \brief Exact integer division, 32bit.
 426  *
 427  * \param a dividend. Function asserts that it is a multiple of divisor
 428  * \param b divisor
 429  *
 430  * \return quotient of division
 431  */
 432 constexpr int32_t exactDiv(int32_t a, int32_t b)
 433 {
 434     return GMX_ASSERT(a % b == 0, "exactDiv called with non-divisible arguments"), a / b;
 435 }
 436
 437 //! Exact integer division, 64bit.
 438 constexpr int64_t exactDiv(int64_t a, int64_t b)
 439 {
 440     return GMX_ASSERT(a % b == 0, "exactDiv called with non-divisible arguments"), a / b;
 441 }
 442
 443 /*! \brief Round float to int
 444  *
 445  * Rounding behavior is round to nearest. Rounding of halfway cases is implemention defined
 446  * (either halway to even or halway away from zero).
 447  */
 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.
 455  */
 456 static inline int roundToInt(float x)
 457 {
 458     return static_cast<int>(rintf(x));
 459 }
 460 //! Round double to int
 461 static inline int roundToInt(double x)
 462 {
 463     return static_cast<int>(rint(x));
 464 }
 465 //! Round float to int64_t
 466 static inline int64_t roundToInt64(float x)
 467 {
 468     return static_cast<int>(rintf(x));
 469 }
 470 //! Round double to int64_t
 471 static inline int64_t roundToInt64(double x)
 472 {
 473     return static_cast<int>(rint(x));
 474 }
 475
 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 */
 481 __host__ __device__
 482 #endif
 483         static inline constexpr bool
 484         isPowerOfTwo(const T v)
 485 {
 486     return (v > 0) && ((v & (v - 1)) == 0);
 487 }
 488
 489 } // namespace gmx
 490
 491
 492 #endif // GMX_MATH_FUNCTIONS_H