src/gromacs/math/tests/functions.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2018,2019 by the GROMACS development team.
   5  * Copyright (c) 2021, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \internal \file
  37  * \brief
  38  * Tests for simple math functions.
  39  *
  40  * \author Erik Lindahl <erik.lindahl@gmail.com>
  41  * \ingroup module_math
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "gromacs/math/functions.h"
  46
  47 #include <cmath>
  48 #include <cstdint>
  49
  50 #include <gtest/gtest.h>
  51
  52 #include "testutils/refdata.h"
  53 #include "testutils/testasserts.h"
  54
  55 namespace
  56 {
  57
  58 TEST(FunctionTest, StaticLog2)
  59 {
  60     gmx::test::TestReferenceData    data;
  61     gmx::test::TestReferenceChecker checker(data.rootChecker());
  62     std::vector<int>                result(11);
  63
  64     // This needs to be expanded manually since it is evaluated at compile time,
  65     // and the compiler chokes if we put it as formal arguments to push_back
  66     result[0]  = gmx::StaticLog2<1>::value;
  67     result[1]  = gmx::StaticLog2<2>::value;
  68     result[2]  = gmx::StaticLog2<3>::value;
  69     result[3]  = gmx::StaticLog2<4>::value;
  70     result[4]  = gmx::StaticLog2<5>::value;
  71     result[5]  = gmx::StaticLog2<6>::value;
  72     result[6]  = gmx::StaticLog2<7>::value;
  73     result[7]  = gmx::StaticLog2<8>::value;
  74     result[8]  = gmx::StaticLog2<0xFFFFFFFF>::value;
  75     result[9]  = gmx::StaticLog2<9876543210>::value; // > 32 bits
  76     result[10] = gmx::StaticLog2<0xFFFFFFFFFFFFFFFFULL>::value;
  77
  78     checker.checkSequence(result.begin(), result.end(), "StaticLog2");
  79 }
  80
  81 TEST(FunctionTest, Log2I32Bit)
  82 {
  83     gmx::test::TestReferenceData    data;
  84     gmx::test::TestReferenceChecker checker(data.rootChecker());
  85     std::vector<int>                result;
  86
  87     for (std::uint32_t i = 1; i <= 0xF; i++)
  88     {
  89         result.push_back(gmx::log2I(i));
  90     }
  91
  92     for (std::uint32_t i = 0; i <= 0xF; i++)
  93     {
  94         result.push_back(gmx::log2I(static_cast<std::uint32_t>(0xFFFFFFF0 + i)));
  95     }
  96     checker.checkSequence(result.begin(), result.end(), "Log2I32Bit");
  97 }
  98
  99 TEST(FunctionTest, Log2I64Bit)
 100 {
 101     gmx::test::TestReferenceData    data;
 102     gmx::test::TestReferenceChecker checker(data.rootChecker());
 103     std::vector<int>                result;
 104
 105     for (std::uint64_t i = 1; i <= 0xF; i++)
 106     {
 107         result.push_back(gmx::log2I(i));
 108     }
 109
 110     for (std::uint64_t i = 0; i <= 0x1F; i++)
 111     {
 112         result.push_back(gmx::log2I(static_cast<std::uint64_t>(0xFFFFFFF0ULL + i)));
 113     }
 114
 115     for (std::uint64_t i = 0; i <= 0xF; i++)
 116     {
 117         result.push_back(gmx::log2I(static_cast<std::uint64_t>(0xFFFFFFFFFFFFFFF0ULL + i)));
 118     }
 119     checker.checkSequence(result.begin(), result.end(), "Log2I64Bit");
 120 }
 121
 122 TEST(FunctionTest, GreatestCommonDivisor)
 123 {
 124     EXPECT_EQ(8, gmx::greatestCommonDivisor(24, 64));
 125     EXPECT_EQ(8, gmx::greatestCommonDivisor(64, 24));
 126     EXPECT_EQ(1, gmx::greatestCommonDivisor(169, 289));
 127 }
 128
 129 TEST(FunctionTest, InvsqrtFloat)
 130 {
 131     gmx::test::TestReferenceData    data;
 132     gmx::test::TestReferenceChecker checker(data.rootChecker());
 133     std::vector<float>              result;
 134
 135     for (float f = 1.0; f < 10.0; f += 1.0)
 136     {
 137         result.push_back(gmx::invsqrt(f));
 138     }
 139     checker.checkSequence(result.begin(), result.end(), "InvsqrtFloat");
 140 }
 141
 142 TEST(FunctionTest, InvsqrtDouble)
 143 {
 144     gmx::test::TestReferenceData    data;
 145     gmx::test::TestReferenceChecker checker(data.rootChecker());
 146     std::vector<double>             result;
 147
 148     for (double f = 1.0; f < 10.0; f += 1.0) // NOLINT(clang-analyzer-security.FloatLoopCounter)
 149     {
 150         result.push_back(gmx::invsqrt(f));
 151     }
 152     checker.checkSequence(result.begin(), result.end(), "InvsqrtDouble");
 153 }
 154
 155 TEST(FunctionTest, InvsqrtInteger)
 156 {
 157     gmx::test::TestReferenceData    data;
 158     gmx::test::TestReferenceChecker checker(data.rootChecker());
 159     std::vector<double>             result;
 160
 161     for (int i = 1; i < 10; i++)
 162     {
 163         result.push_back(gmx::invsqrt(i));
 164     }
 165     checker.checkSequence(result.begin(), result.end(), "InvsqrtInteger");
 166 }
 167
 168 TEST(FunctionTest, InvcbrtFloat)
 169 {
 170     gmx::test::TestReferenceData    data;
 171     gmx::test::TestReferenceChecker checker(data.rootChecker());
 172     std::vector<float>              result;
 173
 174     for (int f : { -5, -4, -3, -2, -1, 1, 2, 3, 4 })
 175     {
 176         result.push_back(gmx::invcbrt(f));
 177     }
 178     checker.checkSequence(result.begin(), result.end(), "InvcbrtFloat");
 179 }
 180
 181 TEST(FunctionTest, InvcbrtDouble)
 182 {
 183     gmx::test::TestReferenceData    data;
 184     gmx::test::TestReferenceChecker checker(data.rootChecker());
 185     std::vector<double>             result;
 186
 187     for (double d : { -5, -4, -3, -2, -1, 1, 2, 3, 4 })
 188     {
 189         result.push_back(gmx::invcbrt(d));
 190     }
 191     checker.checkSequence(result.begin(), result.end(), "InvcbrtDouble");
 192 }
 193
 194 TEST(FunctionTest, InvcbrtInteger)
 195 {
 196     gmx::test::TestReferenceData    data;
 197     gmx::test::TestReferenceChecker checker(data.rootChecker());
 198     std::vector<double>             result;
 199
 200     for (int i : { -5, -4, -3, -2, -1, 1, 2, 3, 4 })
 201     {
 202         result.push_back(gmx::invcbrt(i));
 203     }
 204     checker.checkSequence(result.begin(), result.end(), "InvcbrtInteger");
 205 }
 206
 207
 208 TEST(FunctionTest, SixthrootFloat)
 209 {
 210     gmx::test::TestReferenceData    data;
 211     gmx::test::TestReferenceChecker checker(data.rootChecker());
 212     std::vector<float>              result;
 213
 214     for (float f = 0; f < 10.0; f += 1.0)
 215     {
 216         result.push_back(gmx::sixthroot(f));
 217     }
 218     checker.checkSequence(result.begin(), result.end(), "SixthrootFloat");
 219 }
 220
 221 TEST(FunctionTest, SixthrootDouble)
 222 {
 223     gmx::test::TestReferenceData    data;
 224     gmx::test::TestReferenceChecker checker(data.rootChecker());
 225     std::vector<double>             result;
 226
 227     for (double d = 0; d < 10.0; d += 1.0) // NOLINT(clang-analyzer-security.FloatLoopCounter)
 228     {
 229         result.push_back(gmx::sixthroot(d));
 230     }
 231     checker.checkSequence(result.begin(), result.end(), "SixthrootDouble");
 232 }
 233
 234 TEST(FunctionTest, SixthrootInteger)
 235 {
 236     gmx::test::TestReferenceData    data;
 237     gmx::test::TestReferenceChecker checker(data.rootChecker());
 238     std::vector<double>             result;
 239
 240     result.reserve(10);
 241     for (int i = 0; i < 10; i++)
 242     {
 243         result.push_back(gmx::sixthroot(i));
 244     }
 245     checker.checkSequence(result.begin(), result.end(), "SixthrootInteger");
 246 }
 247
 248 TEST(FunctionTest, InvsixthrootFloat)
 249 {
 250     gmx::test::TestReferenceData    data;
 251     gmx::test::TestReferenceChecker checker(data.rootChecker());
 252     std::vector<float>              result;
 253
 254     for (float f = 1.0; f < 10.0; f += 1.0)
 255     {
 256         result.push_back(gmx::invsixthroot(f));
 257     }
 258     checker.checkSequence(result.begin(), result.end(), "InvsixthrootFloat");
 259 }
 260
 261 TEST(FunctionTest, InvsixthrootDouble)
 262 {
 263     gmx::test::TestReferenceData    data;
 264     gmx::test::TestReferenceChecker checker(data.rootChecker());
 265     std::vector<double>             result;
 266
 267     for (double d = 1.0; d < 10.0; d += 1.0) // NOLINT(clang-analyzer-security.FloatLoopCounter)
 268     {
 269         result.push_back(gmx::invsixthroot(d));
 270     }
 271     checker.checkSequence(result.begin(), result.end(), "InvsixthrootDouble");
 272 }
 273
 274 TEST(FunctionTest, InvsixthrootInteger)
 275 {
 276     gmx::test::TestReferenceData    data;
 277     gmx::test::TestReferenceChecker checker(data.rootChecker());
 278     std::vector<double>             result;
 279
 280     for (int i = 1; i < 10; i++)
 281     {
 282         result.push_back(gmx::invsixthroot(i));
 283     }
 284     checker.checkSequence(result.begin(), result.end(), "InvsixthrootInteger");
 285 }
 286
 287 TEST(FunctionTest, Powers)
 288 {
 289     // These should be remarkably difficult to screw up, but test each
 290     // of them once anyway with integer and fp arguments to catch typos.
 291     EXPECT_EQ(4, gmx::square(2));
 292     EXPECT_EQ(8, gmx::power3(2));
 293     EXPECT_EQ(16, gmx::power4(2));
 294     EXPECT_EQ(32, gmx::power5(2));
 295     EXPECT_EQ(64, gmx::power6(2));
 296     EXPECT_EQ(4096, gmx::power12(2));
 297
 298     EXPECT_EQ(6.25, gmx::square(2.5));
 299     EXPECT_EQ(15.625, gmx::power3(2.5));
 300     EXPECT_EQ(39.0625, gmx::power4(2.5));
 301     EXPECT_EQ(97.65625, gmx::power5(2.5));
 302     EXPECT_EQ(244.140625, gmx::power6(2.5));
 303     EXPECT_EQ(59604.644775390625, gmx::power12(2.5));
 304 }
 305
 306 TEST(FunctionTest, ErfInvFloat)
 307 {
 308     gmx::test::TestReferenceData    data;
 309     gmx::test::TestReferenceChecker checker(data.rootChecker());
 310     std::vector<float>              result;
 311     int                             npoints = 10;
 312
 313     for (int i = 0; i < npoints; i++)
 314     {
 315         float r = float(2 * i - npoints + 1) / float(npoints);
 316
 317         result.push_back(gmx::erfinv(r));
 318     }
 319     checker.checkSequence(result.begin(), result.end(), "ErfInvFloat");
 320 }
 321
 322 TEST(FunctionTest, ErfInvDouble)
 323 {
 324     gmx::test::TestReferenceData    data;
 325     gmx::test::TestReferenceChecker checker(data.rootChecker());
 326     std::vector<double>             result;
 327     int                             npoints = 10;
 328
 329     for (int i = 0; i < npoints; i++)
 330     {
 331         double r = double(2 * i - npoints + 1) / npoints;
 332
 333         result.push_back(gmx::erfinv(r));
 334     }
 335     checker.checkSequence(result.begin(), result.end(), "ErfInvDouble");
 336 }
 337
 338 TEST(FunctionTest, ErfAndErfInvAreInversesFloat)
 339 {
 340     int npoints = 1000;
 341
 342     for (int i = 0; i < npoints; i++)
 343     {
 344         float r = float(2 * i - npoints + 1) / float(npoints);
 345         EXPECT_FLOAT_EQ_TOL(r, std::erf(gmx::erfinv(r)), gmx::test::ulpTolerance(10));
 346     }
 347 }
 348
 349 TEST(FunctionTest, ErfAndErfInvAreInversesDouble)
 350 {
 351     int npoints = 1000;
 352
 353     for (int i = 0; i < npoints; i++)
 354     {
 355         double r = double(2 * i - npoints + 1) / npoints;
 356         EXPECT_DOUBLE_EQ_TOL(r, std::erf(gmx::erfinv(r)), gmx::test::ulpTolerance(10));
 357     }
 358 }
 359
 360 template<typename T>
 361 class FunctionTestIntegerTypes : public ::testing::Test
 362 {
 363 };
 364
 365 typedef ::testing::Types<char, unsigned char, int, unsigned int, long, unsigned long> IntegerTypes;
 366 TYPED_TEST_CASE(FunctionTestIntegerTypes, IntegerTypes);
 367
 368 TYPED_TEST(FunctionTestIntegerTypes, IsPowerOfTwo)
 369 {
 370     if (std::is_signed_v<TypeParam>)
 371     {
 372         EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(std::numeric_limits<TypeParam>::min()));
 373         EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(-16));
 374         EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(-3));
 375         EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(-2));
 376         EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(-1));
 377     }
 378     EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(0));
 379     EXPECT_EQ(true, gmx::isPowerOfTwo<TypeParam>(1));
 380     EXPECT_EQ(true, gmx::isPowerOfTwo<TypeParam>(2));
 381     EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(3));
 382     EXPECT_EQ(true, gmx::isPowerOfTwo<TypeParam>(4));
 383     EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(5));
 384     EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(6));
 385     EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(24));
 386     EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(63));
 387     EXPECT_EQ(true, gmx::isPowerOfTwo<TypeParam>(64));
 388     EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(66));
 389     // Max for any type is always 2^x - 1
 390     EXPECT_EQ(false, gmx::isPowerOfTwo<TypeParam>(std::numeric_limits<TypeParam>::max()));
 391 }
 392
 393 } // namespace