src/gromacs/math/functions.h

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