src/gromacs/math/functions.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  * \brief
  39  * Implements simple math functions
  40  *
  41  * \author Erik Lindahl <erik.lindahl@gmail.com>
  42  * \ingroup module_math
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "functions.h"
  48
  49 #include "config.h"
  50
  51 #include <cstdint>
  52
  53 #include <array>
  54 #include <limits>
  55
  56 #if GMX_NATIVE_WINDOWS
  57 #    include <intrin.h> // _BitScanReverse, _BitScanReverse64
  58 #endif
  59
  60 #include "gromacs/math/utilities.h"
  61 #include "gromacs/utility/gmxassert.h"
  62
  63 namespace gmx
  64 {
  65
  66 unsigned int
  67 log2I(std::uint32_t n)
  68 {
  69     GMX_ASSERT(n > 0, "The behavior of log(0) is undefined");
  70 #if HAVE_BUILTIN_CLZ
  71     // gcc, clang. xor with sign bit should be optimized out
  72     return __builtin_clz(n) ^ 31U;
  73 #elif HAVE_BITSCANREVERSE
  74     // icc, MSVC
  75     {
  76         unsigned long res;
  77         _BitScanReverse(&res, static_cast<unsigned long>(n));
  78         return static_cast<unsigned int>(res);
  79     }
  80 #elif HAVE_CNTLZ4
  81     return 31 - __cntlz4(n);
  82 #else
  83     // http://graphics.stanford.edu/~seander/bithacks.html#IntegerLogLookup
  84
  85     static const std::array<char, 256>
  86     log2TableByte =
  87     {{
  88          0, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
  89          4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4,
  90          5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
  91          5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,
  92          6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  93          6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  94          6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  95          6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6,
  96          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  97          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  98          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
  99          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 100          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 101          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 102          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7,
 103          7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7
 104      }};
 105
 106     unsigned int result;
 107     unsigned int tmp1, tmp2;
 108
 109     if ((tmp1 = n >> 16) != 0)
 110     {
 111         result = ((tmp2 = tmp1 >> 8) != 0) ? 24 + log2TableByte[tmp2] : 16 + log2TableByte[tmp1];
 112     }
 113     else
 114     {
 115         result = ((tmp2 = n >> 8) != 0) ? 8 + log2TableByte[tmp2] : log2TableByte[n];
 116     }
 117     return result;
 118 #endif
 119 }
 120
 121
 122 unsigned int
 123 log2I(std::uint64_t n)
 124 {
 125     GMX_ASSERT(n > 0, "The behavior of log(0) is undefined");
 126 #if HAVE_BUILTIN_CLZLL
 127     // gcc, icc, clang. xor with sign bit should be optimized out
 128     return __builtin_clzll(n) ^ 63U;
 129 #elif HAVE_BITSCANREVERSE64
 130     unsigned long res;
 131     _BitScanReverse64(&res, static_cast<unsigned __int64>(n));
 132     return static_cast<unsigned int>(res);
 133 #elif HAVE_CNTLZ8
 134     return 63 - __cntlz8(n);
 135 #else
 136
 137     // No 64-bit log2 instrinsic available. Solve it by calling our internal
 138     // 32-bit version (which in turn might defer to a software solution)
 139
 140     std::uint32_t high32Bits = static_cast<std::uint32_t>(n>>32);
 141     std::uint32_t result;
 142
 143     if (high32Bits)
 144     {
 145         result = log2I(high32Bits) + 32;
 146     }
 147     else
 148     {
 149         result = log2I(static_cast<std::uint32_t>(n));
 150     }
 151
 152     return result;
 153 #endif
 154 }
 155
 156 unsigned int
 157 log2I(std::int32_t n)
 158 {
 159     GMX_ASSERT(n > 0, "The behavior of log(n) for n<=0 is undefined");
 160     return log2I(static_cast<std::uint32_t>(n));
 161 }
 162
 163 unsigned int
 164 log2I(std::int64_t n)
 165 {
 166     GMX_ASSERT(n > 0, "The behavior of log(n) for n<=0 is undefined");
 167     return log2I(static_cast<std::uint64_t>(n));
 168 }
 169
 170 std::int64_t
 171 greatestCommonDivisor(std::int64_t   p,
 172                       std::int64_t   q)
 173 {
 174     while (q != 0)
 175     {
 176         std::int64_t tmp = q;
 177         q                = p % q;
 178         p                = tmp;
 179     }
 180     return p;
 181 }
 182
 183 double
 184 erfinv(double x)
 185 {
 186     double xabs = std::abs(x);
 187
 188     if (xabs > 1.0)
 189     {
 190         return std::nan("");
 191     }
 192
 193     if (x == 1.0)
 194     {
 195         return std::numeric_limits<double>::infinity();
 196     }
 197
 198     if (x == -1.0)
 199     {
 200         return -std::numeric_limits<double>::infinity();
 201     }
 202
 203     double res;
 204
 205     if (xabs <= 0.7)
 206     {
 207         // Rational approximation in range [0,0.7]
 208         double z = x*x;
 209         double P = (((-0.140543331 * z + 0.914624893) * z - 1.645349621) * z + 0.886226899);
 210         double Q = ((((0.012229801 * z - 0.329097515) * z + 1.442710462) * z - 2.118377725) * z + 1.0);
 211         res = x * P/Q;
 212     }
 213     else
 214     {
 215         // Rational approximation in range [0.7,1)
 216         double z = std::sqrt(-std::log((1.0 - std::abs(x))/2.0));
 217         double P = ((1.641345311 * z + 3.429567803) * z - 1.624906493) * z - 1.970840454;
 218         double Q = (1.637067800 * z + 3.543889200) * z + 1.0;
 219         res = std::copysign(1.0, x) * P/Q;
 220     }
 221
 222     // Double precision requires two N-R iterations
 223     res = res - (std::erf(res) - x)/( (2.0/std::sqrt(M_PI))*std::exp(-res*res));
 224     res = res - (std::erf(res) - x)/( (2.0/std::sqrt(M_PI))*std::exp(-res*res));
 225
 226     return res;
 227 }
 228
 229 float
 230 erfinv(float x)
 231 {
 232     float xabs = std::abs(x);
 233
 234     if (xabs > 1.0F)
 235     {
 236         return std::nan("");
 237     }
 238
 239     if (x == 1.0F)
 240     {
 241         return std::numeric_limits<float>::infinity();
 242     }
 243
 244     if (x == -1.0F)
 245     {
 246         return -std::numeric_limits<float>::infinity();
 247     }
 248
 249     float res;
 250
 251     if (xabs <= 0.7F)
 252     {
 253         // Rational approximation in range [0,0.7]
 254         float z = x*x;
 255         float P = (((-0.140543331F * z + 0.914624893F) * z - 1.645349621F) * z + 0.886226899F);
 256         float Q = ((((0.012229801F * z - 0.329097515F) * z + 1.442710462F) * z - 2.118377725F) * z + 1.0F);
 257         res = x * P/Q;
 258     }
 259     else
 260     {
 261         // Rational approximation in range [0.7,1)
 262         float z = std::sqrt(-std::log((1.0 - std::abs(x))/2.0F));
 263         float P = ((1.641345311F * z + 3.429567803F) * z - 1.624906493F) * z - 1.970840454F;
 264         float Q = (1.637067800F * z + 3.543889200F) * z + 1.0F;
 265         res = std::copysign(1.0F, x) * P/Q;
 266     }
 267
 268     // Single N-R iteration sufficient for single precision
 269     res = res - (std::erf(res) - x)/( (2.0F/std::sqrt(M_PI))*std::exp(-res*res));
 270
 271     return res;
 272 }
 273
 274 } // namespace gmx