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 479