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