src/gromacs/simd/simd_math.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2012,2013,2014, 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 #ifndef GMX_SIMD_SIMD_MATH_H_
  36 #define GMX_SIMD_SIMD_MATH_H_
  37
  38 /*! \libinternal \file
  39  *
  40  * \brief Math functions for SIMD datatypes.
  41  *
  42  * \attention This file is generic for all SIMD architectures, so you cannot
  43  * assume that any of the optional SIMD features (as defined in simd.h) are
  44  * present. In particular, this means you cannot assume support for integers,
  45  * logical operations (neither on floating-point nor integer values), shifts,
  46  * and the architecture might only have SIMD for either float or double.
  47  * Second, to keep this file clean and general, any additions to this file
  48  * must work for all possible SIMD architectures in both single and double
  49  * precision (if they support it), and you cannot make any assumptions about
  50  * SIMD width.
  51  *
  52  * \author Erik Lindahl <erik.lindahl@scilifelab.se>
  53  *
  54  * \inlibraryapi
  55  * \ingroup module_simd
  56  */
  57
  58 #include <math.h>
  59
  60 #include "gromacs/math/utilities.h"
  61 #include "gromacs/simd/simd.h"
  62
  63 /*! \cond libapi */
  64 /*! \addtogroup module_simd */
  65 /*! \{ */
  66
  67 /*! \name Implementation accuracy settings
  68  *  \{
  69  */
  70
  71 /*! \brief We accept lsb errors for 1/sqrt(x) and 1/x, so float target is 22 bits */
  72 #define GMX_SIMD_MATH_TARGET_SINGLE_BITS 22
  73
  74 /*! \brief We accept "double" that has 2x single precision - 44 bits.
  75  *
  76  * This way two Newton-Raphson iterations will suffice in double precision.
  77  */
  78 #define GMX_SIMD_MATH_TARGET_DOUBLE_BITS 44
  79
  80 /*! \} */
  81
  82 #ifdef GMX_SIMD_HAVE_FLOAT
  83
  84 /*! \name Single precision SIMD math functions
  85  *
  86  *  \note In most cases you should use the real-precision functions instead.
  87  *  \{
  88  */
  89
  90 /****************************************
  91  * SINGLE PRECISION SIMD MATH FUNCTIONS *
  92  ****************************************/
  93
  94 /*! \brief SIMD float utility to sum a+b+c+d.
  95  *
  96  * You should normally call the real-precision routine \ref gmx_simd_sum4_r.
  97  *
  98  * \param a term 1 (multiple values)
  99  * \param b term 2 (multiple values)
 100  * \param c term 3 (multiple values)
 101  * \param d term 4 (multiple values)
 102  * \return sum of terms 1-4 (multiple values)
 103  */
 104 static gmx_inline gmx_simd_float_t
 105 gmx_simd_sum4_f(gmx_simd_float_t a, gmx_simd_float_t b,
 106                 gmx_simd_float_t c, gmx_simd_float_t d)
 107 {
 108     return gmx_simd_add_f(gmx_simd_add_f(a, b), gmx_simd_add_f(c, d));
 109 }
 110
 111 /*! \brief Return -a if b is negative, SIMD float.
 112  *
 113  * You should normally call the real-precision routine \ref gmx_simd_xor_sign_r.
 114  *
 115  * \param a Values to set sign for
 116  * \param b Values used to set sign
 117  * \return if b is negative, the sign of a will be changed.
 118  *
 119  * This is equivalent to doing an xor operation on a with the sign bit of b,
 120  * with the exception that negative zero is not considered to be negative
 121  * on architectures where \ref GMX_SIMD_HAVE_LOGICAL is not set.
 122  */
 123 static gmx_inline gmx_simd_float_t
 124 gmx_simd_xor_sign_f(gmx_simd_float_t a, gmx_simd_float_t b)
 125 {
 126 #ifdef GMX_SIMD_HAVE_LOGICAL
 127     return gmx_simd_xor_f(a, gmx_simd_and_f(gmx_simd_set1_f(-0.0), b));
 128 #else
 129     return gmx_simd_blendv_f(a, gmx_simd_fneg_f(a), gmx_simd_cmplt_f(b, gmx_simd_setzero_f()));
 130 #endif
 131 }
 132
 133 /*! \brief Perform one Newton-Raphson iteration to improve 1/sqrt(x) for SIMD float.
 134  *
 135  * This is a low-level routine that should only be used by SIMD math routine
 136  * that evaluates the inverse square root.
 137  *
 138  *  \param lu Approximation of 1/sqrt(x), typically obtained from lookup.
 139  *  \param x  The reference (starting) value x for which we want 1/sqrt(x).
 140  *  \return   An improved approximation with roughly twice as many bits of accuracy.
 141  */
 142 static gmx_inline gmx_simd_float_t
 143 gmx_simd_rsqrt_iter_f(gmx_simd_float_t lu, gmx_simd_float_t x)
 144 {
 145 #    ifdef GMX_SIMD_HAVE_FMA
 146     return gmx_simd_fmadd_f(gmx_simd_fnmadd_f(x, gmx_simd_mul_f(lu, lu), gmx_simd_set1_f(1.0f)), gmx_simd_mul_f(lu, gmx_simd_set1_f(0.5f)), lu);
 147 #    else
 148     return gmx_simd_mul_f(gmx_simd_set1_f(0.5f), gmx_simd_mul_f(gmx_simd_sub_f(gmx_simd_set1_f(3.0f), gmx_simd_mul_f(gmx_simd_mul_f(lu, lu), x)), lu));
 149 #    endif
 150 }
 151
 152 /*! \brief Calculate 1/sqrt(x) for SIMD float.
 153  *
 154  * You should normally call the real-precision routine \ref gmx_simd_invsqrt_r.
 155  *
 156  *  \param x Argument that must be >0. This routine does not check arguments.
 157  *  \return 1/sqrt(x). Result is undefined if your argument was invalid.
 158  */
 159 static gmx_inline gmx_simd_float_t
 160 gmx_simd_invsqrt_f(gmx_simd_float_t x)
 161 {
 162     gmx_simd_float_t lu = gmx_simd_rsqrt_f(x);
 163 #if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
 164     lu = gmx_simd_rsqrt_iter_f(lu, x);
 165 #endif
 166 #if (GMX_SIMD_RSQRT_BITS*2 < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
 167     lu = gmx_simd_rsqrt_iter_f(lu, x);
 168 #endif
 169 #if (GMX_SIMD_RSQRT_BITS*4 < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
 170     lu = gmx_simd_rsqrt_iter_f(lu, x);
 171 #endif
 172     return lu;
 173 }
 174
 175 /*! \brief Calculate 1/sqrt(x) for two SIMD floats.
 176  *
 177  * You should normally call the real-precision routine \ref gmx_simd_invsqrt_pair_r.
 178  *
 179  * \param x0  First set of arguments, x0 must be positive - no argument checking.
 180  * \param x1  Second set of arguments, x1 must be positive - no argument checking.
 181  * \param[out] out0  Result 1/sqrt(x0)
 182  * \param[out] out1  Result 1/sqrt(x1)
 183  *
 184  *  In particular for double precision we can sometimes calculate square root
 185  *  pairs slightly faster by using single precision until the very last step.
 186  */
 187 static gmx_inline void
 188 gmx_simd_invsqrt_pair_f(gmx_simd_float_t x0,    gmx_simd_float_t x1,
 189                         gmx_simd_float_t *out0, gmx_simd_float_t *out1)
 190 {
 191     *out0 = gmx_simd_invsqrt_f(x0);
 192     *out1 = gmx_simd_invsqrt_f(x1);
 193 }
 194
 195 /*! \brief Perform one Newton-Raphson iteration to improve 1/x for SIMD float.
 196  *
 197  * This is a low-level routine that should only be used by SIMD math routine
 198  * that evaluates the reciprocal.
 199  *
 200  *  \param lu Approximation of 1/x, typically obtained from lookup.
 201  *  \param x  The reference (starting) value x for which we want 1/x.
 202  *  \return   An improved approximation with roughly twice as many bits of accuracy.
 203  */
 204 static gmx_inline gmx_simd_float_t
 205 gmx_simd_rcp_iter_f(gmx_simd_float_t lu, gmx_simd_float_t x)
 206 {
 207     return gmx_simd_mul_f(lu, gmx_simd_fnmadd_f(lu, x, gmx_simd_set1_f(2.0f)));
 208 }
 209
 210 /*! \brief Calculate 1/x for SIMD float.
 211  *
 212  * You should normally call the real-precision routine \ref gmx_simd_inv_r.
 213  *
 214  *  \param x Argument that must be nonzero. This routine does not check arguments.
 215  *  \return 1/x. Result is undefined if your argument was invalid.
 216  */
 217 static gmx_inline gmx_simd_float_t
 218 gmx_simd_inv_f(gmx_simd_float_t x)
 219 {
 220     gmx_simd_float_t lu = gmx_simd_rcp_f(x);
 221 #if (GMX_SIMD_RCP_BITS < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
 222     lu = gmx_simd_rcp_iter_f(lu, x);
 223 #endif
 224 #if (GMX_SIMD_RCP_BITS*2 < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
 225     lu = gmx_simd_rcp_iter_f(lu, x);
 226 #endif
 227 #if (GMX_SIMD_RCP_BITS*4 < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
 228     lu = gmx_simd_rcp_iter_f(lu, x);
 229 #endif
 230     return lu;
 231 }
 232
 233 /*! \brief Calculate sqrt(x) correctly for SIMD floats, including argument 0.0.
 234  *
 235  * You should normally call the real-precision routine \ref gmx_simd_sqrt_r.
 236  *
 237  *  \param x Argument that must be >=0.
 238  *  \return sqrt(x). If x=0, the result will correctly be set to 0.
 239  *          The result is undefined if the input value is negative.
 240  */
 241 static gmx_inline gmx_simd_float_t
 242 gmx_simd_sqrt_f(gmx_simd_float_t x)
 243 {
 244     gmx_simd_fbool_t  mask;
 245     gmx_simd_float_t  res;
 246
 247     mask = gmx_simd_cmpeq_f(x, gmx_simd_setzero_f());
 248     res  = gmx_simd_blendnotzero_f(gmx_simd_invsqrt_f(x), mask);
 249     return gmx_simd_mul_f(res, x);
 250 }
 251
 252 /*! \brief SIMD float log(x). This is the natural logarithm.
 253  *
 254  * You should normally call the real-precision routine \ref gmx_simd_log_r.
 255  *
 256  * \param x Argument, should be >0.
 257  * \result The natural logarithm of x. Undefined if argument is invalid.
 258  */
 259 #ifndef gmx_simd_log_f
 260 static gmx_inline gmx_simd_float_t
 261 gmx_simd_log_f(gmx_simd_float_t x)
 262 {
 263     const gmx_simd_float_t  half       = gmx_simd_set1_f(0.5f);
 264     const gmx_simd_float_t  one        = gmx_simd_set1_f(1.0f);
 265     const gmx_simd_float_t  sqrt2      = gmx_simd_set1_f(sqrt(2.0f));
 266     const gmx_simd_float_t  corr       = gmx_simd_set1_f(0.693147180559945286226764f);
 267     const gmx_simd_float_t  CL9        = gmx_simd_set1_f(0.2371599674224853515625f);
 268     const gmx_simd_float_t  CL7        = gmx_simd_set1_f(0.285279005765914916992188f);
 269     const gmx_simd_float_t  CL5        = gmx_simd_set1_f(0.400005519390106201171875f);
 270     const gmx_simd_float_t  CL3        = gmx_simd_set1_f(0.666666567325592041015625f);
 271     const gmx_simd_float_t  CL1        = gmx_simd_set1_f(2.0f);
 272     gmx_simd_float_t        fexp, x2, p;
 273     gmx_simd_fbool_t        mask;
 274
 275     fexp  = gmx_simd_get_exponent_f(x);
 276     x     = gmx_simd_get_mantissa_f(x);
 277
 278     mask  = gmx_simd_cmplt_f(sqrt2, x);
 279     /* Adjust to non-IEEE format for x>sqrt(2): exponent += 1, mantissa *= 0.5 */
 280     fexp  = gmx_simd_add_f(fexp, gmx_simd_blendzero_f(one, mask));
 281     x     = gmx_simd_mul_f(x, gmx_simd_blendv_f(one, half, mask));
 282
 283     x     = gmx_simd_mul_f( gmx_simd_sub_f(x, one), gmx_simd_inv_f( gmx_simd_add_f(x, one) ) );
 284     x2    = gmx_simd_mul_f(x, x);
 285
 286     p     = gmx_simd_fmadd_f(CL9, x2, CL7);
 287     p     = gmx_simd_fmadd_f(p, x2, CL5);
 288     p     = gmx_simd_fmadd_f(p, x2, CL3);
 289     p     = gmx_simd_fmadd_f(p, x2, CL1);
 290     p     = gmx_simd_fmadd_f(p, x, gmx_simd_mul_f(corr, fexp));
 291
 292     return p;
 293 }
 294 #endif
 295
 296 #ifndef gmx_simd_exp2_f
 297 /*! \brief SIMD float 2^x.
 298  *
 299  * You should normally call the real-precision routine \ref gmx_simd_exp2_r.
 300  *
 301  * \param x Argument.
 302  * \result 2^x. Undefined if input argument caused overflow.
 303  */
 304 static gmx_inline gmx_simd_float_t
 305 gmx_simd_exp2_f(gmx_simd_float_t x)
 306 {
 307     /* Lower bound: Disallow numbers that would lead to an IEEE fp exponent reaching +-127. */
 308     const gmx_simd_float_t  arglimit = gmx_simd_set1_f(126.0f);
 309     const gmx_simd_float_t  CC6      = gmx_simd_set1_f(0.0001534581200287996416911311);
 310     const gmx_simd_float_t  CC5      = gmx_simd_set1_f(0.001339993121934088894618990);
 311     const gmx_simd_float_t  CC4      = gmx_simd_set1_f(0.009618488957115180159497841);
 312     const gmx_simd_float_t  CC3      = gmx_simd_set1_f(0.05550328776964726865751735);
 313     const gmx_simd_float_t  CC2      = gmx_simd_set1_f(0.2402264689063408646490722);
 314     const gmx_simd_float_t  CC1      = gmx_simd_set1_f(0.6931472057372680777553816);
 315     const gmx_simd_float_t  one      = gmx_simd_set1_f(1.0f);
 316
 317     gmx_simd_float_t        fexppart;
 318     gmx_simd_float_t        intpart;
 319     gmx_simd_float_t        p;
 320     gmx_simd_fbool_t        valuemask;
 321
 322     fexppart  = gmx_simd_set_exponent_f(x);
 323     intpart   = gmx_simd_round_f(x);
 324     valuemask = gmx_simd_cmple_f(gmx_simd_fabs_f(x), arglimit);
 325     fexppart  = gmx_simd_blendzero_f(fexppart, valuemask);
 326     x         = gmx_simd_sub_f(x, intpart);
 327
 328     p         = gmx_simd_fmadd_f(CC6, x, CC5);
 329     p         = gmx_simd_fmadd_f(p, x, CC4);
 330     p         = gmx_simd_fmadd_f(p, x, CC3);
 331     p         = gmx_simd_fmadd_f(p, x, CC2);
 332     p         = gmx_simd_fmadd_f(p, x, CC1);
 333     p         = gmx_simd_fmadd_f(p, x, one);
 334     x         = gmx_simd_mul_f(p, fexppart);
 335     return x;
 336 }
 337 #endif
 338
 339 #ifndef gmx_simd_exp_f
 340 /*! \brief SIMD float exp(x).
 341  *
 342  * You should normally call the real-precision routine \ref gmx_simd_exp_r.
 343  *
 344  * In addition to scaling the argument for 2^x this routine correctly does
 345  * extended precision arithmetics to improve accuracy.
 346  *
 347  * \param x Argument.
 348  * \result exp(x). Undefined if input argument caused overflow.
 349  */
 350 static gmx_inline gmx_simd_float_t
 351 gmx_simd_exp_f(gmx_simd_float_t x)
 352 {
 353     const gmx_simd_float_t  argscale     = gmx_simd_set1_f(1.44269504088896341f);
 354     /* Lower bound: Disallow numbers that would lead to an IEEE fp exponent reaching +-127. */
 355     const gmx_simd_float_t  arglimit     = gmx_simd_set1_f(126.0f);
 356     const gmx_simd_float_t  invargscale0 = gmx_simd_set1_f(0.693145751953125f);
 357     const gmx_simd_float_t  invargscale1 = gmx_simd_set1_f(1.428606765330187045e-06f);
 358     const gmx_simd_float_t  CC4          = gmx_simd_set1_f(0.00136324646882712841033936f);
 359     const gmx_simd_float_t  CC3          = gmx_simd_set1_f(0.00836596917361021041870117f);
 360     const gmx_simd_float_t  CC2          = gmx_simd_set1_f(0.0416710823774337768554688f);
 361     const gmx_simd_float_t  CC1          = gmx_simd_set1_f(0.166665524244308471679688f);
 362     const gmx_simd_float_t  CC0          = gmx_simd_set1_f(0.499999850988388061523438f);
 363     const gmx_simd_float_t  one          = gmx_simd_set1_f(1.0f);
 364     gmx_simd_float_t        fexppart;
 365     gmx_simd_float_t        intpart;
 366     gmx_simd_float_t        y, p;
 367     gmx_simd_fbool_t        valuemask;
 368
 369     y         = gmx_simd_mul_f(x, argscale);
 370     fexppart  = gmx_simd_set_exponent_f(y);  /* rounds to nearest int internally */
 371     intpart   = gmx_simd_round_f(y);         /* use same rounding algorithm here */
 372     valuemask = gmx_simd_cmple_f(gmx_simd_fabs_f(y), arglimit);
 373     fexppart  = gmx_simd_blendzero_f(fexppart, valuemask);
 374
 375     /* Extended precision arithmetics */
 376     x         = gmx_simd_fnmadd_f(invargscale0, intpart, x);
 377     x         = gmx_simd_fnmadd_f(invargscale1, intpart, x);
 378
 379     p         = gmx_simd_fmadd_f(CC4, x, CC3);
 380     p         = gmx_simd_fmadd_f(p, x, CC2);
 381     p         = gmx_simd_fmadd_f(p, x, CC1);
 382     p         = gmx_simd_fmadd_f(p, x, CC0);
 383     p         = gmx_simd_fmadd_f(gmx_simd_mul_f(x, x), p, x);
 384     p         = gmx_simd_add_f(p, one);
 385     x         = gmx_simd_mul_f(p, fexppart);
 386     return x;
 387 }
 388 #endif
 389
 390 /*! \brief SIMD float erf(x).
 391  *
 392  * You should normally call the real-precision routine \ref gmx_simd_erf_r.
 393  *
 394  * \param x The value to calculate erf(x) for.
 395  * \result erf(x)
 396  *
 397  * This routine achieves very close to full precision, but we do not care about
 398  * the last bit or the subnormal result range.
 399  */
 400 static gmx_inline gmx_simd_float_t
 401 gmx_simd_erf_f(gmx_simd_float_t x)
 402 {
 403     /* Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1] */
 404     const gmx_simd_float_t  CA6      = gmx_simd_set1_f(7.853861353153693e-5f);
 405     const gmx_simd_float_t  CA5      = gmx_simd_set1_f(-8.010193625184903e-4f);
 406     const gmx_simd_float_t  CA4      = gmx_simd_set1_f(5.188327685732524e-3f);
 407     const gmx_simd_float_t  CA3      = gmx_simd_set1_f(-2.685381193529856e-2f);
 408     const gmx_simd_float_t  CA2      = gmx_simd_set1_f(1.128358514861418e-1f);
 409     const gmx_simd_float_t  CA1      = gmx_simd_set1_f(-3.761262582423300e-1f);
 410     const gmx_simd_float_t  CA0      = gmx_simd_set1_f(1.128379165726710f);
 411     /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2] */
 412     const gmx_simd_float_t  CB9      = gmx_simd_set1_f(-0.0018629930017603923f);
 413     const gmx_simd_float_t  CB8      = gmx_simd_set1_f(0.003909821287598495f);
 414     const gmx_simd_float_t  CB7      = gmx_simd_set1_f(-0.0052094582210355615f);
 415     const gmx_simd_float_t  CB6      = gmx_simd_set1_f(0.005685614362160572f);
 416     const gmx_simd_float_t  CB5      = gmx_simd_set1_f(-0.0025367682853477272f);
 417     const gmx_simd_float_t  CB4      = gmx_simd_set1_f(-0.010199799682318782f);
 418     const gmx_simd_float_t  CB3      = gmx_simd_set1_f(0.04369575504816542f);
 419     const gmx_simd_float_t  CB2      = gmx_simd_set1_f(-0.11884063474674492f);
 420     const gmx_simd_float_t  CB1      = gmx_simd_set1_f(0.2732120154030589f);
 421     const gmx_simd_float_t  CB0      = gmx_simd_set1_f(0.42758357702025784f);
 422     /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19] */
 423     const gmx_simd_float_t  CC10     = gmx_simd_set1_f(-0.0445555913112064f);
 424     const gmx_simd_float_t  CC9      = gmx_simd_set1_f(0.21376355144663348f);
 425     const gmx_simd_float_t  CC8      = gmx_simd_set1_f(-0.3473187200259257f);
 426     const gmx_simd_float_t  CC7      = gmx_simd_set1_f(0.016690861551248114f);
 427     const gmx_simd_float_t  CC6      = gmx_simd_set1_f(0.7560973182491192f);
 428     const gmx_simd_float_t  CC5      = gmx_simd_set1_f(-1.2137903600145787f);
 429     const gmx_simd_float_t  CC4      = gmx_simd_set1_f(0.8411872321232948f);
 430     const gmx_simd_float_t  CC3      = gmx_simd_set1_f(-0.08670413896296343f);
 431     const gmx_simd_float_t  CC2      = gmx_simd_set1_f(-0.27124782687240334f);
 432     const gmx_simd_float_t  CC1      = gmx_simd_set1_f(-0.0007502488047806069f);
 433     const gmx_simd_float_t  CC0      = gmx_simd_set1_f(0.5642114853803148f);
 434     const gmx_simd_float_t  one      = gmx_simd_set1_f(1.0f);
 435     const gmx_simd_float_t  two      = gmx_simd_set1_f(2.0f);
 436
 437     gmx_simd_float_t        x2, x4, y;
 438     gmx_simd_float_t        t, t2, w, w2;
 439     gmx_simd_float_t        pA0, pA1, pB0, pB1, pC0, pC1;
 440     gmx_simd_float_t        expmx2;
 441     gmx_simd_float_t        res_erf, res_erfc, res;
 442     gmx_simd_fbool_t        mask;
 443
 444     /* Calculate erf() */
 445     x2   = gmx_simd_mul_f(x, x);
 446     x4   = gmx_simd_mul_f(x2, x2);
 447
 448     pA0  = gmx_simd_fmadd_f(CA6, x4, CA4);
 449     pA1  = gmx_simd_fmadd_f(CA5, x4, CA3);
 450     pA0  = gmx_simd_fmadd_f(pA0, x4, CA2);
 451     pA1  = gmx_simd_fmadd_f(pA1, x4, CA1);
 452     pA0  = gmx_simd_mul_f(pA0, x4);
 453     pA0  = gmx_simd_fmadd_f(pA1, x2, pA0);
 454     /* Constant term must come last for precision reasons */
 455     pA0  = gmx_simd_add_f(pA0, CA0);
 456
 457     res_erf = gmx_simd_mul_f(x, pA0);
 458
 459     /* Calculate erfc */
 460     y       = gmx_simd_fabs_f(x);
 461     t       = gmx_simd_inv_f(y);
 462     w       = gmx_simd_sub_f(t, one);
 463     t2      = gmx_simd_mul_f(t, t);
 464     w2      = gmx_simd_mul_f(w, w);
 465
 466     /* No need for a floating-point sieve here (as in erfc), since erf()
 467      * will never return values that are extremely small for large args.
 468      */
 469     expmx2  = gmx_simd_exp_f( gmx_simd_fneg_f( gmx_simd_mul_f(y, y)));
 470
 471     pB1  = gmx_simd_fmadd_f(CB9, w2, CB7);
 472     pB0  = gmx_simd_fmadd_f(CB8, w2, CB6);
 473     pB1  = gmx_simd_fmadd_f(pB1, w2, CB5);
 474     pB0  = gmx_simd_fmadd_f(pB0, w2, CB4);
 475     pB1  = gmx_simd_fmadd_f(pB1, w2, CB3);
 476     pB0  = gmx_simd_fmadd_f(pB0, w2, CB2);
 477     pB1  = gmx_simd_fmadd_f(pB1, w2, CB1);
 478     pB0  = gmx_simd_fmadd_f(pB0, w2, CB0);
 479     pB0  = gmx_simd_fmadd_f(pB1, w, pB0);
 480
 481     pC0  = gmx_simd_fmadd_f(CC10, t2, CC8);
 482     pC1  = gmx_simd_fmadd_f(CC9, t2, CC7);
 483     pC0  = gmx_simd_fmadd_f(pC0, t2, CC6);
 484     pC1  = gmx_simd_fmadd_f(pC1, t2, CC5);
 485     pC0  = gmx_simd_fmadd_f(pC0, t2, CC4);
 486     pC1  = gmx_simd_fmadd_f(pC1, t2, CC3);
 487     pC0  = gmx_simd_fmadd_f(pC0, t2, CC2);
 488     pC1  = gmx_simd_fmadd_f(pC1, t2, CC1);
 489
 490     pC0  = gmx_simd_fmadd_f(pC0, t2, CC0);
 491     pC0  = gmx_simd_fmadd_f(pC1, t, pC0);
 492     pC0  = gmx_simd_mul_f(pC0, t);
 493
 494     /* SELECT pB0 or pC0 for erfc() */
 495     mask     = gmx_simd_cmplt_f(two, y);
 496     res_erfc = gmx_simd_blendv_f(pB0, pC0, mask);
 497     res_erfc = gmx_simd_mul_f(res_erfc, expmx2);
 498
 499     /* erfc(x<0) = 2-erfc(|x|) */
 500     mask     = gmx_simd_cmplt_f(x, gmx_simd_setzero_f());
 501     res_erfc = gmx_simd_blendv_f(res_erfc, gmx_simd_sub_f(two, res_erfc), mask);
 502
 503     /* Select erf() or erfc() */
 504     mask = gmx_simd_cmplt_f(y, gmx_simd_set1_f(0.75f));
 505     res  = gmx_simd_blendv_f(gmx_simd_sub_f(one, res_erfc), res_erf, mask);
 506
 507     return res;
 508 }
 509
 510 /*! \brief SIMD float erfc(x).
 511  *
 512  * You should normally call the real-precision routine \ref gmx_simd_erfc_r.
 513  *
 514  * \param x The value to calculate erfc(x) for.
 515  * \result erfc(x)
 516  *
 517  * This routine achieves full precision (bar the last bit) over most of the
 518  * input range, but for large arguments where the result is getting close
 519  * to the minimum representable numbers we accept slightly larger errors
 520  * (think results that are in the ballpark of 10^-30 for single precision,
 521  * or 10^-200 for double) since that is not relevant for MD.
 522  */
 523 static gmx_inline gmx_simd_float_t
 524 gmx_simd_erfc_f(gmx_simd_float_t x)
 525 {
 526     /* Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1] */
 527     const gmx_simd_float_t  CA6      = gmx_simd_set1_f(7.853861353153693e-5f);
 528     const gmx_simd_float_t  CA5      = gmx_simd_set1_f(-8.010193625184903e-4f);
 529     const gmx_simd_float_t  CA4      = gmx_simd_set1_f(5.188327685732524e-3f);
 530     const gmx_simd_float_t  CA3      = gmx_simd_set1_f(-2.685381193529856e-2f);
 531     const gmx_simd_float_t  CA2      = gmx_simd_set1_f(1.128358514861418e-1f);
 532     const gmx_simd_float_t  CA1      = gmx_simd_set1_f(-3.761262582423300e-1f);
 533     const gmx_simd_float_t  CA0      = gmx_simd_set1_f(1.128379165726710f);
 534     /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2] */
 535     const gmx_simd_float_t  CB9      = gmx_simd_set1_f(-0.0018629930017603923f);
 536     const gmx_simd_float_t  CB8      = gmx_simd_set1_f(0.003909821287598495f);
 537     const gmx_simd_float_t  CB7      = gmx_simd_set1_f(-0.0052094582210355615f);
 538     const gmx_simd_float_t  CB6      = gmx_simd_set1_f(0.005685614362160572f);
 539     const gmx_simd_float_t  CB5      = gmx_simd_set1_f(-0.0025367682853477272f);
 540     const gmx_simd_float_t  CB4      = gmx_simd_set1_f(-0.010199799682318782f);
 541     const gmx_simd_float_t  CB3      = gmx_simd_set1_f(0.04369575504816542f);
 542     const gmx_simd_float_t  CB2      = gmx_simd_set1_f(-0.11884063474674492f);
 543     const gmx_simd_float_t  CB1      = gmx_simd_set1_f(0.2732120154030589f);
 544     const gmx_simd_float_t  CB0      = gmx_simd_set1_f(0.42758357702025784f);
 545     /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19] */
 546     const gmx_simd_float_t  CC10     = gmx_simd_set1_f(-0.0445555913112064f);
 547     const gmx_simd_float_t  CC9      = gmx_simd_set1_f(0.21376355144663348f);
 548     const gmx_simd_float_t  CC8      = gmx_simd_set1_f(-0.3473187200259257f);
 549     const gmx_simd_float_t  CC7      = gmx_simd_set1_f(0.016690861551248114f);
 550     const gmx_simd_float_t  CC6      = gmx_simd_set1_f(0.7560973182491192f);
 551     const gmx_simd_float_t  CC5      = gmx_simd_set1_f(-1.2137903600145787f);
 552     const gmx_simd_float_t  CC4      = gmx_simd_set1_f(0.8411872321232948f);
 553     const gmx_simd_float_t  CC3      = gmx_simd_set1_f(-0.08670413896296343f);
 554     const gmx_simd_float_t  CC2      = gmx_simd_set1_f(-0.27124782687240334f);
 555     const gmx_simd_float_t  CC1      = gmx_simd_set1_f(-0.0007502488047806069f);
 556     const gmx_simd_float_t  CC0      = gmx_simd_set1_f(0.5642114853803148f);
 557     /* Coefficients for expansion of exp(x) in [0,0.1] */
 558     /* CD0 and CD1 are both 1.0, so no need to declare them separately */
 559     const gmx_simd_float_t  CD2      = gmx_simd_set1_f(0.5000066608081202f);
 560     const gmx_simd_float_t  CD3      = gmx_simd_set1_f(0.1664795422874624f);
 561     const gmx_simd_float_t  CD4      = gmx_simd_set1_f(0.04379839977652482f);
 562     const gmx_simd_float_t  one      = gmx_simd_set1_f(1.0f);
 563     const gmx_simd_float_t  two      = gmx_simd_set1_f(2.0f);
 564
 565     /* We need to use a small trick here, since we cannot assume all SIMD
 566      * architectures support integers, and the flag we want (0xfffff000) would
 567      * evaluate to NaN (i.e., it cannot be expressed as a floating-point num).
 568      * Instead, we represent the flags 0xf0f0f000 and 0x0f0f0000 as valid
 569      * fp numbers, and perform a logical or. Since the expression is constant,
 570      * we can at least hope it is evaluated at compile-time.
 571      */
 572 #ifdef GMX_SIMD_HAVE_LOGICAL
 573     const gmx_simd_float_t  sieve    = gmx_simd_or_f(gmx_simd_set1_f(-5.965323564e+29f), gmx_simd_set1_f(7.05044434e-30f));
 574 #else
 575     const int               isieve   = 0xFFFFF000;
 576     float                   mem[GMX_SIMD_REAL_WIDTH*2];
 577     float *                 pmem = gmx_simd_align_f(mem);
 578     union {
 579         float f; int i;
 580     } conv;
 581     int                     i;
 582 #endif
 583
 584     gmx_simd_float_t        x2, x4, y;
 585     gmx_simd_float_t        q, z, t, t2, w, w2;
 586     gmx_simd_float_t        pA0, pA1, pB0, pB1, pC0, pC1;
 587     gmx_simd_float_t        expmx2, corr;
 588     gmx_simd_float_t        res_erf, res_erfc, res;
 589     gmx_simd_fbool_t        mask;
 590
 591     /* Calculate erf() */
 592     x2     = gmx_simd_mul_f(x, x);
 593     x4     = gmx_simd_mul_f(x2, x2);
 594
 595     pA0  = gmx_simd_fmadd_f(CA6, x4, CA4);
 596     pA1  = gmx_simd_fmadd_f(CA5, x4, CA3);
 597     pA0  = gmx_simd_fmadd_f(pA0, x4, CA2);
 598     pA1  = gmx_simd_fmadd_f(pA1, x4, CA1);
 599     pA1  = gmx_simd_mul_f(pA1, x2);
 600     pA0  = gmx_simd_fmadd_f(pA0, x4, pA1);
 601     /* Constant term must come last for precision reasons */
 602     pA0  = gmx_simd_add_f(pA0, CA0);
 603
 604     res_erf = gmx_simd_mul_f(x, pA0);
 605
 606     /* Calculate erfc */
 607     y       = gmx_simd_fabs_f(x);
 608     t       = gmx_simd_inv_f(y);
 609     w       = gmx_simd_sub_f(t, one);
 610     t2      = gmx_simd_mul_f(t, t);
 611     w2      = gmx_simd_mul_f(w, w);
 612     /*
 613      * We cannot simply calculate exp(-y2) directly in single precision, since
 614      * that will lose a couple of bits of precision due to the multiplication.
 615      * Instead, we introduce y=z+w, where the last 12 bits of precision are in w.
 616      * Then we get exp(-y2) = exp(-z2)*exp((z-y)*(z+y)).
 617      *
 618      * The only drawback with this is that it requires TWO separate exponential
 619      * evaluations, which would be horrible performance-wise. However, the argument
 620      * for the second exp() call is always small, so there we simply use a
 621      * low-order minimax expansion on [0,0.1].
 622      *
 623      * However, this neat idea requires support for logical ops (and) on
 624      * FP numbers, which some vendors decided isn't necessary in their SIMD
 625      * instruction sets (Hi, IBM VSX!). In principle we could use some tricks
 626      * in double, but we still need memory as a backup when that is not available,
 627      * and this case is rare enough that we go directly there...
 628      */
 629 #ifdef GMX_SIMD_HAVE_LOGICAL
 630     z       = gmx_simd_and_f(y, sieve);
 631 #else
 632     gmx_simd_store_f(pmem, y);
 633     for (i = 0; i < GMX_SIMD_FLOAT_WIDTH; i++)
 634     {
 635         conv.f  = pmem[i];
 636         conv.i  = conv.i & isieve;
 637         pmem[i] = conv.f;
 638     }
 639     z = gmx_simd_load_f(pmem);
 640 #endif
 641     q       = gmx_simd_mul_f( gmx_simd_sub_f(z, y), gmx_simd_add_f(z, y) );
 642     corr    = gmx_simd_fmadd_f(CD4, q, CD3);
 643     corr    = gmx_simd_fmadd_f(corr, q, CD2);
 644     corr    = gmx_simd_fmadd_f(corr, q, one);
 645     corr    = gmx_simd_fmadd_f(corr, q, one);
 646
 647     expmx2  = gmx_simd_exp_f( gmx_simd_fneg_f( gmx_simd_mul_f(z, z) ) );
 648     expmx2  = gmx_simd_mul_f(expmx2, corr);
 649
 650     pB1  = gmx_simd_fmadd_f(CB9, w2, CB7);
 651     pB0  = gmx_simd_fmadd_f(CB8, w2, CB6);
 652     pB1  = gmx_simd_fmadd_f(pB1, w2, CB5);
 653     pB0  = gmx_simd_fmadd_f(pB0, w2, CB4);
 654     pB1  = gmx_simd_fmadd_f(pB1, w2, CB3);
 655     pB0  = gmx_simd_fmadd_f(pB0, w2, CB2);
 656     pB1  = gmx_simd_fmadd_f(pB1, w2, CB1);
 657     pB0  = gmx_simd_fmadd_f(pB0, w2, CB0);
 658     pB0  = gmx_simd_fmadd_f(pB1, w, pB0);
 659
 660     pC0  = gmx_simd_fmadd_f(CC10, t2, CC8);
 661     pC1  = gmx_simd_fmadd_f(CC9, t2, CC7);
 662     pC0  = gmx_simd_fmadd_f(pC0, t2, CC6);
 663     pC1  = gmx_simd_fmadd_f(pC1, t2, CC5);
 664     pC0  = gmx_simd_fmadd_f(pC0, t2, CC4);
 665     pC1  = gmx_simd_fmadd_f(pC1, t2, CC3);
 666     pC0  = gmx_simd_fmadd_f(pC0, t2, CC2);
 667     pC1  = gmx_simd_fmadd_f(pC1, t2, CC1);
 668
 669     pC0  = gmx_simd_fmadd_f(pC0, t2, CC0);
 670     pC0  = gmx_simd_fmadd_f(pC1, t, pC0);
 671     pC0  = gmx_simd_mul_f(pC0, t);
 672
 673     /* SELECT pB0 or pC0 for erfc() */
 674     mask     = gmx_simd_cmplt_f(two, y);
 675     res_erfc = gmx_simd_blendv_f(pB0, pC0, mask);
 676     res_erfc = gmx_simd_mul_f(res_erfc, expmx2);
 677
 678     /* erfc(x<0) = 2-erfc(|x|) */
 679     mask     = gmx_simd_cmplt_f(x, gmx_simd_setzero_f());
 680     res_erfc = gmx_simd_blendv_f(res_erfc, gmx_simd_sub_f(two, res_erfc), mask);
 681
 682     /* Select erf() or erfc() */
 683     mask = gmx_simd_cmplt_f(y, gmx_simd_set1_f(0.75f));
 684     res  = gmx_simd_blendv_f(res_erfc, gmx_simd_sub_f(one, res_erf), mask);
 685
 686     return res;
 687 }
 688
 689 /*! \brief SIMD float sin \& cos.
 690  *
 691  * You should normally call the real-precision routine \ref gmx_simd_sincos_r.
 692  *
 693  * \param x The argument to evaluate sin/cos for
 694  * \param[out] sinval Sin(x)
 695  * \param[out] cosval Cos(x)
 696  *
 697  * This version achieves close to machine precision, but for very large
 698  * magnitudes of the argument we inherently begin to lose accuracy due to the
 699  * argument reduction, despite using extended precision arithmetics internally.
 700  */
 701 static gmx_inline void
 702 gmx_simd_sincos_f(gmx_simd_float_t x, gmx_simd_float_t *sinval, gmx_simd_float_t *cosval)
 703 {
 704     /* Constants to subtract Pi/4*x from y while minimizing precision loss */
 705     const gmx_simd_float_t  argred0         = gmx_simd_set1_f(1.5703125);
 706     const gmx_simd_float_t  argred1         = gmx_simd_set1_f(4.83751296997070312500e-04f);
 707     const gmx_simd_float_t  argred2         = gmx_simd_set1_f(7.54953362047672271729e-08f);
 708     const gmx_simd_float_t  argred3         = gmx_simd_set1_f(2.56334406825708960298e-12f);
 709     const gmx_simd_float_t  two_over_pi     = gmx_simd_set1_f(2.0f/M_PI);
 710     const gmx_simd_float_t  const_sin2      = gmx_simd_set1_f(-1.9515295891e-4f);
 711     const gmx_simd_float_t  const_sin1      = gmx_simd_set1_f( 8.3321608736e-3f);
 712     const gmx_simd_float_t  const_sin0      = gmx_simd_set1_f(-1.6666654611e-1f);
 713     const gmx_simd_float_t  const_cos2      = gmx_simd_set1_f( 2.443315711809948e-5f);
 714     const gmx_simd_float_t  const_cos1      = gmx_simd_set1_f(-1.388731625493765e-3f);
 715     const gmx_simd_float_t  const_cos0      = gmx_simd_set1_f( 4.166664568298827e-2f);
 716     const gmx_simd_float_t  half            = gmx_simd_set1_f(0.5f);
 717     const gmx_simd_float_t  one             = gmx_simd_set1_f(1.0f);
 718     gmx_simd_float_t        ssign, csign;
 719     gmx_simd_float_t        x2, y, z, psin, pcos, sss, ccc;
 720     gmx_simd_fbool_t        mask;
 721 #if (defined GMX_SIMD_HAVE_FINT32) && (defined GMX_SIMD_HAVE_FINT32_ARITHMETICS) && (defined GMX_SIMD_HAVE_LOGICAL)
 722     const gmx_simd_fint32_t ione            = gmx_simd_set1_fi(1);
 723     const gmx_simd_fint32_t itwo            = gmx_simd_set1_fi(2);
 724     gmx_simd_fint32_t       iy;
 725
 726     z       = gmx_simd_mul_f(x, two_over_pi);
 727     iy      = gmx_simd_cvt_f2i(z);
 728     y       = gmx_simd_round_f(z);
 729
 730     mask    = gmx_simd_cvt_fib2fb(gmx_simd_cmpeq_fi(gmx_simd_and_fi(iy, ione), gmx_simd_setzero_fi()));
 731     ssign   = gmx_simd_blendzero_f(gmx_simd_set1_f(-0.0f), gmx_simd_cvt_fib2fb(gmx_simd_cmpeq_fi(gmx_simd_and_fi(iy, itwo), itwo)));
 732     csign   = gmx_simd_blendzero_f(gmx_simd_set1_f(-0.0f), gmx_simd_cvt_fib2fb(gmx_simd_cmpeq_fi(gmx_simd_and_fi(gmx_simd_add_fi(iy, ione), itwo), itwo)));
 733 #else
 734     const gmx_simd_float_t  quarter         = gmx_simd_set1_f(0.25f);
 735     const gmx_simd_float_t  minusquarter    = gmx_simd_set1_f(-0.25f);
 736     gmx_simd_float_t        q;
 737     gmx_simd_fbool_t        m1, m2, m3;
 738
 739     /* The most obvious way to find the arguments quadrant in the unit circle
 740      * to calculate the sign is to use integer arithmetic, but that is not
 741      * present in all SIMD implementations. As an alternative, we have devised a
 742      * pure floating-point algorithm that uses truncation for argument reduction
 743      * so that we get a new value 0<=q<1 over the unit circle, and then
 744      * do floating-point comparisons with fractions. This is likely to be
 745      * slightly slower (~10%) due to the longer latencies of floating-point, so
 746      * we only use it when integer SIMD arithmetic is not present.
 747      */
 748     ssign   = x;
 749     x       = gmx_simd_fabs_f(x);
 750     /* It is critical that half-way cases are rounded down */
 751     z       = gmx_simd_fmadd_f(x, two_over_pi, half);
 752     y       = gmx_simd_trunc_f(z);
 753     q       = gmx_simd_mul_f(z, quarter);
 754     q       = gmx_simd_sub_f(q, gmx_simd_trunc_f(q));
 755     /* z now starts at 0.0 for x=-pi/4 (although neg. values cannot occur), and
 756      * then increased by 1.0 as x increases by 2*Pi, when it resets to 0.0.
 757      * This removes the 2*Pi periodicity without using any integer arithmetic.
 758      * First check if y had the value 2 or 3, set csign if true.
 759      */
 760     q       = gmx_simd_sub_f(q, half);
 761     /* If we have logical operations we can work directly on the signbit, which
 762      * saves instructions. Otherwise we need to represent signs as +1.0/-1.0.
 763      * Thus, if you are altering defines to debug alternative code paths, the
 764      * two GMX_SIMD_HAVE_LOGICAL sections in this routine must either both be
 765      * active or inactive - you will get errors if only one is used.
 766      */
 767 #    ifdef GMX_SIMD_HAVE_LOGICAL
 768     ssign   = gmx_simd_and_f(ssign, gmx_simd_set1_f(-0.0f));
 769     csign   = gmx_simd_andnot_f(q, gmx_simd_set1_f(-0.0f));
 770     ssign   = gmx_simd_xor_f(ssign, csign);
 771 #    else
 772     csign   = gmx_simd_xor_sign_f(gmx_simd_set1_f(-1.0f), q);
 773     // ALT: csign = gmx_simd_fneg_f(gmx_simd_copysign(gmx_simd_set1_f(1.0),q));
 774
 775     ssign   = gmx_simd_xor_sign_f(ssign, csign);    /* swap ssign if csign was set. */
 776 #    endif
 777     /* Check if y had value 1 or 3 (remember we subtracted 0.5 from q) */
 778     m1      = gmx_simd_cmplt_f(q, minusquarter);
 779     m2      = gmx_simd_cmple_f(gmx_simd_setzero_f(), q);
 780     m3      = gmx_simd_cmplt_f(q, quarter);
 781     m2      = gmx_simd_and_fb(m2, m3);
 782     mask    = gmx_simd_or_fb(m1, m2);
 783     /* where mask is FALSE, set sign. */
 784     csign   = gmx_simd_xor_sign_f(csign, gmx_simd_blendv_f(gmx_simd_set1_f(-1.0f), one, mask));
 785 #endif
 786     x       = gmx_simd_fnmadd_f(y, argred0, x);
 787     x       = gmx_simd_fnmadd_f(y, argred1, x);
 788     x       = gmx_simd_fnmadd_f(y, argred2, x);
 789     x       = gmx_simd_fnmadd_f(y, argred3, x);
 790     x2      = gmx_simd_mul_f(x, x);
 791
 792     psin    = gmx_simd_fmadd_f(const_sin2, x2, const_sin1);
 793     psin    = gmx_simd_fmadd_f(psin, x2, const_sin0);
 794     psin    = gmx_simd_fmadd_f(psin, gmx_simd_mul_f(x, x2), x);
 795     pcos    = gmx_simd_fmadd_f(const_cos2, x2, const_cos1);
 796     pcos    = gmx_simd_fmadd_f(pcos, x2, const_cos0);
 797     pcos    = gmx_simd_fmsub_f(pcos, x2, half);
 798     pcos    = gmx_simd_fmadd_f(pcos, x2, one);
 799
 800     sss     = gmx_simd_blendv_f(pcos, psin, mask);
 801     ccc     = gmx_simd_blendv_f(psin, pcos, mask);
 802     /* See comment for GMX_SIMD_HAVE_LOGICAL section above. */
 803 #ifdef GMX_SIMD_HAVE_LOGICAL
 804     *sinval = gmx_simd_xor_f(sss, ssign);
 805     *cosval = gmx_simd_xor_f(ccc, csign);
 806 #else
 807     *sinval = gmx_simd_xor_sign_f(sss, ssign);
 808     *cosval = gmx_simd_xor_sign_f(ccc, csign);
 809 #endif
 810 }
 811
 812 /*! \brief SIMD float sin(x).
 813  *
 814  * You should normally call the real-precision routine \ref gmx_simd_sin_r.
 815  *
 816  * \param x The argument to evaluate sin for
 817  * \result Sin(x)
 818  *
 819  * \attention Do NOT call both sin & cos if you need both results, since each of them
 820  * will then call \ref gmx_simd_sincos_r and waste a factor 2 in performance.
 821  */
 822 static gmx_inline gmx_simd_float_t
 823 gmx_simd_sin_f(gmx_simd_float_t x)
 824 {
 825     gmx_simd_float_t s, c;
 826     gmx_simd_sincos_f(x, &s, &c);
 827     return s;
 828 }
 829
 830 /*! \brief SIMD float cos(x).
 831  *
 832  * You should normally call the real-precision routine \ref gmx_simd_cos_r.
 833  *
 834  * \param x The argument to evaluate cos for
 835  * \result Cos(x)
 836  *
 837  * \attention Do NOT call both sin & cos if you need both results, since each of them
 838  * will then call \ref gmx_simd_sincos_r and waste a factor 2 in performance.
 839  */
 840 static gmx_inline gmx_simd_float_t
 841 gmx_simd_cos_f(gmx_simd_float_t x)
 842 {
 843     gmx_simd_float_t s, c;
 844     gmx_simd_sincos_f(x, &s, &c);
 845     return c;
 846 }
 847
 848 /*! \brief SIMD float tan(x).
 849  *
 850  * You should normally call the real-precision routine \ref gmx_simd_tan_r.
 851  *
 852  * \param x The argument to evaluate tan for
 853  * \result Tan(x)
 854  */
 855 static gmx_inline gmx_simd_float_t
 856 gmx_simd_tan_f(gmx_simd_float_t x)
 857 {
 858     const gmx_simd_float_t  argred0         = gmx_simd_set1_f(1.5703125);
 859     const gmx_simd_float_t  argred1         = gmx_simd_set1_f(4.83751296997070312500e-04f);
 860     const gmx_simd_float_t  argred2         = gmx_simd_set1_f(7.54953362047672271729e-08f);
 861     const gmx_simd_float_t  argred3         = gmx_simd_set1_f(2.56334406825708960298e-12f);
 862     const gmx_simd_float_t  two_over_pi     = gmx_simd_set1_f(2.0f/M_PI);
 863     const gmx_simd_float_t  CT6             = gmx_simd_set1_f(0.009498288995810566122993911);
 864     const gmx_simd_float_t  CT5             = gmx_simd_set1_f(0.002895755790837379295226923);
 865     const gmx_simd_float_t  CT4             = gmx_simd_set1_f(0.02460087336161924491836265);
 866     const gmx_simd_float_t  CT3             = gmx_simd_set1_f(0.05334912882656359828045988);
 867     const gmx_simd_float_t  CT2             = gmx_simd_set1_f(0.1333989091464957704418495);
 868     const gmx_simd_float_t  CT1             = gmx_simd_set1_f(0.3333307599244198227797507);
 869
 870     gmx_simd_float_t        x2, p, y, z;
 871     gmx_simd_fbool_t        mask;
 872
 873 #if (defined GMX_SIMD_HAVE_FINT32) && (defined GMX_SIMD_HAVE_FINT32_ARITHMETICS) && (defined GMX_SIMD_HAVE_LOGICAL)
 874     gmx_simd_fint32_t  iy;
 875     gmx_simd_fint32_t  ione = gmx_simd_set1_fi(1);
 876
 877     z       = gmx_simd_mul_f(x, two_over_pi);
 878     iy      = gmx_simd_cvt_f2i(z);
 879     y       = gmx_simd_round_f(z);
 880     mask    = gmx_simd_cvt_fib2fb(gmx_simd_cmpeq_fi(gmx_simd_and_fi(iy, ione), ione));
 881
 882     x       = gmx_simd_fnmadd_f(y, argred0, x);
 883     x       = gmx_simd_fnmadd_f(y, argred1, x);
 884     x       = gmx_simd_fnmadd_f(y, argred2, x);
 885     x       = gmx_simd_fnmadd_f(y, argred3, x);
 886     x       = gmx_simd_xor_f(gmx_simd_blendzero_f(gmx_simd_set1_f(-0.0f), mask), x);
 887 #else
 888     const gmx_simd_float_t  quarter         = gmx_simd_set1_f(0.25f);
 889     const gmx_simd_float_t  half            = gmx_simd_set1_f(0.5f);
 890     const gmx_simd_float_t  threequarter    = gmx_simd_set1_f(0.75f);
 891     gmx_simd_float_t        w, q;
 892     gmx_simd_fbool_t        m1, m2, m3;
 893
 894     w       = gmx_simd_fabs_f(x);
 895     z       = gmx_simd_fmadd_f(w, two_over_pi, half);
 896     y       = gmx_simd_trunc_f(z);
 897     q       = gmx_simd_mul_f(z, quarter);
 898     q       = gmx_simd_sub_f(q, gmx_simd_trunc_f(q));
 899     m1      = gmx_simd_cmple_f(quarter, q);
 900     m2      = gmx_simd_cmplt_f(q, half);
 901     m3      = gmx_simd_cmple_f(threequarter, q);
 902     m1      = gmx_simd_and_fb(m1, m2);
 903     mask    = gmx_simd_or_fb(m1, m3);
 904     w       = gmx_simd_fnmadd_f(y, argred0, w);
 905     w       = gmx_simd_fnmadd_f(y, argred1, w);
 906     w       = gmx_simd_fnmadd_f(y, argred2, w);
 907     w       = gmx_simd_fnmadd_f(y, argred3, w);
 908
 909     w       = gmx_simd_blendv_f(w, gmx_simd_fneg_f(w), mask);
 910     x       = gmx_simd_xor_sign_f(w, x);
 911 #endif
 912     x2      = gmx_simd_mul_f(x, x);
 913     p       = gmx_simd_fmadd_f(CT6, x2, CT5);
 914     p       = gmx_simd_fmadd_f(p, x2, CT4);
 915     p       = gmx_simd_fmadd_f(p, x2, CT3);
 916     p       = gmx_simd_fmadd_f(p, x2, CT2);
 917     p       = gmx_simd_fmadd_f(p, x2, CT1);
 918     p       = gmx_simd_fmadd_f(x2, gmx_simd_mul_f(p, x), x);
 919
 920     p       = gmx_simd_blendv_f( p, gmx_simd_inv_f(p), mask);
 921     return p;
 922 }
 923
 924 /*! \brief SIMD float asin(x).
 925  *
 926  * You should normally call the real-precision routine \ref gmx_simd_asin_r.
 927  *
 928  * \param x The argument to evaluate asin for
 929  * \result Asin(x)
 930  */
 931 static gmx_inline gmx_simd_float_t
 932 gmx_simd_asin_f(gmx_simd_float_t x)
 933 {
 934     const gmx_simd_float_t limitlow   = gmx_simd_set1_f(1e-4f);
 935     const gmx_simd_float_t half       = gmx_simd_set1_f(0.5f);
 936     const gmx_simd_float_t one        = gmx_simd_set1_f(1.0f);
 937     const gmx_simd_float_t halfpi     = gmx_simd_set1_f((float)M_PI/2.0f);
 938     const gmx_simd_float_t CC5        = gmx_simd_set1_f(4.2163199048E-2f);
 939     const gmx_simd_float_t CC4        = gmx_simd_set1_f(2.4181311049E-2f);
 940     const gmx_simd_float_t CC3        = gmx_simd_set1_f(4.5470025998E-2f);
 941     const gmx_simd_float_t CC2        = gmx_simd_set1_f(7.4953002686E-2f);
 942     const gmx_simd_float_t CC1        = gmx_simd_set1_f(1.6666752422E-1f);
 943     gmx_simd_float_t       xabs;
 944     gmx_simd_float_t       z, z1, z2, q, q1, q2;
 945     gmx_simd_float_t       pA, pB;
 946     gmx_simd_fbool_t       mask;
 947
 948     xabs  = gmx_simd_fabs_f(x);
 949     mask  = gmx_simd_cmplt_f(half, xabs);
 950     z1    = gmx_simd_mul_f(half, gmx_simd_sub_f(one, xabs));
 951     q1    = gmx_simd_mul_f(z1, gmx_simd_invsqrt_f(z1));
 952     q1    = gmx_simd_blendnotzero_f(q1, gmx_simd_cmpeq_f(xabs, one));
 953     q2    = xabs;
 954     z2    = gmx_simd_mul_f(q2, q2);
 955     z     = gmx_simd_blendv_f(z2, z1, mask);
 956     q     = gmx_simd_blendv_f(q2, q1, mask);
 957
 958     z2    = gmx_simd_mul_f(z, z);
 959     pA    = gmx_simd_fmadd_f(CC5, z2, CC3);
 960     pB    = gmx_simd_fmadd_f(CC4, z2, CC2);
 961     pA    = gmx_simd_fmadd_f(pA, z2, CC1);
 962     pA    = gmx_simd_mul_f(pA, z);
 963     z     = gmx_simd_fmadd_f(pB, z2, pA);
 964     z     = gmx_simd_fmadd_f(z, q, q);
 965     q2    = gmx_simd_sub_f(halfpi, z);
 966     q2    = gmx_simd_sub_f(q2, z);
 967     z     = gmx_simd_blendv_f(z, q2, mask);
 968
 969     mask  = gmx_simd_cmplt_f(limitlow, xabs);
 970     z     = gmx_simd_blendv_f( xabs, z, mask );
 971     z     = gmx_simd_xor_sign_f(z, x);
 972
 973     return z;
 974 }
 975
 976 /*! \brief SIMD float acos(x).
 977  *
 978  * You should normally call the real-precision routine \ref gmx_simd_acos_r.
 979  *
 980  * \param x The argument to evaluate acos for
 981  * \result Acos(x)
 982  */
 983 static gmx_inline gmx_simd_float_t
 984 gmx_simd_acos_f(gmx_simd_float_t x)
 985 {
 986     const gmx_simd_float_t one       = gmx_simd_set1_f(1.0f);
 987     const gmx_simd_float_t half      = gmx_simd_set1_f(0.5f);
 988     const gmx_simd_float_t pi        = gmx_simd_set1_f((float)M_PI);
 989     const gmx_simd_float_t halfpi    = gmx_simd_set1_f((float)M_PI/2.0f);
 990     gmx_simd_float_t       xabs;
 991     gmx_simd_float_t       z, z1, z2, z3;
 992     gmx_simd_fbool_t       mask1, mask2;
 993
 994     xabs  = gmx_simd_fabs_f(x);
 995     mask1 = gmx_simd_cmplt_f(half, xabs);
 996     mask2 = gmx_simd_cmplt_f(gmx_simd_setzero_f(), x);
 997
 998     z     = gmx_simd_mul_f(half, gmx_simd_sub_f(one, xabs));
 999     z     = gmx_simd_mul_f(z, gmx_simd_invsqrt_f(z));
1000     z     = gmx_simd_blendnotzero_f(z, gmx_simd_cmpeq_f(xabs, one));
1001     z     = gmx_simd_blendv_f(x, z, mask1);
1002     z     = gmx_simd_asin_f(z);
1003
1004     z2    = gmx_simd_add_f(z, z);
1005     z1    = gmx_simd_sub_f(pi, z2);
1006     z3    = gmx_simd_sub_f(halfpi, z);
1007     z     = gmx_simd_blendv_f(z1, z2, mask2);
1008     z     = gmx_simd_blendv_f(z3, z, mask1);
1009
1010     return z;
1011 }
1012
1013 /*! \brief SIMD float asin(x).
1014  *
1015  * You should normally call the real-precision routine \ref gmx_simd_atan_r.
1016  *
1017  * \param x The argument to evaluate atan for
1018  * \result Atan(x), same argument/value range as standard math library.
1019  */
1020 static gmx_inline gmx_simd_float_t
1021 gmx_simd_atan_f(gmx_simd_float_t x)
1022 {
1023     const gmx_simd_float_t halfpi    = gmx_simd_set1_f(M_PI/2);
1024     const gmx_simd_float_t CA17      = gmx_simd_set1_f(0.002823638962581753730774f);
1025     const gmx_simd_float_t CA15      = gmx_simd_set1_f(-0.01595690287649631500244f);
1026     const gmx_simd_float_t CA13      = gmx_simd_set1_f(0.04250498861074447631836f);
1027     const gmx_simd_float_t CA11      = gmx_simd_set1_f(-0.07489009201526641845703f);
1028     const gmx_simd_float_t CA9       = gmx_simd_set1_f(0.1063479334115982055664f);
1029     const gmx_simd_float_t CA7       = gmx_simd_set1_f(-0.1420273631811141967773f);
1030     const gmx_simd_float_t CA5       = gmx_simd_set1_f(0.1999269574880599975585f);
1031     const gmx_simd_float_t CA3       = gmx_simd_set1_f(-0.3333310186862945556640f);
1032     gmx_simd_float_t       x2, x3, x4, pA, pB;
1033     gmx_simd_fbool_t       mask, mask2;
1034
1035     mask  = gmx_simd_cmplt_f(x, gmx_simd_setzero_f());
1036     x     = gmx_simd_fabs_f(x);
1037     mask2 = gmx_simd_cmplt_f(gmx_simd_set1_f(1.0f), x);
1038     x     = gmx_simd_blendv_f(x, gmx_simd_inv_f(x), mask2);
1039
1040     x2    = gmx_simd_mul_f(x, x);
1041     x3    = gmx_simd_mul_f(x2, x);
1042     x4    = gmx_simd_mul_f(x2, x2);
1043     pA    = gmx_simd_fmadd_f(CA17, x4, CA13);
1044     pB    = gmx_simd_fmadd_f(CA15, x4, CA11);
1045     pA    = gmx_simd_fmadd_f(pA, x4, CA9);
1046     pB    = gmx_simd_fmadd_f(pB, x4, CA7);
1047     pA    = gmx_simd_fmadd_f(pA, x4, CA5);
1048     pB    = gmx_simd_fmadd_f(pB, x4, CA3);
1049     pA    = gmx_simd_fmadd_f(pA, x2, pB);
1050     pA    = gmx_simd_fmadd_f(pA, x3, x);
1051
1052     pA    = gmx_simd_blendv_f(pA, gmx_simd_sub_f(halfpi, pA), mask2);
1053     pA    = gmx_simd_blendv_f(pA, gmx_simd_fneg_f(pA), mask);
1054
1055     return pA;
1056 }
1057
1058 /*! \brief SIMD float atan2(y,x).
1059  *
1060  * You should normally call the real-precision routine \ref gmx_simd_atan2_r.
1061  *
1062  * \param y Y component of vector, any quartile
1063  * \param x X component of vector, any quartile
1064  * \result Atan(y,x), same argument/value range as standard math library.
1065  *
1066  * \note This routine should provide correct results for all finite
1067  * non-zero or positive-zero arguments. However, negative zero arguments will
1068  * be treated as positive zero, which means the return value will deviate from
1069  * the standard math library atan2(y,x) for those cases. That should not be
1070  * of any concern in Gromacs, and in particular it will not affect calculations
1071  * of angles from vectors.
1072  */
1073 static gmx_inline gmx_simd_float_t
1074 gmx_simd_atan2_f(gmx_simd_float_t y, gmx_simd_float_t x)
1075 {
1076     const gmx_simd_float_t pi          = gmx_simd_set1_f(M_PI);
1077     const gmx_simd_float_t halfpi      = gmx_simd_set1_f(M_PI/2.0);
1078     gmx_simd_float_t       xinv, p, aoffset;
1079     gmx_simd_fbool_t       mask_x0, mask_y0, mask_xlt0, mask_ylt0;
1080
1081     mask_x0   = gmx_simd_cmpeq_f(x, gmx_simd_setzero_f());
1082     mask_y0   = gmx_simd_cmpeq_f(y, gmx_simd_setzero_f());
1083     mask_xlt0 = gmx_simd_cmplt_f(x, gmx_simd_setzero_f());
1084     mask_ylt0 = gmx_simd_cmplt_f(y, gmx_simd_setzero_f());
1085
1086     aoffset   = gmx_simd_blendzero_f(halfpi, mask_x0);
1087     aoffset   = gmx_simd_blendnotzero_f(aoffset, mask_y0);
1088
1089     aoffset   = gmx_simd_blendv_f(aoffset, pi, mask_xlt0);
1090     aoffset   = gmx_simd_blendv_f(aoffset, gmx_simd_fneg_f(aoffset), mask_ylt0);
1091
1092     xinv      = gmx_simd_blendnotzero_f(gmx_simd_inv_f(x), mask_x0);
1093     p         = gmx_simd_mul_f(y, xinv);
1094     p         = gmx_simd_atan_f(p);
1095     p         = gmx_simd_add_f(p, aoffset);
1096
1097     return p;
1098 }
1099
1100 /*! \brief Calculate the force correction due to PME analytically in SIMD float.
1101  *
1102  * You should normally call the real-precision routine \ref gmx_simd_pmecorrF_r.
1103  *
1104  * \param z2 \f$(r \beta)^2\f$ - see below for details.
1105  * \result Correction factor to coulomb force - see below for details.
1106  *
1107  * This routine is meant to enable analytical evaluation of the
1108  * direct-space PME electrostatic force to avoid tables.
1109  *
1110  * The direct-space potential should be \f$ \mbox{erfc}(\beta r)/r\f$, but there
1111  * are some problems evaluating that:
1112  *
1113  * First, the error function is difficult (read: expensive) to
1114  * approxmiate accurately for intermediate to large arguments, and
1115  * this happens already in ranges of \f$(\beta r)\f$ that occur in simulations.
1116  * Second, we now try to avoid calculating potentials in Gromacs but
1117  * use forces directly.
1118  *
1119  * We can simply things slight by noting that the PME part is really
1120  * a correction to the normal Coulomb force since \f$\mbox{erfc}(z)=1-\mbox{erf}(z)\f$, i.e.
1121  * \f[
1122  * V = \frac{1}{r} - \frac{\mbox{erf}(\beta r)}{r}
1123  * \f]
1124  * The first term we already have from the inverse square root, so
1125  * that we can leave out of this routine.
1126  *
1127  * For pme tolerances of 1e-3 to 1e-8 and cutoffs of 0.5nm to 1.8nm,
1128  * the argument \f$beta r\f$ will be in the range 0.15 to ~4. Use your
1129  * favorite plotting program to realize how well-behaved \f$\frac{\mbox{erf}(z)}{z}\f$ is
1130  * in this range!
1131  *
1132  * We approximate \f$f(z)=\mbox{erf}(z)/z\f$ with a rational minimax polynomial.
1133  * However, it turns out it is more efficient to approximate \f$f(z)/z\f$ and
1134  * then only use even powers. This is another minor optimization, since
1135  * we actually \a want \f$f(z)/z\f$, because it is going to be multiplied by
1136  * the vector between the two atoms to get the vectorial force. The
1137  * fastest flops are the ones we can avoid calculating!
1138  *
1139  * So, here's how it should be used:
1140  *
1141  * 1. Calculate \f$r^2\f$.
1142  * 2. Multiply by \f$\beta^2\f$, so you get \f$z^2=(\beta r)^2\f$.
1143  * 3. Evaluate this routine with \f$z^2\f$ as the argument.
1144  * 4. The return value is the expression:
1145  *
1146  * \f[
1147  *    \frac{2 \exp{-z^2}}{\sqrt{\pi} z^2}-\frac{\mbox{erf}(z)}{z^3}
1148  * \f]
1149  *
1150  * 5. Multiply the entire expression by \f$\beta^3\f$. This will get you
1151  *
1152  *  \f[
1153  *    \frac{2 \beta^3 \exp(-z^2)}{\sqrt{\pi} z^2} - \frac{\beta^3 \mbox{erf}(z)}{z^3}
1154  *  \f]
1155  *
1156  *    or, switching back to \f$r\f$ (since \f$z=r \beta\f$):
1157  *
1158  *  \f[
1159  *    \frac{2 \beta \exp(-r^2 \beta^2)}{\sqrt{\pi} r^2} - \frac{\mbox{erf}(r \beta)}{r^3}
1160  *  \f]
1161  *
1162  *    With a bit of math exercise you should be able to confirm that
1163  *    this is exactly
1164  *
1165  *  \f[
1166  *   \frac{\frac{d}{dr}\left( \frac{\mbox{erf}(\beta r)}{r} \right)}{r}
1167  *  \f]
1168  *
1169  * 6. Add the result to \f$r^{-3}\f$, multiply by the product of the charges,
1170  *    and you have your force (divided by \f$r\f$). A final multiplication
1171  *    with the vector connecting the two particles and you have your
1172  *    vectorial force to add to the particles.
1173  *
1174  * This approximation achieves an accuracy slightly lower than 1e-6; when
1175  * added to \f$1/r\f$ the error will be insignificant.
1176  *
1177  */
1178 static gmx_simd_float_t
1179 gmx_simd_pmecorrF_f(gmx_simd_float_t z2)
1180 {
1181     const gmx_simd_float_t  FN6      = gmx_simd_set1_f(-1.7357322914161492954e-8f);
1182     const gmx_simd_float_t  FN5      = gmx_simd_set1_f(1.4703624142580877519e-6f);
1183     const gmx_simd_float_t  FN4      = gmx_simd_set1_f(-0.000053401640219807709149f);
1184     const gmx_simd_float_t  FN3      = gmx_simd_set1_f(0.0010054721316683106153f);
1185     const gmx_simd_float_t  FN2      = gmx_simd_set1_f(-0.019278317264888380590f);
1186     const gmx_simd_float_t  FN1      = gmx_simd_set1_f(0.069670166153766424023f);
1187     const gmx_simd_float_t  FN0      = gmx_simd_set1_f(-0.75225204789749321333f);
1188
1189     const gmx_simd_float_t  FD4      = gmx_simd_set1_f(0.0011193462567257629232f);
1190     const gmx_simd_float_t  FD3      = gmx_simd_set1_f(0.014866955030185295499f);
1191     const gmx_simd_float_t  FD2      = gmx_simd_set1_f(0.11583842382862377919f);
1192     const gmx_simd_float_t  FD1      = gmx_simd_set1_f(0.50736591960530292870f);
1193     const gmx_simd_float_t  FD0      = gmx_simd_set1_f(1.0f);
1194
1195     gmx_simd_float_t        z4;
1196     gmx_simd_float_t        polyFN0, polyFN1, polyFD0, polyFD1;
1197
1198     z4             = gmx_simd_mul_f(z2, z2);
1199
1200     polyFD0        = gmx_simd_fmadd_f(FD4, z4, FD2);
1201     polyFD1        = gmx_simd_fmadd_f(FD3, z4, FD1);
1202     polyFD0        = gmx_simd_fmadd_f(polyFD0, z4, FD0);
1203     polyFD0        = gmx_simd_fmadd_f(polyFD1, z2, polyFD0);
1204
1205     polyFD0        = gmx_simd_inv_f(polyFD0);
1206
1207     polyFN0        = gmx_simd_fmadd_f(FN6, z4, FN4);
1208     polyFN1        = gmx_simd_fmadd_f(FN5, z4, FN3);
1209     polyFN0        = gmx_simd_fmadd_f(polyFN0, z4, FN2);
1210     polyFN1        = gmx_simd_fmadd_f(polyFN1, z4, FN1);
1211     polyFN0        = gmx_simd_fmadd_f(polyFN0, z4, FN0);
1212     polyFN0        = gmx_simd_fmadd_f(polyFN1, z2, polyFN0);
1213
1214     return gmx_simd_mul_f(polyFN0, polyFD0);
1215 }
1216
1217
1218
1219 /*! \brief Calculate the potential correction due to PME analytically in SIMD float.
1220  *
1221  * You should normally call the real-precision routine \ref gmx_simd_pmecorrV_r.
1222  *
1223  * \param z2 \f$(r \beta)^2\f$ - see below for details.
1224  * \result Correction factor to coulomb potential - see below for details.
1225  *
1226  * See \ref gmx_simd_pmecorrF_f for details about the approximation.
1227  *
1228  * This routine calculates \f$\mbox{erf}(z)/z\f$, although you should provide \f$z^2\f$
1229  * as the input argument.
1230  *
1231  * Here's how it should be used:
1232  *
1233  * 1. Calculate \f$r^2\f$.
1234  * 2. Multiply by \f$\beta^2\f$, so you get \f$z^2=\beta^2*r^2\f$.
1235  * 3. Evaluate this routine with z^2 as the argument.
1236  * 4. The return value is the expression:
1237  *
1238  *  \f[
1239  *   \frac{\mbox{erf}(z)}{z}
1240  *  \f]
1241  *
1242  * 5. Multiply the entire expression by beta and switching back to \f$r\f$ (since \f$z=r \beta\f$):
1243  *
1244  *  \f[
1245  *    \frac{\mbox{erf}(r \beta)}{r}
1246  *  \f]
1247  *
1248  * 6. Subtract the result from \f$1/r\f$, multiply by the product of the charges,
1249  *    and you have your potential.
1250  *
1251  * This approximation achieves an accuracy slightly lower than 1e-6; when
1252  * added to \f$1/r\f$ the error will be insignificant.
1253  */
1254 static gmx_simd_float_t
1255 gmx_simd_pmecorrV_f(gmx_simd_float_t z2)
1256 {
1257     const gmx_simd_float_t  VN6      = gmx_simd_set1_f(1.9296833005951166339e-8f);
1258     const gmx_simd_float_t  VN5      = gmx_simd_set1_f(-1.4213390571557850962e-6f);
1259     const gmx_simd_float_t  VN4      = gmx_simd_set1_f(0.000041603292906656984871f);
1260     const gmx_simd_float_t  VN3      = gmx_simd_set1_f(-0.00013134036773265025626f);
1261     const gmx_simd_float_t  VN2      = gmx_simd_set1_f(0.038657983986041781264f);
1262     const gmx_simd_float_t  VN1      = gmx_simd_set1_f(0.11285044772717598220f);
1263     const gmx_simd_float_t  VN0      = gmx_simd_set1_f(1.1283802385263030286f);
1264
1265     const gmx_simd_float_t  VD3      = gmx_simd_set1_f(0.0066752224023576045451f);
1266     const gmx_simd_float_t  VD2      = gmx_simd_set1_f(0.078647795836373922256f);
1267     const gmx_simd_float_t  VD1      = gmx_simd_set1_f(0.43336185284710920150f);
1268     const gmx_simd_float_t  VD0      = gmx_simd_set1_f(1.0f);
1269
1270     gmx_simd_float_t        z4;
1271     gmx_simd_float_t        polyVN0, polyVN1, polyVD0, polyVD1;
1272
1273     z4             = gmx_simd_mul_f(z2, z2);
1274
1275     polyVD1        = gmx_simd_fmadd_f(VD3, z4, VD1);
1276     polyVD0        = gmx_simd_fmadd_f(VD2, z4, VD0);
1277     polyVD0        = gmx_simd_fmadd_f(polyVD1, z2, polyVD0);
1278
1279     polyVD0        = gmx_simd_inv_f(polyVD0);
1280
1281     polyVN0        = gmx_simd_fmadd_f(VN6, z4, VN4);
1282     polyVN1        = gmx_simd_fmadd_f(VN5, z4, VN3);
1283     polyVN0        = gmx_simd_fmadd_f(polyVN0, z4, VN2);
1284     polyVN1        = gmx_simd_fmadd_f(polyVN1, z4, VN1);
1285     polyVN0        = gmx_simd_fmadd_f(polyVN0, z4, VN0);
1286     polyVN0        = gmx_simd_fmadd_f(polyVN1, z2, polyVN0);
1287
1288     return gmx_simd_mul_f(polyVN0, polyVD0);
1289 }
1290 #endif
1291
1292 /*! \} */
1293
1294 #ifdef GMX_SIMD_HAVE_DOUBLE
1295
1296 /*! \name Double precision SIMD math functions
1297  *
1298  *  \note In most cases you should use the real-precision functions instead.
1299  *  \{
1300  */
1301
1302 /****************************************
1303  * DOUBLE PRECISION SIMD MATH FUNCTIONS *
1304  ****************************************/
1305
1306 /*! \brief SIMD utility function to sum a+b+c+d for SIMD doubles.
1307  *
1308  * \copydetails gmx_simd_sum4_f
1309  */
1310 static gmx_inline gmx_simd_double_t
1311 gmx_simd_sum4_d(gmx_simd_double_t a, gmx_simd_double_t b,
1312                 gmx_simd_double_t c, gmx_simd_double_t d)
1313 {
1314     return gmx_simd_add_d(gmx_simd_add_d(a, b), gmx_simd_add_d(c, d));
1315 }
1316
1317 /*! \brief Return -a if b is negative, SIMD double.
1318  *
1319  * You should normally call the real-precision routine \ref gmx_simd_xor_sign_r.
1320  *
1321  * \param a Values to set sign for
1322  * \param b Values used to set sign
1323  * \return if b is negative, the sign of a will be changed.
1324  *
1325  * This is equivalent to doing an xor operation on a with the sign bit of b,
1326  * with the exception that negative zero is not considered to be negative
1327  * on architectures where \ref GMX_SIMD_HAVE_LOGICAL is not set.
1328  */
1329 static gmx_inline gmx_simd_double_t
1330 gmx_simd_xor_sign_d(gmx_simd_double_t a, gmx_simd_double_t b)
1331 {
1332 #ifdef GMX_SIMD_HAVE_LOGICAL
1333     return gmx_simd_xor_d(a, gmx_simd_and_d(gmx_simd_set1_d(-0.0), b));
1334 #else
1335     return gmx_simd_blendv_d(a, gmx_simd_fneg_d(a), gmx_simd_cmplt_d(b, gmx_simd_setzero_d()));
1336 #endif
1337 }
1338
1339 /*! \brief Perform one Newton-Raphson iteration to improve 1/sqrt(x) for SIMD double.
1340  *
1341  * \copydetails gmx_simd_rsqrt_iter_f
1342  */
1343 static gmx_inline gmx_simd_double_t
1344 gmx_simd_rsqrt_iter_d(gmx_simd_double_t lu, gmx_simd_double_t x)
1345 {
1346 #ifdef GMX_SIMD_HAVE_FMA
1347     return gmx_simd_fmadd_d(gmx_simd_fnmadd_d(x, gmx_simd_mul_d(lu, lu), gmx_simd_set1_d(1.0)), gmx_simd_mul_d(lu, gmx_simd_set1_d(0.5)), lu);
1348 #else
1349     return gmx_simd_mul_d(gmx_simd_set1_d(0.5), gmx_simd_mul_d(gmx_simd_sub_d(gmx_simd_set1_d(3.0), gmx_simd_mul_d(gmx_simd_mul_d(lu, lu), x)), lu));
1350 #endif
1351 }
1352
1353
1354 /*! \brief Calculate 1/sqrt(x) for SIMD double
1355  *
1356  * \copydetails gmx_simd_invsqrt_f
1357  */
1358 static gmx_inline gmx_simd_double_t
1359 gmx_simd_invsqrt_d(gmx_simd_double_t x)
1360 {
1361     gmx_simd_double_t lu = gmx_simd_rsqrt_d(x);
1362 #if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1363     lu = gmx_simd_rsqrt_iter_d(lu, x);
1364 #endif
1365 #if (GMX_SIMD_RSQRT_BITS*2 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1366     lu = gmx_simd_rsqrt_iter_d(lu, x);
1367 #endif
1368 #if (GMX_SIMD_RSQRT_BITS*4 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1369     lu = gmx_simd_rsqrt_iter_d(lu, x);
1370 #endif
1371 #if (GMX_SIMD_RSQRT_BITS*8 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1372     lu = gmx_simd_rsqrt_iter_d(lu, x);
1373 #endif
1374     return lu;
1375 }
1376
1377 /*! \brief Calculate 1/sqrt(x) for two SIMD doubles.
1378  *
1379  * \copydetails gmx_simd_invsqrt_pair_f
1380  */
1381 static gmx_inline void
1382 gmx_simd_invsqrt_pair_d(gmx_simd_double_t x0,    gmx_simd_double_t x1,
1383                         gmx_simd_double_t *out0, gmx_simd_double_t *out1)
1384 {
1385 #if (defined GMX_SIMD_HAVE_FLOAT) && (GMX_SIMD_FLOAT_WIDTH == 2*GMX_SIMD_DOUBLE_WIDTH) && (GMX_SIMD_RSQRT_BITS < 22)
1386     gmx_simd_float_t  xf  = gmx_simd_cvt_dd2f(x0, x1);
1387     gmx_simd_float_t  luf = gmx_simd_rsqrt_f(xf);
1388     gmx_simd_double_t lu0, lu1;
1389     /* Intermediate target is single - mantissa+1 bits */
1390 #if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
1391     luf = gmx_simd_rsqrt_iter_f(luf, xf);
1392 #endif
1393 #if (GMX_SIMD_RSQRT_BITS*2 < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
1394     luf = gmx_simd_rsqrt_iter_f(luf, xf);
1395 #endif
1396 #if (GMX_SIMD_RSQRT_BITS*4 < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
1397     luf = gmx_simd_rsqrt_iter_f(luf, xf);
1398 #endif
1399     gmx_simd_cvt_f2dd(luf, &lu0, &lu1);
1400     /* Last iteration(s) performed in double - if we had 22 bits, this gets us to 44 (~1e-15) */
1401 #if (GMX_SIMD_MATH_TARGET_SINGLE_BITS < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1402     lu0 = gmx_simd_rsqrt_iter_d(lu0, x0);
1403     lu1 = gmx_simd_rsqrt_iter_d(lu1, x1);
1404 #endif
1405 #if (GMX_SIMD_MATH_TARGET_SINGLE_BITS*2 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1406     lu0 = gmx_simd_rsqrt_iter_d(lu0, x0);
1407     lu1 = gmx_simd_rsqrt_iter_d(lu1, x1);
1408 #endif
1409     *out0 = lu0;
1410     *out1 = lu1;
1411 #else
1412     *out0 = gmx_simd_invsqrt_d(x0);
1413     *out1 = gmx_simd_invsqrt_d(x1);
1414 #endif
1415 }
1416
1417 /*! \brief Perform one Newton-Raphson iteration to improve 1/x for SIMD double.
1418  *
1419  * \copydetails gmx_simd_rcp_iter_f
1420  */
1421 static gmx_inline gmx_simd_double_t
1422 gmx_simd_rcp_iter_d(gmx_simd_double_t lu, gmx_simd_double_t x)
1423 {
1424     return gmx_simd_mul_d(lu, gmx_simd_fnmadd_d(lu, x, gmx_simd_set1_d(2.0)));
1425 }
1426
1427 /*! \brief Calculate 1/x for SIMD double.
1428  *
1429  * \copydetails gmx_simd_inv_f
1430  */
1431 static gmx_inline gmx_simd_double_t
1432 gmx_simd_inv_d(gmx_simd_double_t x)
1433 {
1434     gmx_simd_double_t lu = gmx_simd_rcp_d(x);
1435 #if (GMX_SIMD_RCP_BITS < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1436     lu = gmx_simd_rcp_iter_d(lu, x);
1437 #endif
1438 #if (GMX_SIMD_RCP_BITS*2 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1439     lu = gmx_simd_rcp_iter_d(lu, x);
1440 #endif
1441 #if (GMX_SIMD_RCP_BITS*4 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1442     lu = gmx_simd_rcp_iter_d(lu, x);
1443 #endif
1444 #if (GMX_SIMD_RCP_BITS*8 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
1445     lu = gmx_simd_rcp_iter_d(lu, x);
1446 #endif
1447     return lu;
1448 }
1449
1450 /*! \brief Calculate sqrt(x) correctly for SIMD doubles, including argument 0.0.
1451  *
1452  * \copydetails gmx_simd_sqrt_f
1453  */
1454 static gmx_inline gmx_simd_double_t
1455 gmx_simd_sqrt_d(gmx_simd_double_t x)
1456 {
1457     gmx_simd_dbool_t   mask;
1458     gmx_simd_double_t  res;
1459
1460     mask = gmx_simd_cmpeq_d(x, gmx_simd_setzero_d());
1461     res  = gmx_simd_blendnotzero_d(gmx_simd_invsqrt_d(x), mask);
1462     return gmx_simd_mul_d(res, x);
1463 }
1464
1465 /*! \brief SIMD double log(x). This is the natural logarithm.
1466  *
1467  * \copydetails gmx_simd_log_f
1468  */
1469 static gmx_inline gmx_simd_double_t
1470 gmx_simd_log_d(gmx_simd_double_t x)
1471 {
1472     const gmx_simd_double_t  half       = gmx_simd_set1_d(0.5);
1473     const gmx_simd_double_t  one        = gmx_simd_set1_d(1.0);
1474     const gmx_simd_double_t  sqrt2      = gmx_simd_set1_d(sqrt(2.0));
1475     const gmx_simd_double_t  corr       = gmx_simd_set1_d(0.693147180559945286226764);
1476     const gmx_simd_double_t  CL15       = gmx_simd_set1_d(0.148197055177935105296783);
1477     const gmx_simd_double_t  CL13       = gmx_simd_set1_d(0.153108178020442575739679);
1478     const gmx_simd_double_t  CL11       = gmx_simd_set1_d(0.181837339521549679055568);
1479     const gmx_simd_double_t  CL9        = gmx_simd_set1_d(0.22222194152736701733275);
1480     const gmx_simd_double_t  CL7        = gmx_simd_set1_d(0.285714288030134544449368);
1481     const gmx_simd_double_t  CL5        = gmx_simd_set1_d(0.399999999989941956712869);
1482     const gmx_simd_double_t  CL3        = gmx_simd_set1_d(0.666666666666685503450651);
1483     const gmx_simd_double_t  CL1        = gmx_simd_set1_d(2.0);
1484     gmx_simd_double_t        fexp, x2, p;
1485     gmx_simd_dbool_t         mask;
1486
1487     fexp  = gmx_simd_get_exponent_d(x);
1488     x     = gmx_simd_get_mantissa_d(x);
1489
1490     mask  = gmx_simd_cmplt_d(sqrt2, x);
1491     /* Adjust to non-IEEE format for x>sqrt(2): exponent += 1, mantissa *= 0.5 */
1492     fexp  = gmx_simd_add_d(fexp, gmx_simd_blendzero_d(one, mask));
1493     x     = gmx_simd_mul_d(x, gmx_simd_blendv_d(one, half, mask));
1494
1495     x     = gmx_simd_mul_d( gmx_simd_sub_d(x, one), gmx_simd_inv_d( gmx_simd_add_d(x, one) ) );
1496     x2    = gmx_simd_mul_d(x, x);
1497
1498     p     = gmx_simd_fmadd_d(CL15, x2, CL13);
1499     p     = gmx_simd_fmadd_d(p, x2, CL11);
1500     p     = gmx_simd_fmadd_d(p, x2, CL9);
1501     p     = gmx_simd_fmadd_d(p, x2, CL7);
1502     p     = gmx_simd_fmadd_d(p, x2, CL5);
1503     p     = gmx_simd_fmadd_d(p, x2, CL3);
1504     p     = gmx_simd_fmadd_d(p, x2, CL1);
1505     p     = gmx_simd_fmadd_d(p, x, gmx_simd_mul_d(corr, fexp));
1506
1507     return p;
1508 }
1509
1510 /*! \brief SIMD double 2^x.
1511  *
1512  * \copydetails gmx_simd_exp2_f
1513  */
1514 static gmx_inline gmx_simd_double_t
1515 gmx_simd_exp2_d(gmx_simd_double_t x)
1516 {
1517     const gmx_simd_double_t  arglimit      = gmx_simd_set1_d(1022.0);
1518     const gmx_simd_double_t  CE11          = gmx_simd_set1_d(4.435280790452730022081181e-10);
1519     const gmx_simd_double_t  CE10          = gmx_simd_set1_d(7.074105630863314448024247e-09);
1520     const gmx_simd_double_t  CE9           = gmx_simd_set1_d(1.017819803432096698472621e-07);
1521     const gmx_simd_double_t  CE8           = gmx_simd_set1_d(1.321543308956718799557863e-06);
1522     const gmx_simd_double_t  CE7           = gmx_simd_set1_d(0.00001525273348995851746990884);
1523     const gmx_simd_double_t  CE6           = gmx_simd_set1_d(0.0001540353046251466849082632);
1524     const gmx_simd_double_t  CE5           = gmx_simd_set1_d(0.001333355814678995257307880);
1525     const gmx_simd_double_t  CE4           = gmx_simd_set1_d(0.009618129107588335039176502);
1526     const gmx_simd_double_t  CE3           = gmx_simd_set1_d(0.05550410866481992147457793);
1527     const gmx_simd_double_t  CE2           = gmx_simd_set1_d(0.2402265069591015620470894);
1528     const gmx_simd_double_t  CE1           = gmx_simd_set1_d(0.6931471805599453304615075);
1529     const gmx_simd_double_t  one           = gmx_simd_set1_d(1.0);
1530     gmx_simd_double_t        fexppart;
1531     gmx_simd_double_t        intpart;
1532     gmx_simd_double_t        p;
1533     gmx_simd_dbool_t         valuemask;
1534
1535     fexppart  = gmx_simd_set_exponent_d(x);  /* rounds to nearest int internally */
1536     intpart   = gmx_simd_round_d(x);         /* use same rounding mode here */
1537     valuemask = gmx_simd_cmple_d(gmx_simd_fabs_d(x), arglimit);
1538     fexppart  = gmx_simd_blendzero_d(fexppart, valuemask);
1539     x         = gmx_simd_sub_d(x, intpart);
1540
1541     p         = gmx_simd_fmadd_d(CE11, x, CE10);
1542     p         = gmx_simd_fmadd_d(p, x, CE9);
1543     p         = gmx_simd_fmadd_d(p, x, CE8);
1544     p         = gmx_simd_fmadd_d(p, x, CE7);
1545     p         = gmx_simd_fmadd_d(p, x, CE6);
1546     p         = gmx_simd_fmadd_d(p, x, CE5);
1547     p         = gmx_simd_fmadd_d(p, x, CE4);
1548     p         = gmx_simd_fmadd_d(p, x, CE3);
1549     p         = gmx_simd_fmadd_d(p, x, CE2);
1550     p         = gmx_simd_fmadd_d(p, x, CE1);
1551     p         = gmx_simd_fmadd_d(p, x, one);
1552     x         = gmx_simd_mul_d(p, fexppart);
1553     return x;
1554 }
1555
1556 /*! \brief SIMD double exp(x).
1557  *
1558  * \copydetails gmx_simd_exp_f
1559  */
1560 static gmx_inline gmx_simd_double_t
1561 gmx_simd_exp_d(gmx_simd_double_t x)
1562 {
1563     const gmx_simd_double_t  argscale      = gmx_simd_set1_d(1.44269504088896340735992468100);
1564     const gmx_simd_double_t  arglimit      = gmx_simd_set1_d(1022.0);
1565     const gmx_simd_double_t  invargscale0  = gmx_simd_set1_d(0.69314718055966295651160180568695068359375);
1566     const gmx_simd_double_t  invargscale1  = gmx_simd_set1_d(2.8235290563031577122588448175013436025525412068e-13);
1567     const gmx_simd_double_t  CE12          = gmx_simd_set1_d(2.078375306791423699350304e-09);
1568     const gmx_simd_double_t  CE11          = gmx_simd_set1_d(2.518173854179933105218635e-08);
1569     const gmx_simd_double_t  CE10          = gmx_simd_set1_d(2.755842049600488770111608e-07);
1570     const gmx_simd_double_t  CE9           = gmx_simd_set1_d(2.755691815216689746619849e-06);
1571     const gmx_simd_double_t  CE8           = gmx_simd_set1_d(2.480158383706245033920920e-05);
1572     const gmx_simd_double_t  CE7           = gmx_simd_set1_d(0.0001984127043518048611841321);
1573     const gmx_simd_double_t  CE6           = gmx_simd_set1_d(0.001388888889360258341755930);
1574     const gmx_simd_double_t  CE5           = gmx_simd_set1_d(0.008333333332907368102819109);
1575     const gmx_simd_double_t  CE4           = gmx_simd_set1_d(0.04166666666663836745814631);
1576     const gmx_simd_double_t  CE3           = gmx_simd_set1_d(0.1666666666666796929434570);
1577     const gmx_simd_double_t  CE2           = gmx_simd_set1_d(0.5);
1578     const gmx_simd_double_t  one           = gmx_simd_set1_d(1.0);
1579     gmx_simd_double_t        fexppart;
1580     gmx_simd_double_t        intpart;
1581     gmx_simd_double_t        y, p;
1582     gmx_simd_dbool_t         valuemask;
1583
1584     y         = gmx_simd_mul_d(x, argscale);
1585     fexppart  = gmx_simd_set_exponent_d(y);  /* rounds to nearest int internally */
1586     intpart   = gmx_simd_round_d(y);         /* use same rounding mode here */
1587     valuemask = gmx_simd_cmple_d(gmx_simd_fabs_d(y), arglimit);
1588     fexppart  = gmx_simd_blendzero_d(fexppart, valuemask);
1589
1590     /* Extended precision arithmetics */
1591     x         = gmx_simd_fnmadd_d(invargscale0, intpart, x);
1592     x         = gmx_simd_fnmadd_d(invargscale1, intpart, x);
1593
1594     p         = gmx_simd_fmadd_d(CE12, x, CE11);
1595     p         = gmx_simd_fmadd_d(p, x, CE10);
1596     p         = gmx_simd_fmadd_d(p, x, CE9);
1597     p         = gmx_simd_fmadd_d(p, x, CE8);
1598     p         = gmx_simd_fmadd_d(p, x, CE7);
1599     p         = gmx_simd_fmadd_d(p, x, CE6);
1600     p         = gmx_simd_fmadd_d(p, x, CE5);
1601     p         = gmx_simd_fmadd_d(p, x, CE4);
1602     p         = gmx_simd_fmadd_d(p, x, CE3);
1603     p         = gmx_simd_fmadd_d(p, x, CE2);
1604     p         = gmx_simd_fmadd_d(p, gmx_simd_mul_d(x, x), gmx_simd_add_d(x, one));
1605     x         = gmx_simd_mul_d(p, fexppart);
1606     return x;
1607 }
1608
1609 /*! \brief SIMD double erf(x).
1610  *
1611  * \copydetails gmx_simd_erf_f
1612  */
1613 static gmx_inline gmx_simd_double_t
1614 gmx_simd_erf_d(gmx_simd_double_t x)
1615 {
1616     /* Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75] */
1617     const gmx_simd_double_t CAP4      = gmx_simd_set1_d(-0.431780540597889301512e-4);
1618     const gmx_simd_double_t CAP3      = gmx_simd_set1_d(-0.00578562306260059236059);
1619     const gmx_simd_double_t CAP2      = gmx_simd_set1_d(-0.028593586920219752446);
1620     const gmx_simd_double_t CAP1      = gmx_simd_set1_d(-0.315924962948621698209);
1621     const gmx_simd_double_t CAP0      = gmx_simd_set1_d(0.14952975608477029151);
1622
1623     const gmx_simd_double_t CAQ5      = gmx_simd_set1_d(-0.374089300177174709737e-5);
1624     const gmx_simd_double_t CAQ4      = gmx_simd_set1_d(0.00015126584532155383535);
1625     const gmx_simd_double_t CAQ3      = gmx_simd_set1_d(0.00536692680669480725423);
1626     const gmx_simd_double_t CAQ2      = gmx_simd_set1_d(0.0668686825594046122636);
1627     const gmx_simd_double_t CAQ1      = gmx_simd_set1_d(0.402604990869284362773);
1628     /* CAQ0 == 1.0 */
1629     const gmx_simd_double_t CAoffset  = gmx_simd_set1_d(0.9788494110107421875);
1630
1631     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5] */
1632     const gmx_simd_double_t CBP6      = gmx_simd_set1_d(2.49650423685462752497647637088e-10);
1633     const gmx_simd_double_t CBP5      = gmx_simd_set1_d(0.00119770193298159629350136085658);
1634     const gmx_simd_double_t CBP4      = gmx_simd_set1_d(0.0164944422378370965881008942733);
1635     const gmx_simd_double_t CBP3      = gmx_simd_set1_d(0.0984581468691775932063932439252);
1636     const gmx_simd_double_t CBP2      = gmx_simd_set1_d(0.317364595806937763843589437418);
1637     const gmx_simd_double_t CBP1      = gmx_simd_set1_d(0.554167062641455850932670067075);
1638     const gmx_simd_double_t CBP0      = gmx_simd_set1_d(0.427583576155807163756925301060);
1639     const gmx_simd_double_t CBQ7      = gmx_simd_set1_d(0.00212288829699830145976198384930);
1640     const gmx_simd_double_t CBQ6      = gmx_simd_set1_d(0.0334810979522685300554606393425);
1641     const gmx_simd_double_t CBQ5      = gmx_simd_set1_d(0.2361713785181450957579508850717);
1642     const gmx_simd_double_t CBQ4      = gmx_simd_set1_d(0.955364736493055670530981883072);
1643     const gmx_simd_double_t CBQ3      = gmx_simd_set1_d(2.36815675631420037315349279199);
1644     const gmx_simd_double_t CBQ2      = gmx_simd_set1_d(3.55261649184083035537184223542);
1645     const gmx_simd_double_t CBQ1      = gmx_simd_set1_d(2.93501136050160872574376997993);
1646     /* CBQ0 == 1.0 */
1647
1648     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf] */
1649     const gmx_simd_double_t CCP6      = gmx_simd_set1_d(-2.8175401114513378771);
1650     const gmx_simd_double_t CCP5      = gmx_simd_set1_d(-3.22729451764143718517);
1651     const gmx_simd_double_t CCP4      = gmx_simd_set1_d(-2.5518551727311523996);
1652     const gmx_simd_double_t CCP3      = gmx_simd_set1_d(-0.687717681153649930619);
1653     const gmx_simd_double_t CCP2      = gmx_simd_set1_d(-0.212652252872804219852);
1654     const gmx_simd_double_t CCP1      = gmx_simd_set1_d(0.0175389834052493308818);
1655     const gmx_simd_double_t CCP0      = gmx_simd_set1_d(0.00628057170626964891937);
1656
1657     const gmx_simd_double_t CCQ6      = gmx_simd_set1_d(5.48409182238641741584);
1658     const gmx_simd_double_t CCQ5      = gmx_simd_set1_d(13.5064170191802889145);
1659     const gmx_simd_double_t CCQ4      = gmx_simd_set1_d(22.9367376522880577224);
1660     const gmx_simd_double_t CCQ3      = gmx_simd_set1_d(15.930646027911794143);
1661     const gmx_simd_double_t CCQ2      = gmx_simd_set1_d(11.0567237927800161565);
1662     const gmx_simd_double_t CCQ1      = gmx_simd_set1_d(2.79257750980575282228);
1663     /* CCQ0 == 1.0 */
1664     const gmx_simd_double_t CCoffset  = gmx_simd_set1_d(0.5579090118408203125);
1665
1666     const gmx_simd_double_t one       = gmx_simd_set1_d(1.0);
1667     const gmx_simd_double_t two       = gmx_simd_set1_d(2.0);
1668
1669     gmx_simd_double_t       xabs, x2, x4, t, t2, w, w2;
1670     gmx_simd_double_t       PolyAP0, PolyAP1, PolyAQ0, PolyAQ1;
1671     gmx_simd_double_t       PolyBP0, PolyBP1, PolyBQ0, PolyBQ1;
1672     gmx_simd_double_t       PolyCP0, PolyCP1, PolyCQ0, PolyCQ1;
1673     gmx_simd_double_t       res_erf, res_erfcB, res_erfcC, res_erfc, res;
1674     gmx_simd_double_t       expmx2;
1675     gmx_simd_dbool_t        mask;
1676
1677     /* Calculate erf() */
1678     xabs     = gmx_simd_fabs_d(x);
1679     x2       = gmx_simd_mul_d(x, x);
1680     x4       = gmx_simd_mul_d(x2, x2);
1681
1682     PolyAP0  = gmx_simd_mul_d(CAP4, x4);
1683     PolyAP1  = gmx_simd_mul_d(CAP3, x4);
1684     PolyAP0  = gmx_simd_add_d(PolyAP0, CAP2);
1685     PolyAP1  = gmx_simd_add_d(PolyAP1, CAP1);
1686     PolyAP0  = gmx_simd_mul_d(PolyAP0, x4);
1687     PolyAP1  = gmx_simd_mul_d(PolyAP1, x2);
1688     PolyAP0  = gmx_simd_add_d(PolyAP0, CAP0);
1689     PolyAP0  = gmx_simd_add_d(PolyAP0, PolyAP1);
1690
1691     PolyAQ1  = gmx_simd_mul_d(CAQ5, x4);
1692     PolyAQ0  = gmx_simd_mul_d(CAQ4, x4);
1693     PolyAQ1  = gmx_simd_add_d(PolyAQ1, CAQ3);
1694     PolyAQ0  = gmx_simd_add_d(PolyAQ0, CAQ2);
1695     PolyAQ1  = gmx_simd_mul_d(PolyAQ1, x4);
1696     PolyAQ0  = gmx_simd_mul_d(PolyAQ0, x4);
1697     PolyAQ1  = gmx_simd_add_d(PolyAQ1, CAQ1);
1698     PolyAQ0  = gmx_simd_add_d(PolyAQ0, one);
1699     PolyAQ1  = gmx_simd_mul_d(PolyAQ1, x2);
1700     PolyAQ0  = gmx_simd_add_d(PolyAQ0, PolyAQ1);
1701
1702     res_erf  = gmx_simd_mul_d(PolyAP0, gmx_simd_inv_d(PolyAQ0));
1703     res_erf  = gmx_simd_add_d(CAoffset, res_erf);
1704     res_erf  = gmx_simd_mul_d(x, res_erf);
1705
1706     /* Calculate erfc() in range [1,4.5] */
1707     t       = gmx_simd_sub_d(xabs, one);
1708     t2      = gmx_simd_mul_d(t, t);
1709
1710     PolyBP0  = gmx_simd_mul_d(CBP6, t2);
1711     PolyBP1  = gmx_simd_mul_d(CBP5, t2);
1712     PolyBP0  = gmx_simd_add_d(PolyBP0, CBP4);
1713     PolyBP1  = gmx_simd_add_d(PolyBP1, CBP3);
1714     PolyBP0  = gmx_simd_mul_d(PolyBP0, t2);
1715     PolyBP1  = gmx_simd_mul_d(PolyBP1, t2);
1716     PolyBP0  = gmx_simd_add_d(PolyBP0, CBP2);
1717     PolyBP1  = gmx_simd_add_d(PolyBP1, CBP1);
1718     PolyBP0  = gmx_simd_mul_d(PolyBP0, t2);
1719     PolyBP1  = gmx_simd_mul_d(PolyBP1, t);
1720     PolyBP0  = gmx_simd_add_d(PolyBP0, CBP0);
1721     PolyBP0  = gmx_simd_add_d(PolyBP0, PolyBP1);
1722
1723     PolyBQ1 = gmx_simd_mul_d(CBQ7, t2);
1724     PolyBQ0 = gmx_simd_mul_d(CBQ6, t2);
1725     PolyBQ1 = gmx_simd_add_d(PolyBQ1, CBQ5);
1726     PolyBQ0 = gmx_simd_add_d(PolyBQ0, CBQ4);
1727     PolyBQ1 = gmx_simd_mul_d(PolyBQ1, t2);
1728     PolyBQ0 = gmx_simd_mul_d(PolyBQ0, t2);
1729     PolyBQ1 = gmx_simd_add_d(PolyBQ1, CBQ3);
1730     PolyBQ0 = gmx_simd_add_d(PolyBQ0, CBQ2);
1731     PolyBQ1 = gmx_simd_mul_d(PolyBQ1, t2);
1732     PolyBQ0 = gmx_simd_mul_d(PolyBQ0, t2);
1733     PolyBQ1 = gmx_simd_add_d(PolyBQ1, CBQ1);
1734     PolyBQ0 = gmx_simd_add_d(PolyBQ0, one);
1735     PolyBQ1 = gmx_simd_mul_d(PolyBQ1, t);
1736     PolyBQ0 = gmx_simd_add_d(PolyBQ0, PolyBQ1);
1737
1738     res_erfcB = gmx_simd_mul_d(PolyBP0, gmx_simd_inv_d(PolyBQ0));
1739
1740     res_erfcB = gmx_simd_mul_d(res_erfcB, xabs);
1741
1742     /* Calculate erfc() in range [4.5,inf] */
1743     w       = gmx_simd_inv_d(xabs);
1744     w2      = gmx_simd_mul_d(w, w);
1745
1746     PolyCP0  = gmx_simd_mul_d(CCP6, w2);
1747     PolyCP1  = gmx_simd_mul_d(CCP5, w2);
1748     PolyCP0  = gmx_simd_add_d(PolyCP0, CCP4);
1749     PolyCP1  = gmx_simd_add_d(PolyCP1, CCP3);
1750     PolyCP0  = gmx_simd_mul_d(PolyCP0, w2);
1751     PolyCP1  = gmx_simd_mul_d(PolyCP1, w2);
1752     PolyCP0  = gmx_simd_add_d(PolyCP0, CCP2);
1753     PolyCP1  = gmx_simd_add_d(PolyCP1, CCP1);
1754     PolyCP0  = gmx_simd_mul_d(PolyCP0, w2);
1755     PolyCP1  = gmx_simd_mul_d(PolyCP1, w);
1756     PolyCP0  = gmx_simd_add_d(PolyCP0, CCP0);
1757     PolyCP0  = gmx_simd_add_d(PolyCP0, PolyCP1);
1758
1759     PolyCQ0  = gmx_simd_mul_d(CCQ6, w2);
1760     PolyCQ1  = gmx_simd_mul_d(CCQ5, w2);
1761     PolyCQ0  = gmx_simd_add_d(PolyCQ0, CCQ4);
1762     PolyCQ1  = gmx_simd_add_d(PolyCQ1, CCQ3);
1763     PolyCQ0  = gmx_simd_mul_d(PolyCQ0, w2);
1764     PolyCQ1  = gmx_simd_mul_d(PolyCQ1, w2);
1765     PolyCQ0  = gmx_simd_add_d(PolyCQ0, CCQ2);
1766     PolyCQ1  = gmx_simd_add_d(PolyCQ1, CCQ1);
1767     PolyCQ0  = gmx_simd_mul_d(PolyCQ0, w2);
1768     PolyCQ1  = gmx_simd_mul_d(PolyCQ1, w);
1769     PolyCQ0  = gmx_simd_add_d(PolyCQ0, one);
1770     PolyCQ0  = gmx_simd_add_d(PolyCQ0, PolyCQ1);
1771
1772     expmx2   = gmx_simd_exp_d( gmx_simd_fneg_d(x2) );
1773
1774     res_erfcC = gmx_simd_mul_d(PolyCP0, gmx_simd_inv_d(PolyCQ0));
1775     res_erfcC = gmx_simd_add_d(res_erfcC, CCoffset);
1776     res_erfcC = gmx_simd_mul_d(res_erfcC, w);
1777
1778     mask     = gmx_simd_cmplt_d(gmx_simd_set1_d(4.5), xabs);
1779     res_erfc = gmx_simd_blendv_d(res_erfcB, res_erfcC, mask);
1780
1781     res_erfc = gmx_simd_mul_d(res_erfc, expmx2);
1782
1783     /* erfc(x<0) = 2-erfc(|x|) */
1784     mask     = gmx_simd_cmplt_d(x, gmx_simd_setzero_d());
1785     res_erfc = gmx_simd_blendv_d(res_erfc, gmx_simd_sub_d(two, res_erfc), mask);
1786
1787     /* Select erf() or erfc() */
1788     mask = gmx_simd_cmplt_d(xabs, one);
1789     res  = gmx_simd_blendv_d(gmx_simd_sub_d(one, res_erfc), res_erf, mask);
1790
1791     return res;
1792 }
1793
1794 /*! \brief SIMD double erfc(x).
1795  *
1796  * \copydetails gmx_simd_erfc_f
1797  */
1798 static gmx_inline gmx_simd_double_t
1799 gmx_simd_erfc_d(gmx_simd_double_t x)
1800 {
1801     /* Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75] */
1802     const gmx_simd_double_t CAP4      = gmx_simd_set1_d(-0.431780540597889301512e-4);
1803     const gmx_simd_double_t CAP3      = gmx_simd_set1_d(-0.00578562306260059236059);
1804     const gmx_simd_double_t CAP2      = gmx_simd_set1_d(-0.028593586920219752446);
1805     const gmx_simd_double_t CAP1      = gmx_simd_set1_d(-0.315924962948621698209);
1806     const gmx_simd_double_t CAP0      = gmx_simd_set1_d(0.14952975608477029151);
1807
1808     const gmx_simd_double_t CAQ5      = gmx_simd_set1_d(-0.374089300177174709737e-5);
1809     const gmx_simd_double_t CAQ4      = gmx_simd_set1_d(0.00015126584532155383535);
1810     const gmx_simd_double_t CAQ3      = gmx_simd_set1_d(0.00536692680669480725423);
1811     const gmx_simd_double_t CAQ2      = gmx_simd_set1_d(0.0668686825594046122636);
1812     const gmx_simd_double_t CAQ1      = gmx_simd_set1_d(0.402604990869284362773);
1813     /* CAQ0 == 1.0 */
1814     const gmx_simd_double_t CAoffset  = gmx_simd_set1_d(0.9788494110107421875);
1815
1816     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5] */
1817     const gmx_simd_double_t CBP6      = gmx_simd_set1_d(2.49650423685462752497647637088e-10);
1818     const gmx_simd_double_t CBP5      = gmx_simd_set1_d(0.00119770193298159629350136085658);
1819     const gmx_simd_double_t CBP4      = gmx_simd_set1_d(0.0164944422378370965881008942733);
1820     const gmx_simd_double_t CBP3      = gmx_simd_set1_d(0.0984581468691775932063932439252);
1821     const gmx_simd_double_t CBP2      = gmx_simd_set1_d(0.317364595806937763843589437418);
1822     const gmx_simd_double_t CBP1      = gmx_simd_set1_d(0.554167062641455850932670067075);
1823     const gmx_simd_double_t CBP0      = gmx_simd_set1_d(0.427583576155807163756925301060);
1824     const gmx_simd_double_t CBQ7      = gmx_simd_set1_d(0.00212288829699830145976198384930);
1825     const gmx_simd_double_t CBQ6      = gmx_simd_set1_d(0.0334810979522685300554606393425);
1826     const gmx_simd_double_t CBQ5      = gmx_simd_set1_d(0.2361713785181450957579508850717);
1827     const gmx_simd_double_t CBQ4      = gmx_simd_set1_d(0.955364736493055670530981883072);
1828     const gmx_simd_double_t CBQ3      = gmx_simd_set1_d(2.36815675631420037315349279199);
1829     const gmx_simd_double_t CBQ2      = gmx_simd_set1_d(3.55261649184083035537184223542);
1830     const gmx_simd_double_t CBQ1      = gmx_simd_set1_d(2.93501136050160872574376997993);
1831     /* CBQ0 == 1.0 */
1832
1833     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf] */
1834     const gmx_simd_double_t CCP6      = gmx_simd_set1_d(-2.8175401114513378771);
1835     const gmx_simd_double_t CCP5      = gmx_simd_set1_d(-3.22729451764143718517);
1836     const gmx_simd_double_t CCP4      = gmx_simd_set1_d(-2.5518551727311523996);
1837     const gmx_simd_double_t CCP3      = gmx_simd_set1_d(-0.687717681153649930619);
1838     const gmx_simd_double_t CCP2      = gmx_simd_set1_d(-0.212652252872804219852);
1839     const gmx_simd_double_t CCP1      = gmx_simd_set1_d(0.0175389834052493308818);
1840     const gmx_simd_double_t CCP0      = gmx_simd_set1_d(0.00628057170626964891937);
1841
1842     const gmx_simd_double_t CCQ6      = gmx_simd_set1_d(5.48409182238641741584);
1843     const gmx_simd_double_t CCQ5      = gmx_simd_set1_d(13.5064170191802889145);
1844     const gmx_simd_double_t CCQ4      = gmx_simd_set1_d(22.9367376522880577224);
1845     const gmx_simd_double_t CCQ3      = gmx_simd_set1_d(15.930646027911794143);
1846     const gmx_simd_double_t CCQ2      = gmx_simd_set1_d(11.0567237927800161565);
1847     const gmx_simd_double_t CCQ1      = gmx_simd_set1_d(2.79257750980575282228);
1848     /* CCQ0 == 1.0 */
1849     const gmx_simd_double_t CCoffset  = gmx_simd_set1_d(0.5579090118408203125);
1850
1851     const gmx_simd_double_t one       = gmx_simd_set1_d(1.0);
1852     const gmx_simd_double_t two       = gmx_simd_set1_d(2.0);
1853
1854     gmx_simd_double_t       xabs, x2, x4, t, t2, w, w2;
1855     gmx_simd_double_t       PolyAP0, PolyAP1, PolyAQ0, PolyAQ1;
1856     gmx_simd_double_t       PolyBP0, PolyBP1, PolyBQ0, PolyBQ1;
1857     gmx_simd_double_t       PolyCP0, PolyCP1, PolyCQ0, PolyCQ1;
1858     gmx_simd_double_t       res_erf, res_erfcB, res_erfcC, res_erfc, res;
1859     gmx_simd_double_t       expmx2;
1860     gmx_simd_dbool_t        mask;
1861
1862     /* Calculate erf() */
1863     xabs     = gmx_simd_fabs_d(x);
1864     x2       = gmx_simd_mul_d(x, x);
1865     x4       = gmx_simd_mul_d(x2, x2);
1866
1867     PolyAP0  = gmx_simd_mul_d(CAP4, x4);
1868     PolyAP1  = gmx_simd_mul_d(CAP3, x4);
1869     PolyAP0  = gmx_simd_add_d(PolyAP0, CAP2);
1870     PolyAP1  = gmx_simd_add_d(PolyAP1, CAP1);
1871     PolyAP0  = gmx_simd_mul_d(PolyAP0, x4);
1872     PolyAP1  = gmx_simd_mul_d(PolyAP1, x2);
1873     PolyAP0  = gmx_simd_add_d(PolyAP0, CAP0);
1874     PolyAP0  = gmx_simd_add_d(PolyAP0, PolyAP1);
1875
1876     PolyAQ1  = gmx_simd_mul_d(CAQ5, x4);
1877     PolyAQ0  = gmx_simd_mul_d(CAQ4, x4);
1878     PolyAQ1  = gmx_simd_add_d(PolyAQ1, CAQ3);
1879     PolyAQ0  = gmx_simd_add_d(PolyAQ0, CAQ2);
1880     PolyAQ1  = gmx_simd_mul_d(PolyAQ1, x4);
1881     PolyAQ0  = gmx_simd_mul_d(PolyAQ0, x4);
1882     PolyAQ1  = gmx_simd_add_d(PolyAQ1, CAQ1);
1883     PolyAQ0  = gmx_simd_add_d(PolyAQ0, one);
1884     PolyAQ1  = gmx_simd_mul_d(PolyAQ1, x2);
1885     PolyAQ0  = gmx_simd_add_d(PolyAQ0, PolyAQ1);
1886
1887     res_erf  = gmx_simd_mul_d(PolyAP0, gmx_simd_inv_d(PolyAQ0));
1888     res_erf  = gmx_simd_add_d(CAoffset, res_erf);
1889     res_erf  = gmx_simd_mul_d(x, res_erf);
1890
1891     /* Calculate erfc() in range [1,4.5] */
1892     t       = gmx_simd_sub_d(xabs, one);
1893     t2      = gmx_simd_mul_d(t, t);
1894
1895     PolyBP0  = gmx_simd_mul_d(CBP6, t2);
1896     PolyBP1  = gmx_simd_mul_d(CBP5, t2);
1897     PolyBP0  = gmx_simd_add_d(PolyBP0, CBP4);
1898     PolyBP1  = gmx_simd_add_d(PolyBP1, CBP3);
1899     PolyBP0  = gmx_simd_mul_d(PolyBP0, t2);
1900     PolyBP1  = gmx_simd_mul_d(PolyBP1, t2);
1901     PolyBP0  = gmx_simd_add_d(PolyBP0, CBP2);
1902     PolyBP1  = gmx_simd_add_d(PolyBP1, CBP1);
1903     PolyBP0  = gmx_simd_mul_d(PolyBP0, t2);
1904     PolyBP1  = gmx_simd_mul_d(PolyBP1, t);
1905     PolyBP0  = gmx_simd_add_d(PolyBP0, CBP0);
1906     PolyBP0  = gmx_simd_add_d(PolyBP0, PolyBP1);
1907
1908     PolyBQ1 = gmx_simd_mul_d(CBQ7, t2);
1909     PolyBQ0 = gmx_simd_mul_d(CBQ6, t2);
1910     PolyBQ1 = gmx_simd_add_d(PolyBQ1, CBQ5);
1911     PolyBQ0 = gmx_simd_add_d(PolyBQ0, CBQ4);
1912     PolyBQ1 = gmx_simd_mul_d(PolyBQ1, t2);
1913     PolyBQ0 = gmx_simd_mul_d(PolyBQ0, t2);
1914     PolyBQ1 = gmx_simd_add_d(PolyBQ1, CBQ3);
1915     PolyBQ0 = gmx_simd_add_d(PolyBQ0, CBQ2);
1916     PolyBQ1 = gmx_simd_mul_d(PolyBQ1, t2);
1917     PolyBQ0 = gmx_simd_mul_d(PolyBQ0, t2);
1918     PolyBQ1 = gmx_simd_add_d(PolyBQ1, CBQ1);
1919     PolyBQ0 = gmx_simd_add_d(PolyBQ0, one);
1920     PolyBQ1 = gmx_simd_mul_d(PolyBQ1, t);
1921     PolyBQ0 = gmx_simd_add_d(PolyBQ0, PolyBQ1);
1922
1923     res_erfcB = gmx_simd_mul_d(PolyBP0, gmx_simd_inv_d(PolyBQ0));
1924
1925     res_erfcB = gmx_simd_mul_d(res_erfcB, xabs);
1926
1927     /* Calculate erfc() in range [4.5,inf] */
1928     w       = gmx_simd_inv_d(xabs);
1929     w2      = gmx_simd_mul_d(w, w);
1930
1931     PolyCP0  = gmx_simd_mul_d(CCP6, w2);
1932     PolyCP1  = gmx_simd_mul_d(CCP5, w2);
1933     PolyCP0  = gmx_simd_add_d(PolyCP0, CCP4);
1934     PolyCP1  = gmx_simd_add_d(PolyCP1, CCP3);
1935     PolyCP0  = gmx_simd_mul_d(PolyCP0, w2);
1936     PolyCP1  = gmx_simd_mul_d(PolyCP1, w2);
1937     PolyCP0  = gmx_simd_add_d(PolyCP0, CCP2);
1938     PolyCP1  = gmx_simd_add_d(PolyCP1, CCP1);
1939     PolyCP0  = gmx_simd_mul_d(PolyCP0, w2);
1940     PolyCP1  = gmx_simd_mul_d(PolyCP1, w);
1941     PolyCP0  = gmx_simd_add_d(PolyCP0, CCP0);
1942     PolyCP0  = gmx_simd_add_d(PolyCP0, PolyCP1);
1943
1944     PolyCQ0  = gmx_simd_mul_d(CCQ6, w2);
1945     PolyCQ1  = gmx_simd_mul_d(CCQ5, w2);
1946     PolyCQ0  = gmx_simd_add_d(PolyCQ0, CCQ4);
1947     PolyCQ1  = gmx_simd_add_d(PolyCQ1, CCQ3);
1948     PolyCQ0  = gmx_simd_mul_d(PolyCQ0, w2);
1949     PolyCQ1  = gmx_simd_mul_d(PolyCQ1, w2);
1950     PolyCQ0  = gmx_simd_add_d(PolyCQ0, CCQ2);
1951     PolyCQ1  = gmx_simd_add_d(PolyCQ1, CCQ1);
1952     PolyCQ0  = gmx_simd_mul_d(PolyCQ0, w2);
1953     PolyCQ1  = gmx_simd_mul_d(PolyCQ1, w);
1954     PolyCQ0  = gmx_simd_add_d(PolyCQ0, one);
1955     PolyCQ0  = gmx_simd_add_d(PolyCQ0, PolyCQ1);
1956
1957     expmx2   = gmx_simd_exp_d( gmx_simd_fneg_d(x2) );
1958
1959     res_erfcC = gmx_simd_mul_d(PolyCP0, gmx_simd_inv_d(PolyCQ0));
1960     res_erfcC = gmx_simd_add_d(res_erfcC, CCoffset);
1961     res_erfcC = gmx_simd_mul_d(res_erfcC, w);
1962
1963     mask     = gmx_simd_cmplt_d(gmx_simd_set1_d(4.5), xabs);
1964     res_erfc = gmx_simd_blendv_d(res_erfcB, res_erfcC, mask);
1965
1966     res_erfc = gmx_simd_mul_d(res_erfc, expmx2);
1967
1968     /* erfc(x<0) = 2-erfc(|x|) */
1969     mask     = gmx_simd_cmplt_d(x, gmx_simd_setzero_d());
1970     res_erfc = gmx_simd_blendv_d(res_erfc, gmx_simd_sub_d(two, res_erfc), mask);
1971
1972     /* Select erf() or erfc() */
1973     mask = gmx_simd_cmplt_d(xabs, one);
1974     res  = gmx_simd_blendv_d(res_erfc, gmx_simd_sub_d(one, res_erf), mask);
1975
1976     return res;
1977 }
1978
1979 /*! \brief SIMD double sin \& cos.
1980  *
1981  * \copydetails gmx_simd_sincos_f
1982  */
1983 static gmx_inline void
1984 gmx_simd_sincos_d(gmx_simd_double_t x, gmx_simd_double_t *sinval, gmx_simd_double_t *cosval)
1985 {
1986     /* Constants to subtract Pi/4*x from y while minimizing precision loss */
1987     const gmx_simd_double_t  argred0         = gmx_simd_set1_d(2*0.78539816290140151978);
1988     const gmx_simd_double_t  argred1         = gmx_simd_set1_d(2*4.9604678871439933374e-10);
1989     const gmx_simd_double_t  argred2         = gmx_simd_set1_d(2*1.1258708853173288931e-18);
1990     const gmx_simd_double_t  argred3         = gmx_simd_set1_d(2*1.7607799325916000908e-27);
1991     const gmx_simd_double_t  two_over_pi     = gmx_simd_set1_d(2.0/M_PI);
1992     const gmx_simd_double_t  const_sin5      = gmx_simd_set1_d( 1.58938307283228937328511e-10);
1993     const gmx_simd_double_t  const_sin4      = gmx_simd_set1_d(-2.50506943502539773349318e-08);
1994     const gmx_simd_double_t  const_sin3      = gmx_simd_set1_d( 2.75573131776846360512547e-06);
1995     const gmx_simd_double_t  const_sin2      = gmx_simd_set1_d(-0.000198412698278911770864914);
1996     const gmx_simd_double_t  const_sin1      = gmx_simd_set1_d( 0.0083333333333191845961746);
1997     const gmx_simd_double_t  const_sin0      = gmx_simd_set1_d(-0.166666666666666130709393);
1998
1999     const gmx_simd_double_t  const_cos7      = gmx_simd_set1_d(-1.13615350239097429531523e-11);
2000     const gmx_simd_double_t  const_cos6      = gmx_simd_set1_d( 2.08757471207040055479366e-09);
2001     const gmx_simd_double_t  const_cos5      = gmx_simd_set1_d(-2.75573144028847567498567e-07);
2002     const gmx_simd_double_t  const_cos4      = gmx_simd_set1_d( 2.48015872890001867311915e-05);
2003     const gmx_simd_double_t  const_cos3      = gmx_simd_set1_d(-0.00138888888888714019282329);
2004     const gmx_simd_double_t  const_cos2      = gmx_simd_set1_d( 0.0416666666666665519592062);
2005     const gmx_simd_double_t  half            = gmx_simd_set1_d(0.5);
2006     const gmx_simd_double_t  one             = gmx_simd_set1_d(1.0);
2007     gmx_simd_double_t        ssign, csign;
2008     gmx_simd_double_t        x2, y, z, psin, pcos, sss, ccc;
2009     gmx_simd_dbool_t         mask;
2010 #if (defined GMX_SIMD_HAVE_DINT32) && (defined GMX_SIMD_HAVE_DINT32_ARITHMETICS) && (defined GMX_SIMD_HAVE_LOGICAL)
2011     const gmx_simd_dint32_t  ione            = gmx_simd_set1_di(1);
2012     const gmx_simd_dint32_t  itwo            = gmx_simd_set1_di(2);
2013     gmx_simd_dint32_t        iy;
2014
2015     z       = gmx_simd_mul_d(x, two_over_pi);
2016     iy      = gmx_simd_cvt_d2i(z);
2017     y       = gmx_simd_round_d(z);
2018
2019     mask    = gmx_simd_cvt_dib2db(gmx_simd_cmpeq_di(gmx_simd_and_di(iy, ione), gmx_simd_setzero_di()));
2020     ssign   = gmx_simd_blendzero_d(gmx_simd_set1_d(-0.0), gmx_simd_cvt_dib2db(gmx_simd_cmpeq_di(gmx_simd_and_di(iy, itwo), itwo)));
2021     csign   = gmx_simd_blendzero_d(gmx_simd_set1_d(-0.0), gmx_simd_cvt_dib2db(gmx_simd_cmpeq_di(gmx_simd_and_di(gmx_simd_add_di(iy, ione), itwo), itwo)));
2022 #else
2023     const gmx_simd_double_t  quarter         = gmx_simd_set1_d(0.25);
2024     const gmx_simd_double_t  minusquarter    = gmx_simd_set1_d(-0.25);
2025     gmx_simd_double_t        q;
2026     gmx_simd_dbool_t         m1, m2, m3;
2027
2028     /* The most obvious way to find the arguments quadrant in the unit circle
2029      * to calculate the sign is to use integer arithmetic, but that is not
2030      * present in all SIMD implementations. As an alternative, we have devised a
2031      * pure floating-point algorithm that uses truncation for argument reduction
2032      * so that we get a new value 0<=q<1 over the unit circle, and then
2033      * do floating-point comparisons with fractions. This is likely to be
2034      * slightly slower (~10%) due to the longer latencies of floating-point, so
2035      * we only use it when integer SIMD arithmetic is not present.
2036      */
2037     ssign   = x;
2038     x       = gmx_simd_fabs_d(x);
2039     /* It is critical that half-way cases are rounded down */
2040     z       = gmx_simd_fmadd_d(x, two_over_pi, half);
2041     y       = gmx_simd_trunc_d(z);
2042     q       = gmx_simd_mul_d(z, quarter);
2043     q       = gmx_simd_sub_d(q, gmx_simd_trunc_d(q));
2044     /* z now starts at 0.0 for x=-pi/4 (although neg. values cannot occur), and
2045      * then increased by 1.0 as x increases by 2*Pi, when it resets to 0.0.
2046      * This removes the 2*Pi periodicity without using any integer arithmetic.
2047      * First check if y had the value 2 or 3, set csign if true.
2048      */
2049     q       = gmx_simd_sub_d(q, half);
2050     /* If we have logical operations we can work directly on the signbit, which
2051      * saves instructions. Otherwise we need to represent signs as +1.0/-1.0.
2052      * Thus, if you are altering defines to debug alternative code paths, the
2053      * two GMX_SIMD_HAVE_LOGICAL sections in this routine must either both be
2054      * active or inactive - you will get errors if only one is used.
2055      */
2056 #    ifdef GMX_SIMD_HAVE_LOGICAL
2057     ssign   = gmx_simd_and_d(ssign, gmx_simd_set1_d(-0.0));
2058     csign   = gmx_simd_andnot_d(q, gmx_simd_set1_d(-0.0));
2059     ssign   = gmx_simd_xor_d(ssign, csign);
2060 #    else
2061     csign   = gmx_simd_xor_sign_d(gmx_simd_set1_d(-1.0), q);
2062     ssign   = gmx_simd_xor_sign_d(ssign, csign);    /* swap ssign if csign was set. */
2063 #    endif
2064     /* Check if y had value 1 or 3 (remember we subtracted 0.5 from q) */
2065     m1      = gmx_simd_cmplt_d(q, minusquarter);
2066     m2      = gmx_simd_cmple_d(gmx_simd_setzero_d(), q);
2067     m3      = gmx_simd_cmplt_d(q, quarter);
2068     m2      = gmx_simd_and_db(m2, m3);
2069     mask    = gmx_simd_or_db(m1, m2);
2070     /* where mask is FALSE, set sign. */
2071     csign   = gmx_simd_xor_sign_d(csign, gmx_simd_blendv_d(gmx_simd_set1_d(-1.0), one, mask));
2072 #endif
2073     x       = gmx_simd_fnmadd_d(y, argred0, x);
2074     x       = gmx_simd_fnmadd_d(y, argred1, x);
2075     x       = gmx_simd_fnmadd_d(y, argred2, x);
2076     x       = gmx_simd_fnmadd_d(y, argred3, x);
2077     x2      = gmx_simd_mul_d(x, x);
2078
2079     psin    = gmx_simd_fmadd_d(const_sin5, x2, const_sin4);
2080     psin    = gmx_simd_fmadd_d(psin, x2, const_sin3);
2081     psin    = gmx_simd_fmadd_d(psin, x2, const_sin2);
2082     psin    = gmx_simd_fmadd_d(psin, x2, const_sin1);
2083     psin    = gmx_simd_fmadd_d(psin, x2, const_sin0);
2084     psin    = gmx_simd_fmadd_d(psin, gmx_simd_mul_d(x2, x), x);
2085
2086     pcos    = gmx_simd_fmadd_d(const_cos7, x2, const_cos6);
2087     pcos    = gmx_simd_fmadd_d(pcos, x2, const_cos5);
2088     pcos    = gmx_simd_fmadd_d(pcos, x2, const_cos4);
2089     pcos    = gmx_simd_fmadd_d(pcos, x2, const_cos3);
2090     pcos    = gmx_simd_fmadd_d(pcos, x2, const_cos2);
2091     pcos    = gmx_simd_fmsub_d(pcos, x2, half);
2092     pcos    = gmx_simd_fmadd_d(pcos, x2, one);
2093
2094     sss     = gmx_simd_blendv_d(pcos, psin, mask);
2095     ccc     = gmx_simd_blendv_d(psin, pcos, mask);
2096     /* See comment for GMX_SIMD_HAVE_LOGICAL section above. */
2097 #ifdef GMX_SIMD_HAVE_LOGICAL
2098     *sinval = gmx_simd_xor_d(sss, ssign);
2099     *cosval = gmx_simd_xor_d(ccc, csign);
2100 #else
2101     *sinval = gmx_simd_xor_sign_d(sss, ssign);
2102     *cosval = gmx_simd_xor_sign_d(ccc, csign);
2103 #endif
2104 }
2105
2106 /*! \brief SIMD double sin(x).
2107  *
2108  * \copydetails gmx_simd_sin_f
2109  */
2110 static gmx_inline gmx_simd_double_t
2111 gmx_simd_sin_d(gmx_simd_double_t x)
2112 {
2113     gmx_simd_double_t s, c;
2114     gmx_simd_sincos_d(x, &s, &c);
2115     return s;
2116 }
2117
2118 /*! \brief SIMD double cos(x).
2119  *
2120  * \copydetails gmx_simd_cos_f
2121  */
2122 static gmx_inline gmx_simd_double_t
2123 gmx_simd_cos_d(gmx_simd_double_t x)
2124 {
2125     gmx_simd_double_t s, c;
2126     gmx_simd_sincos_d(x, &s, &c);
2127     return c;
2128 }
2129
2130 /*! \brief SIMD double tan(x).
2131  *
2132  * \copydetails gmx_simd_tan_f
2133  */
2134 static gmx_inline gmx_simd_double_t
2135 gmx_simd_tan_d(gmx_simd_double_t x)
2136 {
2137     const gmx_simd_double_t  argred0         = gmx_simd_set1_d(2*0.78539816290140151978);
2138     const gmx_simd_double_t  argred1         = gmx_simd_set1_d(2*4.9604678871439933374e-10);
2139     const gmx_simd_double_t  argred2         = gmx_simd_set1_d(2*1.1258708853173288931e-18);
2140     const gmx_simd_double_t  argred3         = gmx_simd_set1_d(2*1.7607799325916000908e-27);
2141     const gmx_simd_double_t  two_over_pi     = gmx_simd_set1_d(2.0/M_PI);
2142     const gmx_simd_double_t  CT15            = gmx_simd_set1_d(1.01419718511083373224408e-05);
2143     const gmx_simd_double_t  CT14            = gmx_simd_set1_d(-2.59519791585924697698614e-05);
2144     const gmx_simd_double_t  CT13            = gmx_simd_set1_d(5.23388081915899855325186e-05);
2145     const gmx_simd_double_t  CT12            = gmx_simd_set1_d(-3.05033014433946488225616e-05);
2146     const gmx_simd_double_t  CT11            = gmx_simd_set1_d(7.14707504084242744267497e-05);
2147     const gmx_simd_double_t  CT10            = gmx_simd_set1_d(8.09674518280159187045078e-05);
2148     const gmx_simd_double_t  CT9             = gmx_simd_set1_d(0.000244884931879331847054404);
2149     const gmx_simd_double_t  CT8             = gmx_simd_set1_d(0.000588505168743587154904506);
2150     const gmx_simd_double_t  CT7             = gmx_simd_set1_d(0.00145612788922812427978848);
2151     const gmx_simd_double_t  CT6             = gmx_simd_set1_d(0.00359208743836906619142924);
2152     const gmx_simd_double_t  CT5             = gmx_simd_set1_d(0.00886323944362401618113356);
2153     const gmx_simd_double_t  CT4             = gmx_simd_set1_d(0.0218694882853846389592078);
2154     const gmx_simd_double_t  CT3             = gmx_simd_set1_d(0.0539682539781298417636002);
2155     const gmx_simd_double_t  CT2             = gmx_simd_set1_d(0.133333333333125941821962);
2156     const gmx_simd_double_t  CT1             = gmx_simd_set1_d(0.333333333333334980164153);
2157
2158     gmx_simd_double_t        x2, p, y, z;
2159     gmx_simd_dbool_t         mask;
2160
2161 #if (defined GMX_SIMD_HAVE_DINT32) && (defined GMX_SIMD_HAVE_DINT32_ARITHMETICS) && (defined GMX_SIMD_HAVE_LOGICAL)
2162     gmx_simd_dint32_t  iy;
2163     gmx_simd_dint32_t  ione = gmx_simd_set1_di(1);
2164
2165     z       = gmx_simd_mul_d(x, two_over_pi);
2166     iy      = gmx_simd_cvt_d2i(z);
2167     y       = gmx_simd_round_d(z);
2168     mask    = gmx_simd_cvt_dib2db(gmx_simd_cmpeq_di(gmx_simd_and_di(iy, ione), ione));
2169
2170     x       = gmx_simd_fnmadd_d(y, argred0, x);
2171     x       = gmx_simd_fnmadd_d(y, argred1, x);
2172     x       = gmx_simd_fnmadd_d(y, argred2, x);
2173     x       = gmx_simd_fnmadd_d(y, argred3, x);
2174     x       = gmx_simd_xor_d(gmx_simd_blendzero_d(gmx_simd_set1_d(-0.0), mask), x);
2175 #else
2176     const gmx_simd_double_t  quarter         = gmx_simd_set1_d(0.25);
2177     const gmx_simd_double_t  half            = gmx_simd_set1_d(0.5);
2178     const gmx_simd_double_t  threequarter    = gmx_simd_set1_d(0.75);
2179     gmx_simd_double_t        w, q;
2180     gmx_simd_dbool_t         m1, m2, m3;
2181
2182     w       = gmx_simd_fabs_d(x);
2183     z       = gmx_simd_fmadd_d(w, two_over_pi, half);
2184     y       = gmx_simd_trunc_d(z);
2185     q       = gmx_simd_mul_d(z, quarter);
2186     q       = gmx_simd_sub_d(q, gmx_simd_trunc_d(q));
2187     m1      = gmx_simd_cmple_d(quarter, q);
2188     m2      = gmx_simd_cmplt_d(q, half);
2189     m3      = gmx_simd_cmple_d(threequarter, q);
2190     m1      = gmx_simd_and_db(m1, m2);
2191     mask    = gmx_simd_or_db(m1, m3);
2192     w       = gmx_simd_fnmadd_d(y, argred0, w);
2193     w       = gmx_simd_fnmadd_d(y, argred1, w);
2194     w       = gmx_simd_fnmadd_d(y, argred2, w);
2195     w       = gmx_simd_fnmadd_d(y, argred3, w);
2196
2197     w       = gmx_simd_blendv_d(w, gmx_simd_fneg_d(w), mask);
2198     x       = gmx_simd_xor_sign_d(w, x);
2199 #endif
2200     x2      = gmx_simd_mul_d(x, x);
2201     p       = gmx_simd_fmadd_d(CT15, x2, CT14);
2202     p       = gmx_simd_fmadd_d(p, x2, CT13);
2203     p       = gmx_simd_fmadd_d(p, x2, CT12);
2204     p       = gmx_simd_fmadd_d(p, x2, CT11);
2205     p       = gmx_simd_fmadd_d(p, x2, CT10);
2206     p       = gmx_simd_fmadd_d(p, x2, CT9);
2207     p       = gmx_simd_fmadd_d(p, x2, CT8);
2208     p       = gmx_simd_fmadd_d(p, x2, CT7);
2209     p       = gmx_simd_fmadd_d(p, x2, CT6);
2210     p       = gmx_simd_fmadd_d(p, x2, CT5);
2211     p       = gmx_simd_fmadd_d(p, x2, CT4);
2212     p       = gmx_simd_fmadd_d(p, x2, CT3);
2213     p       = gmx_simd_fmadd_d(p, x2, CT2);
2214     p       = gmx_simd_fmadd_d(p, x2, CT1);
2215     p       = gmx_simd_fmadd_d(x2, gmx_simd_mul_d(p, x), x);
2216
2217     p       = gmx_simd_blendv_d( p, gmx_simd_inv_d(p), mask);
2218     return p;
2219 }
2220
2221 /*! \brief SIMD double asin(x).
2222  *
2223  * \copydetails gmx_simd_asin_f
2224  */
2225 static gmx_inline gmx_simd_double_t
2226 gmx_simd_asin_d(gmx_simd_double_t x)
2227 {
2228     /* Same algorithm as cephes library */
2229     const gmx_simd_double_t limit1    = gmx_simd_set1_d(0.625);
2230     const gmx_simd_double_t limit2    = gmx_simd_set1_d(1e-8);
2231     const gmx_simd_double_t one       = gmx_simd_set1_d(1.0);
2232     const gmx_simd_double_t quarterpi = gmx_simd_set1_d(M_PI/4.0);
2233     const gmx_simd_double_t morebits  = gmx_simd_set1_d(6.123233995736765886130e-17);
2234
2235     const gmx_simd_double_t P5        = gmx_simd_set1_d(4.253011369004428248960e-3);
2236     const gmx_simd_double_t P4        = gmx_simd_set1_d(-6.019598008014123785661e-1);
2237     const gmx_simd_double_t P3        = gmx_simd_set1_d(5.444622390564711410273e0);
2238     const gmx_simd_double_t P2        = gmx_simd_set1_d(-1.626247967210700244449e1);
2239     const gmx_simd_double_t P1        = gmx_simd_set1_d(1.956261983317594739197e1);
2240     const gmx_simd_double_t P0        = gmx_simd_set1_d(-8.198089802484824371615e0);
2241
2242     const gmx_simd_double_t Q4        = gmx_simd_set1_d(-1.474091372988853791896e1);
2243     const gmx_simd_double_t Q3        = gmx_simd_set1_d(7.049610280856842141659e1);
2244     const gmx_simd_double_t Q2        = gmx_simd_set1_d(-1.471791292232726029859e2);
2245     const gmx_simd_double_t Q1        = gmx_simd_set1_d(1.395105614657485689735e2);
2246     const gmx_simd_double_t Q0        = gmx_simd_set1_d(-4.918853881490881290097e1);
2247
2248     const gmx_simd_double_t R4        = gmx_simd_set1_d(2.967721961301243206100e-3);
2249     const gmx_simd_double_t R3        = gmx_simd_set1_d(-5.634242780008963776856e-1);
2250     const gmx_simd_double_t R2        = gmx_simd_set1_d(6.968710824104713396794e0);
2251     const gmx_simd_double_t R1        = gmx_simd_set1_d(-2.556901049652824852289e1);
2252     const gmx_simd_double_t R0        = gmx_simd_set1_d(2.853665548261061424989e1);
2253
2254     const gmx_simd_double_t S3        = gmx_simd_set1_d(-2.194779531642920639778e1);
2255     const gmx_simd_double_t S2        = gmx_simd_set1_d(1.470656354026814941758e2);
2256     const gmx_simd_double_t S1        = gmx_simd_set1_d(-3.838770957603691357202e2);
2257     const gmx_simd_double_t S0        = gmx_simd_set1_d(3.424398657913078477438e2);
2258
2259     gmx_simd_double_t       xabs;
2260     gmx_simd_double_t       zz, ww, z, q, w, zz2, ww2;
2261     gmx_simd_double_t       PA, PB;
2262     gmx_simd_double_t       QA, QB;
2263     gmx_simd_double_t       RA, RB;
2264     gmx_simd_double_t       SA, SB;
2265     gmx_simd_double_t       nom, denom;
2266     gmx_simd_dbool_t        mask;
2267
2268     xabs  = gmx_simd_fabs_d(x);
2269
2270     mask  = gmx_simd_cmplt_d(limit1, xabs);
2271
2272     zz    = gmx_simd_sub_d(one, xabs);
2273     ww    = gmx_simd_mul_d(xabs, xabs);
2274     zz2   = gmx_simd_mul_d(zz, zz);
2275     ww2   = gmx_simd_mul_d(ww, ww);
2276
2277     /* R */
2278     RA    = gmx_simd_mul_d(R4, zz2);
2279     RB    = gmx_simd_mul_d(R3, zz2);
2280     RA    = gmx_simd_add_d(RA, R2);
2281     RB    = gmx_simd_add_d(RB, R1);
2282     RA    = gmx_simd_mul_d(RA, zz2);
2283     RB    = gmx_simd_mul_d(RB, zz);
2284     RA    = gmx_simd_add_d(RA, R0);
2285     RA    = gmx_simd_add_d(RA, RB);
2286
2287     /* S, SA = zz2 */
2288     SB    = gmx_simd_mul_d(S3, zz2);
2289     SA    = gmx_simd_add_d(zz2, S2);
2290     SB    = gmx_simd_add_d(SB, S1);
2291     SA    = gmx_simd_mul_d(SA, zz2);
2292     SB    = gmx_simd_mul_d(SB, zz);
2293     SA    = gmx_simd_add_d(SA, S0);
2294     SA    = gmx_simd_add_d(SA, SB);
2295
2296     /* P */
2297     PA    = gmx_simd_mul_d(P5, ww2);
2298     PB    = gmx_simd_mul_d(P4, ww2);
2299     PA    = gmx_simd_add_d(PA, P3);
2300     PB    = gmx_simd_add_d(PB, P2);
2301     PA    = gmx_simd_mul_d(PA, ww2);
2302     PB    = gmx_simd_mul_d(PB, ww2);
2303     PA    = gmx_simd_add_d(PA, P1);
2304     PB    = gmx_simd_add_d(PB, P0);
2305     PA    = gmx_simd_mul_d(PA, ww);
2306     PA    = gmx_simd_add_d(PA, PB);
2307
2308     /* Q, QA = ww2 */
2309     QB    = gmx_simd_mul_d(Q4, ww2);
2310     QA    = gmx_simd_add_d(ww2, Q3);
2311     QB    = gmx_simd_add_d(QB, Q2);
2312     QA    = gmx_simd_mul_d(QA, ww2);
2313     QB    = gmx_simd_mul_d(QB, ww2);
2314     QA    = gmx_simd_add_d(QA, Q1);
2315     QB    = gmx_simd_add_d(QB, Q0);
2316     QA    = gmx_simd_mul_d(QA, ww);
2317     QA    = gmx_simd_add_d(QA, QB);
2318
2319     RA    = gmx_simd_mul_d(RA, zz);
2320     PA    = gmx_simd_mul_d(PA, ww);
2321
2322     nom   = gmx_simd_blendv_d( PA, RA, mask );
2323     denom = gmx_simd_blendv_d( QA, SA, mask );
2324
2325     q     = gmx_simd_mul_d( nom, gmx_simd_inv_d(denom) );
2326
2327     zz    = gmx_simd_add_d(zz, zz);
2328     zz    = gmx_simd_sqrt_d(zz);
2329     z     = gmx_simd_sub_d(quarterpi, zz);
2330     zz    = gmx_simd_mul_d(zz, q);
2331     zz    = gmx_simd_sub_d(zz, morebits);
2332     z     = gmx_simd_sub_d(z, zz);
2333     z     = gmx_simd_add_d(z, quarterpi);
2334
2335     w     = gmx_simd_mul_d(xabs, q);
2336     w     = gmx_simd_add_d(w, xabs);
2337
2338     z     = gmx_simd_blendv_d( w, z, mask );
2339
2340     mask  = gmx_simd_cmplt_d(limit2, xabs);
2341     z     = gmx_simd_blendv_d( xabs, z, mask );
2342
2343     z = gmx_simd_xor_sign_d(z, x);
2344
2345     return z;
2346 }
2347
2348 /*! \brief SIMD double acos(x).
2349  *
2350  * \copydetails gmx_simd_acos_f
2351  */
2352 static gmx_inline gmx_simd_double_t
2353 gmx_simd_acos_d(gmx_simd_double_t x)
2354 {
2355     const gmx_simd_double_t one        = gmx_simd_set1_d(1.0);
2356     const gmx_simd_double_t half       = gmx_simd_set1_d(0.5);
2357     const gmx_simd_double_t quarterpi0 = gmx_simd_set1_d(7.85398163397448309616e-1);
2358     const gmx_simd_double_t quarterpi1 = gmx_simd_set1_d(6.123233995736765886130e-17);
2359
2360     gmx_simd_dbool_t        mask1;
2361     gmx_simd_double_t       z, z1, z2;
2362
2363     mask1 = gmx_simd_cmplt_d(half, x);
2364     z1    = gmx_simd_mul_d(half, gmx_simd_sub_d(one, x));
2365     z1    = gmx_simd_sqrt_d(z1);
2366     z     = gmx_simd_blendv_d( x, z1, mask1 );
2367
2368     z     = gmx_simd_asin_d(z);
2369
2370     z1    = gmx_simd_add_d(z, z);
2371
2372     z2    = gmx_simd_sub_d(quarterpi0, z);
2373     z2    = gmx_simd_add_d(z2, quarterpi1);
2374     z2    = gmx_simd_add_d(z2, quarterpi0);
2375
2376     z     = gmx_simd_blendv_d(z2, z1, mask1);
2377
2378     return z;
2379 }
2380
2381 /*! \brief SIMD double atan(x).
2382  *
2383  * \copydetails gmx_simd_atan_f
2384  */
2385 static gmx_inline gmx_simd_double_t
2386 gmx_simd_atan_d(gmx_simd_double_t x)
2387 {
2388     /* Same algorithm as cephes library */
2389     const gmx_simd_double_t limit1    = gmx_simd_set1_d(0.66);
2390     const gmx_simd_double_t limit2    = gmx_simd_set1_d(2.41421356237309504880);
2391     const gmx_simd_double_t quarterpi = gmx_simd_set1_d(M_PI/4.0);
2392     const gmx_simd_double_t halfpi    = gmx_simd_set1_d(M_PI/2.0);
2393     const gmx_simd_double_t mone      = gmx_simd_set1_d(-1.0);
2394     const gmx_simd_double_t morebits1 = gmx_simd_set1_d(0.5*6.123233995736765886130E-17);
2395     const gmx_simd_double_t morebits2 = gmx_simd_set1_d(6.123233995736765886130E-17);
2396
2397     const gmx_simd_double_t P4        = gmx_simd_set1_d(-8.750608600031904122785E-1);
2398     const gmx_simd_double_t P3        = gmx_simd_set1_d(-1.615753718733365076637E1);
2399     const gmx_simd_double_t P2        = gmx_simd_set1_d(-7.500855792314704667340E1);
2400     const gmx_simd_double_t P1        = gmx_simd_set1_d(-1.228866684490136173410E2);
2401     const gmx_simd_double_t P0        = gmx_simd_set1_d(-6.485021904942025371773E1);
2402
2403     const gmx_simd_double_t Q4        = gmx_simd_set1_d(2.485846490142306297962E1);
2404     const gmx_simd_double_t Q3        = gmx_simd_set1_d(1.650270098316988542046E2);
2405     const gmx_simd_double_t Q2        = gmx_simd_set1_d(4.328810604912902668951E2);
2406     const gmx_simd_double_t Q1        = gmx_simd_set1_d(4.853903996359136964868E2);
2407     const gmx_simd_double_t Q0        = gmx_simd_set1_d(1.945506571482613964425E2);
2408
2409     gmx_simd_double_t       y, xabs, t1, t2;
2410     gmx_simd_double_t       z, z2;
2411     gmx_simd_double_t       P_A, P_B, Q_A, Q_B;
2412     gmx_simd_dbool_t        mask1, mask2;
2413
2414     xabs   = gmx_simd_fabs_d(x);
2415
2416     mask1  = gmx_simd_cmplt_d(limit1, xabs);
2417     mask2  = gmx_simd_cmplt_d(limit2, xabs);
2418
2419     t1     = gmx_simd_mul_d(gmx_simd_add_d(xabs, mone), gmx_simd_inv_d(gmx_simd_sub_d(xabs, mone)));
2420     t2     = gmx_simd_mul_d(mone, gmx_simd_inv_d(xabs));
2421
2422     y      = gmx_simd_blendzero_d(quarterpi, mask1);
2423     y      = gmx_simd_blendv_d(y, halfpi, mask2);
2424     xabs   = gmx_simd_blendv_d(xabs, t1, mask1);
2425     xabs   = gmx_simd_blendv_d(xabs, t2, mask2);
2426
2427     z      = gmx_simd_mul_d(xabs, xabs);
2428     z2     = gmx_simd_mul_d(z, z);
2429
2430     P_A    = gmx_simd_mul_d(P4, z2);
2431     P_B    = gmx_simd_mul_d(P3, z2);
2432     P_A    = gmx_simd_add_d(P_A, P2);
2433     P_B    = gmx_simd_add_d(P_B, P1);
2434     P_A    = gmx_simd_mul_d(P_A, z2);
2435     P_B    = gmx_simd_mul_d(P_B, z);
2436     P_A    = gmx_simd_add_d(P_A, P0);
2437     P_A    = gmx_simd_add_d(P_A, P_B);
2438
2439     /* Q_A = z2 */
2440     Q_B    = gmx_simd_mul_d(Q4, z2);
2441     Q_A    = gmx_simd_add_d(z2, Q3);
2442     Q_B    = gmx_simd_add_d(Q_B, Q2);
2443     Q_A    = gmx_simd_mul_d(Q_A, z2);
2444     Q_B    = gmx_simd_mul_d(Q_B, z2);
2445     Q_A    = gmx_simd_add_d(Q_A, Q1);
2446     Q_B    = gmx_simd_add_d(Q_B, Q0);
2447     Q_A    = gmx_simd_mul_d(Q_A, z);
2448     Q_A    = gmx_simd_add_d(Q_A, Q_B);
2449
2450     z      = gmx_simd_mul_d(z, P_A);
2451     z      = gmx_simd_mul_d(z, gmx_simd_inv_d(Q_A));
2452     z      = gmx_simd_mul_d(z, xabs);
2453     z      = gmx_simd_add_d(z, xabs);
2454
2455     t1     = gmx_simd_blendzero_d(morebits1, mask1);
2456     t1     = gmx_simd_blendv_d(t1, morebits2, mask2);
2457
2458     z      = gmx_simd_add_d(z, t1);
2459     y      = gmx_simd_add_d(y, z);
2460
2461     y      = gmx_simd_xor_sign_d(y, x);
2462
2463     return y;
2464 }
2465
2466 /*! \brief SIMD double atan2(y,x).
2467  *
2468  * \copydetails gmx_simd_atan2_f
2469  */
2470 static gmx_inline gmx_simd_double_t
2471 gmx_simd_atan2_d(gmx_simd_double_t y, gmx_simd_double_t x)
2472 {
2473     const gmx_simd_double_t pi          = gmx_simd_set1_d(M_PI);
2474     const gmx_simd_double_t halfpi      = gmx_simd_set1_d(M_PI/2.0);
2475     gmx_simd_double_t       xinv, p, aoffset;
2476     gmx_simd_dbool_t        mask_x0, mask_y0, mask_xlt0, mask_ylt0;
2477
2478     mask_x0   = gmx_simd_cmpeq_d(x, gmx_simd_setzero_d());
2479     mask_y0   = gmx_simd_cmpeq_d(y, gmx_simd_setzero_d());
2480     mask_xlt0 = gmx_simd_cmplt_d(x, gmx_simd_setzero_d());
2481     mask_ylt0 = gmx_simd_cmplt_d(y, gmx_simd_setzero_d());
2482
2483     aoffset   = gmx_simd_blendzero_d(halfpi, mask_x0);
2484     aoffset   = gmx_simd_blendnotzero_d(aoffset, mask_y0);
2485
2486     aoffset   = gmx_simd_blendv_d(aoffset, pi, mask_xlt0);
2487     aoffset   = gmx_simd_blendv_d(aoffset, gmx_simd_fneg_d(aoffset), mask_ylt0);
2488
2489     xinv      = gmx_simd_blendnotzero_d(gmx_simd_inv_d(x), mask_x0);
2490     p         = gmx_simd_mul_d(y, xinv);
2491     p         = gmx_simd_atan_d(p);
2492     p         = gmx_simd_add_d(p, aoffset);
2493
2494     return p;
2495 }
2496
2497
2498 /*! \brief Calculate the force correction due to PME analytically for SIMD double.
2499  *
2500  * \copydetails gmx_simd_pmecorrF_f
2501  */
2502 static gmx_simd_double_t
2503 gmx_simd_pmecorrF_d(gmx_simd_double_t z2)
2504 {
2505     const gmx_simd_double_t  FN10     = gmx_simd_set1_d(-8.0072854618360083154e-14);
2506     const gmx_simd_double_t  FN9      = gmx_simd_set1_d(1.1859116242260148027e-11);
2507     const gmx_simd_double_t  FN8      = gmx_simd_set1_d(-8.1490406329798423616e-10);
2508     const gmx_simd_double_t  FN7      = gmx_simd_set1_d(3.4404793543907847655e-8);
2509     const gmx_simd_double_t  FN6      = gmx_simd_set1_d(-9.9471420832602741006e-7);
2510     const gmx_simd_double_t  FN5      = gmx_simd_set1_d(0.000020740315999115847456);
2511     const gmx_simd_double_t  FN4      = gmx_simd_set1_d(-0.00031991745139313364005);
2512     const gmx_simd_double_t  FN3      = gmx_simd_set1_d(0.0035074449373659008203);
2513     const gmx_simd_double_t  FN2      = gmx_simd_set1_d(-0.031750380176100813405);
2514     const gmx_simd_double_t  FN1      = gmx_simd_set1_d(0.13884101728898463426);
2515     const gmx_simd_double_t  FN0      = gmx_simd_set1_d(-0.75225277815249618847);
2516
2517     const gmx_simd_double_t  FD5      = gmx_simd_set1_d(0.000016009278224355026701);
2518     const gmx_simd_double_t  FD4      = gmx_simd_set1_d(0.00051055686934806966046);
2519     const gmx_simd_double_t  FD3      = gmx_simd_set1_d(0.0081803507497974289008);
2520     const gmx_simd_double_t  FD2      = gmx_simd_set1_d(0.077181146026670287235);
2521     const gmx_simd_double_t  FD1      = gmx_simd_set1_d(0.41543303143712535988);
2522     const gmx_simd_double_t  FD0      = gmx_simd_set1_d(1.0);
2523
2524     gmx_simd_double_t        z4;
2525     gmx_simd_double_t        polyFN0, polyFN1, polyFD0, polyFD1;
2526
2527     z4             = gmx_simd_mul_d(z2, z2);
2528
2529     polyFD1        = gmx_simd_fmadd_d(FD5, z4, FD3);
2530     polyFD1        = gmx_simd_fmadd_d(polyFD1, z4, FD1);
2531     polyFD1        = gmx_simd_mul_d(polyFD1, z2);
2532     polyFD0        = gmx_simd_fmadd_d(FD4, z4, FD2);
2533     polyFD0        = gmx_simd_fmadd_d(polyFD0, z4, FD0);
2534     polyFD0        = gmx_simd_add_d(polyFD0, polyFD1);
2535
2536     polyFD0        = gmx_simd_inv_d(polyFD0);
2537
2538     polyFN0        = gmx_simd_fmadd_d(FN10, z4, FN8);
2539     polyFN0        = gmx_simd_fmadd_d(polyFN0, z4, FN6);
2540     polyFN0        = gmx_simd_fmadd_d(polyFN0, z4, FN4);
2541     polyFN0        = gmx_simd_fmadd_d(polyFN0, z4, FN2);
2542     polyFN0        = gmx_simd_fmadd_d(polyFN0, z4, FN0);
2543     polyFN1        = gmx_simd_fmadd_d(FN9, z4, FN7);
2544     polyFN1        = gmx_simd_fmadd_d(polyFN1, z4, FN5);
2545     polyFN1        = gmx_simd_fmadd_d(polyFN1, z4, FN3);
2546     polyFN1        = gmx_simd_fmadd_d(polyFN1, z4, FN1);
2547     polyFN0        = gmx_simd_fmadd_d(polyFN1, z2, polyFN0);
2548
2549
2550     return gmx_simd_mul_d(polyFN0, polyFD0);
2551 }
2552
2553
2554
2555 /*! \brief Calculate the potential correction due to PME analytically for SIMD double.
2556  *
2557  * \copydetails gmx_simd_pmecorrV_f
2558  */
2559 static gmx_simd_double_t
2560 gmx_simd_pmecorrV_d(gmx_simd_double_t z2)
2561 {
2562     const gmx_simd_double_t  VN9      = gmx_simd_set1_d(-9.3723776169321855475e-13);
2563     const gmx_simd_double_t  VN8      = gmx_simd_set1_d(1.2280156762674215741e-10);
2564     const gmx_simd_double_t  VN7      = gmx_simd_set1_d(-7.3562157912251309487e-9);
2565     const gmx_simd_double_t  VN6      = gmx_simd_set1_d(2.6215886208032517509e-7);
2566     const gmx_simd_double_t  VN5      = gmx_simd_set1_d(-4.9532491651265819499e-6);
2567     const gmx_simd_double_t  VN4      = gmx_simd_set1_d(0.00025907400778966060389);
2568     const gmx_simd_double_t  VN3      = gmx_simd_set1_d(0.0010585044856156469792);
2569     const gmx_simd_double_t  VN2      = gmx_simd_set1_d(0.045247661136833092885);
2570     const gmx_simd_double_t  VN1      = gmx_simd_set1_d(0.11643931522926034421);
2571     const gmx_simd_double_t  VN0      = gmx_simd_set1_d(1.1283791671726767970);
2572
2573     const gmx_simd_double_t  VD5      = gmx_simd_set1_d(0.000021784709867336150342);
2574     const gmx_simd_double_t  VD4      = gmx_simd_set1_d(0.00064293662010911388448);
2575     const gmx_simd_double_t  VD3      = gmx_simd_set1_d(0.0096311444822588683504);
2576     const gmx_simd_double_t  VD2      = gmx_simd_set1_d(0.085608012351550627051);
2577     const gmx_simd_double_t  VD1      = gmx_simd_set1_d(0.43652499166614811084);
2578     const gmx_simd_double_t  VD0      = gmx_simd_set1_d(1.0);
2579
2580     gmx_simd_double_t        z4;
2581     gmx_simd_double_t        polyVN0, polyVN1, polyVD0, polyVD1;
2582
2583     z4             = gmx_simd_mul_d(z2, z2);
2584
2585     polyVD1        = gmx_simd_fmadd_d(VD5, z4, VD3);
2586     polyVD0        = gmx_simd_fmadd_d(VD4, z4, VD2);
2587     polyVD1        = gmx_simd_fmadd_d(polyVD1, z4, VD1);
2588     polyVD0        = gmx_simd_fmadd_d(polyVD0, z4, VD0);
2589     polyVD0        = gmx_simd_fmadd_d(polyVD1, z2, polyVD0);
2590
2591     polyVD0        = gmx_simd_inv_d(polyVD0);
2592
2593     polyVN1        = gmx_simd_fmadd_d(VN9, z4, VN7);
2594     polyVN0        = gmx_simd_fmadd_d(VN8, z4, VN6);
2595     polyVN1        = gmx_simd_fmadd_d(polyVN1, z4, VN5);
2596     polyVN0        = gmx_simd_fmadd_d(polyVN0, z4, VN4);
2597     polyVN1        = gmx_simd_fmadd_d(polyVN1, z4, VN3);
2598     polyVN0        = gmx_simd_fmadd_d(polyVN0, z4, VN2);
2599     polyVN1        = gmx_simd_fmadd_d(polyVN1, z4, VN1);
2600     polyVN0        = gmx_simd_fmadd_d(polyVN0, z4, VN0);
2601     polyVN0        = gmx_simd_fmadd_d(polyVN1, z2, polyVN0);
2602
2603     return gmx_simd_mul_d(polyVN0, polyVD0);
2604 }
2605
2606 /*! \} */
2607
2608 #endif
2609
2610
2611 /*! \name SIMD4 math functions
2612  *
2613  * \note Only a subset of the math functions are implemented for SIMD4.
2614  *  \{
2615  */
2616
2617
2618 #ifdef GMX_SIMD4_HAVE_FLOAT
2619
2620 /*************************************************************************
2621  * SINGLE PRECISION SIMD4 MATH FUNCTIONS - JUST A SMALL SUBSET SUPPORTED *
2622  *************************************************************************/
2623
2624 /*! \brief SIMD4 utility function to sum a+b+c+d for SIMD4 floats.
2625  *
2626  * \copydetails gmx_simd_sum4_f
2627  */
2628 static gmx_inline gmx_simd4_float_t
2629 gmx_simd4_sum4_f(gmx_simd4_float_t a, gmx_simd4_float_t b,
2630                  gmx_simd4_float_t c, gmx_simd4_float_t d)
2631 {
2632     return gmx_simd4_add_f(gmx_simd4_add_f(a, b), gmx_simd4_add_f(c, d));
2633 }
2634
2635 /*! \brief Perform one Newton-Raphson iteration to improve 1/sqrt(x) for SIMD4 float.
2636  *
2637  * \copydetails gmx_simd_rsqrt_iter_f
2638  */
2639 static gmx_inline gmx_simd4_float_t
2640 gmx_simd4_rsqrt_iter_f(gmx_simd4_float_t lu, gmx_simd4_float_t x)
2641 {
2642 #    ifdef GMX_SIMD_HAVE_FMA
2643     return gmx_simd4_fmadd_f(gmx_simd4_fnmadd_f(x, gmx_simd4_mul_f(lu, lu), gmx_simd4_set1_f(1.0f)), gmx_simd4_mul_f(lu, gmx_simd4_set1_f(0.5f)), lu);
2644 #    else
2645     return gmx_simd4_mul_f(gmx_simd4_set1_f(0.5f), gmx_simd4_mul_f(gmx_simd4_sub_f(gmx_simd4_set1_f(3.0f), gmx_simd4_mul_f(gmx_simd4_mul_f(lu, lu), x)), lu));
2646 #    endif
2647 }
2648
2649 /*! \brief Calculate 1/sqrt(x) for SIMD4 float.
2650  *
2651  * \copydetails gmx_simd_invsqrt_f
2652  */
2653 static gmx_inline gmx_simd4_float_t
2654 gmx_simd4_invsqrt_f(gmx_simd4_float_t x)
2655 {
2656     gmx_simd4_float_t lu = gmx_simd4_rsqrt_f(x);
2657 #if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
2658     lu = gmx_simd4_rsqrt_iter_f(lu, x);
2659 #endif
2660 #if (GMX_SIMD_RSQRT_BITS*2 < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
2661     lu = gmx_simd4_rsqrt_iter_f(lu, x);
2662 #endif
2663 #if (GMX_SIMD_RSQRT_BITS*4 < GMX_SIMD_MATH_TARGET_SINGLE_BITS)
2664     lu = gmx_simd4_rsqrt_iter_f(lu, x);
2665 #endif
2666     return lu;
2667 }
2668
2669 #endif /* GMX_SIMD4_HAVE_FLOAT */
2670
2671
2672
2673 #ifdef GMX_SIMD4_HAVE_DOUBLE
2674 /*************************************************************************
2675  * DOUBLE PRECISION SIMD4 MATH FUNCTIONS - JUST A SMALL SUBSET SUPPORTED *
2676  *************************************************************************/
2677
2678
2679 /*! \brief SIMD4 utility function to sum a+b+c+d for SIMD4 doubles.
2680  *
2681  * \copydetails gmx_simd_sum4_f
2682  */
2683 static gmx_inline gmx_simd4_double_t
2684 gmx_simd4_sum4_d(gmx_simd4_double_t a, gmx_simd4_double_t b,
2685                  gmx_simd4_double_t c, gmx_simd4_double_t d)
2686 {
2687     return gmx_simd4_add_d(gmx_simd4_add_d(a, b), gmx_simd4_add_d(c, d));
2688 }
2689
2690 /*! \brief Perform one Newton-Raphson iteration to improve 1/sqrt(x) for SIMD4 double.
2691  *
2692  * \copydetails gmx_simd_rsqrt_iter_f
2693  */
2694 static gmx_inline gmx_simd4_double_t
2695 gmx_simd4_rsqrt_iter_d(gmx_simd4_double_t lu, gmx_simd4_double_t x)
2696 {
2697 #ifdef GMX_SIMD_HAVE_FMA
2698     return gmx_simd4_fmadd_d(gmx_simd4_fnmadd_d(x, gmx_simd4_mul_d(lu, lu), gmx_simd4_set1_d(1.0)), gmx_simd4_mul_d(lu, gmx_simd4_set1_d(0.5)), lu);
2699 #else
2700     return gmx_simd4_mul_d(gmx_simd4_set1_d(0.5), gmx_simd4_mul_d(gmx_simd4_sub_d(gmx_simd4_set1_d(3.0), gmx_simd4_mul_d(gmx_simd4_mul_d(lu, lu), x)), lu));
2701 #endif
2702 }
2703
2704 /*! \brief Calculate 1/sqrt(x) for SIMD4 double.
2705  *
2706  * \copydetails gmx_simd_invsqrt_f
2707  */
2708 static gmx_inline gmx_simd4_double_t
2709 gmx_simd4_invsqrt_d(gmx_simd4_double_t x)
2710 {
2711     gmx_simd4_double_t lu = gmx_simd4_rsqrt_d(x);
2712 #if (GMX_SIMD_RSQRT_BITS < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
2713     lu = gmx_simd4_rsqrt_iter_d(lu, x);
2714 #endif
2715 #if (GMX_SIMD_RSQRT_BITS*2 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
2716     lu = gmx_simd4_rsqrt_iter_d(lu, x);
2717 #endif
2718 #if (GMX_SIMD_RSQRT_BITS*4 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
2719     lu = gmx_simd4_rsqrt_iter_d(lu, x);
2720 #endif
2721 #if (GMX_SIMD_RSQRT_BITS*8 < GMX_SIMD_MATH_TARGET_DOUBLE_BITS)
2722     lu = gmx_simd4_rsqrt_iter_d(lu, x);
2723 #endif
2724     return lu;
2725 }
2726 #endif /* GMX_SIMD4_HAVE_DOUBLE */
2727
2728 /*! \} */
2729
2730
2731 /* Set defines based on default Gromacs precision */
2732 #ifdef GMX_DOUBLE
2733 /* Documentation in single branch below */
2734 #    define gmx_simd_sum4_r           gmx_simd_sum4_d
2735 #    define gmx_simd_xor_sign_r       gmx_simd_xor_sign_d
2736 #    define gmx_simd_invsqrt_r        gmx_simd_invsqrt_d
2737 #    define gmx_simd_invsqrt_pair_r   gmx_simd_invsqrt_pair_d
2738 #    define gmx_simd_sqrt_r           gmx_simd_sqrt_d
2739 #    define gmx_simd_inv_r            gmx_simd_inv_d
2740 #    define gmx_simd_log_r            gmx_simd_log_d
2741 #    define gmx_simd_exp2_r           gmx_simd_exp2_d
2742 #    define gmx_simd_exp_r            gmx_simd_exp_d
2743 #    define gmx_simd_erf_r            gmx_simd_erf_d
2744 #    define gmx_simd_erfc_r           gmx_simd_erfc_d
2745 #    define gmx_simd_sincos_r         gmx_simd_sincos_d
2746 #    define gmx_simd_sin_r            gmx_simd_sin_d
2747 #    define gmx_simd_cos_r            gmx_simd_cos_d
2748 #    define gmx_simd_tan_r            gmx_simd_tan_d
2749 #    define gmx_simd_asin_r           gmx_simd_asin_d
2750 #    define gmx_simd_acos_r           gmx_simd_acos_d
2751 #    define gmx_simd_atan_r           gmx_simd_atan_d
2752 #    define gmx_simd_atan2_r          gmx_simd_atan2_d
2753 #    define gmx_simd_pmecorrF_r       gmx_simd_pmecorrF_d
2754 #    define gmx_simd_pmecorrV_r       gmx_simd_pmecorrV_d
2755 #    define gmx_simd4_sum4_r          gmx_simd4_sum4_d
2756 #    define gmx_simd4_invsqrt_r       gmx_simd4_invsqrt_d
2757
2758 #else /* GMX_DOUBLE */
2759
2760 /*! \name Real-precision SIMD math functions
2761  *
2762  *  These are the ones you should typically call in Gromacs.
2763  * \{
2764  */
2765
2766 /*! \brief SIMD utility function to sum a+b+c+d for SIMD reals.
2767  *
2768  * \copydetails gmx_simd_sum4_f
2769  */
2770 #    define gmx_simd_sum4_r           gmx_simd_sum4_f
2771
2772 /*! \brief Return -a if b is negative, SIMD real.
2773  *
2774  * \copydetails gmx_simd_xor_sign_f
2775  */
2776 #    define gmx_simd_xor_sign_r       gmx_simd_xor_sign_f
2777
2778 /*! \brief Calculate 1/sqrt(x) for SIMD real.
2779  *
2780  * \copydetails gmx_simd_invsqrt_f
2781  */
2782 #    define gmx_simd_invsqrt_r        gmx_simd_invsqrt_f
2783
2784 /*! \brief Calculate 1/sqrt(x) for two SIMD reals.
2785  *
2786  * \copydetails gmx_simd_invsqrt_pair_f
2787  */
2788 #    define gmx_simd_invsqrt_pair_r   gmx_simd_invsqrt_pair_f
2789
2790 /*! \brief Calculate sqrt(x) correctly for SIMD real, including argument 0.0.
2791  *
2792  * \copydetails gmx_simd_sqrt_f
2793  */
2794 #    define gmx_simd_sqrt_r           gmx_simd_sqrt_f
2795
2796 /*! \brief Calculate 1/x for SIMD real.
2797  *
2798  * \copydetails gmx_simd_inv_f
2799  */
2800 #    define gmx_simd_inv_r            gmx_simd_inv_f
2801
2802 /*! \brief SIMD real log(x). This is the natural logarithm.
2803  *
2804  * \copydetails gmx_simd_log_f
2805  */
2806 #    define gmx_simd_log_r            gmx_simd_log_f
2807
2808 /*! \brief SIMD real 2^x.
2809  *
2810  * \copydetails gmx_simd_exp2_f
2811  */
2812 #    define gmx_simd_exp2_r           gmx_simd_exp2_f
2813
2814 /*! \brief SIMD real e^x.
2815  *
2816  * \copydetails gmx_simd_exp_f
2817  */
2818 #    define gmx_simd_exp_r            gmx_simd_exp_f
2819
2820 /*! \brief SIMD real erf(x).
2821  *
2822  * \copydetails gmx_simd_erf_f
2823  */
2824 #    define gmx_simd_erf_r            gmx_simd_erf_f
2825
2826 /*! \brief SIMD real erfc(x).
2827  *
2828  * \copydetails gmx_simd_erfc_f
2829  */
2830 #    define gmx_simd_erfc_r           gmx_simd_erfc_f
2831
2832 /*! \brief SIMD real sin \& cos.
2833  *
2834  * \copydetails gmx_simd_sincos_f
2835  */
2836 #    define gmx_simd_sincos_r         gmx_simd_sincos_f
2837
2838 /*! \brief SIMD real sin(x).
2839  *
2840  * \copydetails gmx_simd_sin_f
2841  */
2842 #    define gmx_simd_sin_r            gmx_simd_sin_f
2843
2844 /*! \brief SIMD real cos(x).
2845  *
2846  * \copydetails gmx_simd_cos_f
2847  */
2848 #    define gmx_simd_cos_r            gmx_simd_cos_f
2849
2850 /*! \brief SIMD real tan(x).
2851  *
2852  * \copydetails gmx_simd_tan_f
2853  */
2854 #    define gmx_simd_tan_r            gmx_simd_tan_f
2855
2856 /*! \brief SIMD real asin(x).
2857  *
2858  * \copydetails gmx_simd_asin_f
2859  */
2860 #    define gmx_simd_asin_r           gmx_simd_asin_f
2861
2862 /*! \brief SIMD real acos(x).
2863  *
2864  * \copydetails gmx_simd_acos_f
2865  */
2866 #    define gmx_simd_acos_r           gmx_simd_acos_f
2867
2868 /*! \brief SIMD real atan(x).
2869  *
2870  * \copydetails gmx_simd_atan_f
2871  */
2872 #    define gmx_simd_atan_r           gmx_simd_atan_f
2873
2874 /*! \brief SIMD real atan2(y,x).
2875  *
2876  * \copydetails gmx_simd_atan2_f
2877  */
2878 #    define gmx_simd_atan2_r          gmx_simd_atan2_f
2879
2880 /*! \brief SIMD Analytic PME force correction.
2881  *
2882  * \copydetails gmx_simd_pmecorrF_f
2883  */
2884 #    define gmx_simd_pmecorrF_r       gmx_simd_pmecorrF_f
2885
2886 /*! \brief SIMD Analytic PME potential correction.
2887  *
2888  * \copydetails gmx_simd_pmecorrV_f
2889  */
2890 #    define gmx_simd_pmecorrV_r       gmx_simd_pmecorrV_f
2891
2892 /*! \}
2893  * \name SIMD4 math functions
2894  * \{
2895  */
2896
2897 /*! \brief SIMD4 utility function to sum a+b+c+d for SIMD4 reals.
2898  *
2899  * \copydetails gmx_simd_sum4_f
2900  */
2901 #    define gmx_simd4_sum4_r          gmx_simd4_sum4_f
2902
2903 /*! \brief Calculate 1/sqrt(x) for SIMD4 real.
2904  *
2905  * \copydetails gmx_simd_invsqrt_f
2906  */
2907 #    define gmx_simd4_invsqrt_r       gmx_simd4_invsqrt_f
2908
2909 /*! \} */
2910
2911 #endif /* GMX_DOUBLE */
2912
2913 /*! \} */
2914 /*! \endcond */
2915
2916 #endif /* GMX_SIMD_SIMD_MATH_H_ */