src/gromacs/simd/tests/simd_math.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2017,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 #include "gmxpre.h"
  36
  37 #include "gromacs/simd/simd_math.h"
  38
  39 #include "config.h"
  40
  41 #include <cmath>
  42 #include <cstdint>
  43
  44 #include <map>
  45 #include <utility>
  46 #include <vector>
  47
  48 #include "gromacs/math/utilities.h"
  49 #include "gromacs/options/basicoptions.h"
  50 #include "gromacs/simd/simd.h"
  51
  52 #include "testutils/refdata.h"
  53 #include "testutils/testasserts.h"
  54
  55 #include "simd.h"
  56
  57 #if GMX_SIMD
  58
  59 namespace gmx
  60 {
  61 namespace test
  62 {
  63
  64 /*! \cond internal */
  65 /*! \addtogroup module_simd */
  66 /*! \{ */
  67
  68 #    if GMX_SIMD_HAVE_REAL
  69
  70 class SimdMathTest : public SimdTest
  71 {
  72 public:
  73     /*! \brief Type for half-open intervals specifying test ranges */
  74     typedef std::pair<real, real> Range;
  75
  76     /*! \brief Control what is considered matching values
  77      *
  78      * Normal simply means that we request the values to be equal
  79      * to within the specified tolerance.
  80      * However, there are also two more cases that are special:
  81      *
  82      * - Even if we only care about normal (i.e., not denormal) values, some math
  83      *   libraries might clamp the value to zero, which means our SIMD output
  84      *   might not match their values. By using MatchRule::Dtz, we will consider
  85      *   all values both from the reference and test functions that are within the
  86      *   requested ulp tolerance of a denormal number to be equivalent to 0.0.
  87      * - For some older architectures without fused multiply-add units (e.g. x86 SSE2),
  88      *   we might end up clamping the results to zero just before reaching
  89      *   denormal output, since the intermediate results e.g. in polynomial
  90      *   approximations can be smaller than the final one. We often simply don't
  91      *   care about those values, and then one can use
  92      *   MatchRule::ReferenceOrZero to allow the test value to either match
  93      *   the reference or be zero.
  94      */
  95     enum class MatchRule
  96     {
  97         Normal, //!< Match function values
  98         Dtz, //!< Match function values after setting denormals to zero both in test and reference
  99         ReferenceOrZero, //!< Test values can either match reference or be zero
 100     };
 101
 102     const std::map<MatchRule, std::string> matchRuleNames_ = {
 103         { MatchRule::Normal, "Test should match reference." },
 104         { MatchRule::Dtz, "Test should match reference, with denormals treated as 0.0." },
 105         { MatchRule::ReferenceOrZero, "Test should match reference or 0.0." }
 106     };
 107
 108     /*! \brief Settings used for simd math function comparisons */
 109     struct CompareSettings
 110     {
 111         Range        range;     //!< Range over which to test function
 112         std::int64_t ulpTol;    //!< Ulp tolerance
 113         real         absTol;    //!< Absolute tolerance
 114         MatchRule    matchRule; //!< Decide what we consider a match
 115     };
 116
 117     ::testing::AssertionResult compareSimdMathFunction(const char* refFuncExpr,
 118                                                        const char* simdFuncExpr,
 119                                                        const char* compareSettingsExpr,
 120                                                        real        refFunc(real x),
 121                                                        SimdReal gmx_simdcall  simdFunc(SimdReal x),
 122                                                        const CompareSettings& compareSettings);
 123
 124     /*! \brief Generate test point vector
 125      *
 126      *  \param range  The test interval, half open. Upper limit is not included.
 127      *                Pass by value, since we need to modify in method anyway.
 128      *  \param points Number of points to generate. This might be increased
 129      *                slightly to account both for extra special values like 0.0
 130      *                and the SIMD width.
 131      *
 132      * This routine generates a vector with test points separated by constant
 133      * multiplicative factors, based on the range and number of points in the
 134      * class. If the range includes both negative and positive values, points
 135      * will be generated separately for the negative/positive intervals down
 136      * to the smallest real number that can be represented, and we also include
 137      * 0.0 explicitly.
 138      *
 139      * This is highly useful for large test ranges. For example, with a linear
 140      * 1000-point division of the range (1,1e10) the first three values to test
 141      * would be 1, 10000000.999, and 20000000.998, etc. For large values we would
 142      * commonly hit the point where adding the small delta has no effect due to
 143      * limited numerical precision.
 144      * When we instead use this routine, the values will be 1, 1.0239, 1.0471, etc.
 145      * This will spread the entropy over all bits in the IEEE754 representation,
 146      * and be a much better test of all potential input values.
 147      *
 148      *  \note We do not use the static variable s_nPoints in the parent class
 149      *        to avoid altering any value the user has set on the command line; since
 150      *        it's a static member, changing it would have permanent effect.
 151      */
 152     std::vector<real> generateTestPoints(Range range, std::size_t points);
 153
 154     /*! \brief Test routine for the test point vector generation
 155      */
 156     void generateTestPointsTest();
 157 };
 158
 159 /*! \brief Test approximate equality of SIMD vs reference version of a function.
 160  *
 161  * This macro takes vanilla C and SIMD flavors of a function and tests it with
 162  * the number of points, range, and tolerances specified by the test fixture class.
 163  *
 164  * The third option controls the range, tolerances, and match settings.
 165  */
 166 #        define GMX_EXPECT_SIMD_FUNC_NEAR(refFunc, tstFunc, compareSettings) \
 167             EXPECT_PRED_FORMAT3(compareSimdMathFunction, refFunc, tstFunc, compareSettings)
 168
 169 std::vector<real> SimdMathTest::generateTestPoints(Range inputRange, std::size_t inputPoints)
 170 {
 171
 172     std::vector<real> testPoints;
 173     testPoints.reserve(inputPoints);
 174
 175     GMX_RELEASE_ASSERT(inputRange.first < inputRange.second,
 176                        "The start of the interval must come before the end");
 177
 178     std::vector<Range> testRanges;
 179
 180     if (inputRange.first < 0 && inputRange.second > 0)
 181     {
 182         testRanges.emplace_back(Range({ inputRange.first, -std::numeric_limits<real>::min() }));
 183         testRanges.emplace_back(Range({ 0.0, inputRange.second }));
 184     }
 185     else
 186     {
 187         if (inputRange.second == 0)
 188         {
 189             inputRange.second = -std::numeric_limits<real>::min();
 190             inputRange.first  = std::min(inputRange.first, inputRange.second);
 191         }
 192         testRanges.push_back(inputRange);
 193     }
 194
 195     for (Range& range : testRanges)
 196     {
 197         std::size_t points = inputPoints / testRanges.size();
 198
 199         // The value 0 is special, and can only occur at the start of
 200         // the interval after the corrections outside this loop.
 201         // Add it explicitly, and adjust the interval to continue
 202         // at the first valid non-zero positive number.
 203         if (range.first == 0)
 204         {
 205             testPoints.push_back(0.0);
 206             range.first = std::numeric_limits<real>::min();
 207             points--; // Used one point
 208         }
 209
 210         union {
 211             real                                                                               r;
 212             std::conditional<sizeof(real) == sizeof(double), std::int64_t, std::int32_t>::type i;
 213         } low, high, x;
 214
 215         low.r  = range.first;
 216         high.r = range.second;
 217
 218         // IEEE754 floating-point numbers have the cool property that for any range of
 219         // constant sign, for all non-zero numbers a constant (i.e., linear) difference
 220         // in the bitwise representation corresponds to a constant multiplicative factor.
 221         //
 222         // Divide the ulp difference evenly
 223         std::int64_t ulpDiff = high.i - low.i;
 224         // dividend and divisor must both be signed types
 225         std::int64_t ulpDelta    = ulpDiff / static_cast<std::int64_t>(points);
 226         std::int64_t minUlpDelta = (ulpDiff > 0) ? 1 : -1;
 227
 228         if (ulpDelta == 0)
 229         {
 230             // Very short interval or very many points caused round-to-zero.
 231             // Select the smallest possible change, which is one ulp (with correct sign)
 232             ulpDelta = minUlpDelta;
 233             points   = std::abs(ulpDiff);
 234         }
 235
 236         x.r = low.r;
 237         // Use an index-based loop to avoid floating-point comparisons with
 238         // values that might have overflowed. Save one point for the very last
 239         // bitwise value that is part of the interval
 240         for (std::size_t i = 0; i < points - 1; i++)
 241         {
 242             testPoints.push_back(x.r);
 243             x.i += ulpDelta;
 244         }
 245
 246         // Make sure we test the very last point that is inside the interval
 247         x.r = high.r;
 248         x.i -= minUlpDelta;
 249         testPoints.push_back(x.r);
 250     }
 251     return testPoints;
 252 }
 253
 254 /*! \brief Implementation routine to compare SIMD vs reference functions.
 255  *
 256  * \param refFuncExpr         Description of reference function expression
 257  * \param simdFuncExpr        Description of SIMD function expression
 258  * \param compareSettingsExpr Description of compareSettings
 259  * \param refFunc             Reference math function pointer
 260  * \param simdFunc            SIMD math function pointer
 261  * \param compareSettings     Structure with the range, tolerances, and
 262  *                            matching rules to use for the comparison.
 263  *
 264  * \note You should not never call this function directly,  but use the
 265  *       macro GMX_EXPECT_SIMD_FUNC_NEAR(refFunc,tstFunc,matchRule) instead.
 266  */
 267 ::testing::AssertionResult SimdMathTest::compareSimdMathFunction(const char* refFuncExpr,
 268                                                                  const char* simdFuncExpr,
 269                                                                  const char gmx_unused* compareSettingsExpr,
 270                                                                  real refFunc(real x),
 271                                                                  SimdReal gmx_simdcall simdFunc(SimdReal x),
 272                                                                  const CompareSettings& compareSettings)
 273 {
 274     std::vector<real> vx(GMX_SIMD_REAL_WIDTH);
 275     std::vector<real> vref(GMX_SIMD_REAL_WIDTH);
 276     std::vector<real> vtst(GMX_SIMD_REAL_WIDTH);
 277     real              absDiff;
 278     std::int64_t      ulpDiff;
 279     std::int64_t      maxUlpDiff = 0;
 280     real              maxUlpDiffPos;
 281     real              refValMaxUlpDiff, simdValMaxUlpDiff;
 282     const int         niter = s_nPoints / GMX_SIMD_REAL_WIDTH;
 283
 284     union {
 285         real                                                                               r;
 286         std::conditional<sizeof(real) == sizeof(double), std::int64_t, std::int32_t>::type i;
 287     } conv0, conv1;
 288
 289     // Allow zero-size intervals - nothing to test means we succeeded at it
 290     if (compareSettings.range.first == compareSettings.range.second)
 291     {
 292         ::testing::AssertionSuccess();
 293     }
 294
 295     // Calculate the tolerance limit to use for denormals - we want
 296     // values that are within the ulp tolerance of denormals to be considered matching
 297     conv0.r = std::numeric_limits<real>::min();
 298     conv0.i += compareSettings.ulpTol - 1; // min() itself is not denormal, but one ulp larger
 299     const real denormalLimit = conv0.r;
 300
 301     // We want to test as many diverse bit combinations as possible over the range requested,
 302     // and in particular do it evenly spaced in bit-space.
 303     // Due to the way IEEE754 floating-point is represented, that means we should have a
 304     // constant multiplicative factor between adjacent values. This gets a bit complicated
 305     // when we have both positive and negative values, so we offload the generation of the
 306     // specific testing values to a separate routine
 307     std::vector<real> testPoints = generateTestPoints(compareSettings.range, s_nPoints);
 308
 309     size_t pointIndex = 0;
 310
 311     for (int iter = 0; iter < niter; iter++)
 312     {
 313         for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 314         {
 315             vx[i]   = testPoints[pointIndex];
 316             vref[i] = refFunc(vx[i]);
 317             // If we reach the end of the points, stop increasing index so we pad with
 318             // extra copies of the last element up to the SIMD width
 319             if (pointIndex + 1 < testPoints.size())
 320             {
 321                 pointIndex++;
 322             }
 323         }
 324         vtst = simdReal2Vector(simdFunc(vector2SimdReal(vx)));
 325
 326         bool absOk = true, signOk = true;
 327         for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 328         {
 329             if (compareSettings.matchRule == MatchRule::Dtz && std::abs(vref[i]) <= denormalLimit
 330                 && std::abs(vtst[i]) <= denormalLimit)
 331             {
 332                 continue;
 333             }
 334
 335             if (compareSettings.matchRule == MatchRule::ReferenceOrZero && vtst[i] == 0.0)
 336             {
 337                 // If we accept 0.0 for the test function, we can continue to the next loop iteration.
 338                 continue;
 339             }
 340
 341             absDiff = std::abs(vref[i] - vtst[i]);
 342             absOk   = absOk && (absDiff < compareSettings.absTol);
 343             signOk  = signOk && ((vref[i] >= 0 && vtst[i] >= 0) || (vref[i] <= 0 && vtst[i] <= 0));
 344
 345             if (absDiff >= compareSettings.absTol)
 346             {
 347                 /* We replicate the trivial ulp differences comparison here rather than
 348                  * calling the lower-level routine for comparing them, since this enables
 349                  * us to run through the entire test range and report the largest deviation
 350                  * without lots of extra glue routines.
 351                  */
 352                 conv0.r = vref[i];
 353                 conv1.r = vtst[i];
 354                 ulpDiff = llabs(conv0.i - conv1.i);
 355                 if (ulpDiff > maxUlpDiff)
 356                 {
 357                     maxUlpDiff        = ulpDiff;
 358                     maxUlpDiffPos     = vx[i];
 359                     refValMaxUlpDiff  = vref[i];
 360                     simdValMaxUlpDiff = vtst[i];
 361                 }
 362             }
 363         }
 364         if ((!absOk) && (!signOk))
 365         {
 366             return ::testing::AssertionFailure()
 367                    << "Failing SIMD math function comparison due to sign differences." << std::endl
 368                    << "Reference function: " << refFuncExpr << std::endl
 369                    << "Simd function:      " << simdFuncExpr << std::endl
 370                    << "Test range is ( " << compareSettings.range.first << " , "
 371                    << compareSettings.range.second << " ) " << std::endl
 372                    << "Match rule: " << matchRuleNames_.at(compareSettings.matchRule) << std::endl
 373                    << "First sign difference around x=" << std::setprecision(20)
 374                    << ::testing::PrintToString(vx) << std::endl
 375                    << "Ref values:  " << std::setprecision(20) << ::testing::PrintToString(vref)
 376                    << std::endl
 377                    << "SIMD values: " << std::setprecision(20) << ::testing::PrintToString(vtst)
 378                    << std::endl;
 379         }
 380     }
 381
 382     GMX_RELEASE_ASSERT(compareSettings.ulpTol >= 0, "Invalid ulp value.");
 383     if (maxUlpDiff <= compareSettings.ulpTol)
 384     {
 385         return ::testing::AssertionSuccess();
 386     }
 387     else
 388     {
 389         return ::testing::AssertionFailure()
 390                << "Failing SIMD math function ulp comparison between " << refFuncExpr << " and "
 391                << simdFuncExpr << std::endl
 392                << "Requested ulp tolerance: " << compareSettings.ulpTol << std::endl
 393                << "Requested abs tolerance: " << compareSettings.absTol << std::endl
 394                << "Match rule: " << matchRuleNames_.at(compareSettings.matchRule) << std::endl
 395                << "Largest Ulp difference occurs for x=" << std::setprecision(20) << maxUlpDiffPos
 396                << std::endl
 397                << "Ref  values: " << std::setprecision(20) << refValMaxUlpDiff << std::endl
 398                << "SIMD values: " << std::setprecision(20) << simdValMaxUlpDiff << std::endl
 399                << "Ulp diff.:   " << std::setprecision(20) << maxUlpDiff << std::endl;
 400     }
 401 }
 402
 403 // Actual routine to generate a small set of test points in current precision. This will
 404 // be called by either the double or single precision test fixture, since we need different
 405 // test names to compare to the right reference data.
 406 void SimdMathTest::generateTestPointsTest()
 407 {
 408     int                             points(10);
 409     gmx::test::TestReferenceData    data;
 410     gmx::test::TestReferenceChecker checker(data.rootChecker());
 411
 412     std::vector<real> result;
 413
 414     result = generateTestPoints(Range(-1e10, -1), points);
 415     checker.checkSequence(result.begin(), result.end(), "Test points for interval [-1e10,-1[");
 416
 417     result = generateTestPoints(Range(-1e10, -1e-10), points);
 418     checker.checkSequence(result.begin(), result.end(), "Test points for interval [-1e10, -1e-10[");
 419
 420     result = generateTestPoints(Range(1, 1e10), points);
 421     checker.checkSequence(result.begin(), result.end(), "Test points for interval [1, 1e10[");
 422
 423     result = generateTestPoints(Range(1e-10, 1e10), points);
 424     checker.checkSequence(result.begin(), result.end(), "Test points for interval [1e-10, 1e10[");
 425
 426     result = generateTestPoints(Range(-1e10, 1e-10), points);
 427     checker.checkSequence(result.begin(), result.end(), "Test points for interval [-1e10, 1e-10[");
 428
 429     result = generateTestPoints(Range(-1e-10, 1e-10), points);
 430     checker.checkSequence(result.begin(), result.end(), "Test points for interval [-1e-10, 1e-10[");
 431
 432     result = generateTestPoints(Range(-1e-10, 1e10), points);
 433     checker.checkSequence(result.begin(), result.end(), "Test points for interval [-1e-10, 1e10[");
 434
 435     result = generateTestPoints(Range(-1e10, 1e10), points);
 436     checker.checkSequence(result.begin(), result.end(), "Test points for interval [-1e10, 1e10[");
 437
 438     result = generateTestPoints(Range(-1000, 0), points);
 439     checker.checkSequence(result.begin(), result.end(), "Test points for interval [-1000, 0[");
 440
 441     result = generateTestPoints(Range(0, 1000), points);
 442     checker.checkSequence(result.begin(), result.end(), "Test points for interval [0, 1000[");
 443 }
 444
 445 /*! \} */
 446 /*! \endcond */
 447
 448
 449 // Actual math function tests below
 450
 451
 452 namespace
 453 {
 454
 455 /*! \cond internal */
 456 /*! \addtogroup module_simd */
 457 /*! \{ */
 458
 459 // Reference data is selected based on test name, so make the test name precision-dependent
 460 #        if GMX_DOUBLE
 461 TEST_F(SimdMathTest, generateTestPointsDouble)
 462 {
 463     generateTestPointsTest();
 464 }
 465 #        else
 466 TEST_F(SimdMathTest, generateTestPointsFloat)
 467 {
 468     generateTestPointsTest();
 469 }
 470 #        endif
 471
 472 TEST_F(SimdMathTest, copysign)
 473 {
 474     GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(-c0, c1, c2),
 475                             copysign(setSimdRealFrom3R(c0, c1, c2), setSimdRealFrom3R(-c3, c4, 0)));
 476     GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(-c0, c1, c2),
 477                             copysign(setSimdRealFrom3R(-c0, -c1, -c2), setSimdRealFrom3R(-c3, c4, 0)));
 478     GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(c0, -c1, c2),
 479                             copysign(setSimdRealFrom3R(c0, c1, c2), setSimdRealFrom3R(c3, -c4, 0)));
 480     GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(c0, -c1, c2),
 481                             copysign(setSimdRealFrom3R(-c0, -c1, -c2), setSimdRealFrom3R(c3, -c4, 0)));
 482 }
 483
 484 /*! \brief Function wrapper to evaluate reference 1/sqrt(x) */
 485 real refInvsqrt(real x)
 486 {
 487     return 1.0 / std::sqrt(x);
 488 }
 489
 490 TEST_F(SimdMathTest, invsqrt)
 491 {
 492     const real      low  = std::numeric_limits<float>::min();
 493     const real      high = std::numeric_limits<float>::max();
 494     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
 495
 496     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, invsqrt, settings);
 497 }
 498
 499 TEST_F(SimdMathTest, maskzInvsqrt)
 500 {
 501     SimdReal x   = setSimdRealFrom3R(c1, 0.0, c2);
 502     SimdBool m   = (setZero() < x);
 503     SimdReal ref = setSimdRealFrom3R(1.0 / std::sqrt(c1), 0.0, 1.0 / std::sqrt(c2));
 504     GMX_EXPECT_SIMD_REAL_NEAR(ref, maskzInvsqrt(x, m));
 505 }
 506
 507 /*! \brief Function wrapper to return first result when testing \ref invsqrtPair */
 508 SimdReal gmx_simdcall tstInvsqrtPair0(SimdReal x)
 509 {
 510     SimdReal r0, r1;
 511     invsqrtPair(x, x, &r0, &r1);
 512     return r0;
 513 }
 514
 515 /*! \brief Function wrapper to return second result when testing \ref invsqrtPair */
 516 SimdReal gmx_simdcall tstInvsqrtPair1(SimdReal x)
 517 {
 518     SimdReal r0, r1;
 519     invsqrtPair(x, x, &r0, &r1);
 520     return r1;
 521 }
 522
 523 TEST_F(SimdMathTest, invsqrtPair)
 524 {
 525     const real low  = std::numeric_limits<float>::min();
 526     const real high = std::numeric_limits<float>::max();
 527
 528     // Accuracy conversions lose a bit of accuracy compared to all-double,
 529     // so increase the tolerance to 4*ulpTol_
 530     CompareSettings settings{ Range(low, high), 4 * ulpTol_, absTol_, MatchRule::Normal };
 531
 532     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, tstInvsqrtPair0, settings);
 533     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, tstInvsqrtPair1, settings);
 534 }
 535
 536 /*! \brief Function wrapper to evaluate reference sqrt(x) */
 537 real refSqrt(real x)
 538 {
 539     return std::sqrt(x);
 540 }
 541
 542 TEST_F(SimdMathTest, sqrt)
 543 {
 544     // Since the first lookup step is sometimes performed in single precision,
 545     // our SIMD sqrt can only handle single-precision input values, even when
 546     // compiled in double precision.
 547
 548     const real      minFloat     = std::numeric_limits<float>::min();
 549     const real      minSafeFloat = minFloat * 10;
 550     const real      maxSafeFloat = std::numeric_limits<float>::max() * 0.1;
 551     CompareSettings settings;
 552     // The accuracy conversions lose a bit of extra accuracy compared to
 553     // doing the iterations in all-double.
 554     setUlpTol(4 * ulpTol_);
 555
 556     // First test that 0.0 and a few other values work
 557     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(0, std::sqrt(c1), std::sqrt(c2)),
 558                               sqrt(setSimdRealFrom3R(0, c1, c2)));
 559
 560 #        if GMX_DOUBLE
 561     // As mentioned above, we cannot guarantee that very small double precision
 562     // input values (below std::numeric_limits<float>::min()) are handled correctly,
 563     // so our implementation will clamp it to zero. In this range we allow either
 564     // the correct value or zero, but it's important that it does not result in NaN or Inf values.
 565     //
 566     // This test range must not be called for single precision, since if we try to divide
 567     // the interval (0.0, low( in npoints we will try to multiply by factors so small that
 568     // they end up being flushed to zero, and the loop would never end.
 569     settings = { Range(0.0, minFloat), ulpTol_, absTol_, MatchRule::ReferenceOrZero };
 570     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrt, settings);
 571 #        endif
 572
 573     // Next range: Just about minFloat the lookup should always work, but the results
 574     // might be a bit fragile due to issues with the N-R iterations being flushed to zero
 575     // for denormals. We can probably relax the latter in double precision, but since we
 576     // anyway cannot handle numbers that cannot be represented in single it's not worth
 577     // worrying too much about whether we have zero or an exact values around 10^-38....
 578     settings = { Range(minFloat, minSafeFloat), ulpTol_, absTol_, MatchRule::ReferenceOrZero };
 579     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrt, settings);
 580
 581     settings = { Range(minSafeFloat, maxSafeFloat), ulpTol_, absTol_, MatchRule::Normal };
 582     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrt, settings);
 583 }
 584
 585 TEST_F(SimdMathTest, sqrtUnsafe)
 586 {
 587     const real minSafeFloat = std::numeric_limits<float>::min() * 10;
 588     const real maxSafeFloat = std::numeric_limits<float>::max() * 0.1;
 589
 590     // The accuracy conversions lose a bit of extra accuracy compared to
 591     // doing the iterations in all-double, so we use 4*ulpTol_
 592     setUlpTol(4 * ulpTol_);
 593
 594     CompareSettings settings{ Range(minSafeFloat, maxSafeFloat), 4 * ulpTol_, absTol_, MatchRule::Normal };
 595     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrt<MathOptimization::Unsafe>, settings);
 596 }
 597
 598 /*! \brief Function wrapper to evaluate reference 1/x */
 599 real refInv(real x)
 600 {
 601     return 1.0 / x;
 602 }
 603
 604 TEST_F(SimdMathTest, inv)
 605 {
 606     // Since the first lookup step is sometimes performed in single precision,
 607     // our SIMD 1/x can only handle single-precision input values, even when
 608     // compiled in double precision.
 609
 610     // Relevant threshold points
 611     const real minSafeFloat = std::numeric_limits<float>::min()
 612                               * 10; // X value guaranteed not to result in Inf intermediates for 1/x calc.
 613     const real maxSafeFloat = std::numeric_limits<float>::max()
 614                               * 0.1; // X value guaranteed not to result in DTZ intermediates for 1/x calc.
 615     // Scale highest value by 1-eps, since we will do some arithmetics on this value
 616     const real maxFloat =
 617             std::numeric_limits<float>::max() * (1.0 - std::numeric_limits<float>::epsilon());
 618     CompareSettings settings;
 619
 620     // Danger zone where intermediates might be flushed to zero and produce 1/x==0.0
 621     settings = { Range(-maxFloat, -maxSafeFloat), ulpTol_, absTol_, MatchRule::ReferenceOrZero };
 622     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv, settings);
 623
 624     // Normal checks for x < 0
 625     settings = { Range(-maxSafeFloat, -minSafeFloat), ulpTol_, absTol_, MatchRule::Normal };
 626     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv, settings);
 627
 628     // We do not care about the small range -minSafeFloat < x < +minSafeFloat where the result can be +/- Inf, since we don't require strict IEEE754.
 629
 630     // Normal checks for x > 0
 631     settings = { Range(minSafeFloat, maxSafeFloat), ulpTol_, absTol_, MatchRule::Normal };
 632     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv, settings);
 633
 634     // Danger zone where intermediates might be flushed to zero and produce 1/x==0.0
 635     settings = { Range(maxSafeFloat, maxFloat), ulpTol_, absTol_, MatchRule::ReferenceOrZero };
 636     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv, settings);
 637 }
 638
 639 TEST_F(SimdMathTest, maskzInv)
 640 {
 641     SimdReal x   = setSimdRealFrom3R(c1, 0.0, c2);
 642     SimdBool m   = (setZero() < x);
 643     SimdReal ref = setSimdRealFrom3R(1.0 / c1, 0.0, 1.0 / c2);
 644     GMX_EXPECT_SIMD_REAL_NEAR(ref, maskzInv(x, m));
 645 }
 646
 647 TEST_F(SimdMathTest, cbrt)
 648 {
 649     const real low  = -std::numeric_limits<real>::max();
 650     const real high = std::numeric_limits<real>::max();
 651
 652     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
 653     GMX_EXPECT_SIMD_FUNC_NEAR(std::cbrt, cbrt, settings);
 654 }
 655
 656 /*! \brief Function wrapper to evaluate reference 1/cbrt(x) */
 657 real refInvCbrt(real x)
 658 {
 659     return 1.0 / std::cbrt(x);
 660 }
 661
 662 TEST_F(SimdMathTest, invcbrt)
 663 {
 664     // Negative values first
 665     real low  = -std::numeric_limits<real>::max();
 666     real high = -std::numeric_limits<real>::min();
 667
 668     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
 669     GMX_EXPECT_SIMD_FUNC_NEAR(refInvCbrt, invcbrt, settings);
 670
 671     // Positive values
 672     low      = std::numeric_limits<real>::min();
 673     high     = std::numeric_limits<real>::max();
 674     settings = { Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
 675     GMX_EXPECT_SIMD_FUNC_NEAR(refInvCbrt, invcbrt, settings);
 676 }
 677
 678 TEST_F(SimdMathTest, log2)
 679 {
 680     const real low  = std::numeric_limits<real>::min();
 681     const real high = std::numeric_limits<real>::max();
 682
 683     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
 684     GMX_EXPECT_SIMD_FUNC_NEAR(std::log2, log2, settings);
 685 }
 686
 687 TEST_F(SimdMathTest, log)
 688 {
 689     const real low  = std::numeric_limits<real>::min();
 690     const real high = std::numeric_limits<real>::max();
 691
 692     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
 693     GMX_EXPECT_SIMD_FUNC_NEAR(std::log, log, settings);
 694 }
 695
 696 TEST_F(SimdMathTest, exp2)
 697 {
 698     // Relevant threshold points
 699     constexpr real lowestReal = -std::numeric_limits<real>::max();
 700     constexpr real lowestRealThatProducesNormal =
 701             std::numeric_limits<real>::min_exponent
 702             - 1; // adding the significant corresponds to one more unit in exponent
 703     constexpr real lowestRealThatProducesDenormal =
 704             lowestRealThatProducesNormal
 705             - std::numeric_limits<real>::digits; // digits refer to bits in significand, so 24/53 for float/double
 706     constexpr real highestRealThatProducesNormal =
 707             std::numeric_limits<real>::max_exponent
 708             - 1; // adding the significant corresponds to one more unit in exponent
 709     CompareSettings settings;
 710
 711     // Below subnormal range all results should be zero (so, match the reference)
 712     settings = { Range(lowestReal, lowestRealThatProducesDenormal), ulpTol_, absTol_, MatchRule::Normal };
 713     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2, settings);
 714
 715     // Subnormal range, require matching, but DTZ is fine
 716     settings = { Range(lowestRealThatProducesDenormal, lowestRealThatProducesNormal), ulpTol_,
 717                  absTol_, MatchRule::Dtz };
 718     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2, settings);
 719
 720     // Normal range, standard result expected
 721     settings = { Range(lowestRealThatProducesNormal, highestRealThatProducesNormal), ulpTol_,
 722                  absTol_, MatchRule::Normal };
 723     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2, settings);
 724 }
 725
 726 TEST_F(SimdMathTest, exp2Unsafe)
 727 {
 728     // The unsafe version is only defined in the normal range
 729     constexpr real lowestRealThatProducesNormal =
 730             std::numeric_limits<real>::min_exponent
 731             - 1; // adding the significant corresponds to one more unit in exponent
 732     constexpr real highestRealThatProducesNormal =
 733             std::numeric_limits<real>::max_exponent
 734             - 1; // adding the significant corresponds to one more unit in exponent
 735
 736     CompareSettings settings{ Range(lowestRealThatProducesNormal, highestRealThatProducesNormal),
 737                               ulpTol_, absTol_, MatchRule::Normal };
 738     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2<MathOptimization::Unsafe>, settings);
 739 }
 740
 741 TEST_F(SimdMathTest, exp)
 742 {
 743     // Relevant threshold points. See the exp2 test for more details about the values; these are
 744     // simply scaled by log(2) due to the difference between exp2 and exp.
 745     const real lowestReal = -std::numeric_limits<real>::max();
 746     // In theory the smallest value should be (min_exponent-1)*log(2), but rounding after the multiplication will cause this
 747     // value to be a single ulp too low. This might cause failed tests on CPUs that use different DTZ modes for SIMD vs.
 748     // non-SIMD arithmetics (ARM v7), so multiply by (1.0-eps) to increase it by a single ulp.
 749     const real lowestRealThatProducesNormal = (std::numeric_limits<real>::min_exponent - 1)
 750                                               * std::log(2.0)
 751                                               * (1 - std::numeric_limits<real>::epsilon());
 752     const real lowestRealThatProducesDenormal =
 753             lowestRealThatProducesNormal - std::numeric_limits<real>::digits * std::log(2.0);
 754     const real highestRealThatProducesNormal =
 755             (std::numeric_limits<real>::max_exponent - 1) * std::log(2.0);
 756     CompareSettings settings;
 757
 758     // Below subnormal range all results should be zero (so, match the reference)
 759     settings = { Range(lowestReal, lowestRealThatProducesDenormal), ulpTol_, absTol_, MatchRule::Normal };
 760     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, exp, settings);
 761
 762     // Subnormal range, require matching, but DTZ is fine
 763     settings = { Range(lowestRealThatProducesDenormal, lowestRealThatProducesNormal), ulpTol_,
 764                  absTol_, MatchRule::Dtz };
 765     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, exp, settings);
 766
 767     // Normal range, standard result expected
 768     settings = { Range(lowestRealThatProducesNormal, highestRealThatProducesNormal), ulpTol_,
 769                  absTol_, MatchRule::Normal };
 770     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, exp, settings);
 771 }
 772
 773 TEST_F(SimdMathTest, expUnsafe)
 774 {
 775     // See test of exp() for comments about test ranges
 776     const real lowestRealThatProducesNormal = (std::numeric_limits<real>::min_exponent - 1)
 777                                               * std::log(2.0)
 778                                               * (1 - std::numeric_limits<real>::epsilon());
 779     const real highestRealThatProducesNormal =
 780             (std::numeric_limits<real>::max_exponent - 1) * std::log(2.0);
 781
 782     CompareSettings settings{ Range(lowestRealThatProducesNormal, highestRealThatProducesNormal),
 783                               ulpTol_, absTol_, MatchRule::Normal };
 784     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, exp<MathOptimization::Unsafe>, settings);
 785 }
 786
 787 TEST_F(SimdMathTest, pow)
 788 {
 789     // We already test the log2/exp2 components of pow() extensively above, and it's a very
 790     // simple single-line function, so here we just test a handful of values to catch typos
 791     // and then some special values.
 792
 793     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::pow(c0, c3), std::pow(c1, c4), std::pow(c2, c5)),
 794                               pow(rSimd_c0c1c2, rSimd_c3c4c5));
 795
 796     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::pow(c0, -c3), std::pow(c1, -c0), std::pow(c2, -c4)),
 797                               pow(rSimd_c0c1c2, rSimd_m3m0m4));
 798
 799     // 0^0 = 1 , 0^c1=0, -c1^0=1
 800     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(1.0, 0.0, 1.0),
 801                               pow(setSimdRealFrom3R(0, 0.0, -c1), setSimdRealFrom3R(0.0, c1, 0.0)));
 802 }
 803
 804 TEST_F(SimdMathTest, powUnsafe)
 805 {
 806     // We already test the log2/exp2 components of pow() extensively above, and it's a very
 807     // simple single-line function, so here we just test a handful of values to catch typos
 808     // and then some special values.
 809
 810     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::pow(c0, c3), std::pow(c1, c4), std::pow(c2, c5)),
 811                               pow<MathOptimization::Unsafe>(rSimd_c0c1c2, rSimd_c3c4c5));
 812
 813     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::pow(c0, -c3), std::pow(c1, -c0), std::pow(c2, -c4)),
 814                               pow<MathOptimization::Unsafe>(rSimd_c0c1c2, rSimd_m3m0m4));
 815 }
 816
 817 /*! \brief Function wrapper for erf(x), with argument/return in default Gromacs precision.
 818  *
 819  * \note Single-precision erf() in some libraries can be slightly lower precision
 820  * than the SIMD flavor, so we use a cast to force double precision for reference.
 821  */
 822 real refErf(real x)
 823 {
 824     return std::erf(static_cast<double>(x));
 825 }
 826
 827 TEST_F(SimdMathTest, erf)
 828 {
 829     CompareSettings settings{ Range(-9, 9), ulpTol_, std::numeric_limits<real>::min(), MatchRule::Normal };
 830     GMX_EXPECT_SIMD_FUNC_NEAR(refErf, erf, settings);
 831 }
 832
 833 /*! \brief Function wrapper for erfc(x), with argument/return in default Gromacs precision.
 834  *
 835  * \note Single-precision erfc() in some libraries can be slightly lower precision
 836  * than the SIMD flavor, so we use a cast to force double precision for reference.
 837  */
 838 real refErfc(real x)
 839 {
 840     return std::erfc(static_cast<double>(x));
 841 }
 842
 843 TEST_F(SimdMathTest, erfc)
 844 {
 845     // Our erfc algorithm has 4 ulp accuracy, so relax tolerance a bit to 4*ulpTol
 846     CompareSettings settings{ Range(-9, 9), 4 * ulpTol_, std::numeric_limits<real>::min(),
 847                               MatchRule::Normal };
 848     GMX_EXPECT_SIMD_FUNC_NEAR(refErfc, erfc, settings);
 849 }
 850
 851 TEST_F(SimdMathTest, sin)
 852 {
 853     CompareSettings settings{ Range(-8 * M_PI, 8 * M_PI), ulpTol_, absTol_, MatchRule::Normal };
 854     GMX_EXPECT_SIMD_FUNC_NEAR(std::sin, sin, settings);
 855
 856     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 857     settings = { Range(-10000, 10000), 2 * ulpTol_, absTol_, MatchRule::Normal };
 858     GMX_EXPECT_SIMD_FUNC_NEAR(std::sin, sin, settings);
 859 }
 860
 861 TEST_F(SimdMathTest, cos)
 862 {
 863     CompareSettings settings{ Range(-8 * M_PI, 8 * M_PI), ulpTol_, absTol_, MatchRule::Normal };
 864     GMX_EXPECT_SIMD_FUNC_NEAR(std::cos, cos, settings);
 865
 866     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 867     settings = { Range(-10000, 10000), 2 * ulpTol_, absTol_, MatchRule::Normal };
 868     GMX_EXPECT_SIMD_FUNC_NEAR(std::cos, cos, settings);
 869 }
 870
 871 TEST_F(SimdMathTest, tan)
 872 {
 873     // Tan(x) is a little sensitive due to the division in the algorithm.
 874     // Rather than using lots of extra FP operations, we accept the algorithm
 875     // presently only achieves a ~3 ulp error and use the medium tolerance.
 876     CompareSettings settings{ Range(-8 * M_PI, 8 * M_PI), ulpTol_, absTol_, MatchRule::Normal };
 877     GMX_EXPECT_SIMD_FUNC_NEAR(std::tan, tan, settings);
 878
 879     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 880     settings = { Range(-10000, 10000), 2 * ulpTol_, absTol_, MatchRule::Normal };
 881     GMX_EXPECT_SIMD_FUNC_NEAR(std::tan, tan, settings);
 882 }
 883
 884 TEST_F(SimdMathTest, asin)
 885 {
 886     // Our present asin(x) algorithm achieves 2-3 ulp accuracy
 887     CompareSettings settings{ Range(-1, 1), ulpTol_, absTol_, MatchRule::Normal };
 888     GMX_EXPECT_SIMD_FUNC_NEAR(std::asin, asin, settings);
 889 }
 890
 891 TEST_F(SimdMathTest, acos)
 892 {
 893     // Our present acos(x) algorithm achieves 2-3 ulp accuracy
 894     CompareSettings settings{ Range(-1, 1), ulpTol_, absTol_, MatchRule::Normal };
 895     GMX_EXPECT_SIMD_FUNC_NEAR(std::acos, acos, settings);
 896 }
 897
 898 TEST_F(SimdMathTest, atan)
 899 {
 900     // Our present atan(x) algorithm achieves 1 ulp accuracy
 901     CompareSettings settings{ Range(-10000, 10000), ulpTol_, absTol_, MatchRule::Normal };
 902     GMX_EXPECT_SIMD_FUNC_NEAR(std::atan, atan, settings);
 903 }
 904
 905 TEST_F(SimdMathTest, atan2)
 906 {
 907     // test each quadrant
 908     GMX_EXPECT_SIMD_REAL_NEAR(
 909             setSimdRealFrom3R(std::atan2(c0, c3), std::atan2(c1, c4), std::atan2(c2, c5)),
 910             atan2(rSimd_c0c1c2, rSimd_c3c4c5));
 911     GMX_EXPECT_SIMD_REAL_NEAR(
 912             setSimdRealFrom3R(std::atan2(-c0, c3), std::atan2(-c1, c4), std::atan2(-c2, c5)),
 913             atan2(rSimd_m0m1m2, rSimd_c3c4c5));
 914     GMX_EXPECT_SIMD_REAL_NEAR(
 915             setSimdRealFrom3R(std::atan2(-c0, -c3), std::atan2(-c1, -c0), std::atan2(-c2, -c4)),
 916             atan2(rSimd_m0m1m2, rSimd_m3m0m4));
 917     GMX_EXPECT_SIMD_REAL_NEAR(
 918             setSimdRealFrom3R(std::atan2(c0, -c3), std::atan2(c1, -c0), std::atan2(c2, -c4)),
 919             atan2(rSimd_c0c1c2, rSimd_m3m0m4));
 920
 921     // cases important for calculating angles
 922     // values on coordinate axes
 923     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(0, c0), std::atan2(0, c1), std::atan2(0, c2)),
 924                               atan2(setZero(), rSimd_c0c1c2));
 925     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(c0, 0), std::atan2(c1, 0), std::atan2(c2, 0)),
 926                               atan2(rSimd_c0c1c2, setZero()));
 927     GMX_EXPECT_SIMD_REAL_NEAR(
 928             setSimdRealFrom3R(std::atan2(0, -c0), std::atan2(0, -c1), std::atan2(0, -c2)),
 929             atan2(setZero(), rSimd_m0m1m2));
 930     GMX_EXPECT_SIMD_REAL_NEAR(
 931             setSimdRealFrom3R(std::atan2(-c0, 0), std::atan2(-c1, 0), std::atan2(-c2, 0)),
 932             atan2(rSimd_m0m1m2, setZero()));
 933     // degenerate value (origin) should return 0.0. At least IBM xlc 13.1.5 gets the reference
 934     // value wrong (-nan) at -O3 optimization, so we compare to the correct value (0.0) instead.
 935     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(0.0), atan2(setSimdRealFrom3R(0.0, 0.0, 0.0), setZero()));
 936 }
 937
 938 /*! \brief Evaluate reference version of PME force correction. */
 939 real refPmeForceCorrection(real x)
 940 {
 941     if (x != 0)
 942     {
 943         real y = std::sqrt(x);
 944         return 2 * std::exp(-x) / (std::sqrt(M_PI) * x) - std::erf(static_cast<double>(y)) / (x * y);
 945     }
 946     else
 947     {
 948         return -4 / (3 * std::sqrt(M_PI));
 949     }
 950 }
 951
 952 // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
 953 TEST_F(SimdMathTest, pmeForceCorrection)
 954 {
 955     // Pme correction relative accuracy only needs to be ~1e-6 accuracy single, 1e-10 double
 956     const std::int64_t ulpTol = (GMX_DOUBLE ? 5e-10 : 5e-6) / GMX_REAL_EPS;
 957
 958     CompareSettings settings{ Range(0.15, 4), ulpTol, GMX_REAL_EPS, MatchRule::Normal };
 959     GMX_EXPECT_SIMD_FUNC_NEAR(refPmeForceCorrection, pmeForceCorrection, settings);
 960 }
 961
 962 /*! \brief Evaluate reference version of PME potential correction. */
 963 real refPmePotentialCorrection(real x)
 964 {
 965     real y = std::sqrt(x);
 966     return std::erf(static_cast<double>(y)) / y;
 967 }
 968
 969 // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
 970 TEST_F(SimdMathTest, pmePotentialCorrection)
 971 {
 972     // Pme correction relative accuracy only needs to be ~1e-6 accuracy single, 1e-10 double
 973     const std::int64_t ulpTol = (GMX_DOUBLE ? 5e-10 : 5e-6) / GMX_REAL_EPS;
 974
 975     CompareSettings settings{ Range(0.15, 4), ulpTol, GMX_REAL_EPS, MatchRule::Normal };
 976     GMX_EXPECT_SIMD_FUNC_NEAR(refPmePotentialCorrection, pmePotentialCorrection, settings);
 977 }
 978
 979 // Functions that only target single accuracy, even for double SIMD data
 980
 981 TEST_F(SimdMathTest, invsqrtSingleAccuracy)
 982 {
 983     // Here we always use float limits, since the lookup is not defined for numbers that
 984     // cannot be represented in single precision.
 985     const real low  = std::numeric_limits<float>::min();
 986     const real high = std::numeric_limits<float>::max();
 987     /* Increase the allowed error by the difference between the actual precision and single */
 988     setUlpTolSingleAccuracy(ulpTol_);
 989
 990     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
 991     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, invsqrtSingleAccuracy, settings);
 992 }
 993
 994 /*! \brief Function wrapper to return first result when testing \ref invsqrtPairSingleAccuracy */
 995 SimdReal gmx_simdcall tst_invsqrt_SingleAccuracy_pair0(SimdReal x)
 996 {
 997     SimdReal r0, r1;
 998     invsqrtPairSingleAccuracy(x, x, &r0, &r1);
 999     return r0;
1000 }
1001
1002 /*! \brief Function wrapper to return second result when testing \ref invsqrtPairSingleAccuracy */
1003 SimdReal gmx_simdcall tst_invsqrt_SingleAccuracy_pair1(SimdReal x)
1004 {
1005     SimdReal r0, r1;
1006     invsqrtPairSingleAccuracy(x, x, &r0, &r1);
1007     return r1;
1008 }
1009
1010 TEST_F(SimdMathTest, invsqrtPairSingleAccuracy)
1011 {
1012     // Float limits since lookup is always performed in single
1013     const real low  = std::numeric_limits<float>::min();
1014     const real high = std::numeric_limits<float>::max();
1015     /* Increase the allowed error by the difference between the actual precision and single */
1016     setUlpTolSingleAccuracy(ulpTol_);
1017
1018     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
1019     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, tst_invsqrt_SingleAccuracy_pair0, settings);
1020     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, tst_invsqrt_SingleAccuracy_pair1, settings);
1021 }
1022
1023 TEST_F(SimdMathTest, sqrtSingleAccuracy)
1024 {
1025     // Since the first lookup step is sometimes performed in single precision,
1026     // our SIMD sqrt can only handle single-precision input values, even when
1027     // compiled in double precision - thus we use single precision limits here.
1028
1029     // Scale lowest value by 1+eps, since we will do some arithmetics on this value
1030     const real low = std::numeric_limits<float>::min() * (1.0 + std::numeric_limits<float>::epsilon());
1031     const real      high = std::numeric_limits<float>::max();
1032     CompareSettings settings;
1033
1034     // Increase the allowed error by the difference between the actual precision and single
1035     setUlpTolSingleAccuracy(ulpTol_);
1036
1037     // First test that 0.0 and a few other values works
1038     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(0, std::sqrt(c0), std::sqrt(c1)),
1039                               sqrtSingleAccuracy(setSimdRealFrom3R(0, c0, c1)));
1040
1041 #        if GMX_DOUBLE
1042     // As mentioned above, we cannot guarantee that very small double precision
1043     // input values (below std::numeric_limits<float>::min()) are handled correctly,
1044     // so our implementation will clamp it to zero. In this range we allow either
1045     // the correct value or zero, but it's important that it does not result in NaN or Inf values.
1046     //
1047     // This test range must not be called for single precision, since if we try to divide
1048     // the interval (0.0, low( in npoints we will try to multiply by factors so small that
1049     // they end up being flushed to zero, and the loop would never end.
1050     settings = { Range(0.0, low), ulpTol_, absTol_, MatchRule::ReferenceOrZero };
1051     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrtSingleAccuracy, settings);
1052 #        endif
1053
1054     settings = { Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
1055     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrtSingleAccuracy, settings);
1056 }
1057
1058 TEST_F(SimdMathTest, sqrtSingleAccuracyUnsafe)
1059 {
1060     // Test the full range, but stick to float limits since lookup is done in single.
1061     const real low  = std::numeric_limits<float>::min();
1062     const real high = std::numeric_limits<float>::max();
1063
1064     /* Increase the allowed error by the difference between the actual precision and single */
1065     setUlpTolSingleAccuracy(ulpTol_);
1066
1067     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
1068     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrtSingleAccuracy<MathOptimization::Unsafe>, settings);
1069 }
1070
1071 TEST_F(SimdMathTest, invSingleAccuracy)
1072 {
1073     // Since the first lookup step is sometimes performed in single precision,
1074     // our SIMD 1/x can only handle single-precision input values, even when
1075     // compiled in double precision.
1076
1077     // Relevant threshold points
1078     const real minSafeFloat = std::numeric_limits<float>::min()
1079                               * 10; // X value guaranteed not to result in Inf intermediates for 1/x calc.
1080     const real maxSafeFloat = std::numeric_limits<float>::max()
1081                               * 0.1; // X value guaranteed not to result in DTZ intermediates for 1/x calc.
1082     // Scale highest value by 1-eps, since we will do some arithmetics on this value
1083     const real maxFloat =
1084             std::numeric_limits<float>::max() * (1.0 - std::numeric_limits<float>::epsilon());
1085     CompareSettings settings;
1086
1087     // Increase the allowed error by the difference between the actual precision and single
1088     setUlpTolSingleAccuracy(ulpTol_);
1089
1090     // Danger zone where intermediates might be flushed to zero and produce 1/x==0.0
1091     settings = { Range(-maxFloat, -maxSafeFloat), ulpTol_, absTol_, MatchRule::ReferenceOrZero };
1092     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv, settings);
1093
1094     // Normal checks for x < 0
1095     settings = { Range(-maxSafeFloat, -minSafeFloat), ulpTol_, absTol_, MatchRule::Normal };
1096     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv, settings);
1097
1098     // We do not care about the small range -minSafeFloat < x < +minSafeFloat where the result can be +/- Inf, since we don't require strict IEEE754.
1099
1100     // Normal checks for x > 0
1101     settings = { Range(minSafeFloat, maxSafeFloat), ulpTol_, absTol_, MatchRule::Normal };
1102     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv, settings);
1103
1104     // Danger zone where intermediates might be flushed to zero and produce 1/x==0.0
1105     settings = { Range(maxSafeFloat, maxFloat), ulpTol_, absTol_, MatchRule::ReferenceOrZero };
1106     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv, settings);
1107 }
1108
1109 TEST_F(SimdMathTest, cbrtSingleAccuracy)
1110 {
1111     const real low  = -std::numeric_limits<real>::max();
1112     const real high = std::numeric_limits<real>::max();
1113
1114     // Increase the allowed error by the difference between the actual precision and single
1115     setUlpTolSingleAccuracy(ulpTol_);
1116
1117     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
1118     GMX_EXPECT_SIMD_FUNC_NEAR(std::cbrt, cbrtSingleAccuracy, settings);
1119 }
1120
1121 TEST_F(SimdMathTest, invcbrtSingleAccuracy)
1122 {
1123     // Increase the allowed error by the difference between the actual precision and single
1124     setUlpTolSingleAccuracy(ulpTol_);
1125
1126     // Negative values first
1127     real low  = -std::numeric_limits<real>::max();
1128     real high = -std::numeric_limits<real>::min();
1129
1130     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
1131     GMX_EXPECT_SIMD_FUNC_NEAR(refInvCbrt, invcbrtSingleAccuracy, settings);
1132
1133     // Positive values
1134     low      = std::numeric_limits<real>::min();
1135     high     = std::numeric_limits<real>::max();
1136     settings = { Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
1137     GMX_EXPECT_SIMD_FUNC_NEAR(refInvCbrt, invcbrtSingleAccuracy, settings);
1138 }
1139
1140 TEST_F(SimdMathTest, log2SingleAccuracy)
1141 {
1142     const real low  = std::numeric_limits<real>::min();
1143     const real high = std::numeric_limits<real>::max();
1144
1145     // Increase the allowed error by the difference between the actual precision and single
1146     setUlpTolSingleAccuracy(ulpTol_);
1147
1148     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
1149     GMX_EXPECT_SIMD_FUNC_NEAR(std::log2, log2SingleAccuracy, settings);
1150 }
1151
1152 TEST_F(SimdMathTest, logSingleAccuracy)
1153 {
1154     const real low  = std::numeric_limits<real>::min();
1155     const real high = std::numeric_limits<real>::max();
1156
1157     // Increase the allowed error by the difference between the actual precision and single
1158     setUlpTolSingleAccuracy(ulpTol_);
1159
1160     CompareSettings settings{ Range(low, high), ulpTol_, absTol_, MatchRule::Normal };
1161     GMX_EXPECT_SIMD_FUNC_NEAR(std::log, logSingleAccuracy, settings);
1162 }
1163
1164 TEST_F(SimdMathTest, exp2SingleAccuracy)
1165 {
1166     // Relevant threshold points - float limits since we only target single accuracy
1167     constexpr real lowestReal = -std::numeric_limits<real>::max();
1168     constexpr real lowestRealThatProducesNormal =
1169             std::numeric_limits<real>::min_exponent
1170             - 1; // adding the significant corresponds to one more unit in exponent
1171     constexpr real lowestRealThatProducesDenormal =
1172             lowestRealThatProducesNormal
1173             - std::numeric_limits<real>::digits; // digits refer to bits in significand, so 24/53 for float/double
1174     constexpr real highestRealThatProducesNormal =
1175             std::numeric_limits<real>::max_exponent
1176             - 1; // adding the significant corresponds to one more unit in exponent
1177     CompareSettings settings;
1178
1179     // Increase the allowed error by the difference between the actual precision and single
1180     setUlpTolSingleAccuracy(ulpTol_);
1181
1182     // Below subnormal range all results should be zero
1183     settings = { Range(lowestReal, lowestRealThatProducesDenormal), ulpTol_, absTol_, MatchRule::Normal };
1184     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2SingleAccuracy, settings);
1185
1186     // Subnormal range, require matching, but DTZ is fine
1187     settings = { Range(lowestRealThatProducesDenormal, lowestRealThatProducesNormal), ulpTol_,
1188                  absTol_, MatchRule::Dtz };
1189     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2SingleAccuracy, settings);
1190
1191     // Normal range, standard result expected
1192     settings = { Range(lowestRealThatProducesNormal, highestRealThatProducesNormal), ulpTol_,
1193                  absTol_, MatchRule::Normal };
1194     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2SingleAccuracy, settings);
1195 }
1196
1197 TEST_F(SimdMathTest, exp2SingleAccuracyUnsafe)
1198 {
1199     // The unsafe version is only defined in the normal range
1200     constexpr real lowestRealThatProducesNormal =
1201             std::numeric_limits<real>::min_exponent
1202             - 1; // adding the significant corresponds to one more unit in exponent
1203     constexpr real highestRealThatProducesNormal =
1204             std::numeric_limits<real>::max_exponent
1205             - 1; // adding the significant corresponds to one more unit in exponent
1206
1207     /* Increase the allowed error by the difference between the actual precision and single */
1208     setUlpTolSingleAccuracy(ulpTol_);
1209
1210     CompareSettings settings{ Range(lowestRealThatProducesNormal, highestRealThatProducesNormal),
1211                               ulpTol_, absTol_, MatchRule::Normal };
1212     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2SingleAccuracy<MathOptimization::Unsafe>, settings);
1213 }
1214
1215 TEST_F(SimdMathTest, expSingleAccuracy)
1216 {
1217     // See threshold point comments in normal exp() test
1218     const real lowestReal = -std::numeric_limits<real>::max();
1219     // In theory the smallest value should be (min_exponent-1)*log(2), but rounding after the multiplication will cause this
1220     // value to be a single ulp too low. This might cause failed tests on CPUs that use different DTZ modes for SIMD vs.
1221     // non-SIMD arithmetics (ARM v7), so multiply by (1.0-eps) to increase it by a single ulp.
1222     const real lowestRealThatProducesNormal = (std::numeric_limits<real>::min_exponent - 1)
1223                                               * std::log(2.0)
1224                                               * (1.0 - std::numeric_limits<real>::epsilon());
1225     const real lowestRealThatProducesDenormal =
1226             lowestRealThatProducesNormal - std::numeric_limits<real>::digits * std::log(2.0);
1227     const real highestRealThatProducesNormal =
1228             (std::numeric_limits<real>::max_exponent - 1) * std::log(2.0);
1229     CompareSettings settings;
1230
1231     // Increase the allowed error by the difference between the actual precision and single
1232     setUlpTolSingleAccuracy(ulpTol_);
1233
1234     // Below subnormal range all results should be zero (so, match the reference)
1235     settings = { Range(lowestReal, lowestRealThatProducesDenormal), ulpTol_, absTol_, MatchRule::Normal };
1236     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, expSingleAccuracy, settings);
1237
1238     // Subnormal range, require matching, but DTZ is fine
1239     settings = { Range(lowestRealThatProducesDenormal, lowestRealThatProducesNormal), ulpTol_,
1240                  absTol_, MatchRule::Dtz };
1241     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, expSingleAccuracy, settings);
1242
1243     // Normal range, standard result expected
1244     settings = { Range(lowestRealThatProducesNormal, highestRealThatProducesNormal), ulpTol_,
1245                  absTol_, MatchRule::Normal };
1246     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, expSingleAccuracy, settings);
1247 }
1248
1249 TEST_F(SimdMathTest, expSingleAccuracyUnsafe)
1250 {
1251     // See test of exp() for comments about test ranges
1252     const real lowestRealThatProducesNormal = (std::numeric_limits<real>::min_exponent - 1)
1253                                               * std::log(2.0)
1254                                               * (1 - std::numeric_limits<real>::epsilon());
1255     const real highestRealThatProducesNormal =
1256             (std::numeric_limits<real>::max_exponent - 1) * std::log(2.0);
1257
1258     // Increase the allowed error by the difference between the actual precision and single
1259     setUlpTolSingleAccuracy(ulpTol_);
1260
1261     CompareSettings settings{ Range(lowestRealThatProducesNormal, highestRealThatProducesNormal),
1262                               ulpTol_, absTol_, MatchRule::Normal };
1263     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, expSingleAccuracy<MathOptimization::Unsafe>, settings);
1264 }
1265
1266 TEST_F(SimdMathTest, powSingleAccuracy)
1267 {
1268     // We already test the log2/exp2 components of pow() extensively above, and it's a very
1269     // simple single-line function, so here we just test a handful of values to catch typos
1270     // and then some special values.
1271
1272     // Increase the allowed error by the difference between the actual precision and single
1273     setUlpTolSingleAccuracy(ulpTol_);
1274
1275     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::pow(c0, c3), std::pow(c1, c4), std::pow(c2, c5)),
1276                               powSingleAccuracy(rSimd_c0c1c2, rSimd_c3c4c5));
1277
1278     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::pow(c0, -c3), std::pow(c1, -c0), std::pow(c2, -c4)),
1279                               powSingleAccuracy(rSimd_c0c1c2, rSimd_m3m0m4));
1280
1281     // 0^0 = 1 , 0^c1=0, -c1^0=1
1282     GMX_EXPECT_SIMD_REAL_NEAR(
1283             setSimdRealFrom3R(1.0, 0.0, 1.0),
1284             powSingleAccuracy(setSimdRealFrom3R(0, 0.0, -c1), setSimdRealFrom3R(0.0, c1, 0.0)));
1285 }
1286
1287 TEST_F(SimdMathTest, powSingleAccuracyUnsafe)
1288 {
1289     // We already test the log2/exp2 components of pow() extensively above, and it's a very
1290     // simple single-line function, so here we just test a handful of values to catch typos
1291     // and then some special values.
1292
1293     // Increase the allowed error by the difference between the actual precision and single
1294     setUlpTolSingleAccuracy(ulpTol_);
1295
1296     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::pow(c0, c3), std::pow(c1, c4), std::pow(c2, c5)),
1297                               powSingleAccuracy<MathOptimization::Unsafe>(rSimd_c0c1c2, rSimd_c3c4c5));
1298
1299     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::pow(c0, -c3), std::pow(c1, -c0), std::pow(c2, -c4)),
1300                               powSingleAccuracy<MathOptimization::Unsafe>(rSimd_c0c1c2, rSimd_m3m0m4));
1301 }
1302
1303 TEST_F(SimdMathTest, erfSingleAccuracy)
1304 {
1305     // Increase the allowed error by the difference between the actual precision and single
1306     setUlpTolSingleAccuracy(ulpTol_);
1307
1308     CompareSettings settings{ Range(-9, 9), ulpTol_, GMX_REAL_MIN, MatchRule::Normal };
1309     GMX_EXPECT_SIMD_FUNC_NEAR(refErf, erfSingleAccuracy, settings);
1310 }
1311
1312 TEST_F(SimdMathTest, erfcSingleAccuracy)
1313 {
1314     // Increase the allowed error by the difference between the actual precision and single
1315     setUlpTolSingleAccuracy(ulpTol_);
1316
1317     // Our erfc algorithm has 4 ulp accuracy, so relax tolerance a bit
1318     CompareSettings settings{ Range(-9, 9), 4 * ulpTol_, GMX_REAL_MIN, MatchRule::Normal };
1319     GMX_EXPECT_SIMD_FUNC_NEAR(refErfc, erfcSingleAccuracy, settings);
1320 }
1321
1322
1323 TEST_F(SimdMathTest, sinSingleAccuracy)
1324 {
1325     /* Increase the allowed error by the difference between the actual precision and single */
1326     setUlpTolSingleAccuracy(ulpTol_);
1327
1328     CompareSettings settings{ Range(-8 * M_PI, 8 * M_PI), ulpTol_, absTol_, MatchRule::Normal };
1329     GMX_EXPECT_SIMD_FUNC_NEAR(std::sin, sinSingleAccuracy, settings);
1330
1331     // Range reduction leads to accuracy loss, so we might want higher tolerance here
1332     settings = { Range(-10000, 10000), 2 * ulpTol_, absTol_, MatchRule::Normal };
1333     GMX_EXPECT_SIMD_FUNC_NEAR(std::sin, sinSingleAccuracy, settings);
1334 }
1335
1336 TEST_F(SimdMathTest, cosSingleAccuracy)
1337 {
1338     /* Increase the allowed error by the difference between the actual precision and single */
1339     setUlpTolSingleAccuracy(ulpTol_);
1340
1341     CompareSettings settings{ Range(-8 * M_PI, 8 * M_PI), ulpTol_, absTol_, MatchRule::Normal };
1342     GMX_EXPECT_SIMD_FUNC_NEAR(std::cos, cosSingleAccuracy, settings);
1343
1344     // Range reduction leads to accuracy loss, so we might want higher tolerance here
1345     settings = { Range(-10000, 10000), ulpTol_, absTol_, MatchRule::Normal };
1346     GMX_EXPECT_SIMD_FUNC_NEAR(std::cos, cosSingleAccuracy, settings);
1347 }
1348
1349 TEST_F(SimdMathTest, tanSingleAccuracy)
1350 {
1351     /* Increase the allowed error by the difference between the actual precision and single */
1352     setUlpTolSingleAccuracy(ulpTol_);
1353
1354     // Tan(x) is a little sensitive due to the division in the algorithm.
1355     // Rather than using lots of extra FP operations, we accept the algorithm
1356     // presently only achieves a ~3 ulp error and use the medium tolerance.
1357     CompareSettings settings{ Range(-8 * M_PI, 8 * M_PI), ulpTol_, absTol_, MatchRule::Normal };
1358     GMX_EXPECT_SIMD_FUNC_NEAR(std::tan, tanSingleAccuracy, settings);
1359
1360     // Range reduction leads to accuracy loss, so we might want higher tolerance here
1361     settings = { Range(-10000, 10000), ulpTol_, absTol_, MatchRule::Normal };
1362     GMX_EXPECT_SIMD_FUNC_NEAR(std::tan, tanSingleAccuracy, settings);
1363 }
1364
1365 TEST_F(SimdMathTest, asinSingleAccuracy)
1366 {
1367     /* Increase the allowed error by the difference between the actual precision and single */
1368     setUlpTolSingleAccuracy(ulpTol_);
1369
1370     // Our present asin(x) algorithm achieves 2-3 ulp accuracy
1371     CompareSettings settings{ Range(-1, 1), ulpTol_, absTol_, MatchRule::Normal };
1372     GMX_EXPECT_SIMD_FUNC_NEAR(std::asin, asinSingleAccuracy, settings);
1373 }
1374
1375 TEST_F(SimdMathTest, acosSingleAccuracy)
1376 {
1377     /* Increase the allowed error by the difference between the actual precision and single */
1378     setUlpTolSingleAccuracy(ulpTol_);
1379
1380     // Our present acos(x) algorithm achieves 2-3 ulp accuracy
1381     CompareSettings settings{ Range(-1, 1), ulpTol_, absTol_, MatchRule::Normal };
1382     GMX_EXPECT_SIMD_FUNC_NEAR(std::acos, acosSingleAccuracy, settings);
1383 }
1384
1385 TEST_F(SimdMathTest, atanSingleAccuracy)
1386 {
1387     /* Increase the allowed error by the difference between the actual precision and single */
1388     setUlpTolSingleAccuracy(ulpTol_);
1389
1390     // Our present atan(x) algorithm achieves 1 ulp accuracy
1391     CompareSettings settings{ Range(-10000, 10000), ulpTol_, absTol_, MatchRule::Normal };
1392     GMX_EXPECT_SIMD_FUNC_NEAR(std::atan, atanSingleAccuracy, settings);
1393 }
1394
1395 TEST_F(SimdMathTest, atan2SingleAccuracy)
1396 {
1397     /* Increase the allowed error by the difference between the actual precision and single */
1398     setUlpTolSingleAccuracy(ulpTol_);
1399
1400     // test each quadrant
1401     GMX_EXPECT_SIMD_REAL_NEAR(
1402             setSimdRealFrom3R(std::atan2(c0, c3), std::atan2(c1, c4), std::atan2(c2, c5)),
1403             atan2SingleAccuracy(rSimd_c0c1c2, rSimd_c3c4c5));
1404     GMX_EXPECT_SIMD_REAL_NEAR(
1405             setSimdRealFrom3R(std::atan2(-c0, c3), std::atan2(-c1, c4), std::atan2(-c2, c5)),
1406             atan2SingleAccuracy(rSimd_m0m1m2, rSimd_c3c4c5));
1407     GMX_EXPECT_SIMD_REAL_NEAR(
1408             setSimdRealFrom3R(std::atan2(-c0, -c3), std::atan2(-c1, -c0), std::atan2(-c2, -c4)),
1409             atan2SingleAccuracy(rSimd_m0m1m2, rSimd_m3m0m4));
1410     GMX_EXPECT_SIMD_REAL_NEAR(
1411             setSimdRealFrom3R(std::atan2(c0, -c3), std::atan2(c1, -c0), std::atan2(c2, -c4)),
1412             atan2SingleAccuracy(rSimd_c0c1c2, rSimd_m3m0m4));
1413     // cases important for calculating angles
1414     // values on coordinate axes
1415     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(0, c0), std::atan2(0, c1), std::atan2(0, c2)),
1416                               atan2SingleAccuracy(setZero(), rSimd_c0c1c2));
1417     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(c0, 0), std::atan2(c1, 0), std::atan2(c2, 0)),
1418                               atan2SingleAccuracy(rSimd_c0c1c2, setZero()));
1419     GMX_EXPECT_SIMD_REAL_NEAR(
1420             setSimdRealFrom3R(std::atan2(0, -c0), std::atan2(0, -c1), std::atan2(0, -c2)),
1421             atan2SingleAccuracy(setZero(), rSimd_m0m1m2));
1422     GMX_EXPECT_SIMD_REAL_NEAR(
1423             setSimdRealFrom3R(std::atan2(-c0, 0), std::atan2(-c1, 0), std::atan2(-c2, 0)),
1424             atan2SingleAccuracy(rSimd_m0m1m2, setZero()));
1425
1426     // degenerate value (origin) should return 0.0. At least IBM xlc 13.1.5 gets the reference
1427     // value wrong (-nan) at -O3 optimization, so we compare to the correct value (0.0) instead.
1428     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(0.0),
1429                               atan2SingleAccuracy(setSimdRealFrom3R(0.0, 0.0, 0.0), setZero()));
1430 }
1431
1432 TEST_F(SimdMathTest, pmeForceCorrectionSingleAccuracy)
1433 {
1434     // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
1435     // Pme correction only needs to be ~1e-6 accuracy single.
1436     // Then increase the allowed error by the difference between the actual precision and single.
1437     setUlpTolSingleAccuracy(std::int64_t(5e-6 / GMX_FLOAT_EPS));
1438
1439     CompareSettings settings{ Range(0.15, 4), ulpTol_, GMX_FLOAT_EPS, MatchRule::Normal };
1440     GMX_EXPECT_SIMD_FUNC_NEAR(refPmeForceCorrection, pmeForceCorrectionSingleAccuracy, settings);
1441 }
1442
1443 TEST_F(SimdMathTest, pmePotentialCorrectionSingleAccuracy)
1444 {
1445     // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
1446     // Pme correction only needs to be ~1e-6 accuracy single.
1447     // Then increase the allowed error by the difference between the actual precision and single.
1448     setUlpTolSingleAccuracy(std::int64_t(5e-6 / GMX_FLOAT_EPS));
1449
1450     CompareSettings settings{ Range(0.15, 4), ulpTol_, GMX_FLOAT_EPS, MatchRule::Normal };
1451     GMX_EXPECT_SIMD_FUNC_NEAR(refPmePotentialCorrection, pmePotentialCorrectionSingleAccuracy, settings);
1452 }
1453
1454 } // namespace
1455
1456 #    endif // GMX_SIMD_HAVE_REAL
1457
1458 /*! \} */
1459 /*! \endcond */
1460
1461 } // namespace test
1462 } // namespace gmx
1463
1464 #endif // GMX_SIMD