src/gromacs/math/functions.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016, 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/real.h"
  50
  51 namespace gmx
  52 {
  53
  54 /*! \brief Evaluate log2(n) for integer n statically at compile time.
  55  *
  56  * Use as staticLog2<n>::value, where n must be a positive integer.
  57  * Negative n will be reinterpreted as the corresponding unsigned integer,
  58  * and you will get a compile-time error if n==0.
  59  * The calculation is done by recursively dividing n by 2 (until it is 1),
  60  * and incrementing the result by 1 in each step.
  61  *
  62  * \tparam n Value to recursively calculate log2(n) for
  63  */
  64 template<std::uint64_t n>
  65 struct StaticLog2
  66 {
  67     static const int value = StaticLog2<n/2>::value+1; //!< Variable value used for recursive static calculation of Log2(int)
  68 };
  69
  70 /*! \brief Specialization of StaticLog2<n> for n==1.
  71  *
  72  *  This specialization provides the final value in the recursion; never
  73  *  call it directly, but use StaticLog2<n>::value.
  74  */
  75 template<>
  76 struct StaticLog2<1>
  77 {
  78     static const int value = 0; //!< Base value for recursive static calculation of Log2(int)
  79 };
  80
  81 /*! \brief Specialization of StaticLog2<n> for n==0.
  82  *
  83  *  This specialization should never actually be used since log2(0) is
  84  *  negative infinity. However, since Log2() is often used to calculate the number
  85  *  of bits needed for a number, we might be using the value 0 with a conditional
  86  *  statement around the logarithm. Depending on the compiler the expansion of
  87  *  the template can occur before the conditional statement, so to avoid infinite
  88  *  recursion we need a specialization for the case n==0.
  89  */
  90 template<>
  91 struct StaticLog2<0>
  92 {
  93     static const int value = -1; //!< Base value for recursive static calculation of Log2(int)
  94 };
  95
  96
  97 /*! \brief Compute floor of logarithm to base 2, 32 bit signed argument
  98  *
  99  *  \param x 32-bit signed argument
 100  *
 101  *  \return log2(x)
 102  *
 103  *  \note This version of the overloaded function will assert that x is
 104  *        not negative.
 105  */
 106 unsigned int
 107 log2I(std::int32_t x);
 108
 109 /*! \brief Compute floor of logarithm to base 2, 64 bit signed argument
 110  *
 111  *  \param x 64-bit signed argument
 112  *
 113  *  \return log2(x)
 114  *
 115  *  \note This version of the overloaded function will assert that x is
 116  *        not negative.
 117  */
 118 unsigned int
 119 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
 131 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
 143 log2I(std::uint64_t x);
 144
 145 /*! \brief Find greatest common divisor of two numbers
 146  *
 147  *  \param p First number, positive
 148  *  \param q Second number, positive
 149  *
 150  * \return Greatest common divisor of p and q
 151  */
 152 std::int64_t
 153 greatestCommonDivisor(std::int64_t p, std::int64_t q);
 154
 155
 156 /*! \brief Calculate 1.0/sqrt(x) in single precision
 157  *
 158  * \param  x  Positive value to calculate inverse square root for
 159  *
 160  * For now this is implemented with std::sqrt(x) since gcc seems to do a
 161  * decent job optimizing it. However, we might decide to use instrinsics
 162  * or compiler-specific functions in the future.
 163  *
 164  * \return 1.0/sqrt(x)
 165  */
 166 static inline float
 167 invsqrt(float x)
 168 {
 169     return 1.0f/std::sqrt(x);
 170 }
 171
 172 /*! \brief Calculate 1.0/sqrt(x) in double precision, but single range
 173  *
 174  * \param  x  Positive value to calculate inverse square root for, must be
 175  *            in the input domain valid for single precision.
 176  *
 177  * For now this is implemented with std::sqrt(x). However, we might
 178  * decide to use instrinsics or compiler-specific functions in the future, and
 179  * then we want to have the freedom to do the first step in single precision.
 180  *
 181  * \return 1.0/sqrt(x)
 182  */
 183 static inline double
 184 invsqrt(double x)
 185 {
 186     return 1.0/std::sqrt(x);
 187 }
 188
 189 /*! \brief Calculate 1.0/sqrt(x) for integer x in double precision.
 190  *
 191  * \param  x  Positive value to calculate inverse square root for.
 192  *
 193  * \return 1.0/sqrt(x)
 194  */
 195 static inline double
 196 invsqrt(int x)
 197 {
 198     return invsqrt(static_cast<double>(x));
 199 }
 200
 201 /*! \brief Calculate inverse cube root of x in single precision
 202  *
 203  *  \param  x  Argument
 204  *
 205  *  \return x^(-1/3)
 206  *
 207  *  This routine is typically faster than using std::pow().
 208  */
 209 static inline float
 210 invcbrt(float x)
 211 {
 212     return 1.0f/std::cbrt(x);
 213 }
 214
 215 /*! \brief Calculate inverse sixth root of x in double precision
 216  *
 217  *  \param  x  Argument
 218  *
 219  *  \return x^(-1/3)
 220  *
 221  *  This routine is typically faster than using std::pow().
 222  */
 223 static inline double
 224 invcbrt(double x)
 225 {
 226     return 1.0/std::cbrt(x);
 227 }
 228
 229 /*! \brief Calculate inverse sixth root of integer x in double precision
 230  *
 231  *  \param  x  Argument
 232  *
 233  *  \return x^(-1/3)
 234  *
 235  *  This routine is typically faster than using std::pow().
 236  */
 237 static inline double
 238 invcbrt(int x)
 239 {
 240     return 1.0/std::cbrt(x);
 241 }
 242
 243 /*! \brief Calculate sixth root of x in single precision.
 244  *
 245  *  \param  x  Argument, must be greater than or equal to zero.
 246  *
 247  *  \return x^(1/6)
 248  *
 249  *  This routine is typically faster than using std::pow().
 250  */
 251 static inline float
 252 sixthroot(float x)
 253 {
 254     return std::sqrt(std::cbrt(x));
 255 }
 256
 257 /*! \brief Calculate sixth root of x in double precision.
 258  *
 259  *  \param  x  Argument, must be greater than or equal to zero.
 260  *
 261  *  \return x^(1/6)
 262  *
 263  *  This routine is typically faster than using std::pow().
 264  */
 265 static inline double
 266 sixthroot(double x)
 267 {
 268     return std::sqrt(std::cbrt(x));
 269 }
 270
 271 /*! \brief Calculate sixth root of integer x, return double.
 272  *
 273  *  \param  x  Argument, must be greater than or equal to zero.
 274  *
 275  *  \return x^(1/6)
 276  *
 277  *  This routine is typically faster than using std::pow().
 278  */
 279 static inline double
 280 sixthroot(int x)
 281 {
 282     return std::sqrt(std::cbrt(x));
 283 }
 284
 285 /*! \brief Calculate inverse sixth root of x in single precision
 286  *
 287  *  \param  x  Argument, must be greater than zero.
 288  *
 289  *  \return x^(-1/6)
 290  *
 291  *  This routine is typically faster than using std::pow().
 292  */
 293 static inline float
 294 invsixthroot(float x)
 295 {
 296     return invsqrt(std::cbrt(x));
 297 }
 298
 299 /*! \brief Calculate inverse sixth root of 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
 308 invsixthroot(double x)
 309 {
 310     return invsqrt(std::cbrt(x));
 311 }
 312
 313 /*! \brief Calculate inverse sixth root of integer x in double precision
 314  *
 315  *  \param  x  Argument, must be greater than zero.
 316  *
 317  *  \return x^(-1/6)
 318  *
 319  *  This routine is typically faster than using std::pow().
 320  */
 321 static inline double
 322 invsixthroot(int x)
 323 {
 324     return invsqrt(std::cbrt(x));
 325 }
 326
 327 /*! \brief calculate x^2
 328  *
 329  *  \tparam T  Type of argument and return value
 330  *  \param  x  argument
 331  *
 332  *  \return x^2
 333  */
 334 template <typename T>
 335 T
 336 square(T x)
 337 {
 338     return x*x;
 339 }
 340
 341 /*! \brief calculate x^3
 342  *
 343  *  \tparam T  Type of argument and return value
 344  *  \param  x  argument
 345  *
 346  *  \return x^3
 347  */
 348 template <typename T>
 349 T
 350 power3(T x)
 351 {
 352     return x*square(x);
 353 }
 354
 355 /*! \brief calculate x^4
 356  *
 357  *  \tparam T  Type of argument and return value
 358  *  \param  x  argument
 359  *
 360  *  \return x^4
 361  */
 362 template <typename T>
 363 T
 364 power4(T x)
 365 {
 366     return square(square(x));
 367 }
 368
 369 /*! \brief calculate x^5
 370  *
 371  *  \tparam T  Type of argument and return value
 372  *  \param  x  argument
 373  *
 374  *  \return x^5
 375  */
 376 template <typename T>
 377 T
 378 power5(T x)
 379 {
 380     return x*power4(x);
 381 }
 382
 383 /*! \brief calculate x^6
 384  *
 385  *  \tparam T  Type of argument and return value
 386  *  \param  x  argument
 387  *
 388  *  \return x^6
 389  */
 390 template <typename T>
 391 T
 392 power6(T x)
 393 {
 394     return square(power3(x));
 395 }
 396
 397 /*! \brief calculate x^12
 398  *
 399  *  \tparam T  Type of argument and return value
 400  *  \param  x  argument
 401  *
 402  *  \return x^12
 403  */
 404 template <typename T>
 405 T
 406 power12(T x)
 407 {
 408     return square(power6(x));
 409 }
 410
 411 /*! \brief Maclaurin series for sinh(x)/x.
 412  *
 413  * Used for NH chains and MTTK pressure control.
 414  * Here, we compute it to 10th order, which might be an overkill.
 415  * 8th is probably enough, but it's not very much more expensive.
 416  */
 417 static inline real series_sinhx(real x)
 418 {
 419     real x2 = x*x;
 420     return (1 + (x2/6.0)*(1 + (x2/20.0)*(1 + (x2/42.0)*(1 + (x2/72.0)*(1 + (x2/110.0))))));
 421 }
 422
 423 /*! \brief Inverse error function, double precision.
 424  *
 425  *  \param x Argument, should be in the range -1.0 < x < 1.0
 426  *
 427  *  \return The inverse of the error function if the argument is inside the
 428  *          range, +/- infinity if it is exactly 1.0 or -1.0, and NaN otherwise.
 429  */
 430 double
 431 erfinv(double x);
 432
 433 /*! \brief Inverse error function, single precision.
 434  *
 435  *  \param x Argument, should be in the range -1.0 < x < 1.0
 436  *
 437  *  \return The inverse of the error function if the argument is inside the
 438  *          range, +/- infinity if it is exactly 1.0 or -1.0, and NaN otherwise.
 439  */
 440 float
 441 erfinv(float x);
 442
 443 } // namespace gmx
 444
 445
 446 #endif // GMX_MATH_FUNCTIONS_H