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
  69                              + 1; //!< Variable value used for recursive static calculation of Log2(int)
  70 };
  71
  72 /*! \brief Specialization of StaticLog2<n> for n==1.
  73  *
  74  *  This specialization provides the final value in the recursion; never
  75  *  call it directly, but use StaticLog2<n>::value.
  76  */
  77 template<>
  78 struct StaticLog2<1>
  79 {
  80     static const int value = 0; //!< Base value for recursive static calculation of Log2(int)
  81 };
  82
  83 /*! \brief Specialization of StaticLog2<n> for n==0.
  84  *
  85  *  This specialization should never actually be used since log2(0) is
  86  *  negative infinity. However, since Log2() is often used to calculate the number
  87  *  of bits needed for a number, we might be using the value 0 with a conditional
  88  *  statement around the logarithm. Depending on the compiler the expansion of
  89  *  the template can occur before the conditional statement, so to avoid infinite
  90  *  recursion we need a specialization for the case n==0.
  91  */
  92 template<>
  93 struct StaticLog2<0>
  94 {
  95     static const int value = -1; //!< Base value for recursive static calculation of Log2(int)
  96 };
  97
  98
  99 /*! \brief Compute floor of logarithm to base 2, 32 bit signed argument
 100  *
 101  *  \param x 32-bit signed argument
 102  *
 103  *  \return log2(x)
 104  *
 105  *  \note This version of the overloaded function will assert that x is
 106  *        not negative.
 107  */
 108 unsigned int 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 log2I(std::int64_t x);
 120
 121 /*! \brief Compute floor of logarithm to base 2, 32 bit unsigned argument
 122  *
 123  *  \param x 32-bit unsigned argument
 124  *
 125  *  \return log2(x)
 126  *
 127  *  \note This version of the overloaded function uses unsigned arguments to
 128  *        be able to handle arguments using all 32 bits.
 129  */
 130 unsigned int log2I(std::uint32_t x);
 131
 132 /*! \brief Compute floor of logarithm to base 2, 64 bit unsigned argument
 133  *
 134  *  \param x 64-bit unsigned argument
 135  *
 136  *  \return log2(x)
 137  *
 138  *  \note This version of the overloaded function uses unsigned arguments to
 139  *        be able to handle arguments using all 64 bits.
 140  */
 141 unsigned int log2I(std::uint64_t x);
 142
 143 /*! \brief Find greatest common divisor of two numbers
 144  *
 145  *  \param p First number, positive
 146  *  \param q Second number, positive
 147  *
 148  * \return Greatest common divisor of p and q
 149  */
 150 std::int64_t greatestCommonDivisor(std::int64_t p, std::int64_t q);
 151
 152
 153 /*! \brief Calculate 1.0/sqrt(x) in single precision
 154  *
 155  * \param  x  Positive value to calculate inverse square root for
 156  *
 157  * For now this is implemented with std::sqrt(x) since gcc seems to do a
 158  * decent job optimizing it. However, we might decide to use instrinsics
 159  * or compiler-specific functions in the future.
 160  *
 161  * \return 1.0/sqrt(x)
 162  */
 163 static inline float invsqrt(float x)
 164 {
 165     return 1.0F / std::sqrt(x);
 166 }
 167
 168 /*! \brief Calculate 1.0/sqrt(x) in double precision, but single range
 169  *
 170  * \param  x  Positive value to calculate inverse square root for, must be
 171  *            in the input domain valid for single precision.
 172  *
 173  * For now this is implemented with std::sqrt(x). However, we might
 174  * decide to use instrinsics or compiler-specific functions in the future, and
 175  * then we want to have the freedom to do the first step in single precision.
 176  *
 177  * \return 1.0/sqrt(x)
 178  */
 179 static inline double invsqrt(double x)
 180 {
 181     return 1.0 / std::sqrt(x);
 182 }
 183
 184 /*! \brief Calculate 1.0/sqrt(x) for integer x in double precision.
 185  *
 186  * \param  x  Positive value to calculate inverse square root for.
 187  *
 188  * \return 1.0/sqrt(x)
 189  */
 190 static inline double invsqrt(int x)
 191 {
 192     return invsqrt(static_cast<double>(x));
 193 }
 194
 195 /*! \brief Calculate inverse cube root of x in single precision
 196  *
 197  *  \param  x  Argument
 198  *
 199  *  \return x^(-1/3)
 200  *
 201  *  This routine is typically faster than using std::pow().
 202  */
 203 static inline float invcbrt(float x)
 204 {
 205     return 1.0F / std::cbrt(x);
 206 }
 207
 208 /*! \brief Calculate inverse sixth root of x in double precision
 209  *
 210  *  \param  x  Argument
 211  *
 212  *  \return x^(-1/3)
 213  *
 214  *  This routine is typically faster than using std::pow().
 215  */
 216 static inline double invcbrt(double x)
 217 {
 218     return 1.0 / std::cbrt(x);
 219 }
 220
 221 /*! \brief Calculate inverse sixth root of integer x in double precision
 222  *
 223  *  \param  x  Argument
 224  *
 225  *  \return x^(-1/3)
 226  *
 227  *  This routine is typically faster than using std::pow().
 228  */
 229 static inline double invcbrt(int x)
 230 {
 231     return 1.0 / std::cbrt(x);
 232 }
 233
 234 /*! \brief Calculate sixth root of x in single precision.
 235  *
 236  *  \param  x  Argument, must be greater than or equal to zero.
 237  *
 238  *  \return x^(1/6)
 239  *
 240  *  This routine is typically faster than using std::pow().
 241  */
 242 static inline float sixthroot(float x)
 243 {
 244     return std::sqrt(std::cbrt(x));
 245 }
 246
 247 /*! \brief Calculate sixth root of x in double precision.
 248  *
 249  *  \param  x  Argument, must be greater than or equal to zero.
 250  *
 251  *  \return x^(1/6)
 252  *
 253  *  This routine is typically faster than using std::pow().
 254  */
 255 static inline double sixthroot(double x)
 256 {
 257     return std::sqrt(std::cbrt(x));
 258 }
 259
 260 /*! \brief Calculate sixth root of integer x, return double.
 261  *
 262  *  \param  x  Argument, must be greater than or equal to zero.
 263  *
 264  *  \return x^(1/6)
 265  *
 266  *  This routine is typically faster than using std::pow().
 267  */
 268 static inline double sixthroot(int x)
 269 {
 270     return std::sqrt(std::cbrt(x));
 271 }
 272
 273 /*! \brief Calculate inverse sixth root of x in single precision
 274  *
 275  *  \param  x  Argument, must be greater than zero.
 276  *
 277  *  \return x^(-1/6)
 278  *
 279  *  This routine is typically faster than using std::pow().
 280  */
 281 static inline float invsixthroot(float x)
 282 {
 283     return invsqrt(std::cbrt(x));
 284 }
 285
 286 /*! \brief Calculate inverse sixth root of x in double 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 double invsixthroot(double x)
 295 {
 296     return invsqrt(std::cbrt(x));
 297 }
 298
 299 /*! \brief Calculate inverse sixth root of integer x in double precision
 300  *
 301  *  \param  x  Argument, must be greater than zero.
 302  *
 303  *  \return x^(-1/6)
 304  *
 305  *  This routine is typically faster than using std::pow().
 306  */
 307 static inline double invsixthroot(int x)
 308 {
 309     return invsqrt(std::cbrt(x));
 310 }
 311
 312 /*! \brief calculate x^2
 313  *
 314  *  \tparam T  Type of argument and return value
 315  *  \param  x  argument
 316  *
 317  *  \return x^2
 318  */
 319 template<typename T>
 320 T square(T x)
 321 {
 322     return x * x;
 323 }
 324
 325 /*! \brief calculate x^3
 326  *
 327  *  \tparam T  Type of argument and return value
 328  *  \param  x  argument
 329  *
 330  *  \return x^3
 331  */
 332 template<typename T>
 333 T power3(T x)
 334 {
 335     return x * square(x);
 336 }
 337
 338 /*! \brief calculate x^4
 339  *
 340  *  \tparam T  Type of argument and return value
 341  *  \param  x  argument
 342  *
 343  *  \return x^4
 344  */
 345 template<typename T>
 346 T power4(T x)
 347 {
 348     return square(square(x));
 349 }
 350
 351 /*! \brief calculate x^5
 352  *
 353  *  \tparam T  Type of argument and return value
 354  *  \param  x  argument
 355  *
 356  *  \return x^5
 357  */
 358 template<typename T>
 359 T power5(T x)
 360 {
 361     return x * power4(x);
 362 }
 363
 364 /*! \brief calculate x^6
 365  *
 366  *  \tparam T  Type of argument and return value
 367  *  \param  x  argument
 368  *
 369  *  \return x^6
 370  */
 371 template<typename T>
 372 T power6(T x)
 373 {
 374     return square(power3(x));
 375 }
 376
 377 /*! \brief calculate x^12
 378  *
 379  *  \tparam T  Type of argument and return value
 380  *  \param  x  argument
 381  *
 382  *  \return x^12
 383  */
 384 template<typename T>
 385 T power12(T x)
 386 {
 387     return square(power6(x));
 388 }
 389
 390 /*! \brief Maclaurin series for sinh(x)/x.
 391  *
 392  * Used for NH chains and MTTK pressure control.
 393  * Here, we compute it to 10th order, which might be an overkill.
 394  * 8th is probably enough, but it's not very much more expensive.
 395  */
 396 static inline real series_sinhx(real x)
 397 {
 398     real x2 = x * x;
 399     return (1
 400             + (x2 / 6.0_real)
 401                       * (1
 402                          + (x2 / 20.0_real)
 403                                    * (1 + (x2 / 42.0_real) * (1 + (x2 / 72.0_real) * (1 + (x2 / 110.0_real))))));
 404 }
 405
 406 /*! \brief Inverse error function, double precision.
 407  *
 408  *  \param x Argument, should be in the range -1.0 < x < 1.0
 409  *
 410  *  \return The inverse of the error function if the argument is inside the
 411  *          range, +/- infinity if it is exactly 1.0 or -1.0, and NaN otherwise.
 412  */
 413 double erfinv(double x);
 414
 415 /*! \brief Inverse error function, single precision.
 416  *
 417  *  \param x Argument, should be in the range -1.0 < x < 1.0
 418  *
 419  *  \return The inverse of the error function if the argument is inside the
 420  *          range, +/- infinity if it is exactly 1.0 or -1.0, and NaN otherwise.
 421  */
 422 float erfinv(float x);
 423
 424 /*! \brief Exact integer division, 32bit.
 425  *
 426  * \param a dividend. Function asserts that it is a multiple of divisor
 427  * \param b divisor
 428  *
 429  * \return quotient of division
 430  */
 431 constexpr int32_t exactDiv(int32_t a, int32_t b)
 432 {
 433     return GMX_ASSERT(a % b == 0, "exactDiv called with non-divisible arguments"), a / b;
 434 }
 435
 436 //! Exact integer division, 64bit.
 437 constexpr int64_t exactDiv(int64_t a, int64_t b)
 438 {
 439     return GMX_ASSERT(a % b == 0, "exactDiv called with non-divisible arguments"), a / b;
 440 }
 441
 442 /*! \brief Round float to int
 443  *
 444  * Rounding behavior is round to nearest. Rounding of halfway cases is implemention defined
 445  * (either halway to even or halway away from zero).
 446  */
 447 /* Implementation details: It is assumed that FE_TONEAREST is default and not changed by anyone.
 448  * Currently the implementation is using rint(f) because 1) on all known HW that is faster than
 449  * lround and 2) some compilers (e.g. clang (#22944) and icc) don't optimize (l)lrint(f) well.
 450  * GCC(>=4.7) optimizes (l)lrint(f) well but with "-fno-math-errno -funsafe-math-optimizations"
 451  * rint(f) is optimized as well. This avoids using intrinsics.
 452  * rint(f) followed by float/double to int/int64 conversion produces the same result as directly
 453  * rounding to int/int64.
 454  */
 455 static inline int roundToInt(float x)
 456 {
 457     return static_cast<int>(rintf(x));
 458 }
 459 //! Round double to int
 460 static inline int roundToInt(double x)
 461 {
 462     return static_cast<int>(rint(x));
 463 }
 464 //! Round float to int64_t
 465 static inline int64_t roundToInt64(float x)
 466 {
 467     return static_cast<int>(rintf(x));
 468 }
 469 //! Round double to int64_t
 470 static inline int64_t roundToInt64(double x)
 471 {
 472     return static_cast<int>(rint(x));
 473 }
 474
 475 } // namespace gmx
 476
 477
 478 #endif // GMX_MATH_FUNCTIONS_H