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, 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 <vector>
  45
  46 #include "gromacs/math/utilities.h"
  47 #include "gromacs/options/basicoptions.h"
  48 #include "gromacs/simd/simd.h"
  49
  50 #include "simd.h"
  51
  52 #if GMX_SIMD
  53
  54 namespace gmx
  55 {
  56 namespace test
  57 {
  58
  59 /*! \cond internal */
  60 /*! \addtogroup module_simd */
  61 /*! \{ */
  62
  63 #if GMX_SIMD_HAVE_REAL
  64
  65 class SimdMathTest : public SimdTest
  66 {
  67     public:
  68         ::testing::AssertionResult
  69                             compareSimdMathFunction(const char *  refFuncExpr,
  70                                                     const char *  simdFuncExpr,
  71                                                     const char *  denormalsToZeroExpr,
  72                                                     real          refFunc(real x),
  73                                                     SimdReal      gmx_simdcall simdFunc(SimdReal x),
  74                                                     bool          denormalsToZero);
  75 };
  76
  77 /*! \brief Test approximate equality of SIMD vs reference version of a function.
  78  *
  79  * This macro takes vanilla C and SIMD flavors of a function and tests it with
  80  * the number of points, range, and tolerances specified by the test fixture class.
  81  */
  82 #define GMX_EXPECT_SIMD_FUNC_NEAR(refFunc, tstFunc) \
  83     EXPECT_PRED_FORMAT3(compareSimdMathFunction, refFunc, tstFunc, false)
  84
  85 /*! \brief Test approximate equality of SIMD vs reference function, denormals can be zero.
  86  *
  87  * This macro takes vanilla C and SIMD flavors of a function and tests it with
  88  * the number of points, range, and tolerances specified by the test fixture class.
  89  *
  90  * This version of the function will also return success if the test function
  91  * returns zero where the reference function returns a denormal value.
  92  */
  93 #define GMX_EXPECT_SIMD_FUNC_NEAR_DTZ(refFunc, tstFunc) \
  94     EXPECT_PRED_FORMAT3(compareSimdMathFunction, refFunc, tstFunc, true)
  95
  96 /*! \brief Implementation routine to compare SIMD vs reference functions.
  97  *
  98  * \param refFuncExpr         Description of reference function expression
  99  * \param simdFuncExpr        Description of SIMD function expression
 100  * \param denormalsToZeroExpr Description of denormal-to-zero setting
 101  * \param refFunc             Reference math function pointer
 102  * \param simdFunc            SIMD math function pointer
 103  * \param denormalsToZero     If true, the function will consider denormal
 104  *                            values equivalent to 0.0.
 105  *
 106  * The function will be tested with the range and tolerances specified in
 107  * the SimdBaseTest class. You should not never call this function directly,
 108  * but use the macro GMX_EXPECT_SIMD_FUNC_NEAR(refFunc,tstFunc) instead.
 109  */
 110 ::testing::AssertionResult
 111 SimdMathTest::compareSimdMathFunction(const char              * refFuncExpr,
 112                                       const char              * simdFuncExpr,
 113                                       const char              * denormalsToZeroExpr,
 114                                       real                      refFunc(real x),
 115                                       SimdReal     gmx_simdcall simdFunc(SimdReal x),
 116                                       bool                      denormalsToZero)
 117 {
 118     std::vector<real>            vx(GMX_SIMD_REAL_WIDTH);
 119     std::vector<real>            vref(GMX_SIMD_REAL_WIDTH);
 120     std::vector<real>            vtst(GMX_SIMD_REAL_WIDTH);
 121     real                         dx, absDiff;
 122     std::int64_t                 ulpDiff, maxUlpDiff;
 123     real                         maxUlpDiffPos;
 124     real                         refValMaxUlpDiff, simdValMaxUlpDiff;
 125     const int                    niter   = s_nPoints/GMX_SIMD_REAL_WIDTH;
 126     int                          npoints = niter*GMX_SIMD_REAL_WIDTH;
 127 #    if GMX_DOUBLE
 128     union {
 129         double r; std::int64_t i;
 130     } conv0, conv1;
 131 #    else
 132     union {
 133         float  r; std::int32_t i;
 134     } conv0, conv1;
 135 #    endif
 136
 137     maxUlpDiff = 0;
 138     dx         = (range_.second-range_.first)/npoints;
 139
 140     for (int iter = 0; iter < niter; iter++)
 141     {
 142         for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 143         {
 144             vx[i]   = range_.first+dx*(iter*GMX_SIMD_REAL_WIDTH+i);
 145             vref[i] = refFunc(vx[i]);
 146         }
 147         vtst  = simdReal2Vector(simdFunc(vector2SimdReal(vx)));
 148
 149         bool absOk = true, signOk = true;
 150         for (int i = 0; i < GMX_SIMD_REAL_WIDTH; i++)
 151         {
 152             if (denormalsToZero)
 153             {
 154                 // Clamp denormal values to zero if requested
 155                 if (std::abs(vref[i]) <= GMX_REAL_MIN)
 156                 {
 157                     vref[i] = 0.0;
 158                 }
 159                 if (std::abs(vtst[i]) <= GMX_REAL_MIN)
 160                 {
 161                     vtst[i] = 0.0;
 162                 }
 163             }
 164
 165             absDiff = std::abs(vref[i]-vtst[i]);
 166             absOk   = absOk  && ( absDiff < absTol_ );
 167             signOk  = signOk && ( (vref[i] >= 0 && vtst[i] >= 0) ||
 168                                   (vref[i] <= 0 && vtst[i] <= 0));
 169
 170             if (absDiff >= absTol_)
 171             {
 172                 /* We replicate the trivial ulp differences comparison here rather than
 173                  * calling the lower-level routine for comparing them, since this enables
 174                  * us to run through the entire test range and report the largest deviation
 175                  * without lots of extra glue routines.
 176                  */
 177                 conv0.r           = vref[i];
 178                 conv1.r           = vtst[i];
 179                 ulpDiff           = llabs(conv0.i-conv1.i);
 180                 if (ulpDiff > maxUlpDiff)
 181                 {
 182                     maxUlpDiff        = ulpDiff;
 183                     maxUlpDiffPos     = vx[i];
 184                     refValMaxUlpDiff  = vref[i];
 185                     simdValMaxUlpDiff = vtst[i];
 186                 }
 187             }
 188         }
 189         if ( (!absOk) && (!signOk) )
 190         {
 191             return ::testing::AssertionFailure()
 192                    << "Failing SIMD math function comparison due to sign differences." << std::endl
 193                    << "Reference function: " << refFuncExpr << std::endl
 194                    << "Simd function:      " << simdFuncExpr << std::endl
 195                    << "Test range is ( " << range_.first << " , " << range_.second << " ) " << std::endl
 196                    << "Denormals can be 0: " << denormalsToZeroExpr << std::endl
 197                    << "First sign difference around x=" << std::setprecision(20) << ::testing::PrintToString(vx) << std::endl
 198                    << "Ref values:  " << std::setprecision(20) << ::testing::PrintToString(vref) << std::endl
 199                    << "SIMD values: " << std::setprecision(20) << ::testing::PrintToString(vtst) << std::endl;
 200         }
 201     }
 202     GMX_RELEASE_ASSERT(ulpTol_ >= 0, "Invalid ulp value.");
 203     if (maxUlpDiff <= ulpTol_)
 204     {
 205         return ::testing::AssertionSuccess();
 206     }
 207     else
 208     {
 209         return ::testing::AssertionFailure()
 210                << "Failing SIMD math function ulp comparison between " << refFuncExpr << " and " << simdFuncExpr << std::endl
 211                << "Requested ulp tolerance: " << ulpTol_ << std::endl
 212                << "Requested abs tolerance: " << absTol_ << std::endl
 213                << "Denormals can be 0: " << denormalsToZeroExpr << std::endl
 214                << "Largest Ulp difference occurs for x=" << std::setprecision(20) << maxUlpDiffPos << std::endl
 215                << "Ref  values: " << std::setprecision(20) << refValMaxUlpDiff << std::endl
 216                << "SIMD values: " << std::setprecision(20) << simdValMaxUlpDiff << std::endl
 217                << "Ulp diff.:   " << std::setprecision(20) << maxUlpDiff << std::endl;
 218     }
 219 }
 220
 221 /*! \} */
 222 /*! \endcond */
 223
 224
 225 // Actual math function tests below
 226
 227
 228 namespace
 229 {
 230
 231 /*! \cond internal */
 232 /*! \addtogroup module_simd */
 233 /*! \{ */
 234
 235 TEST_F(SimdMathTest, copysign)
 236 {
 237     GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(-c0,  c1, c2), copysign(setSimdRealFrom3R( c0,  c1,  c2), setSimdRealFrom3R(-c3,  c4, 0)));
 238     GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R(-c0,  c1, c2), copysign(setSimdRealFrom3R(-c0, -c1, -c2), setSimdRealFrom3R(-c3,  c4, 0)));
 239     GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R( c0, -c1, c2), copysign(setSimdRealFrom3R( c0,  c1,  c2), setSimdRealFrom3R( c3, -c4, 0)));
 240     GMX_EXPECT_SIMD_REAL_EQ(setSimdRealFrom3R( c0, -c1, c2), copysign(setSimdRealFrom3R(-c0, -c1, -c2), setSimdRealFrom3R( c3, -c4, 0)));
 241 }
 242
 243 /*! \brief Function wrapper to evaluate reference 1/sqrt(x) */
 244 real
 245 refInvsqrt(real x)
 246 {
 247     return 1.0/std::sqrt(x);
 248 }
 249
 250 TEST_F(SimdMathTest, invsqrt)
 251 {
 252     setRange(1.01*GMX_FLOAT_MIN, GMX_FLOAT_MAX);
 253     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, invsqrt);
 254 }
 255
 256 TEST_F(SimdMathTest, maskzInvsqrt)
 257 {
 258     SimdReal x   = setSimdRealFrom3R(c1, 0.0, c2);
 259     SimdBool m   = (setZero() < x);
 260     SimdReal ref = setSimdRealFrom3R(1.0/std::sqrt(c1), 0.0, 1.0/std::sqrt(c2));
 261     GMX_EXPECT_SIMD_REAL_NEAR(ref, maskzInvsqrt(x, m));
 262 }
 263
 264 /*! \brief Function wrapper to return first result when testing \ref invsqrtPair */
 265 SimdReal gmx_simdcall
 266 tstInvsqrtPair0(SimdReal x)
 267 {
 268     SimdReal r0, r1;
 269     invsqrtPair(x, x, &r0, &r1);
 270     return r0;
 271 }
 272
 273 /*! \brief Function wrapper to return second result when testing \ref invsqrtPair */
 274 SimdReal gmx_simdcall
 275 tstInvsqrtPair1(SimdReal x)
 276 {
 277     SimdReal r0, r1;
 278     invsqrtPair(x, x, &r0, &r1);
 279     return r1;
 280 }
 281
 282 TEST_F(SimdMathTest, invsqrtPair)
 283 {
 284     setRange(1.01*GMX_FLOAT_MIN, GMX_FLOAT_MAX);
 285     // The accuracy conversions lose a bit of extra accuracy compared to
 286     // doing the iterations in all-double.
 287     setUlpTol(4*ulpTol_);
 288
 289     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, tstInvsqrtPair0);
 290     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, tstInvsqrtPair1);
 291 }
 292
 293 /*! \brief Function wrapper to evaluate reference sqrt(x) */
 294 real
 295 refSqrt(real x)
 296 {
 297     return std::sqrt(x);
 298 }
 299
 300 /*! \brief Dummy function returning 0.0 to test function ranges that should be zero */
 301 gmx_unused real refZero(real gmx_unused x)
 302 {
 303     return 0.0;
 304 }
 305
 306
 307 TEST_F(SimdMathTest, sqrt)
 308 {
 309     // The accuracy conversions lose a bit of extra accuracy compared to
 310     // doing the iterations in all-double.
 311     setUlpTol(4*ulpTol_);
 312
 313     // First test that 0.0 and a few other values works
 314     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(0, std::sqrt(c1), std::sqrt(c2)), sqrt(setSimdRealFrom3R(0, c1, c2)));
 315
 316     // Values smaller-than-or-equal to GMX_FLOAT_MIN will be clamped to 0.0,
 317     // so only test larger values
 318     setRange(1.01*GMX_FLOAT_MIN, GMX_FLOAT_MAX);
 319     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrt);
 320
 321 #if GMX_DOUBLE
 322     // Make sure that values smaller than GMX_FLOAT_MIN lead to result 0.0
 323     setRange(0.0, 0.99*GMX_FLOAT_MIN);
 324     GMX_EXPECT_SIMD_FUNC_NEAR(refZero, sqrt);
 325 #endif
 326 }
 327
 328 TEST_F(SimdMathTest, sqrtUnsafe)
 329 {
 330     // The accuracy conversions lose a bit of extra accuracy compared to
 331     // doing the iterations in all-double.
 332     setUlpTol(4*ulpTol_);
 333
 334     setRange(GMX_FLOAT_MIN, GMX_FLOAT_MAX);
 335     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrt<MathOptimization::Unsafe>);
 336 }
 337
 338 /*! \brief Function wrapper to evaluate reference 1/x */
 339 real refInv(real x)
 340 {
 341     return 1.0/x;
 342 }
 343
 344 TEST_F(SimdMathTest, inv)
 345 {
 346     // test <0
 347     setRange(-1e10, -1e-10);
 348     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv);
 349     setRange(1e-10, 1e10);
 350     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, inv);
 351 }
 352
 353 TEST_F(SimdMathTest, maskzInv)
 354 {
 355     SimdReal x   = setSimdRealFrom3R(c1, 0.0, c2);
 356     SimdBool m   = (setZero() < x);
 357     SimdReal ref = setSimdRealFrom3R(1.0/c1, 0.0, 1.0/c2);
 358     GMX_EXPECT_SIMD_REAL_NEAR(ref, maskzInv(x, m));
 359 }
 360
 361 TEST_F(SimdMathTest, log)
 362 {
 363     setRange(1e-30, 1e30);
 364     GMX_EXPECT_SIMD_FUNC_NEAR(std::log, log);
 365 }
 366
 367 TEST_F(SimdMathTest, exp2)
 368 {
 369     // Test normal/denormal/zero range separately to make errors clearer
 370
 371     // First test the range where we get normalized (non-denormal) results,
 372     // since we don't require denormal results to be reproduced correctly.
 373 #if GMX_DOUBLE
 374     setRange(-1022, 1023);
 375 #else
 376     setRange(-126, 127);
 377 #endif
 378     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2);
 379
 380     // Some implementations might have denormal support, in which case they
 381     // support an extended range, adding roughly the number of bits in the
 382     // mantissa to the smallest allowed arg (1023+52 in double, 127+23 single).
 383     // In this range we allow the value to be either correct (denormal) or 0.0
 384 #if GMX_DOUBLE
 385     setRange(-1075, -1022);
 386 #else
 387     setRange(-150, -126);
 388 #endif
 389     GMX_EXPECT_SIMD_FUNC_NEAR_DTZ(std::exp2, exp2);
 390
 391     // For arguments smaller than the subnormal the result should be zero
 392     // both in the reference and our implementations.
 393 #if GMX_DOUBLE
 394     setRange(-1000000.0, -1075.0);
 395 #else
 396     setRange(-100000.0, -150.0);
 397 #endif
 398     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2);
 399
 400     // Test a few very negative values, including values so small that they
 401     // will start to cause inf values in the polynomial interpolations
 402     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::exp2(-GMX_FLOAT_MAX), std::exp2(-0.1*GMX_REAL_MAX), std::exp2(-GMX_REAL_MAX)),
 403                               exp2(setSimdRealFrom3R(-GMX_FLOAT_MAX, -0.1*GMX_REAL_MAX, -GMX_REAL_MAX)));
 404 }
 405
 406 TEST_F(SimdMathTest, exp2Unsafe)
 407 {
 408     // The unsafe version is only defined in this range
 409 #if GMX_DOUBLE
 410     setRange(-1022, 1023);
 411 #else
 412     setRange(-126, 127);
 413 #endif
 414     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2<MathOptimization::Unsafe>);
 415 }
 416
 417 TEST_F(SimdMathTest, exp)
 418 {
 419     // Test normal/denormal/zero range separately to make errors clearer
 420
 421     // First test the range where we get normalized (non-denormal) results,
 422     // since we don't require denormal results to be reproduced correctly.
 423     //
 424     // For very small arguments that would produce results close to the
 425     // smallest representable value, some of the intermediate values might
 426     // trigger flush-to-zero denormals without FMA operations,
 427     // e.g. for the icc compiler. Since we never use such values in Gromacs, we
 428     // shrink the range a bit in that case instead of requiring the compiler to
 429     // handle denormals (which might reduce performance).
 430 #if GMX_DOUBLE
 431 #if GMX_SIMD_HAVE_FMA
 432     setRange(-708.3, 709.1);
 433 #else
 434     setRange(-690, 709.1);
 435 #endif
 436 #else
 437 #if GMX_SIMD_HAVE_FMA
 438     setRange(-87.3, 88.0);
 439 #else
 440     setRange(-80, 88.0);
 441 #endif
 442 #endif
 443     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, exp);
 444
 445     // Some implementations might have denormal support, in which case they
 446     // support an extended range, adding roughly the number of bits in the
 447     // mantissa to the smallest allowed arg (1023+52 in double, 127+23 single).
 448     // Then multiply with ln(2) to get our limit for exp().
 449     // In this range we allow the value to be either correct (denormal) or 0.0
 450 #if GMX_DOUBLE
 451     setRange(-746.0, -708.4);
 452 #else
 453     setRange(-104.0, -87.3);
 454 #endif
 455     GMX_EXPECT_SIMD_FUNC_NEAR_DTZ(std::exp, exp);
 456
 457     // For arguments smaller than the subnormal the result should be zero
 458     // both in the reference and our implementations.
 459 #if GMX_DOUBLE
 460     setRange(-1000000.0, -746.0);
 461 #else
 462     setRange(-100000.0, -104.0);
 463 #endif
 464     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, exp);
 465
 466     // Test a few very negative values, including values so small that they
 467     // will start to cause inf values in the polynomial interpolations
 468     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::exp(-GMX_FLOAT_MAX), std::exp(-0.1*GMX_REAL_MAX), std::exp(-GMX_REAL_MAX)),
 469                               exp(setSimdRealFrom3R(-GMX_FLOAT_MAX, -0.1*GMX_REAL_MAX, -GMX_REAL_MAX)));
 470 }
 471
 472 TEST_F(SimdMathTest, expUnsafe)
 473 {
 474 #if GMX_DOUBLE
 475 #if GMX_SIMD_HAVE_FMA
 476     setRange(-708.3, 709.1);
 477 #else
 478     setRange(-690, 709.1);
 479 #endif
 480 #else
 481 #if GMX_SIMD_HAVE_FMA
 482     setRange(-87.3, 88.0);
 483 #else
 484     setRange(-80, 88.0);
 485 #endif
 486 #endif
 487     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, exp<MathOptimization::Unsafe>);
 488 }
 489
 490 /*! \brief Function wrapper for erf(x), with argument/return in default Gromacs precision.
 491  *
 492  * \note Single-precision erf() in some libraries can be slightly lower precision
 493  * than the SIMD flavor, so we use a cast to force double precision for reference.
 494  */
 495 real
 496 refErf(real x)
 497 {
 498     return std::erf(static_cast<double>(x));
 499 }
 500
 501 TEST_F(SimdMathTest, erf)
 502 {
 503     setRange(-9, 9);
 504     setAbsTol(GMX_REAL_MIN);
 505     GMX_EXPECT_SIMD_FUNC_NEAR(refErf, erf);
 506 }
 507
 508 /*! \brief Function wrapper for erfc(x), with argument/return in default Gromacs precision.
 509  *
 510  * \note Single-precision erfc() in some libraries can be slightly lower precision
 511  * than the SIMD flavor, so we use a cast to force double precision for reference.
 512  */
 513 real
 514 refErfc(real x)
 515 {
 516     return std::erfc(static_cast<double>(x));
 517 }
 518
 519 TEST_F(SimdMathTest, erfc)
 520 {
 521     setRange(-9, 9);
 522     setAbsTol(GMX_REAL_MIN);
 523     // Our erfc algorithm has 4 ulp accuracy, so relax defaultTol a bit
 524     setUlpTol(4*ulpTol_);
 525     GMX_EXPECT_SIMD_FUNC_NEAR(refErfc, erfc);
 526 }
 527
 528 TEST_F(SimdMathTest, sin)
 529 {
 530     setRange(-8*M_PI, 8*M_PI);
 531     GMX_EXPECT_SIMD_FUNC_NEAR(std::sin, sin);
 532     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 533     setRange(-10000, 10000);
 534     setUlpTol(2*ulpTol_);
 535     GMX_EXPECT_SIMD_FUNC_NEAR(std::sin, sin);
 536 }
 537
 538 TEST_F(SimdMathTest, cos)
 539 {
 540     setRange(-8*M_PI, 8*M_PI);
 541     GMX_EXPECT_SIMD_FUNC_NEAR(std::cos, cos);
 542     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 543     setRange(-10000, 10000);
 544     setUlpTol(2*ulpTol_);
 545     GMX_EXPECT_SIMD_FUNC_NEAR(std::cos, cos);
 546 }
 547
 548 TEST_F(SimdMathTest, tan)
 549 {
 550     // Tan(x) is a little sensitive due to the division in the algorithm.
 551     // Rather than using lots of extra FP operations, we accept the algorithm
 552     // presently only achieves a ~3 ulp error and use the medium tolerance.
 553     setRange(-8*M_PI, 8*M_PI);
 554     GMX_EXPECT_SIMD_FUNC_NEAR(std::tan, tan);
 555     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 556     setRange(-10000, 10000);
 557     setUlpTol(2*ulpTol_);
 558     GMX_EXPECT_SIMD_FUNC_NEAR(std::tan, tan);
 559 }
 560
 561 TEST_F(SimdMathTest, asin)
 562 {
 563     // Our present asin(x) algorithm achieves 2-3 ulp accuracy
 564     setRange(-1, 1);
 565     GMX_EXPECT_SIMD_FUNC_NEAR(std::asin, asin);
 566 }
 567
 568 TEST_F(SimdMathTest, acos)
 569 {
 570     // Our present acos(x) algorithm achieves 2-3 ulp accuracy
 571     setRange(-1, 1);
 572     GMX_EXPECT_SIMD_FUNC_NEAR(std::acos, acos);
 573 }
 574
 575 TEST_F(SimdMathTest, atan)
 576 {
 577     // Our present atan(x) algorithm achieves 1 ulp accuracy
 578     setRange(-10000, 10000);
 579     GMX_EXPECT_SIMD_FUNC_NEAR(std::atan, atan);
 580 }
 581
 582 TEST_F(SimdMathTest, atan2)
 583 {
 584     // test each quadrant
 585     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(c0, c3), std::atan2(c1, c4), std::atan2(c2, c5)),
 586                               atan2(rSimd_c0c1c2, rSimd_c3c4c5));
 587     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(-c0, c3), std::atan2(-c1, c4), std::atan2(-c2, c5)),
 588                               atan2(rSimd_m0m1m2, rSimd_c3c4c5));
 589     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(-c0, -c3), std::atan2(-c1, -c0), std::atan2(-c2, -c4)),
 590                               atan2(rSimd_m0m1m2, rSimd_m3m0m4));
 591     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(c0, -c3), std::atan2(c1, -c0), std::atan2(c2, -c4)),
 592                               atan2(rSimd_c0c1c2, rSimd_m3m0m4));
 593
 594     // cases important for calculating angles
 595     // values on coordinate axes
 596     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(0, c0), std::atan2(0, c1), std::atan2(0, c2)),
 597                               atan2(setZero(), rSimd_c0c1c2));
 598     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(c0, 0), std::atan2(c1, 0), std::atan2(c2, 0)),
 599                               atan2(rSimd_c0c1c2, setZero()));
 600     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(0, -c0), std::atan2(0, -c1), std::atan2(0, -c2)),
 601                               atan2(setZero(), rSimd_m0m1m2));
 602     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(-c0, 0), std::atan2(-c1, 0), std::atan2(-c2, 0)),
 603                               atan2(rSimd_m0m1m2, setZero()));
 604     // degenerate value (origin) should return 0.0. At least IBM xlc 13.1.5 gets the reference
 605     // value wrong (-nan) at -O3 optimization, so we compare to the correct value (0.0) instead.
 606     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(0.0), atan2(setSimdRealFrom3R(0.0, 0.0, 0.0), setZero()));
 607 }
 608
 609 /*! \brief Evaluate reference version of PME force correction. */
 610 real
 611 refPmeForceCorrection(real x)
 612 {
 613     if (x != 0)
 614     {
 615         real y = std::sqrt(x);
 616         return 2*std::exp(-x)/(std::sqrt(M_PI)*x) - std::erf(static_cast<double>(y))/(x*y);
 617     }
 618     else
 619     {
 620         return -4/(3*std::sqrt(M_PI));
 621     }
 622 }
 623
 624 // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
 625 TEST_F(SimdMathTest, pmeForceCorrection)
 626 {
 627     // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
 628 #if GMX_DOUBLE
 629     setUlpTol(std::int64_t(5e-10/GMX_REAL_EPS));
 630 #else
 631     setUlpTol(std::int64_t(5e-6/GMX_REAL_EPS));
 632 #endif
 633
 634     setRange(0.15, 4);
 635     setAbsTol(GMX_REAL_EPS);
 636     GMX_EXPECT_SIMD_FUNC_NEAR(refPmeForceCorrection, pmeForceCorrection);
 637 }
 638
 639 /*! \brief Evaluate reference version of PME potential correction. */
 640 real
 641 refPmePotentialCorrection(real x)
 642 {
 643     real y = std::sqrt(x);
 644     return std::erf(static_cast<double>(y))/y;
 645 }
 646
 647 // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
 648 TEST_F(SimdMathTest, pmePotentialCorrection)
 649 {
 650     // Pme correction only needs to be ~1e-6 accuracy single, 1e-10 double
 651 #if GMX_DOUBLE
 652     setUlpTol(std::int64_t(5e-10/GMX_REAL_EPS));
 653 #else
 654     setUlpTol(std::int64_t(5e-6/GMX_REAL_EPS));
 655 #endif
 656     setRange(0.15, 4);
 657     setAbsTol(GMX_REAL_EPS);
 658     GMX_EXPECT_SIMD_FUNC_NEAR(refPmePotentialCorrection, pmePotentialCorrection);
 659 }
 660
 661 // Functions that only target single accuracy, even for double SIMD data
 662
 663 TEST_F(SimdMathTest, invsqrtSingleAccuracy)
 664 {
 665     /* Increase the allowed error by the difference between the actual precision and single */
 666     setUlpTolSingleAccuracy(ulpTol_);
 667
 668     setRange(1.01*GMX_FLOAT_MIN, GMX_FLOAT_MAX);
 669     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, invsqrtSingleAccuracy);
 670 }
 671
 672 /*! \brief Function wrapper to return first result when testing \ref invsqrtPairSingleAccuracy */
 673 SimdReal gmx_simdcall
 674 tst_invsqrt_SingleAccuracy_pair0(SimdReal x)
 675 {
 676     SimdReal r0, r1;
 677     invsqrtPairSingleAccuracy(x, x, &r0, &r1);
 678     return r0;
 679 }
 680
 681 /*! \brief Function wrapper to return second result when testing \ref invsqrtPairSingleAccuracy */
 682 SimdReal gmx_simdcall
 683 tst_invsqrt_SingleAccuracy_pair1(SimdReal x)
 684 {
 685     SimdReal r0, r1;
 686     invsqrtPairSingleAccuracy(x, x, &r0, &r1);
 687     return r1;
 688 }
 689
 690 TEST_F(SimdMathTest, invsqrtPairSingleAccuracy)
 691 {
 692     /* Increase the allowed error by the difference between the actual precision and single */
 693     setUlpTolSingleAccuracy(ulpTol_);
 694
 695     setRange(1.01*GMX_FLOAT_MIN, GMX_FLOAT_MAX);
 696     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, tst_invsqrt_SingleAccuracy_pair0);
 697     GMX_EXPECT_SIMD_FUNC_NEAR(refInvsqrt, tst_invsqrt_SingleAccuracy_pair1);
 698 }
 699
 700 TEST_F(SimdMathTest, sqrtSingleAccuracy)
 701 {
 702     /* Increase the allowed error by the difference between the actual precision and single */
 703     setUlpTolSingleAccuracy(ulpTol_);
 704
 705     // First test that 0.0 and a few other values works
 706     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(0, std::sqrt(c0), std::sqrt(c1)), sqrtSingleAccuracy(setSimdRealFrom3R(0, c0, c1)));
 707
 708     // Values smaller-than-or-equal to GMX_FLOAT_MIN will be clamped to 0.0,
 709     // so only test larger values
 710     setRange(1.01*GMX_FLOAT_MIN, GMX_FLOAT_MAX);
 711     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrtSingleAccuracy);
 712
 713 #if GMX_DOUBLE
 714     // Make sure that values smaller than GMX_FLOAT_MIN lead to result 0.0
 715     setRange(0.0, 0.99*GMX_FLOAT_MIN);
 716     GMX_EXPECT_SIMD_FUNC_NEAR(refZero, sqrtSingleAccuracy);
 717 #endif
 718 }
 719
 720 TEST_F(SimdMathTest, sqrtSingleAccuracyUnsafe)
 721 {
 722     /* Increase the allowed error by the difference between the actual precision and single */
 723     setUlpTolSingleAccuracy(ulpTol_);
 724
 725     // Test the full range
 726     setRange(GMX_FLOAT_MIN, GMX_FLOAT_MAX);
 727     GMX_EXPECT_SIMD_FUNC_NEAR(refSqrt, sqrtSingleAccuracy<MathOptimization::Unsafe>);
 728 }
 729
 730 TEST_F(SimdMathTest, invSingleAccuracy)
 731 {
 732     /* Increase the allowed error by the difference between the actual precision and single */
 733     setUlpTolSingleAccuracy(ulpTol_);
 734
 735     // test <0
 736     setRange(-1e10, -1e-10);
 737     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, invSingleAccuracy);
 738     setRange(1e-10, 1e10);
 739     GMX_EXPECT_SIMD_FUNC_NEAR(refInv, invSingleAccuracy);
 740 }
 741
 742 TEST_F(SimdMathTest, logSingleAccuracy)
 743 {
 744     /* Increase the allowed error by the difference between the actual precision and single */
 745     setUlpTolSingleAccuracy(ulpTol_);
 746
 747     setRange(1e-30, 1e30);
 748     GMX_EXPECT_SIMD_FUNC_NEAR(std::log, logSingleAccuracy);
 749 }
 750
 751 TEST_F(SimdMathTest, exp2SingleAccuracy)
 752 {
 753     /* Increase the allowed error by the difference between the actual precision and single */
 754     setUlpTolSingleAccuracy(ulpTol_);
 755
 756 #if GMX_DOUBLE
 757     setRange(-1022.49, 1023.49);
 758 #else
 759     setRange(-126, 127.49);
 760 #endif
 761     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2SingleAccuracy);
 762
 763     // Test a range that should be zero both for reference and simd versions.
 764     // Some implementations might have subnormal support, in which case they
 765     // support an extended range, adding roughly the number of bits in the
 766     // mantissa to the smallest allowed arg (1023+52 in double, 127+23 single).
 767 #if GMX_DOUBLE
 768     setRange(-1000000.0, -1075.0);
 769 #else
 770     setRange(-100000.0, -150.0);
 771 #endif
 772     GMX_EXPECT_SIMD_FUNC_NEAR_DTZ(std::exp2, exp2SingleAccuracy);
 773
 774     // Test a few very negative values
 775     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::exp2(-GMX_FLOAT_MAX), std::exp2(-0.1*GMX_REAL_MAX), std::exp2(-GMX_REAL_MAX)),
 776                               exp2SingleAccuracy(setSimdRealFrom3R(-GMX_FLOAT_MAX, -0.1*GMX_REAL_MAX, -GMX_REAL_MAX)));
 777 }
 778
 779 TEST_F(SimdMathTest, exp2SingleAccuracyUnsafe)
 780 {
 781     /* Increase the allowed error by the difference between the actual precision and single */
 782     setUlpTolSingleAccuracy(ulpTol_);
 783
 784 #if GMX_DOUBLE
 785     setRange(-1022.49, 1023.49);
 786 #else
 787     setRange(-126, 127.49);
 788 #endif
 789     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp2, exp2SingleAccuracy<MathOptimization::Unsafe>);
 790 }
 791
 792 TEST_F(SimdMathTest, expSingleAccuracy)
 793 {
 794     /* Increase the allowed error by the difference between the actual precision and single */
 795     setUlpTolSingleAccuracy(ulpTol_);
 796
 797 #if GMX_DOUBLE
 798     setRange(-708.7, 709.4);
 799 #else
 800     setRange(-87.3, 88.3);
 801 #endif
 802     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, expSingleAccuracy);
 803
 804     // Test a range that should be zero both for reference and simd versions.
 805     // Some implementations might have subnormal support, in which case they
 806     // support an extended range, adding roughly the number of bits in the
 807     // mantissa to the smallest result exponent (1023+52 in double, 127+23 single).
 808     // Then multiply with ln(2) to get our limit for exp().
 809 #if GMX_DOUBLE
 810     setRange(-1000000.0, -746.0);
 811 #else
 812     setRange(-100000.0, -104.0);
 813 #endif
 814     GMX_EXPECT_SIMD_FUNC_NEAR_DTZ(std::exp, expSingleAccuracy);
 815
 816     // Test a few very negative values, including values so small that they
 817     // will start to cause inf values in the polynomial interpolations
 818     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::exp(-GMX_FLOAT_MAX), std::exp(-0.1*GMX_REAL_MAX), std::exp(-GMX_REAL_MAX)),
 819                               expSingleAccuracy(setSimdRealFrom3R(-GMX_FLOAT_MAX, -0.1*GMX_REAL_MAX, -GMX_REAL_MAX)));
 820 }
 821
 822 TEST_F(SimdMathTest, expSingleAccuracyUnsafe)
 823 {
 824     /* Increase the allowed error by the difference between the actual precision and single */
 825     setUlpTolSingleAccuracy(ulpTol_);
 826
 827 #if GMX_DOUBLE
 828     setRange(-708.7, 709.4);
 829 #else
 830     setRange(-87.3, 88.3);
 831 #endif
 832     GMX_EXPECT_SIMD_FUNC_NEAR(std::exp, expSingleAccuracy<MathOptimization::Unsafe>);
 833 }
 834
 835 TEST_F(SimdMathTest, erfSingleAccuracy)
 836 {
 837     /* Increase the allowed error by the difference between the actual precision and single */
 838     setUlpTolSingleAccuracy(ulpTol_);
 839
 840     setRange(-9, 9);
 841     setAbsTol(GMX_REAL_MIN);
 842     GMX_EXPECT_SIMD_FUNC_NEAR(refErf, erfSingleAccuracy);
 843 }
 844
 845 TEST_F(SimdMathTest, erfcSingleAccuracy)
 846 {
 847     /* Increase the allowed error by the difference between the actual precision and single */
 848     setUlpTolSingleAccuracy(ulpTol_);
 849
 850     setRange(-9, 9);
 851     setAbsTol(GMX_REAL_MIN);
 852     // Our erfc algorithm has 4 ulp accuracy, so relax tolerance a bit
 853     setUlpTol(4*ulpTol_);
 854     GMX_EXPECT_SIMD_FUNC_NEAR(refErfc, erfcSingleAccuracy);
 855 }
 856
 857
 858 TEST_F(SimdMathTest, sinSingleAccuracy)
 859 {
 860     /* Increase the allowed error by the difference between the actual precision and single */
 861     setUlpTolSingleAccuracy(ulpTol_);
 862
 863     setRange(-8*M_PI, 8*M_PI);
 864     GMX_EXPECT_SIMD_FUNC_NEAR(std::sin, sinSingleAccuracy);
 865     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 866     setRange(-10000, 10000);
 867     setUlpTol(2*ulpTol_);
 868     GMX_EXPECT_SIMD_FUNC_NEAR(std::sin, sinSingleAccuracy);
 869 }
 870
 871 TEST_F(SimdMathTest, cosSingleAccuracy)
 872 {
 873     /* Increase the allowed error by the difference between the actual precision and single */
 874     setUlpTolSingleAccuracy(ulpTol_);
 875
 876     setRange(-8*M_PI, 8*M_PI);
 877     GMX_EXPECT_SIMD_FUNC_NEAR(std::cos, cosSingleAccuracy);
 878     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 879     setRange(-10000, 10000);
 880     GMX_EXPECT_SIMD_FUNC_NEAR(std::cos, cosSingleAccuracy);
 881 }
 882
 883 TEST_F(SimdMathTest, tanSingleAccuracy)
 884 {
 885     /* Increase the allowed error by the difference between the actual precision and single */
 886     setUlpTolSingleAccuracy(ulpTol_);
 887
 888     // Tan(x) is a little sensitive due to the division in the algorithm.
 889     // Rather than using lots of extra FP operations, we accept the algorithm
 890     // presently only achieves a ~3 ulp error and use the medium tolerance.
 891     setRange(-8*M_PI, 8*M_PI);
 892     GMX_EXPECT_SIMD_FUNC_NEAR(std::tan, tanSingleAccuracy);
 893     // Range reduction leads to accuracy loss, so we might want higher tolerance here
 894     setRange(-10000, 10000);
 895     GMX_EXPECT_SIMD_FUNC_NEAR(std::tan, tanSingleAccuracy);
 896 }
 897
 898 TEST_F(SimdMathTest, asinSingleAccuracy)
 899 {
 900     /* Increase the allowed error by the difference between the actual precision and single */
 901     setUlpTolSingleAccuracy(ulpTol_);
 902
 903     // Our present asin(x) algorithm achieves 2-3 ulp accuracy
 904     setRange(-1, 1);
 905     GMX_EXPECT_SIMD_FUNC_NEAR(std::asin, asinSingleAccuracy);
 906 }
 907
 908 TEST_F(SimdMathTest, acosSingleAccuracy)
 909 {
 910     /* Increase the allowed error by the difference between the actual precision and single */
 911     setUlpTolSingleAccuracy(ulpTol_);
 912
 913     // Our present acos(x) algorithm achieves 2-3 ulp accuracy
 914     setRange(-1, 1);
 915     GMX_EXPECT_SIMD_FUNC_NEAR(std::acos, acosSingleAccuracy);
 916 }
 917
 918 TEST_F(SimdMathTest, atanSingleAccuracy)
 919 {
 920     /* Increase the allowed error by the difference between the actual precision and single */
 921     setUlpTolSingleAccuracy(ulpTol_);
 922
 923     // Our present atan(x) algorithm achieves 1 ulp accuracy
 924     setRange(-10000, 10000);
 925     GMX_EXPECT_SIMD_FUNC_NEAR(std::atan, atanSingleAccuracy);
 926 }
 927
 928 TEST_F(SimdMathTest, atan2SingleAccuracy)
 929 {
 930     /* Increase the allowed error by the difference between the actual precision and single */
 931     setUlpTolSingleAccuracy(ulpTol_);
 932
 933     // test each quadrant
 934     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(c0, c3), std::atan2(c1, c4), std::atan2(c2, c5)),
 935                               atan2SingleAccuracy(rSimd_c0c1c2, rSimd_c3c4c5));
 936     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(-c0, c3), std::atan2(-c1, c4), std::atan2(-c2, c5)),
 937                               atan2SingleAccuracy(rSimd_m0m1m2, rSimd_c3c4c5));
 938     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(-c0, -c3), std::atan2(-c1, -c0), std::atan2(-c2, -c4)),
 939                               atan2SingleAccuracy(rSimd_m0m1m2, rSimd_m3m0m4));
 940     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(c0, -c3), std::atan2(c1, -c0), std::atan2(c2, -c4)),
 941                               atan2SingleAccuracy(rSimd_c0c1c2, rSimd_m3m0m4));
 942     // cases important for calculating angles
 943     // values on coordinate axes
 944     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(0, c0), std::atan2(0, c1), std::atan2(0, c2)),
 945                               atan2SingleAccuracy(setZero(), rSimd_c0c1c2));
 946     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(c0, 0), std::atan2(c1, 0), std::atan2(c2, 0)),
 947                               atan2SingleAccuracy(rSimd_c0c1c2, setZero()));
 948     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(0, -c0), std::atan2(0, -c1), std::atan2(0, -c2)),
 949                               atan2SingleAccuracy(setZero(), rSimd_m0m1m2));
 950     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom3R(std::atan2(-c0, 0), std::atan2(-c1, 0), std::atan2(-c2, 0)),
 951                               atan2SingleAccuracy(rSimd_m0m1m2, setZero()));
 952
 953     // degenerate value (origin) should return 0.0. At least IBM xlc 13.1.5 gets the reference
 954     // value wrong (-nan) at -O3 optimization, so we compare to the correct value (0.0) instead.
 955     GMX_EXPECT_SIMD_REAL_NEAR(setSimdRealFrom1R(0.0), atan2SingleAccuracy(setSimdRealFrom3R(0.0, 0.0, 0.0), setZero()));
 956 }
 957
 958 TEST_F(SimdMathTest, pmeForceCorrectionSingleAccuracy)
 959 {
 960     // The PME corrections will be added to ~1/r2, so absolute tolerance of EPS is fine.
 961     // Pme correction only needs to be ~1e-6 accuracy single.
 962     // Then increase the allowed error by the difference between the actual precision and single.
 963     setUlpTolSingleAccuracy(std::int64_t(5e-6/GMX_FLOAT_EPS));
 964
 965     setRange(0.15, 4);
 966     setAbsTol(GMX_FLOAT_EPS);
 967     GMX_EXPECT_SIMD_FUNC_NEAR(refPmeForceCorrection, pmeForceCorrectionSingleAccuracy);
 968 }
 969
 970 TEST_F(SimdMathTest, pmePotentialCorrectionSingleAccuracy)
 971 {
 972     // The PME corrections will be added to ~1/r, so absolute tolerance of EPS is fine.
 973     // Pme correction only needs to be ~1e-6 accuracy single.
 974     // Then increase the allowed error by the difference between the actual precision and single.
 975     setUlpTolSingleAccuracy(std::int64_t(5e-6/GMX_FLOAT_EPS));
 976
 977     setRange(0.15, 4);
 978     setAbsTol(GMX_FLOAT_EPS);
 979     GMX_EXPECT_SIMD_FUNC_NEAR(refPmePotentialCorrection, pmePotentialCorrectionSingleAccuracy);
 980 }
 981
 982 }      // namespace
 983
 984 #endif // GMX_SIMD_HAVE_REAL
 985
 986 /*! \} */
 987 /*! \endcond */
 988
 989 }      // namespace
 990 }      // namespace
 991
 992 #endif // GMX_SIMD