src/gromacs/legacyheaders/gmx_math_x86_avx_256_double.h

   1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  * This file is part of GROMACS.
   5  * Copyright (c) 2012-
   6  *
   7  * Written by the Gromacs development team under coordination of
   8  * David van der Spoel, Berk Hess, and Erik Lindahl.
   9  *
  10  * This library is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2
  13  * of the License, or (at your option) any later version.
  14  *
  15  * To help us fund GROMACS development, we humbly ask that you cite
  16  * the research papers on the package. Check out http://www.gromacs.org
  17  *
  18  * And Hey:
  19  * Gnomes, ROck Monsters And Chili Sauce
  20  */
  21 #ifndef _gmx_math_x86_avx_256_double_h_
  22 #define _gmx_math_x86_avx_256_double_h_
  23
  24 #include <math.h>
  25
  26 #include "gmx_x86_avx_256.h"
  27
  28 #ifndef M_PI
  29 #  define M_PI 3.14159265358979323846264338327950288
  30 #endif
  31
  32
  33 /************************
  34  *                      *
  35  * Simple math routines *
  36  *                      *
  37  ************************/
  38
  39 /* 1.0/sqrt(x), 256 bit wide */
  40 static gmx_inline __m256d
  41 gmx_mm256_invsqrt_pd(__m256d x)
  42 {
  43     const __m256d half  = _mm256_set1_pd(0.5);
  44     const __m256d three = _mm256_set1_pd(3.0);
  45
  46     /* Lookup instruction only exists in single precision, convert back and forth... */
  47     __m256d lu = _mm256_cvtps_pd(_mm_rsqrt_ps( _mm256_cvtpd_ps(x)));
  48
  49     lu = _mm256_mul_pd(half, _mm256_mul_pd(_mm256_sub_pd(three, _mm256_mul_pd(_mm256_mul_pd(lu, lu), x)), lu));
  50     return _mm256_mul_pd(half, _mm256_mul_pd(_mm256_sub_pd(three, _mm256_mul_pd(_mm256_mul_pd(lu, lu), x)), lu));
  51 }
  52
  53 /* 1.0/sqrt(x), done for a pair of arguments to improve throughput */
  54 static void
  55 gmx_mm256_invsqrt_pair_pd(__m256d x1, __m256d x2, __m256d *invsqrt1, __m256d *invsqrt2)
  56 {
  57     const __m256d half   = _mm256_set1_pd(0.5);
  58     const __m256d three  = _mm256_set1_pd(3.0);
  59     const __m256  halff  = _mm256_set1_ps(0.5f);
  60     const __m256  threef = _mm256_set1_ps(3.0f);
  61
  62     __m256        xf, luf;
  63     __m256d       lu1, lu2;
  64
  65     /* Do first N-R step in float for 2x throughput */
  66     xf  = _mm256_insertf128_ps(_mm256_castps128_ps256(_mm256_cvtpd_ps(x1)), _mm256_cvtpd_ps(x2), 0x1);
  67     luf = _mm256_rsqrt_ps(xf);
  68
  69     luf = _mm256_mul_ps(halff, _mm256_mul_ps(_mm256_sub_ps(threef, _mm256_mul_ps(_mm256_mul_ps(luf, luf), xf)), luf));
  70
  71     lu2 = _mm256_cvtps_pd(_mm256_extractf128_ps(luf, 0x1));
  72     lu1 = _mm256_cvtps_pd(_mm256_castps256_ps128(luf));
  73
  74     *invsqrt1 = _mm256_mul_pd(half, _mm256_mul_pd(_mm256_sub_pd(three, _mm256_mul_pd(_mm256_mul_pd(lu1, lu1), x1)), lu1));
  75     *invsqrt2 = _mm256_mul_pd(half, _mm256_mul_pd(_mm256_sub_pd(three, _mm256_mul_pd(_mm256_mul_pd(lu2, lu2), x2)), lu2));
  76 }
  77
  78 /* 1.0/sqrt(x), 128 bit wide */
  79 static gmx_inline __m128d
  80 gmx_mm_invsqrt_pd(__m128d x)
  81 {
  82     const __m128d half  = _mm_set1_pd(0.5);
  83     const __m128d three = _mm_set1_pd(3.0);
  84
  85     /* Lookup instruction only exists in single precision, convert back and forth... */
  86     __m128d lu = _mm_cvtps_pd(_mm_rsqrt_ps( _mm_cvtpd_ps(x)));
  87
  88     lu = _mm_mul_pd(half, _mm_mul_pd(_mm_sub_pd(three, _mm_mul_pd(_mm_mul_pd(lu, lu), x)), lu));
  89     return _mm_mul_pd(half, _mm_mul_pd(_mm_sub_pd(three, _mm_mul_pd(_mm_mul_pd(lu, lu), x)), lu));
  90 }
  91
  92 /* 1.0/sqrt(x), done for two pairs to improve throughput */
  93 static void
  94 gmx_mm_invsqrt_pair_pd(__m128d x1, __m128d x2, __m128d *invsqrt1, __m128d *invsqrt2)
  95 {
  96     const __m128d half   = _mm_set1_pd(0.5);
  97     const __m128d three  = _mm_set1_pd(3.0);
  98     const __m128  halff  = _mm_set1_ps(0.5f);
  99     const __m128  threef = _mm_set1_ps(3.0f);
 100
 101     __m128        xf, luf;
 102     __m128d       lu1, lu2;
 103
 104     /* Do first N-R step in float for 2x throughput */
 105     xf  = _mm_shuffle_ps(_mm_cvtpd_ps(x1), _mm_cvtpd_ps(x2), _MM_SHUFFLE(1, 0, 1, 0));
 106     luf = _mm_rsqrt_ps(xf);
 107     luf = _mm_mul_ps(halff, _mm_mul_ps(_mm_sub_ps(threef, _mm_mul_ps(_mm_mul_ps(luf, luf), xf)), luf));
 108
 109     lu2 = _mm_cvtps_pd(_mm_shuffle_ps(luf, luf, _MM_SHUFFLE(3, 2, 3, 2)));
 110     lu1 = _mm_cvtps_pd(luf);
 111
 112     *invsqrt1 = _mm_mul_pd(half, _mm_mul_pd(_mm_sub_pd(three, _mm_mul_pd(_mm_mul_pd(lu1, lu1), x1)), lu1));
 113     *invsqrt2 = _mm_mul_pd(half, _mm_mul_pd(_mm_sub_pd(three, _mm_mul_pd(_mm_mul_pd(lu2, lu2), x2)), lu2));
 114 }
 115
 116 /* sqrt(x) (256 bit)- Do NOT use this (but rather invsqrt) if you actually need 1.0/sqrt(x) */
 117 static gmx_inline __m256d
 118 gmx_mm256_sqrt_pd(__m256d x)
 119 {
 120     __m256d mask;
 121     __m256d res;
 122
 123     mask = _mm256_cmp_pd(x, _mm256_setzero_pd(), _CMP_EQ_OQ);
 124     res  = _mm256_andnot_pd(mask, gmx_mm256_invsqrt_pd(x));
 125
 126     res  = _mm256_mul_pd(x, res);
 127
 128     return res;
 129 }
 130
 131 /* sqrt(x) (128 bit) - Do NOT use this (but rather invsqrt) if you actually need 1.0/sqrt(x) */
 132 static gmx_inline __m128d
 133 gmx_mm_sqrt_pd(__m128d x)
 134 {
 135     __m128d mask;
 136     __m128d res;
 137
 138     mask = _mm_cmpeq_pd(x, _mm_setzero_pd());
 139     res  = _mm_andnot_pd(mask, gmx_mm_invsqrt_pd(x));
 140
 141     res  = _mm_mul_pd(x, res);
 142
 143     return res;
 144 }
 145
 146
 147 /* 1.0/x, 256 bit wide */
 148 static gmx_inline __m256d
 149 gmx_mm256_inv_pd(__m256d x)
 150 {
 151     const __m256d two  = _mm256_set1_pd(2.0);
 152
 153     /* Lookup instruction only exists in single precision, convert back and forth... */
 154     __m256d lu = _mm256_cvtps_pd(_mm_rcp_ps( _mm256_cvtpd_ps(x)));
 155
 156     /* Perform two N-R steps for double precision */
 157     lu         = _mm256_mul_pd(lu, _mm256_sub_pd(two, _mm256_mul_pd(x, lu)));
 158     return _mm256_mul_pd(lu, _mm256_sub_pd(two, _mm256_mul_pd(x, lu)));
 159 }
 160
 161 /* 1.0/x, 128 bit */
 162 static gmx_inline __m128d
 163 gmx_mm_inv_pd(__m128d x)
 164 {
 165     const __m128d two  = _mm_set1_pd(2.0);
 166
 167     /* Lookup instruction only exists in single precision, convert back and forth... */
 168     __m128d lu = _mm_cvtps_pd(_mm_rcp_ps( _mm_cvtpd_ps(x)));
 169
 170     /* Perform two N-R steps for double precision */
 171     lu         = _mm_mul_pd(lu, _mm_sub_pd(two, _mm_mul_pd(x, lu)));
 172     return _mm_mul_pd(lu, _mm_sub_pd(two, _mm_mul_pd(x, lu)));
 173 }
 174
 175
 176 static gmx_inline __m256d
 177 gmx_mm256_abs_pd(__m256d x)
 178 {
 179     const __m256d signmask  = _mm256_castsi256_pd( _mm256_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF,
 180                                                                     0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
 181
 182     return _mm256_and_pd(x, signmask);
 183 }
 184
 185 static gmx_inline __m128d
 186 gmx_mm_abs_pd(__m128d x)
 187 {
 188     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
 189
 190     return _mm_and_pd(x, signmask);
 191 }
 192
 193
 194 /*
 195  * 2^x function, 256 bit
 196  *
 197  * The 2^w term is calculated from a (6,0)-th order (no denominator) Minimax polynomia on the interval
 198  * [-0.5,0.5].
 199  *
 200  * The approximation on [-0.5,0.5] is a rational Padé approximation, 1+2*P(x^2)/(Q(x^2)-P(x^2)),
 201  * according to the same algorithm as used in the Cephes/netlib math routines.
 202  */
 203 static __m256d
 204 gmx_mm256_exp2_pd(__m256d x)
 205 {
 206     /* Lower bound: We do not allow numbers that would lead to an IEEE fp representation exponent smaller than -126. */
 207     const __m256d arglimit = _mm256_set1_pd(1022.0);
 208     const __m128i expbase  = _mm_set1_epi32(1023);
 209
 210     const __m256d P2       = _mm256_set1_pd(2.30933477057345225087e-2);
 211     const __m256d P1       = _mm256_set1_pd(2.02020656693165307700e1);
 212     const __m256d P0       = _mm256_set1_pd(1.51390680115615096133e3);
 213     /* Q2 == 1.0 */
 214     const __m256d Q1       = _mm256_set1_pd(2.33184211722314911771e2);
 215     const __m256d Q0       = _mm256_set1_pd(4.36821166879210612817e3);
 216     const __m256d one      = _mm256_set1_pd(1.0);
 217     const __m256d two      = _mm256_set1_pd(2.0);
 218
 219     __m256d       valuemask;
 220     __m256i       iexppart;
 221     __m128i       iexppart128a, iexppart128b;
 222     __m256d       fexppart;
 223     __m256d       intpart;
 224     __m256d       z, z2;
 225     __m256d       PolyP, PolyQ;
 226
 227     iexppart128a  = _mm256_cvtpd_epi32(x);
 228     intpart       = _mm256_round_pd(x, _MM_FROUND_TO_NEAREST_INT);
 229
 230     /* Add exponent bias */
 231     iexppart128a   = _mm_add_epi32(iexppart128a, expbase);
 232
 233     /* We now want to shift the exponent 52 positions left, but to achieve this we need
 234      * to separate the 128-bit register data into two registers (4x64-bit > 128bit)
 235      * shift them, and then merge into a single __m256d.
 236      * Elements 0/1 should end up in iexppart128a, and 2/3 in iexppart128b.
 237      * It doesnt matter what we put in the 2nd/4th position, since that data will be
 238      * shifted out and replaced with zeros.
 239      */
 240     iexppart128b   = _mm_shuffle_epi32(iexppart128a, _MM_SHUFFLE(3, 3, 2, 2));
 241     iexppart128a   = _mm_shuffle_epi32(iexppart128a, _MM_SHUFFLE(1, 1, 0, 0));
 242
 243     iexppart128b   = _mm_slli_epi64(iexppart128b, 52);
 244     iexppart128a   = _mm_slli_epi64(iexppart128a, 52);
 245
 246     iexppart  = _mm256_castsi128_si256(iexppart128a);
 247     iexppart  = _mm256_insertf128_si256(iexppart, iexppart128b, 0x1);
 248
 249     valuemask = _mm256_cmp_pd(arglimit, gmx_mm256_abs_pd(x), _CMP_GE_OQ);
 250     fexppart  = _mm256_and_pd(valuemask, _mm256_castsi256_pd(iexppart));
 251
 252     z         = _mm256_sub_pd(x, intpart);
 253
 254     z2        = _mm256_mul_pd(z, z);
 255
 256     PolyP     = _mm256_mul_pd(P2, z2);
 257     PolyP     = _mm256_add_pd(PolyP, P1);
 258     PolyQ     = _mm256_add_pd(z2, Q1);
 259     PolyP     = _mm256_mul_pd(PolyP, z2);
 260     PolyQ     = _mm256_mul_pd(PolyQ, z2);
 261     PolyP     = _mm256_add_pd(PolyP, P0);
 262     PolyQ     = _mm256_add_pd(PolyQ, Q0);
 263     PolyP     = _mm256_mul_pd(PolyP, z);
 264
 265     z         = _mm256_mul_pd(PolyP, gmx_mm256_inv_pd(_mm256_sub_pd(PolyQ, PolyP)));
 266     z         = _mm256_add_pd(one, _mm256_mul_pd(two, z));
 267
 268     z         = _mm256_mul_pd(z, fexppart);
 269
 270     return z;
 271 }
 272
 273 /* 2^x, 128 bit */
 274 static __m128d
 275 gmx_mm_exp2_pd(__m128d x)
 276 {
 277     /* Lower bound: We do not allow numbers that would lead to an IEEE fp representation exponent smaller than -126. */
 278     const __m128d arglimit = _mm_set1_pd(1022.0);
 279     const __m128i expbase  = _mm_set1_epi32(1023);
 280
 281     const __m128d P2       = _mm_set1_pd(2.30933477057345225087e-2);
 282     const __m128d P1       = _mm_set1_pd(2.02020656693165307700e1);
 283     const __m128d P0       = _mm_set1_pd(1.51390680115615096133e3);
 284     /* Q2 == 1.0 */
 285     const __m128d Q1       = _mm_set1_pd(2.33184211722314911771e2);
 286     const __m128d Q0       = _mm_set1_pd(4.36821166879210612817e3);
 287     const __m128d one      = _mm_set1_pd(1.0);
 288     const __m128d two      = _mm_set1_pd(2.0);
 289
 290     __m128d       valuemask;
 291     __m128i       iexppart;
 292     __m128d       fexppart;
 293     __m128d       intpart;
 294     __m128d       z, z2;
 295     __m128d       PolyP, PolyQ;
 296
 297     iexppart  = _mm_cvtpd_epi32(x);
 298     intpart   = _mm_round_pd(x, _MM_FROUND_TO_NEAREST_INT);
 299
 300     /* The two lowest elements of iexppart now contains 32-bit numbers with a correctly biased exponent.
 301      * To be able to shift it into the exponent for a double precision number we first need to
 302      * shuffle so that the lower half contains the first element, and the upper half the second.
 303      * This should really be done as a zero-extension, but since the next instructions will shift
 304      * the registers left by 52 bits it doesn't matter what we put there - it will be shifted out.
 305      * (thus we just use element 2 from iexppart).
 306      */
 307     iexppart  = _mm_shuffle_epi32(iexppart, _MM_SHUFFLE(2, 1, 2, 0));
 308
 309     /* Do the shift operation on the 64-bit registers */
 310     iexppart  = _mm_add_epi32(iexppart, expbase);
 311     iexppart  = _mm_slli_epi64(iexppart, 52);
 312
 313     valuemask = _mm_cmpge_pd(arglimit, gmx_mm_abs_pd(x));
 314     fexppart  = _mm_and_pd(valuemask, gmx_mm_castsi128_pd(iexppart));
 315
 316     z         = _mm_sub_pd(x, intpart);
 317     z2        = _mm_mul_pd(z, z);
 318
 319     PolyP     = _mm_mul_pd(P2, z2);
 320     PolyP     = _mm_add_pd(PolyP, P1);
 321     PolyQ     = _mm_add_pd(z2, Q1);
 322     PolyP     = _mm_mul_pd(PolyP, z2);
 323     PolyQ     = _mm_mul_pd(PolyQ, z2);
 324     PolyP     = _mm_add_pd(PolyP, P0);
 325     PolyQ     = _mm_add_pd(PolyQ, Q0);
 326     PolyP     = _mm_mul_pd(PolyP, z);
 327
 328     z         = _mm_mul_pd(PolyP, gmx_mm_inv_pd(_mm_sub_pd(PolyQ, PolyP)));
 329     z         = _mm_add_pd(one, _mm_mul_pd(two, z));
 330
 331     z         = _mm_mul_pd(z, fexppart);
 332
 333     return z;
 334 }
 335
 336
 337 /* Exponential function, 256 bit. This could be calculated from 2^x as Exp(x)=2^(y),
 338  * where y=log2(e)*x, but there will then be a small rounding error since we lose
 339  * some precision due to the multiplication. This will then be magnified a lot by
 340  * the exponential.
 341  *
 342  * Instead, we calculate the fractional part directly as a Padé approximation of
 343  * Exp(z) on [-0.5,0.5]. We use extended precision arithmetics to calculate the fraction
 344  * remaining after 2^y, which avoids the precision-loss.
 345  */
 346 static __m256d
 347 gmx_mm256_exp_pd(__m256d exparg)
 348 {
 349     const __m256d argscale = _mm256_set1_pd(1.4426950408889634073599);
 350     /* Lower bound: We do not allow numbers that would lead to an IEEE fp representation exponent smaller than -126. */
 351     const __m256d arglimit = _mm256_set1_pd(1022.0);
 352     const __m128i expbase  = _mm_set1_epi32(1023);
 353
 354     const __m256d invargscale0  = _mm256_set1_pd(6.93145751953125e-1);
 355     const __m256d invargscale1  = _mm256_set1_pd(1.42860682030941723212e-6);
 356
 357     const __m256d P2       = _mm256_set1_pd(1.26177193074810590878e-4);
 358     const __m256d P1       = _mm256_set1_pd(3.02994407707441961300e-2);
 359     /* P0 == 1.0 */
 360     const __m256d Q3       = _mm256_set1_pd(3.00198505138664455042E-6);
 361     const __m256d Q2       = _mm256_set1_pd(2.52448340349684104192E-3);
 362     const __m256d Q1       = _mm256_set1_pd(2.27265548208155028766E-1);
 363     /* Q0 == 2.0 */
 364     const __m256d one      = _mm256_set1_pd(1.0);
 365     const __m256d two      = _mm256_set1_pd(2.0);
 366
 367     __m256d       valuemask;
 368     __m256i       iexppart;
 369     __m128i       iexppart128a, iexppart128b;
 370     __m256d       fexppart;
 371     __m256d       intpart;
 372     __m256d       x, z, z2;
 373     __m256d       PolyP, PolyQ;
 374
 375     x             = _mm256_mul_pd(exparg, argscale);
 376
 377     iexppart128a  = _mm256_cvtpd_epi32(x);
 378     intpart       = _mm256_round_pd(x, _MM_FROUND_TO_NEAREST_INT);
 379
 380     /* Add exponent bias */
 381     iexppart128a   = _mm_add_epi32(iexppart128a, expbase);
 382
 383     /* We now want to shift the exponent 52 positions left, but to achieve this we need
 384      * to separate the 128-bit register data into two registers (4x64-bit > 128bit)
 385      * shift them, and then merge into a single __m256d.
 386      * Elements 0/1 should end up in iexppart128a, and 2/3 in iexppart128b.
 387      * It doesnt matter what we put in the 2nd/4th position, since that data will be
 388      * shifted out and replaced with zeros.
 389      */
 390     iexppart128b   = _mm_shuffle_epi32(iexppart128a, _MM_SHUFFLE(3, 3, 2, 2));
 391     iexppart128a   = _mm_shuffle_epi32(iexppart128a, _MM_SHUFFLE(1, 1, 0, 0));
 392
 393     iexppart128b   = _mm_slli_epi64(iexppart128b, 52);
 394     iexppart128a   = _mm_slli_epi64(iexppart128a, 52);
 395
 396     iexppart  = _mm256_castsi128_si256(iexppart128a);
 397     iexppart  = _mm256_insertf128_si256(iexppart, iexppart128b, 0x1);
 398
 399     valuemask = _mm256_cmp_pd(arglimit, gmx_mm256_abs_pd(x), _CMP_GE_OQ);
 400     fexppart  = _mm256_and_pd(valuemask, _mm256_castsi256_pd(iexppart));
 401
 402     z         = _mm256_sub_pd(exparg, _mm256_mul_pd(invargscale0, intpart));
 403     z         = _mm256_sub_pd(z, _mm256_mul_pd(invargscale1, intpart));
 404
 405     z2        = _mm256_mul_pd(z, z);
 406
 407     PolyQ     = _mm256_mul_pd(Q3, z2);
 408     PolyQ     = _mm256_add_pd(PolyQ, Q2);
 409     PolyP     = _mm256_mul_pd(P2, z2);
 410     PolyQ     = _mm256_mul_pd(PolyQ, z2);
 411     PolyP     = _mm256_add_pd(PolyP, P1);
 412     PolyQ     = _mm256_add_pd(PolyQ, Q1);
 413     PolyP     = _mm256_mul_pd(PolyP, z2);
 414     PolyQ     = _mm256_mul_pd(PolyQ, z2);
 415     PolyP     = _mm256_add_pd(PolyP, one);
 416     PolyQ     = _mm256_add_pd(PolyQ, two);
 417
 418     PolyP     = _mm256_mul_pd(PolyP, z);
 419
 420     z         = _mm256_mul_pd(PolyP, gmx_mm256_inv_pd(_mm256_sub_pd(PolyQ, PolyP)));
 421     z         = _mm256_add_pd(one, _mm256_mul_pd(two, z));
 422
 423     z         = _mm256_mul_pd(z, fexppart);
 424
 425     return z;
 426 }
 427
 428 /* exp(), 128 bit */
 429 static __m128d
 430 gmx_mm_exp_pd(__m128d exparg)
 431 {
 432     const __m128d argscale = _mm_set1_pd(1.4426950408889634073599);
 433     /* Lower bound: We do not allow numbers that would lead to an IEEE fp representation exponent smaller than -126. */
 434     const __m128d arglimit = _mm_set1_pd(1022.0);
 435     const __m128i expbase  = _mm_set1_epi32(1023);
 436
 437     const __m128d invargscale0  = _mm_set1_pd(6.93145751953125e-1);
 438     const __m128d invargscale1  = _mm_set1_pd(1.42860682030941723212e-6);
 439
 440     const __m128d P2       = _mm_set1_pd(1.26177193074810590878e-4);
 441     const __m128d P1       = _mm_set1_pd(3.02994407707441961300e-2);
 442     /* P0 == 1.0 */
 443     const __m128d Q3       = _mm_set1_pd(3.00198505138664455042E-6);
 444     const __m128d Q2       = _mm_set1_pd(2.52448340349684104192E-3);
 445     const __m128d Q1       = _mm_set1_pd(2.27265548208155028766E-1);
 446     /* Q0 == 2.0 */
 447     const __m128d one      = _mm_set1_pd(1.0);
 448     const __m128d two      = _mm_set1_pd(2.0);
 449
 450     __m128d       valuemask;
 451     __m128i       iexppart;
 452     __m128d       fexppart;
 453     __m128d       intpart;
 454     __m128d       x, z, z2;
 455     __m128d       PolyP, PolyQ;
 456
 457     x             = _mm_mul_pd(exparg, argscale);
 458
 459     iexppart  = _mm_cvtpd_epi32(x);
 460     intpart   = _mm_round_pd(x, _MM_FROUND_TO_NEAREST_INT);
 461
 462     /* The two lowest elements of iexppart now contains 32-bit numbers with a correctly biased exponent.
 463      * To be able to shift it into the exponent for a double precision number we first need to
 464      * shuffle so that the lower half contains the first element, and the upper half the second.
 465      * This should really be done as a zero-extension, but since the next instructions will shift
 466      * the registers left by 52 bits it doesn't matter what we put there - it will be shifted out.
 467      * (thus we just use element 2 from iexppart).
 468      */
 469     iexppart  = _mm_shuffle_epi32(iexppart, _MM_SHUFFLE(2, 1, 2, 0));
 470
 471     /* Do the shift operation on the 64-bit registers */
 472     iexppart  = _mm_add_epi32(iexppart, expbase);
 473     iexppart  = _mm_slli_epi64(iexppart, 52);
 474
 475     valuemask = _mm_cmpge_pd(arglimit, gmx_mm_abs_pd(x));
 476     fexppart  = _mm_and_pd(valuemask, gmx_mm_castsi128_pd(iexppart));
 477
 478     z         = _mm_sub_pd(exparg, _mm_mul_pd(invargscale0, intpart));
 479     z         = _mm_sub_pd(z, _mm_mul_pd(invargscale1, intpart));
 480
 481     z2        = _mm_mul_pd(z, z);
 482
 483     PolyQ     = _mm_mul_pd(Q3, z2);
 484     PolyQ     = _mm_add_pd(PolyQ, Q2);
 485     PolyP     = _mm_mul_pd(P2, z2);
 486     PolyQ     = _mm_mul_pd(PolyQ, z2);
 487     PolyP     = _mm_add_pd(PolyP, P1);
 488     PolyQ     = _mm_add_pd(PolyQ, Q1);
 489     PolyP     = _mm_mul_pd(PolyP, z2);
 490     PolyQ     = _mm_mul_pd(PolyQ, z2);
 491     PolyP     = _mm_add_pd(PolyP, one);
 492     PolyQ     = _mm_add_pd(PolyQ, two);
 493
 494     PolyP     = _mm_mul_pd(PolyP, z);
 495
 496     z         = _mm_mul_pd(PolyP, gmx_mm_inv_pd(_mm_sub_pd(PolyQ, PolyP)));
 497     z         = _mm_add_pd(one, _mm_mul_pd(two, z));
 498
 499     z         = _mm_mul_pd(z, fexppart);
 500
 501     return z;
 502 }
 503
 504
 505 static __m256d
 506 gmx_mm256_log_pd(__m256d x)
 507 {
 508     /* Same algorithm as cephes library */
 509     const __m256d expmask    = _mm256_castsi256_pd( _mm256_set_epi32(0x7FF00000, 0x00000000, 0x7FF00000, 0x00000000,
 510                                                                      0x7FF00000, 0x00000000, 0x7FF00000, 0x00000000) );
 511
 512     const __m128i expbase_m1 = _mm_set1_epi32(1023-1); /* We want non-IEEE format */
 513
 514     const __m256d half       = _mm256_set1_pd(0.5);
 515     const __m256d one        = _mm256_set1_pd(1.0);
 516     const __m256d two        = _mm256_set1_pd(2.0);
 517     const __m256d invsq2     = _mm256_set1_pd(1.0/sqrt(2.0));
 518
 519     const __m256d corr1      = _mm256_set1_pd(-2.121944400546905827679e-4);
 520     const __m256d corr2      = _mm256_set1_pd(0.693359375);
 521
 522     const __m256d P5         = _mm256_set1_pd(1.01875663804580931796e-4);
 523     const __m256d P4         = _mm256_set1_pd(4.97494994976747001425e-1);
 524     const __m256d P3         = _mm256_set1_pd(4.70579119878881725854e0);
 525     const __m256d P2         = _mm256_set1_pd(1.44989225341610930846e1);
 526     const __m256d P1         = _mm256_set1_pd(1.79368678507819816313e1);
 527     const __m256d P0         = _mm256_set1_pd(7.70838733755885391666e0);
 528
 529     const __m256d Q4         = _mm256_set1_pd(1.12873587189167450590e1);
 530     const __m256d Q3         = _mm256_set1_pd(4.52279145837532221105e1);
 531     const __m256d Q2         = _mm256_set1_pd(8.29875266912776603211e1);
 532     const __m256d Q1         = _mm256_set1_pd(7.11544750618563894466e1);
 533     const __m256d Q0         = _mm256_set1_pd(2.31251620126765340583e1);
 534
 535     const __m256d R2         = _mm256_set1_pd(-7.89580278884799154124e-1);
 536     const __m256d R1         = _mm256_set1_pd(1.63866645699558079767e1);
 537     const __m256d R0         = _mm256_set1_pd(-6.41409952958715622951e1);
 538
 539     const __m256d S2         = _mm256_set1_pd(-3.56722798256324312549E1);
 540     const __m256d S1         = _mm256_set1_pd(3.12093766372244180303E2);
 541     const __m256d S0         = _mm256_set1_pd(-7.69691943550460008604E2);
 542
 543     __m256d       fexp;
 544     __m256i       iexp;
 545     __m128i       iexp128a, iexp128b;
 546
 547     __m256d       mask1, mask2;
 548     __m256d       corr, t1, t2, q;
 549     __m256d       zA, yA, xA, zB, yB, xB, z;
 550     __m256d       polyR, polyS;
 551     __m256d       polyP1, polyP2, polyQ1, polyQ2;
 552
 553     /* Separate x into exponent and mantissa, with a mantissa in the range [0.5..1[ (not IEEE754 standard!) */
 554     fexp     = _mm256_and_pd(x, expmask);
 555
 556     iexp     = _mm256_castpd_si256(fexp);
 557     iexp128b = _mm256_extractf128_si256(iexp, 0x1);
 558     iexp128a = _mm256_castsi256_si128(iexp);
 559
 560     iexp128a  = _mm_srli_epi64(iexp128a, 52);
 561     iexp128b  = _mm_srli_epi64(iexp128b, 52);
 562     /* Merge into a single register */
 563     iexp128a  = _mm_shuffle_epi32(iexp128a, _MM_SHUFFLE(1, 1, 2, 0));
 564     iexp128b  = _mm_shuffle_epi32(iexp128b, _MM_SHUFFLE(2, 0, 1, 1));
 565     iexp128a  = _mm_or_si128(iexp128a, iexp128b);
 566     iexp128a  = _mm_sub_epi32(iexp128a, expbase_m1);
 567
 568     fexp      = _mm256_cvtepi32_pd(iexp128a);
 569
 570     x         = _mm256_andnot_pd(expmask, x); /* Get mantissa */
 571     x         = _mm256_or_pd(x, one);
 572     x         = _mm256_mul_pd(x, half);
 573
 574     mask1     = _mm256_cmp_pd(gmx_mm256_abs_pd(fexp), two, _CMP_GT_OQ);
 575     mask2     = _mm256_cmp_pd(x, invsq2, _CMP_LT_OQ);
 576
 577     fexp      = _mm256_sub_pd(fexp, _mm256_and_pd(mask2, one));
 578
 579     /* If mask1 is set ('A') */
 580     zA     = _mm256_sub_pd(x, half);
 581     t1     = _mm256_blendv_pd( zA, x, mask2 );
 582     zA     = _mm256_sub_pd(t1, half);
 583     t2     = _mm256_blendv_pd( x, zA, mask2 );
 584     yA     = _mm256_mul_pd(half, _mm256_add_pd(t2, one));
 585
 586     xA     = _mm256_mul_pd(zA, gmx_mm256_inv_pd(yA));
 587     zA     = _mm256_mul_pd(xA, xA);
 588
 589     /* EVALUATE POLY */
 590     polyR  = _mm256_mul_pd(R2, zA);
 591     polyR  = _mm256_add_pd(polyR, R1);
 592     polyR  = _mm256_mul_pd(polyR, zA);
 593     polyR  = _mm256_add_pd(polyR, R0);
 594
 595     polyS  = _mm256_add_pd(zA, S2);
 596     polyS  = _mm256_mul_pd(polyS, zA);
 597     polyS  = _mm256_add_pd(polyS, S1);
 598     polyS  = _mm256_mul_pd(polyS, zA);
 599     polyS  = _mm256_add_pd(polyS, S0);
 600
 601     q      = _mm256_mul_pd(polyR, gmx_mm256_inv_pd(polyS));
 602     zA     = _mm256_mul_pd(_mm256_mul_pd(xA, zA), q);
 603
 604     zA     = _mm256_add_pd(zA, _mm256_mul_pd(corr1, fexp));
 605     zA     = _mm256_add_pd(zA, xA);
 606     zA     = _mm256_add_pd(zA, _mm256_mul_pd(corr2, fexp));
 607
 608     /* If mask1 is not set ('B') */
 609     corr   = _mm256_and_pd(mask2, x);
 610     xB     = _mm256_add_pd(x, corr);
 611     xB     = _mm256_sub_pd(xB, one);
 612     zB     = _mm256_mul_pd(xB, xB);
 613
 614     polyP1 = _mm256_mul_pd(P5, zB);
 615     polyP2 = _mm256_mul_pd(P4, zB);
 616     polyP1 = _mm256_add_pd(polyP1, P3);
 617     polyP2 = _mm256_add_pd(polyP2, P2);
 618     polyP1 = _mm256_mul_pd(polyP1, zB);
 619     polyP2 = _mm256_mul_pd(polyP2, zB);
 620     polyP1 = _mm256_add_pd(polyP1, P1);
 621     polyP2 = _mm256_add_pd(polyP2, P0);
 622     polyP1 = _mm256_mul_pd(polyP1, xB);
 623     polyP1 = _mm256_add_pd(polyP1, polyP2);
 624
 625     polyQ2 = _mm256_mul_pd(Q4, zB);
 626     polyQ1 = _mm256_add_pd(zB, Q3);
 627     polyQ2 = _mm256_add_pd(polyQ2, Q2);
 628     polyQ1 = _mm256_mul_pd(polyQ1, zB);
 629     polyQ2 = _mm256_mul_pd(polyQ2, zB);
 630     polyQ1 = _mm256_add_pd(polyQ1, Q1);
 631     polyQ2 = _mm256_add_pd(polyQ2, Q0);
 632     polyQ1 = _mm256_mul_pd(polyQ1, xB);
 633     polyQ1 = _mm256_add_pd(polyQ1, polyQ2);
 634
 635     fexp   = _mm256_and_pd(fexp, _mm256_cmp_pd(fexp, _mm256_setzero_pd(), _CMP_NEQ_OQ));
 636
 637     q      = _mm256_mul_pd(polyP1, gmx_mm256_inv_pd(polyQ1));
 638     yB     = _mm256_mul_pd(_mm256_mul_pd(xB, zB), q);
 639
 640     yB     = _mm256_add_pd(yB, _mm256_mul_pd(corr1, fexp));
 641     yB     = _mm256_sub_pd(yB, _mm256_mul_pd(half, zB));
 642     zB     = _mm256_add_pd(xB, yB);
 643     zB     = _mm256_add_pd(zB, _mm256_mul_pd(corr2, fexp));
 644
 645     z      = _mm256_blendv_pd( zB, zA, mask1 );
 646
 647     return z;
 648 }
 649
 650 static __m128d
 651 gmx_mm_log_pd(__m128d x)
 652 {
 653     /* Same algorithm as cephes library */
 654     const __m128d expmask    = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FF00000, 0x00000000, 0x7FF00000, 0x00000000) );
 655
 656     const __m128i expbase_m1 = _mm_set1_epi32(1023-1); /* We want non-IEEE format */
 657
 658     const __m128d half       = _mm_set1_pd(0.5);
 659     const __m128d one        = _mm_set1_pd(1.0);
 660     const __m128d two        = _mm_set1_pd(2.0);
 661     const __m128d invsq2     = _mm_set1_pd(1.0/sqrt(2.0));
 662
 663     const __m128d corr1      = _mm_set1_pd(-2.121944400546905827679e-4);
 664     const __m128d corr2      = _mm_set1_pd(0.693359375);
 665
 666     const __m128d P5         = _mm_set1_pd(1.01875663804580931796e-4);
 667     const __m128d P4         = _mm_set1_pd(4.97494994976747001425e-1);
 668     const __m128d P3         = _mm_set1_pd(4.70579119878881725854e0);
 669     const __m128d P2         = _mm_set1_pd(1.44989225341610930846e1);
 670     const __m128d P1         = _mm_set1_pd(1.79368678507819816313e1);
 671     const __m128d P0         = _mm_set1_pd(7.70838733755885391666e0);
 672
 673     const __m128d Q4         = _mm_set1_pd(1.12873587189167450590e1);
 674     const __m128d Q3         = _mm_set1_pd(4.52279145837532221105e1);
 675     const __m128d Q2         = _mm_set1_pd(8.29875266912776603211e1);
 676     const __m128d Q1         = _mm_set1_pd(7.11544750618563894466e1);
 677     const __m128d Q0         = _mm_set1_pd(2.31251620126765340583e1);
 678
 679     const __m128d R2         = _mm_set1_pd(-7.89580278884799154124e-1);
 680     const __m128d R1         = _mm_set1_pd(1.63866645699558079767e1);
 681     const __m128d R0         = _mm_set1_pd(-6.41409952958715622951e1);
 682
 683     const __m128d S2         = _mm_set1_pd(-3.56722798256324312549E1);
 684     const __m128d S1         = _mm_set1_pd(3.12093766372244180303E2);
 685     const __m128d S0         = _mm_set1_pd(-7.69691943550460008604E2);
 686
 687     __m128d       fexp;
 688     __m128i       iexp;
 689
 690     __m128d       mask1, mask2;
 691     __m128d       corr, t1, t2, q;
 692     __m128d       zA, yA, xA, zB, yB, xB, z;
 693     __m128d       polyR, polyS;
 694     __m128d       polyP1, polyP2, polyQ1, polyQ2;
 695
 696     /* Separate x into exponent and mantissa, with a mantissa in the range [0.5..1[ (not IEEE754 standard!) */
 697     fexp   = _mm_and_pd(x, expmask);
 698     iexp   = gmx_mm_castpd_si128(fexp);
 699     iexp   = _mm_srli_epi64(iexp, 52);
 700     iexp   = _mm_sub_epi32(iexp, expbase_m1);
 701     iexp   = _mm_shuffle_epi32(iexp, _MM_SHUFFLE(1, 1, 2, 0) );
 702     fexp   = _mm_cvtepi32_pd(iexp);
 703
 704     x      = _mm_andnot_pd(expmask, x);
 705     x      = _mm_or_pd(x, one);
 706     x      = _mm_mul_pd(x, half);
 707
 708     mask1     = _mm_cmpgt_pd(gmx_mm_abs_pd(fexp), two);
 709     mask2     = _mm_cmplt_pd(x, invsq2);
 710
 711     fexp   = _mm_sub_pd(fexp, _mm_and_pd(mask2, one));
 712
 713     /* If mask1 is set ('A') */
 714     zA     = _mm_sub_pd(x, half);
 715     t1     = _mm_blendv_pd( zA, x, mask2 );
 716     zA     = _mm_sub_pd(t1, half);
 717     t2     = _mm_blendv_pd( x, zA, mask2 );
 718     yA     = _mm_mul_pd(half, _mm_add_pd(t2, one));
 719
 720     xA     = _mm_mul_pd(zA, gmx_mm_inv_pd(yA));
 721     zA     = _mm_mul_pd(xA, xA);
 722
 723     /* EVALUATE POLY */
 724     polyR  = _mm_mul_pd(R2, zA);
 725     polyR  = _mm_add_pd(polyR, R1);
 726     polyR  = _mm_mul_pd(polyR, zA);
 727     polyR  = _mm_add_pd(polyR, R0);
 728
 729     polyS  = _mm_add_pd(zA, S2);
 730     polyS  = _mm_mul_pd(polyS, zA);
 731     polyS  = _mm_add_pd(polyS, S1);
 732     polyS  = _mm_mul_pd(polyS, zA);
 733     polyS  = _mm_add_pd(polyS, S0);
 734
 735     q      = _mm_mul_pd(polyR, gmx_mm_inv_pd(polyS));
 736     zA     = _mm_mul_pd(_mm_mul_pd(xA, zA), q);
 737
 738     zA     = _mm_add_pd(zA, _mm_mul_pd(corr1, fexp));
 739     zA     = _mm_add_pd(zA, xA);
 740     zA     = _mm_add_pd(zA, _mm_mul_pd(corr2, fexp));
 741
 742     /* If mask1 is not set ('B') */
 743     corr   = _mm_and_pd(mask2, x);
 744     xB     = _mm_add_pd(x, corr);
 745     xB     = _mm_sub_pd(xB, one);
 746     zB     = _mm_mul_pd(xB, xB);
 747
 748     polyP1 = _mm_mul_pd(P5, zB);
 749     polyP2 = _mm_mul_pd(P4, zB);
 750     polyP1 = _mm_add_pd(polyP1, P3);
 751     polyP2 = _mm_add_pd(polyP2, P2);
 752     polyP1 = _mm_mul_pd(polyP1, zB);
 753     polyP2 = _mm_mul_pd(polyP2, zB);
 754     polyP1 = _mm_add_pd(polyP1, P1);
 755     polyP2 = _mm_add_pd(polyP2, P0);
 756     polyP1 = _mm_mul_pd(polyP1, xB);
 757     polyP1 = _mm_add_pd(polyP1, polyP2);
 758
 759     polyQ2 = _mm_mul_pd(Q4, zB);
 760     polyQ1 = _mm_add_pd(zB, Q3);
 761     polyQ2 = _mm_add_pd(polyQ2, Q2);
 762     polyQ1 = _mm_mul_pd(polyQ1, zB);
 763     polyQ2 = _mm_mul_pd(polyQ2, zB);
 764     polyQ1 = _mm_add_pd(polyQ1, Q1);
 765     polyQ2 = _mm_add_pd(polyQ2, Q0);
 766     polyQ1 = _mm_mul_pd(polyQ1, xB);
 767     polyQ1 = _mm_add_pd(polyQ1, polyQ2);
 768
 769     fexp   = _mm_and_pd(fexp, _mm_cmpneq_pd(fexp, _mm_setzero_pd()));
 770
 771     q      = _mm_mul_pd(polyP1, gmx_mm_inv_pd(polyQ1));
 772     yB     = _mm_mul_pd(_mm_mul_pd(xB, zB), q);
 773
 774     yB     = _mm_add_pd(yB, _mm_mul_pd(corr1, fexp));
 775     yB     = _mm_sub_pd(yB, _mm_mul_pd(half, zB));
 776     zB     = _mm_add_pd(xB, yB);
 777     zB     = _mm_add_pd(zB, _mm_mul_pd(corr2, fexp));
 778
 779     z      = _mm_blendv_pd( zB, zA, mask1 );
 780
 781     return z;
 782 }
 783
 784
 785 static __m256d
 786 gmx_mm256_erf_pd(__m256d x)
 787 {
 788     /* Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75] */
 789     const __m256d CAP4      = _mm256_set1_pd(-0.431780540597889301512e-4);
 790     const __m256d CAP3      = _mm256_set1_pd(-0.00578562306260059236059);
 791     const __m256d CAP2      = _mm256_set1_pd(-0.028593586920219752446);
 792     const __m256d CAP1      = _mm256_set1_pd(-0.315924962948621698209);
 793     const __m256d CAP0      = _mm256_set1_pd(0.14952975608477029151);
 794
 795     const __m256d CAQ5      = _mm256_set1_pd(-0.374089300177174709737e-5);
 796     const __m256d CAQ4      = _mm256_set1_pd(0.00015126584532155383535);
 797     const __m256d CAQ3      = _mm256_set1_pd(0.00536692680669480725423);
 798     const __m256d CAQ2      = _mm256_set1_pd(0.0668686825594046122636);
 799     const __m256d CAQ1      = _mm256_set1_pd(0.402604990869284362773);
 800     /* CAQ0 == 1.0 */
 801     const __m256d CAoffset  = _mm256_set1_pd(0.9788494110107421875);
 802
 803     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5] */
 804     const __m256d CBP6      = _mm256_set1_pd(2.49650423685462752497647637088e-10);
 805     const __m256d CBP5      = _mm256_set1_pd(0.00119770193298159629350136085658);
 806     const __m256d CBP4      = _mm256_set1_pd(0.0164944422378370965881008942733);
 807     const __m256d CBP3      = _mm256_set1_pd(0.0984581468691775932063932439252);
 808     const __m256d CBP2      = _mm256_set1_pd(0.317364595806937763843589437418);
 809     const __m256d CBP1      = _mm256_set1_pd(0.554167062641455850932670067075);
 810     const __m256d CBP0      = _mm256_set1_pd(0.427583576155807163756925301060);
 811     const __m256d CBQ7      = _mm256_set1_pd(0.00212288829699830145976198384930);
 812     const __m256d CBQ6      = _mm256_set1_pd(0.0334810979522685300554606393425);
 813     const __m256d CBQ5      = _mm256_set1_pd(0.2361713785181450957579508850717);
 814     const __m256d CBQ4      = _mm256_set1_pd(0.955364736493055670530981883072);
 815     const __m256d CBQ3      = _mm256_set1_pd(2.36815675631420037315349279199);
 816     const __m256d CBQ2      = _mm256_set1_pd(3.55261649184083035537184223542);
 817     const __m256d CBQ1      = _mm256_set1_pd(2.93501136050160872574376997993);
 818     /* CBQ0 == 1.0 */
 819
 820     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf] */
 821     const __m256d CCP6      = _mm256_set1_pd(-2.8175401114513378771);
 822     const __m256d CCP5      = _mm256_set1_pd(-3.22729451764143718517);
 823     const __m256d CCP4      = _mm256_set1_pd(-2.5518551727311523996);
 824     const __m256d CCP3      = _mm256_set1_pd(-0.687717681153649930619);
 825     const __m256d CCP2      = _mm256_set1_pd(-0.212652252872804219852);
 826     const __m256d CCP1      = _mm256_set1_pd(0.0175389834052493308818);
 827     const __m256d CCP0      = _mm256_set1_pd(0.00628057170626964891937);
 828
 829     const __m256d CCQ6      = _mm256_set1_pd(5.48409182238641741584);
 830     const __m256d CCQ5      = _mm256_set1_pd(13.5064170191802889145);
 831     const __m256d CCQ4      = _mm256_set1_pd(22.9367376522880577224);
 832     const __m256d CCQ3      = _mm256_set1_pd(15.930646027911794143);
 833     const __m256d CCQ2      = _mm256_set1_pd(11.0567237927800161565);
 834     const __m256d CCQ1      = _mm256_set1_pd(2.79257750980575282228);
 835     /* CCQ0 == 1.0 */
 836     const __m256d CCoffset  = _mm256_set1_pd(0.5579090118408203125);
 837
 838     const __m256d one       = _mm256_set1_pd(1.0);
 839     const __m256d two       = _mm256_set1_pd(2.0);
 840
 841     const __m256d signbit   = _mm256_castsi256_pd( _mm256_set_epi32(0x80000000, 0x00000000, 0x80000000, 0x00000000,
 842                                                                     0x80000000, 0x00000000, 0x80000000, 0x00000000) );
 843
 844     __m256d xabs, x2, x4, t, t2, w, w2;
 845     __m256d PolyAP0, PolyAP1, PolyAQ0, PolyAQ1;
 846     __m256d PolyBP0, PolyBP1, PolyBQ0, PolyBQ1;
 847     __m256d PolyCP0, PolyCP1, PolyCQ0, PolyCQ1;
 848     __m256d res_erf, res_erfcB, res_erfcC, res_erfc, res;
 849     __m256d mask, expmx2;
 850
 851     /* Calculate erf() */
 852     xabs     = gmx_mm256_abs_pd(x);
 853     x2       = _mm256_mul_pd(x, x);
 854     x4       = _mm256_mul_pd(x2, x2);
 855
 856     PolyAP0  = _mm256_mul_pd(CAP4, x4);
 857     PolyAP1  = _mm256_mul_pd(CAP3, x4);
 858     PolyAP0  = _mm256_add_pd(PolyAP0, CAP2);
 859     PolyAP1  = _mm256_add_pd(PolyAP1, CAP1);
 860     PolyAP0  = _mm256_mul_pd(PolyAP0, x4);
 861     PolyAP1  = _mm256_mul_pd(PolyAP1, x2);
 862     PolyAP0  = _mm256_add_pd(PolyAP0, CAP0);
 863     PolyAP0  = _mm256_add_pd(PolyAP0, PolyAP1);
 864
 865     PolyAQ1  = _mm256_mul_pd(CAQ5, x4);
 866     PolyAQ0  = _mm256_mul_pd(CAQ4, x4);
 867     PolyAQ1  = _mm256_add_pd(PolyAQ1, CAQ3);
 868     PolyAQ0  = _mm256_add_pd(PolyAQ0, CAQ2);
 869     PolyAQ1  = _mm256_mul_pd(PolyAQ1, x4);
 870     PolyAQ0  = _mm256_mul_pd(PolyAQ0, x4);
 871     PolyAQ1  = _mm256_add_pd(PolyAQ1, CAQ1);
 872     PolyAQ0  = _mm256_add_pd(PolyAQ0, one);
 873     PolyAQ1  = _mm256_mul_pd(PolyAQ1, x2);
 874     PolyAQ0  = _mm256_add_pd(PolyAQ0, PolyAQ1);
 875
 876     res_erf  = _mm256_mul_pd(PolyAP0, gmx_mm256_inv_pd(PolyAQ0));
 877     res_erf  = _mm256_add_pd(CAoffset, res_erf);
 878     res_erf  = _mm256_mul_pd(x, res_erf);
 879
 880     /* Calculate erfc() in range [1,4.5] */
 881     t       = _mm256_sub_pd(xabs, one);
 882     t2      = _mm256_mul_pd(t, t);
 883
 884     PolyBP0  = _mm256_mul_pd(CBP6, t2);
 885     PolyBP1  = _mm256_mul_pd(CBP5, t2);
 886     PolyBP0  = _mm256_add_pd(PolyBP0, CBP4);
 887     PolyBP1  = _mm256_add_pd(PolyBP1, CBP3);
 888     PolyBP0  = _mm256_mul_pd(PolyBP0, t2);
 889     PolyBP1  = _mm256_mul_pd(PolyBP1, t2);
 890     PolyBP0  = _mm256_add_pd(PolyBP0, CBP2);
 891     PolyBP1  = _mm256_add_pd(PolyBP1, CBP1);
 892     PolyBP0  = _mm256_mul_pd(PolyBP0, t2);
 893     PolyBP1  = _mm256_mul_pd(PolyBP1, t);
 894     PolyBP0  = _mm256_add_pd(PolyBP0, CBP0);
 895     PolyBP0  = _mm256_add_pd(PolyBP0, PolyBP1);
 896
 897     PolyBQ1 = _mm256_mul_pd(CBQ7, t2);
 898     PolyBQ0 = _mm256_mul_pd(CBQ6, t2);
 899     PolyBQ1 = _mm256_add_pd(PolyBQ1, CBQ5);
 900     PolyBQ0 = _mm256_add_pd(PolyBQ0, CBQ4);
 901     PolyBQ1 = _mm256_mul_pd(PolyBQ1, t2);
 902     PolyBQ0 = _mm256_mul_pd(PolyBQ0, t2);
 903     PolyBQ1 = _mm256_add_pd(PolyBQ1, CBQ3);
 904     PolyBQ0 = _mm256_add_pd(PolyBQ0, CBQ2);
 905     PolyBQ1 = _mm256_mul_pd(PolyBQ1, t2);
 906     PolyBQ0 = _mm256_mul_pd(PolyBQ0, t2);
 907     PolyBQ1 = _mm256_add_pd(PolyBQ1, CBQ1);
 908     PolyBQ0 = _mm256_add_pd(PolyBQ0, one);
 909     PolyBQ1 = _mm256_mul_pd(PolyBQ1, t);
 910     PolyBQ0 = _mm256_add_pd(PolyBQ0, PolyBQ1);
 911
 912     res_erfcB = _mm256_mul_pd(PolyBP0, gmx_mm256_inv_pd(PolyBQ0));
 913
 914     res_erfcB = _mm256_mul_pd(res_erfcB, xabs);
 915
 916     /* Calculate erfc() in range [4.5,inf] */
 917     w       = gmx_mm256_inv_pd(xabs);
 918     w2      = _mm256_mul_pd(w, w);
 919
 920     PolyCP0  = _mm256_mul_pd(CCP6, w2);
 921     PolyCP1  = _mm256_mul_pd(CCP5, w2);
 922     PolyCP0  = _mm256_add_pd(PolyCP0, CCP4);
 923     PolyCP1  = _mm256_add_pd(PolyCP1, CCP3);
 924     PolyCP0  = _mm256_mul_pd(PolyCP0, w2);
 925     PolyCP1  = _mm256_mul_pd(PolyCP1, w2);
 926     PolyCP0  = _mm256_add_pd(PolyCP0, CCP2);
 927     PolyCP1  = _mm256_add_pd(PolyCP1, CCP1);
 928     PolyCP0  = _mm256_mul_pd(PolyCP0, w2);
 929     PolyCP1  = _mm256_mul_pd(PolyCP1, w);
 930     PolyCP0  = _mm256_add_pd(PolyCP0, CCP0);
 931     PolyCP0  = _mm256_add_pd(PolyCP0, PolyCP1);
 932
 933     PolyCQ0  = _mm256_mul_pd(CCQ6, w2);
 934     PolyCQ1  = _mm256_mul_pd(CCQ5, w2);
 935     PolyCQ0  = _mm256_add_pd(PolyCQ0, CCQ4);
 936     PolyCQ1  = _mm256_add_pd(PolyCQ1, CCQ3);
 937     PolyCQ0  = _mm256_mul_pd(PolyCQ0, w2);
 938     PolyCQ1  = _mm256_mul_pd(PolyCQ1, w2);
 939     PolyCQ0  = _mm256_add_pd(PolyCQ0, CCQ2);
 940     PolyCQ1  = _mm256_add_pd(PolyCQ1, CCQ1);
 941     PolyCQ0  = _mm256_mul_pd(PolyCQ0, w2);
 942     PolyCQ1  = _mm256_mul_pd(PolyCQ1, w);
 943     PolyCQ0  = _mm256_add_pd(PolyCQ0, one);
 944     PolyCQ0  = _mm256_add_pd(PolyCQ0, PolyCQ1);
 945
 946     expmx2   = gmx_mm256_exp_pd( _mm256_or_pd(signbit, x2) );
 947
 948     res_erfcC = _mm256_mul_pd(PolyCP0, gmx_mm256_inv_pd(PolyCQ0));
 949     res_erfcC = _mm256_add_pd(res_erfcC, CCoffset);
 950     res_erfcC = _mm256_mul_pd(res_erfcC, w);
 951
 952     mask     = _mm256_cmp_pd(xabs, _mm256_set1_pd(4.5), _CMP_GT_OQ);
 953     res_erfc = _mm256_blendv_pd(res_erfcB, res_erfcC, mask);
 954
 955     res_erfc = _mm256_mul_pd(res_erfc, expmx2);
 956
 957     /* erfc(x<0) = 2-erfc(|x|) */
 958     mask     = _mm256_cmp_pd(x, _mm256_setzero_pd(), _CMP_LT_OQ);
 959     res_erfc = _mm256_blendv_pd(res_erfc, _mm256_sub_pd(two, res_erfc), mask);
 960
 961     /* Select erf() or erfc() */
 962     mask = _mm256_cmp_pd(xabs, one, _CMP_LT_OQ);
 963     res  = _mm256_blendv_pd(_mm256_sub_pd(one, res_erfc), res_erf, mask);
 964
 965     return res;
 966 }
 967
 968 static __m128d
 969 gmx_mm_erf_pd(__m128d x)
 970 {
 971     /* Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75] */
 972     const __m128d CAP4      = _mm_set1_pd(-0.431780540597889301512e-4);
 973     const __m128d CAP3      = _mm_set1_pd(-0.00578562306260059236059);
 974     const __m128d CAP2      = _mm_set1_pd(-0.028593586920219752446);
 975     const __m128d CAP1      = _mm_set1_pd(-0.315924962948621698209);
 976     const __m128d CAP0      = _mm_set1_pd(0.14952975608477029151);
 977
 978     const __m128d CAQ5      = _mm_set1_pd(-0.374089300177174709737e-5);
 979     const __m128d CAQ4      = _mm_set1_pd(0.00015126584532155383535);
 980     const __m128d CAQ3      = _mm_set1_pd(0.00536692680669480725423);
 981     const __m128d CAQ2      = _mm_set1_pd(0.0668686825594046122636);
 982     const __m128d CAQ1      = _mm_set1_pd(0.402604990869284362773);
 983     /* CAQ0 == 1.0 */
 984     const __m128d CAoffset  = _mm_set1_pd(0.9788494110107421875);
 985
 986     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5] */
 987     const __m128d CBP6      = _mm_set1_pd(2.49650423685462752497647637088e-10);
 988     const __m128d CBP5      = _mm_set1_pd(0.00119770193298159629350136085658);
 989     const __m128d CBP4      = _mm_set1_pd(0.0164944422378370965881008942733);
 990     const __m128d CBP3      = _mm_set1_pd(0.0984581468691775932063932439252);
 991     const __m128d CBP2      = _mm_set1_pd(0.317364595806937763843589437418);
 992     const __m128d CBP1      = _mm_set1_pd(0.554167062641455850932670067075);
 993     const __m128d CBP0      = _mm_set1_pd(0.427583576155807163756925301060);
 994     const __m128d CBQ7      = _mm_set1_pd(0.00212288829699830145976198384930);
 995     const __m128d CBQ6      = _mm_set1_pd(0.0334810979522685300554606393425);
 996     const __m128d CBQ5      = _mm_set1_pd(0.2361713785181450957579508850717);
 997     const __m128d CBQ4      = _mm_set1_pd(0.955364736493055670530981883072);
 998     const __m128d CBQ3      = _mm_set1_pd(2.36815675631420037315349279199);
 999     const __m128d CBQ2      = _mm_set1_pd(3.55261649184083035537184223542);
1000     const __m128d CBQ1      = _mm_set1_pd(2.93501136050160872574376997993);
1001     /* CBQ0 == 1.0 */
1002
1003     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf] */
1004     const __m128d CCP6      = _mm_set1_pd(-2.8175401114513378771);
1005     const __m128d CCP5      = _mm_set1_pd(-3.22729451764143718517);
1006     const __m128d CCP4      = _mm_set1_pd(-2.5518551727311523996);
1007     const __m128d CCP3      = _mm_set1_pd(-0.687717681153649930619);
1008     const __m128d CCP2      = _mm_set1_pd(-0.212652252872804219852);
1009     const __m128d CCP1      = _mm_set1_pd(0.0175389834052493308818);
1010     const __m128d CCP0      = _mm_set1_pd(0.00628057170626964891937);
1011
1012     const __m128d CCQ6      = _mm_set1_pd(5.48409182238641741584);
1013     const __m128d CCQ5      = _mm_set1_pd(13.5064170191802889145);
1014     const __m128d CCQ4      = _mm_set1_pd(22.9367376522880577224);
1015     const __m128d CCQ3      = _mm_set1_pd(15.930646027911794143);
1016     const __m128d CCQ2      = _mm_set1_pd(11.0567237927800161565);
1017     const __m128d CCQ1      = _mm_set1_pd(2.79257750980575282228);
1018     /* CCQ0 == 1.0 */
1019     const __m128d CCoffset  = _mm_set1_pd(0.5579090118408203125);
1020
1021     const __m128d one       = _mm_set1_pd(1.0);
1022     const __m128d two       = _mm_set1_pd(2.0);
1023
1024     const __m128d signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000, 0x00000000, 0x80000000, 0x00000000) );
1025
1026     __m128d       xabs, x2, x4, t, t2, w, w2;
1027     __m128d       PolyAP0, PolyAP1, PolyAQ0, PolyAQ1;
1028     __m128d       PolyBP0, PolyBP1, PolyBQ0, PolyBQ1;
1029     __m128d       PolyCP0, PolyCP1, PolyCQ0, PolyCQ1;
1030     __m128d       res_erf, res_erfcB, res_erfcC, res_erfc, res;
1031     __m128d       mask, expmx2;
1032
1033     /* Calculate erf() */
1034     xabs     = gmx_mm_abs_pd(x);
1035     x2       = _mm_mul_pd(x, x);
1036     x4       = _mm_mul_pd(x2, x2);
1037
1038     PolyAP0  = _mm_mul_pd(CAP4, x4);
1039     PolyAP1  = _mm_mul_pd(CAP3, x4);
1040     PolyAP0  = _mm_add_pd(PolyAP0, CAP2);
1041     PolyAP1  = _mm_add_pd(PolyAP1, CAP1);
1042     PolyAP0  = _mm_mul_pd(PolyAP0, x4);
1043     PolyAP1  = _mm_mul_pd(PolyAP1, x2);
1044     PolyAP0  = _mm_add_pd(PolyAP0, CAP0);
1045     PolyAP0  = _mm_add_pd(PolyAP0, PolyAP1);
1046
1047     PolyAQ1  = _mm_mul_pd(CAQ5, x4);
1048     PolyAQ0  = _mm_mul_pd(CAQ4, x4);
1049     PolyAQ1  = _mm_add_pd(PolyAQ1, CAQ3);
1050     PolyAQ0  = _mm_add_pd(PolyAQ0, CAQ2);
1051     PolyAQ1  = _mm_mul_pd(PolyAQ1, x4);
1052     PolyAQ0  = _mm_mul_pd(PolyAQ0, x4);
1053     PolyAQ1  = _mm_add_pd(PolyAQ1, CAQ1);
1054     PolyAQ0  = _mm_add_pd(PolyAQ0, one);
1055     PolyAQ1  = _mm_mul_pd(PolyAQ1, x2);
1056     PolyAQ0  = _mm_add_pd(PolyAQ0, PolyAQ1);
1057
1058     res_erf  = _mm_mul_pd(PolyAP0, gmx_mm_inv_pd(PolyAQ0));
1059     res_erf  = _mm_add_pd(CAoffset, res_erf);
1060     res_erf  = _mm_mul_pd(x, res_erf);
1061
1062     /* Calculate erfc() in range [1,4.5] */
1063     t       = _mm_sub_pd(xabs, one);
1064     t2      = _mm_mul_pd(t, t);
1065
1066     PolyBP0  = _mm_mul_pd(CBP6, t2);
1067     PolyBP1  = _mm_mul_pd(CBP5, t2);
1068     PolyBP0  = _mm_add_pd(PolyBP0, CBP4);
1069     PolyBP1  = _mm_add_pd(PolyBP1, CBP3);
1070     PolyBP0  = _mm_mul_pd(PolyBP0, t2);
1071     PolyBP1  = _mm_mul_pd(PolyBP1, t2);
1072     PolyBP0  = _mm_add_pd(PolyBP0, CBP2);
1073     PolyBP1  = _mm_add_pd(PolyBP1, CBP1);
1074     PolyBP0  = _mm_mul_pd(PolyBP0, t2);
1075     PolyBP1  = _mm_mul_pd(PolyBP1, t);
1076     PolyBP0  = _mm_add_pd(PolyBP0, CBP0);
1077     PolyBP0  = _mm_add_pd(PolyBP0, PolyBP1);
1078
1079     PolyBQ1 = _mm_mul_pd(CBQ7, t2);
1080     PolyBQ0 = _mm_mul_pd(CBQ6, t2);
1081     PolyBQ1 = _mm_add_pd(PolyBQ1, CBQ5);
1082     PolyBQ0 = _mm_add_pd(PolyBQ0, CBQ4);
1083     PolyBQ1 = _mm_mul_pd(PolyBQ1, t2);
1084     PolyBQ0 = _mm_mul_pd(PolyBQ0, t2);
1085     PolyBQ1 = _mm_add_pd(PolyBQ1, CBQ3);
1086     PolyBQ0 = _mm_add_pd(PolyBQ0, CBQ2);
1087     PolyBQ1 = _mm_mul_pd(PolyBQ1, t2);
1088     PolyBQ0 = _mm_mul_pd(PolyBQ0, t2);
1089     PolyBQ1 = _mm_add_pd(PolyBQ1, CBQ1);
1090     PolyBQ0 = _mm_add_pd(PolyBQ0, one);
1091     PolyBQ1 = _mm_mul_pd(PolyBQ1, t);
1092     PolyBQ0 = _mm_add_pd(PolyBQ0, PolyBQ1);
1093
1094     res_erfcB = _mm_mul_pd(PolyBP0, gmx_mm_inv_pd(PolyBQ0));
1095
1096     res_erfcB = _mm_mul_pd(res_erfcB, xabs);
1097
1098     /* Calculate erfc() in range [4.5,inf] */
1099     w       = gmx_mm_inv_pd(xabs);
1100     w2      = _mm_mul_pd(w, w);
1101
1102     PolyCP0  = _mm_mul_pd(CCP6, w2);
1103     PolyCP1  = _mm_mul_pd(CCP5, w2);
1104     PolyCP0  = _mm_add_pd(PolyCP0, CCP4);
1105     PolyCP1  = _mm_add_pd(PolyCP1, CCP3);
1106     PolyCP0  = _mm_mul_pd(PolyCP0, w2);
1107     PolyCP1  = _mm_mul_pd(PolyCP1, w2);
1108     PolyCP0  = _mm_add_pd(PolyCP0, CCP2);
1109     PolyCP1  = _mm_add_pd(PolyCP1, CCP1);
1110     PolyCP0  = _mm_mul_pd(PolyCP0, w2);
1111     PolyCP1  = _mm_mul_pd(PolyCP1, w);
1112     PolyCP0  = _mm_add_pd(PolyCP0, CCP0);
1113     PolyCP0  = _mm_add_pd(PolyCP0, PolyCP1);
1114
1115     PolyCQ0  = _mm_mul_pd(CCQ6, w2);
1116     PolyCQ1  = _mm_mul_pd(CCQ5, w2);
1117     PolyCQ0  = _mm_add_pd(PolyCQ0, CCQ4);
1118     PolyCQ1  = _mm_add_pd(PolyCQ1, CCQ3);
1119     PolyCQ0  = _mm_mul_pd(PolyCQ0, w2);
1120     PolyCQ1  = _mm_mul_pd(PolyCQ1, w2);
1121     PolyCQ0  = _mm_add_pd(PolyCQ0, CCQ2);
1122     PolyCQ1  = _mm_add_pd(PolyCQ1, CCQ1);
1123     PolyCQ0  = _mm_mul_pd(PolyCQ0, w2);
1124     PolyCQ1  = _mm_mul_pd(PolyCQ1, w);
1125     PolyCQ0  = _mm_add_pd(PolyCQ0, one);
1126     PolyCQ0  = _mm_add_pd(PolyCQ0, PolyCQ1);
1127
1128     expmx2   = gmx_mm_exp_pd( _mm_or_pd(signbit, x2) );
1129
1130     res_erfcC = _mm_mul_pd(PolyCP0, gmx_mm_inv_pd(PolyCQ0));
1131     res_erfcC = _mm_add_pd(res_erfcC, CCoffset);
1132     res_erfcC = _mm_mul_pd(res_erfcC, w);
1133
1134     mask     = _mm_cmpgt_pd(xabs, _mm_set1_pd(4.5));
1135     res_erfc = _mm_blendv_pd(res_erfcB, res_erfcC, mask);
1136
1137     res_erfc = _mm_mul_pd(res_erfc, expmx2);
1138
1139     /* erfc(x<0) = 2-erfc(|x|) */
1140     mask     = _mm_cmplt_pd(x, _mm_setzero_pd());
1141     res_erfc = _mm_blendv_pd(res_erfc, _mm_sub_pd(two, res_erfc), mask);
1142
1143     /* Select erf() or erfc() */
1144     mask = _mm_cmplt_pd(xabs, one);
1145     res  = _mm_blendv_pd(_mm_sub_pd(one, res_erfc), res_erf, mask);
1146
1147     return res;
1148 }
1149
1150
1151 static __m256d
1152 gmx_mm256_erfc_pd(__m256d x)
1153 {
1154     /* Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75] */
1155     const __m256d CAP4      = _mm256_set1_pd(-0.431780540597889301512e-4);
1156     const __m256d CAP3      = _mm256_set1_pd(-0.00578562306260059236059);
1157     const __m256d CAP2      = _mm256_set1_pd(-0.028593586920219752446);
1158     const __m256d CAP1      = _mm256_set1_pd(-0.315924962948621698209);
1159     const __m256d CAP0      = _mm256_set1_pd(0.14952975608477029151);
1160
1161     const __m256d CAQ5      = _mm256_set1_pd(-0.374089300177174709737e-5);
1162     const __m256d CAQ4      = _mm256_set1_pd(0.00015126584532155383535);
1163     const __m256d CAQ3      = _mm256_set1_pd(0.00536692680669480725423);
1164     const __m256d CAQ2      = _mm256_set1_pd(0.0668686825594046122636);
1165     const __m256d CAQ1      = _mm256_set1_pd(0.402604990869284362773);
1166     /* CAQ0 == 1.0 */
1167     const __m256d CAoffset  = _mm256_set1_pd(0.9788494110107421875);
1168
1169     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5] */
1170     const __m256d CBP6      = _mm256_set1_pd(2.49650423685462752497647637088e-10);
1171     const __m256d CBP5      = _mm256_set1_pd(0.00119770193298159629350136085658);
1172     const __m256d CBP4      = _mm256_set1_pd(0.0164944422378370965881008942733);
1173     const __m256d CBP3      = _mm256_set1_pd(0.0984581468691775932063932439252);
1174     const __m256d CBP2      = _mm256_set1_pd(0.317364595806937763843589437418);
1175     const __m256d CBP1      = _mm256_set1_pd(0.554167062641455850932670067075);
1176     const __m256d CBP0      = _mm256_set1_pd(0.427583576155807163756925301060);
1177     const __m256d CBQ7      = _mm256_set1_pd(0.00212288829699830145976198384930);
1178     const __m256d CBQ6      = _mm256_set1_pd(0.0334810979522685300554606393425);
1179     const __m256d CBQ5      = _mm256_set1_pd(0.2361713785181450957579508850717);
1180     const __m256d CBQ4      = _mm256_set1_pd(0.955364736493055670530981883072);
1181     const __m256d CBQ3      = _mm256_set1_pd(2.36815675631420037315349279199);
1182     const __m256d CBQ2      = _mm256_set1_pd(3.55261649184083035537184223542);
1183     const __m256d CBQ1      = _mm256_set1_pd(2.93501136050160872574376997993);
1184     /* CBQ0 == 1.0 */
1185
1186     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf] */
1187     const __m256d CCP6      = _mm256_set1_pd(-2.8175401114513378771);
1188     const __m256d CCP5      = _mm256_set1_pd(-3.22729451764143718517);
1189     const __m256d CCP4      = _mm256_set1_pd(-2.5518551727311523996);
1190     const __m256d CCP3      = _mm256_set1_pd(-0.687717681153649930619);
1191     const __m256d CCP2      = _mm256_set1_pd(-0.212652252872804219852);
1192     const __m256d CCP1      = _mm256_set1_pd(0.0175389834052493308818);
1193     const __m256d CCP0      = _mm256_set1_pd(0.00628057170626964891937);
1194
1195     const __m256d CCQ6      = _mm256_set1_pd(5.48409182238641741584);
1196     const __m256d CCQ5      = _mm256_set1_pd(13.5064170191802889145);
1197     const __m256d CCQ4      = _mm256_set1_pd(22.9367376522880577224);
1198     const __m256d CCQ3      = _mm256_set1_pd(15.930646027911794143);
1199     const __m256d CCQ2      = _mm256_set1_pd(11.0567237927800161565);
1200     const __m256d CCQ1      = _mm256_set1_pd(2.79257750980575282228);
1201     /* CCQ0 == 1.0 */
1202     const __m256d CCoffset  = _mm256_set1_pd(0.5579090118408203125);
1203
1204     const __m256d one       = _mm256_set1_pd(1.0);
1205     const __m256d two       = _mm256_set1_pd(2.0);
1206
1207     const __m256d signbit   = _mm256_castsi256_pd( _mm256_set_epi32(0x80000000, 0x00000000, 0x80000000, 0x00000000,
1208                                                                     0x80000000, 0x00000000, 0x80000000, 0x00000000) );
1209
1210     __m256d xabs, x2, x4, t, t2, w, w2;
1211     __m256d PolyAP0, PolyAP1, PolyAQ0, PolyAQ1;
1212     __m256d PolyBP0, PolyBP1, PolyBQ0, PolyBQ1;
1213     __m256d PolyCP0, PolyCP1, PolyCQ0, PolyCQ1;
1214     __m256d res_erf, res_erfcB, res_erfcC, res_erfc, res;
1215     __m256d mask, expmx2;
1216
1217     /* Calculate erf() */
1218     xabs     = gmx_mm256_abs_pd(x);
1219     x2       = _mm256_mul_pd(x, x);
1220     x4       = _mm256_mul_pd(x2, x2);
1221
1222     PolyAP0  = _mm256_mul_pd(CAP4, x4);
1223     PolyAP1  = _mm256_mul_pd(CAP3, x4);
1224     PolyAP0  = _mm256_add_pd(PolyAP0, CAP2);
1225     PolyAP1  = _mm256_add_pd(PolyAP1, CAP1);
1226     PolyAP0  = _mm256_mul_pd(PolyAP0, x4);
1227     PolyAP1  = _mm256_mul_pd(PolyAP1, x2);
1228     PolyAP0  = _mm256_add_pd(PolyAP0, CAP0);
1229     PolyAP0  = _mm256_add_pd(PolyAP0, PolyAP1);
1230
1231     PolyAQ1  = _mm256_mul_pd(CAQ5, x4);
1232     PolyAQ0  = _mm256_mul_pd(CAQ4, x4);
1233     PolyAQ1  = _mm256_add_pd(PolyAQ1, CAQ3);
1234     PolyAQ0  = _mm256_add_pd(PolyAQ0, CAQ2);
1235     PolyAQ1  = _mm256_mul_pd(PolyAQ1, x4);
1236     PolyAQ0  = _mm256_mul_pd(PolyAQ0, x4);
1237     PolyAQ1  = _mm256_add_pd(PolyAQ1, CAQ1);
1238     PolyAQ0  = _mm256_add_pd(PolyAQ0, one);
1239     PolyAQ1  = _mm256_mul_pd(PolyAQ1, x2);
1240     PolyAQ0  = _mm256_add_pd(PolyAQ0, PolyAQ1);
1241
1242     res_erf  = _mm256_mul_pd(PolyAP0, gmx_mm256_inv_pd(PolyAQ0));
1243     res_erf  = _mm256_add_pd(CAoffset, res_erf);
1244     res_erf  = _mm256_mul_pd(x, res_erf);
1245
1246     /* Calculate erfc() in range [1,4.5] */
1247     t       = _mm256_sub_pd(xabs, one);
1248     t2      = _mm256_mul_pd(t, t);
1249
1250     PolyBP0  = _mm256_mul_pd(CBP6, t2);
1251     PolyBP1  = _mm256_mul_pd(CBP5, t2);
1252     PolyBP0  = _mm256_add_pd(PolyBP0, CBP4);
1253     PolyBP1  = _mm256_add_pd(PolyBP1, CBP3);
1254     PolyBP0  = _mm256_mul_pd(PolyBP0, t2);
1255     PolyBP1  = _mm256_mul_pd(PolyBP1, t2);
1256     PolyBP0  = _mm256_add_pd(PolyBP0, CBP2);
1257     PolyBP1  = _mm256_add_pd(PolyBP1, CBP1);
1258     PolyBP0  = _mm256_mul_pd(PolyBP0, t2);
1259     PolyBP1  = _mm256_mul_pd(PolyBP1, t);
1260     PolyBP0  = _mm256_add_pd(PolyBP0, CBP0);
1261     PolyBP0  = _mm256_add_pd(PolyBP0, PolyBP1);
1262
1263     PolyBQ1 = _mm256_mul_pd(CBQ7, t2);
1264     PolyBQ0 = _mm256_mul_pd(CBQ6, t2);
1265     PolyBQ1 = _mm256_add_pd(PolyBQ1, CBQ5);
1266     PolyBQ0 = _mm256_add_pd(PolyBQ0, CBQ4);
1267     PolyBQ1 = _mm256_mul_pd(PolyBQ1, t2);
1268     PolyBQ0 = _mm256_mul_pd(PolyBQ0, t2);
1269     PolyBQ1 = _mm256_add_pd(PolyBQ1, CBQ3);
1270     PolyBQ0 = _mm256_add_pd(PolyBQ0, CBQ2);
1271     PolyBQ1 = _mm256_mul_pd(PolyBQ1, t2);
1272     PolyBQ0 = _mm256_mul_pd(PolyBQ0, t2);
1273     PolyBQ1 = _mm256_add_pd(PolyBQ1, CBQ1);
1274     PolyBQ0 = _mm256_add_pd(PolyBQ0, one);
1275     PolyBQ1 = _mm256_mul_pd(PolyBQ1, t);
1276     PolyBQ0 = _mm256_add_pd(PolyBQ0, PolyBQ1);
1277
1278     res_erfcB = _mm256_mul_pd(PolyBP0, gmx_mm256_inv_pd(PolyBQ0));
1279
1280     res_erfcB = _mm256_mul_pd(res_erfcB, xabs);
1281
1282     /* Calculate erfc() in range [4.5,inf] */
1283     w       = gmx_mm256_inv_pd(xabs);
1284     w2      = _mm256_mul_pd(w, w);
1285
1286     PolyCP0  = _mm256_mul_pd(CCP6, w2);
1287     PolyCP1  = _mm256_mul_pd(CCP5, w2);
1288     PolyCP0  = _mm256_add_pd(PolyCP0, CCP4);
1289     PolyCP1  = _mm256_add_pd(PolyCP1, CCP3);
1290     PolyCP0  = _mm256_mul_pd(PolyCP0, w2);
1291     PolyCP1  = _mm256_mul_pd(PolyCP1, w2);
1292     PolyCP0  = _mm256_add_pd(PolyCP0, CCP2);
1293     PolyCP1  = _mm256_add_pd(PolyCP1, CCP1);
1294     PolyCP0  = _mm256_mul_pd(PolyCP0, w2);
1295     PolyCP1  = _mm256_mul_pd(PolyCP1, w);
1296     PolyCP0  = _mm256_add_pd(PolyCP0, CCP0);
1297     PolyCP0  = _mm256_add_pd(PolyCP0, PolyCP1);
1298
1299     PolyCQ0  = _mm256_mul_pd(CCQ6, w2);
1300     PolyCQ1  = _mm256_mul_pd(CCQ5, w2);
1301     PolyCQ0  = _mm256_add_pd(PolyCQ0, CCQ4);
1302     PolyCQ1  = _mm256_add_pd(PolyCQ1, CCQ3);
1303     PolyCQ0  = _mm256_mul_pd(PolyCQ0, w2);
1304     PolyCQ1  = _mm256_mul_pd(PolyCQ1, w2);
1305     PolyCQ0  = _mm256_add_pd(PolyCQ0, CCQ2);
1306     PolyCQ1  = _mm256_add_pd(PolyCQ1, CCQ1);
1307     PolyCQ0  = _mm256_mul_pd(PolyCQ0, w2);
1308     PolyCQ1  = _mm256_mul_pd(PolyCQ1, w);
1309     PolyCQ0  = _mm256_add_pd(PolyCQ0, one);
1310     PolyCQ0  = _mm256_add_pd(PolyCQ0, PolyCQ1);
1311
1312     expmx2   = gmx_mm256_exp_pd( _mm256_or_pd(signbit, x2) );
1313
1314     res_erfcC = _mm256_mul_pd(PolyCP0, gmx_mm256_inv_pd(PolyCQ0));
1315     res_erfcC = _mm256_add_pd(res_erfcC, CCoffset);
1316     res_erfcC = _mm256_mul_pd(res_erfcC, w);
1317
1318     mask     = _mm256_cmp_pd(xabs, _mm256_set1_pd(4.5), _CMP_GT_OQ);
1319     res_erfc = _mm256_blendv_pd(res_erfcB, res_erfcC, mask);
1320
1321     res_erfc = _mm256_mul_pd(res_erfc, expmx2);
1322
1323     /* erfc(x<0) = 2-erfc(|x|) */
1324     mask     = _mm256_cmp_pd(x, _mm256_setzero_pd(), _CMP_LT_OQ);
1325     res_erfc = _mm256_blendv_pd(res_erfc, _mm256_sub_pd(two, res_erfc), mask);
1326
1327     /* Select erf() or erfc() */
1328     mask = _mm256_cmp_pd(xabs, one, _CMP_LT_OQ);
1329     res  = _mm256_blendv_pd(res_erfc, _mm256_sub_pd(one, res_erf), mask);
1330
1331     return res;
1332 }
1333
1334
1335 static __m128d
1336 gmx_mm_erfc_pd(__m128d x)
1337 {
1338     /* Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75] */
1339     const __m128d CAP4      = _mm_set1_pd(-0.431780540597889301512e-4);
1340     const __m128d CAP3      = _mm_set1_pd(-0.00578562306260059236059);
1341     const __m128d CAP2      = _mm_set1_pd(-0.028593586920219752446);
1342     const __m128d CAP1      = _mm_set1_pd(-0.315924962948621698209);
1343     const __m128d CAP0      = _mm_set1_pd(0.14952975608477029151);
1344
1345     const __m128d CAQ5      = _mm_set1_pd(-0.374089300177174709737e-5);
1346     const __m128d CAQ4      = _mm_set1_pd(0.00015126584532155383535);
1347     const __m128d CAQ3      = _mm_set1_pd(0.00536692680669480725423);
1348     const __m128d CAQ2      = _mm_set1_pd(0.0668686825594046122636);
1349     const __m128d CAQ1      = _mm_set1_pd(0.402604990869284362773);
1350     /* CAQ0 == 1.0 */
1351     const __m128d CAoffset  = _mm_set1_pd(0.9788494110107421875);
1352
1353     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5] */
1354     const __m128d CBP6      = _mm_set1_pd(2.49650423685462752497647637088e-10);
1355     const __m128d CBP5      = _mm_set1_pd(0.00119770193298159629350136085658);
1356     const __m128d CBP4      = _mm_set1_pd(0.0164944422378370965881008942733);
1357     const __m128d CBP3      = _mm_set1_pd(0.0984581468691775932063932439252);
1358     const __m128d CBP2      = _mm_set1_pd(0.317364595806937763843589437418);
1359     const __m128d CBP1      = _mm_set1_pd(0.554167062641455850932670067075);
1360     const __m128d CBP0      = _mm_set1_pd(0.427583576155807163756925301060);
1361     const __m128d CBQ7      = _mm_set1_pd(0.00212288829699830145976198384930);
1362     const __m128d CBQ6      = _mm_set1_pd(0.0334810979522685300554606393425);
1363     const __m128d CBQ5      = _mm_set1_pd(0.2361713785181450957579508850717);
1364     const __m128d CBQ4      = _mm_set1_pd(0.955364736493055670530981883072);
1365     const __m128d CBQ3      = _mm_set1_pd(2.36815675631420037315349279199);
1366     const __m128d CBQ2      = _mm_set1_pd(3.55261649184083035537184223542);
1367     const __m128d CBQ1      = _mm_set1_pd(2.93501136050160872574376997993);
1368     /* CBQ0 == 1.0 */
1369
1370     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf] */
1371     const __m128d CCP6      = _mm_set1_pd(-2.8175401114513378771);
1372     const __m128d CCP5      = _mm_set1_pd(-3.22729451764143718517);
1373     const __m128d CCP4      = _mm_set1_pd(-2.5518551727311523996);
1374     const __m128d CCP3      = _mm_set1_pd(-0.687717681153649930619);
1375     const __m128d CCP2      = _mm_set1_pd(-0.212652252872804219852);
1376     const __m128d CCP1      = _mm_set1_pd(0.0175389834052493308818);
1377     const __m128d CCP0      = _mm_set1_pd(0.00628057170626964891937);
1378
1379     const __m128d CCQ6      = _mm_set1_pd(5.48409182238641741584);
1380     const __m128d CCQ5      = _mm_set1_pd(13.5064170191802889145);
1381     const __m128d CCQ4      = _mm_set1_pd(22.9367376522880577224);
1382     const __m128d CCQ3      = _mm_set1_pd(15.930646027911794143);
1383     const __m128d CCQ2      = _mm_set1_pd(11.0567237927800161565);
1384     const __m128d CCQ1      = _mm_set1_pd(2.79257750980575282228);
1385     /* CCQ0 == 1.0 */
1386     const __m128d CCoffset  = _mm_set1_pd(0.5579090118408203125);
1387
1388     const __m128d one       = _mm_set1_pd(1.0);
1389     const __m128d two       = _mm_set1_pd(2.0);
1390
1391     const __m128d signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000, 0x00000000, 0x80000000, 0x00000000) );
1392
1393     __m128d       xabs, x2, x4, t, t2, w, w2;
1394     __m128d       PolyAP0, PolyAP1, PolyAQ0, PolyAQ1;
1395     __m128d       PolyBP0, PolyBP1, PolyBQ0, PolyBQ1;
1396     __m128d       PolyCP0, PolyCP1, PolyCQ0, PolyCQ1;
1397     __m128d       res_erf, res_erfcB, res_erfcC, res_erfc, res;
1398     __m128d       mask, expmx2;
1399
1400     /* Calculate erf() */
1401     xabs     = gmx_mm_abs_pd(x);
1402     x2       = _mm_mul_pd(x, x);
1403     x4       = _mm_mul_pd(x2, x2);
1404
1405     PolyAP0  = _mm_mul_pd(CAP4, x4);
1406     PolyAP1  = _mm_mul_pd(CAP3, x4);
1407     PolyAP0  = _mm_add_pd(PolyAP0, CAP2);
1408     PolyAP1  = _mm_add_pd(PolyAP1, CAP1);
1409     PolyAP0  = _mm_mul_pd(PolyAP0, x4);
1410     PolyAP1  = _mm_mul_pd(PolyAP1, x2);
1411     PolyAP0  = _mm_add_pd(PolyAP0, CAP0);
1412     PolyAP0  = _mm_add_pd(PolyAP0, PolyAP1);
1413
1414     PolyAQ1  = _mm_mul_pd(CAQ5, x4);
1415     PolyAQ0  = _mm_mul_pd(CAQ4, x4);
1416     PolyAQ1  = _mm_add_pd(PolyAQ1, CAQ3);
1417     PolyAQ0  = _mm_add_pd(PolyAQ0, CAQ2);
1418     PolyAQ1  = _mm_mul_pd(PolyAQ1, x4);
1419     PolyAQ0  = _mm_mul_pd(PolyAQ0, x4);
1420     PolyAQ1  = _mm_add_pd(PolyAQ1, CAQ1);
1421     PolyAQ0  = _mm_add_pd(PolyAQ0, one);
1422     PolyAQ1  = _mm_mul_pd(PolyAQ1, x2);
1423     PolyAQ0  = _mm_add_pd(PolyAQ0, PolyAQ1);
1424
1425     res_erf  = _mm_mul_pd(PolyAP0, gmx_mm_inv_pd(PolyAQ0));
1426     res_erf  = _mm_add_pd(CAoffset, res_erf);
1427     res_erf  = _mm_mul_pd(x, res_erf);
1428
1429     /* Calculate erfc() in range [1,4.5] */
1430     t       = _mm_sub_pd(xabs, one);
1431     t2      = _mm_mul_pd(t, t);
1432
1433     PolyBP0  = _mm_mul_pd(CBP6, t2);
1434     PolyBP1  = _mm_mul_pd(CBP5, t2);
1435     PolyBP0  = _mm_add_pd(PolyBP0, CBP4);
1436     PolyBP1  = _mm_add_pd(PolyBP1, CBP3);
1437     PolyBP0  = _mm_mul_pd(PolyBP0, t2);
1438     PolyBP1  = _mm_mul_pd(PolyBP1, t2);
1439     PolyBP0  = _mm_add_pd(PolyBP0, CBP2);
1440     PolyBP1  = _mm_add_pd(PolyBP1, CBP1);
1441     PolyBP0  = _mm_mul_pd(PolyBP0, t2);
1442     PolyBP1  = _mm_mul_pd(PolyBP1, t);
1443     PolyBP0  = _mm_add_pd(PolyBP0, CBP0);
1444     PolyBP0  = _mm_add_pd(PolyBP0, PolyBP1);
1445
1446     PolyBQ1 = _mm_mul_pd(CBQ7, t2);
1447     PolyBQ0 = _mm_mul_pd(CBQ6, t2);
1448     PolyBQ1 = _mm_add_pd(PolyBQ1, CBQ5);
1449     PolyBQ0 = _mm_add_pd(PolyBQ0, CBQ4);
1450     PolyBQ1 = _mm_mul_pd(PolyBQ1, t2);
1451     PolyBQ0 = _mm_mul_pd(PolyBQ0, t2);
1452     PolyBQ1 = _mm_add_pd(PolyBQ1, CBQ3);
1453     PolyBQ0 = _mm_add_pd(PolyBQ0, CBQ2);
1454     PolyBQ1 = _mm_mul_pd(PolyBQ1, t2);
1455     PolyBQ0 = _mm_mul_pd(PolyBQ0, t2);
1456     PolyBQ1 = _mm_add_pd(PolyBQ1, CBQ1);
1457     PolyBQ0 = _mm_add_pd(PolyBQ0, one);
1458     PolyBQ1 = _mm_mul_pd(PolyBQ1, t);
1459     PolyBQ0 = _mm_add_pd(PolyBQ0, PolyBQ1);
1460
1461     res_erfcB = _mm_mul_pd(PolyBP0, gmx_mm_inv_pd(PolyBQ0));
1462
1463     res_erfcB = _mm_mul_pd(res_erfcB, xabs);
1464
1465     /* Calculate erfc() in range [4.5,inf] */
1466     w       = gmx_mm_inv_pd(xabs);
1467     w2      = _mm_mul_pd(w, w);
1468
1469     PolyCP0  = _mm_mul_pd(CCP6, w2);
1470     PolyCP1  = _mm_mul_pd(CCP5, w2);
1471     PolyCP0  = _mm_add_pd(PolyCP0, CCP4);
1472     PolyCP1  = _mm_add_pd(PolyCP1, CCP3);
1473     PolyCP0  = _mm_mul_pd(PolyCP0, w2);
1474     PolyCP1  = _mm_mul_pd(PolyCP1, w2);
1475     PolyCP0  = _mm_add_pd(PolyCP0, CCP2);
1476     PolyCP1  = _mm_add_pd(PolyCP1, CCP1);
1477     PolyCP0  = _mm_mul_pd(PolyCP0, w2);
1478     PolyCP1  = _mm_mul_pd(PolyCP1, w);
1479     PolyCP0  = _mm_add_pd(PolyCP0, CCP0);
1480     PolyCP0  = _mm_add_pd(PolyCP0, PolyCP1);
1481
1482     PolyCQ0  = _mm_mul_pd(CCQ6, w2);
1483     PolyCQ1  = _mm_mul_pd(CCQ5, w2);
1484     PolyCQ0  = _mm_add_pd(PolyCQ0, CCQ4);
1485     PolyCQ1  = _mm_add_pd(PolyCQ1, CCQ3);
1486     PolyCQ0  = _mm_mul_pd(PolyCQ0, w2);
1487     PolyCQ1  = _mm_mul_pd(PolyCQ1, w2);
1488     PolyCQ0  = _mm_add_pd(PolyCQ0, CCQ2);
1489     PolyCQ1  = _mm_add_pd(PolyCQ1, CCQ1);
1490     PolyCQ0  = _mm_mul_pd(PolyCQ0, w2);
1491     PolyCQ1  = _mm_mul_pd(PolyCQ1, w);
1492     PolyCQ0  = _mm_add_pd(PolyCQ0, one);
1493     PolyCQ0  = _mm_add_pd(PolyCQ0, PolyCQ1);
1494
1495     expmx2   = gmx_mm_exp_pd( _mm_or_pd(signbit, x2) );
1496
1497     res_erfcC = _mm_mul_pd(PolyCP0, gmx_mm_inv_pd(PolyCQ0));
1498     res_erfcC = _mm_add_pd(res_erfcC, CCoffset);
1499     res_erfcC = _mm_mul_pd(res_erfcC, w);
1500
1501     mask     = _mm_cmpgt_pd(xabs, _mm_set1_pd(4.5));
1502     res_erfc = _mm_blendv_pd(res_erfcB, res_erfcC, mask);
1503
1504     res_erfc = _mm_mul_pd(res_erfc, expmx2);
1505
1506     /* erfc(x<0) = 2-erfc(|x|) */
1507     mask     = _mm_cmplt_pd(x, _mm_setzero_pd());
1508     res_erfc = _mm_blendv_pd(res_erfc, _mm_sub_pd(two, res_erfc), mask);
1509
1510     /* Select erf() or erfc() */
1511     mask = _mm_cmplt_pd(xabs, one);
1512     res  = _mm_blendv_pd(res_erfc, _mm_sub_pd(one, res_erf), mask);
1513
1514     return res;
1515 }
1516
1517
1518 /* Calculate the force correction due to PME analytically.
1519  *
1520  * This routine is meant to enable analytical evaluation of the
1521  * direct-space PME electrostatic force to avoid tables.
1522  *
1523  * The direct-space potential should be Erfc(beta*r)/r, but there
1524  * are some problems evaluating that:
1525  *
1526  * First, the error function is difficult (read: expensive) to
1527  * approxmiate accurately for intermediate to large arguments, and
1528  * this happens already in ranges of beta*r that occur in simulations.
1529  * Second, we now try to avoid calculating potentials in Gromacs but
1530  * use forces directly.
1531  *
1532  * We can simply things slight by noting that the PME part is really
1533  * a correction to the normal Coulomb force since Erfc(z)=1-Erf(z), i.e.
1534  *
1535  * V= 1/r - Erf(beta*r)/r
1536  *
1537  * The first term we already have from the inverse square root, so
1538  * that we can leave out of this routine.
1539  *
1540  * For pme tolerances of 1e-3 to 1e-8 and cutoffs of 0.5nm to 1.8nm,
1541  * the argument beta*r will be in the range 0.15 to ~4. Use your
1542  * favorite plotting program to realize how well-behaved Erf(z)/z is
1543  * in this range!
1544  *
1545  * We approximate f(z)=erf(z)/z with a rational minimax polynomial.
1546  * However, it turns out it is more efficient to approximate f(z)/z and
1547  * then only use even powers. This is another minor optimization, since
1548  * we actually WANT f(z)/z, because it is going to be multiplied by
1549  * the vector between the two atoms to get the vectorial force. The
1550  * fastest flops are the ones we can avoid calculating!
1551  *
1552  * So, here's how it should be used:
1553  *
1554  * 1. Calculate r^2.
1555  * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
1556  * 3. Evaluate this routine with z^2 as the argument.
1557  * 4. The return value is the expression:
1558  *
1559  *
1560  *       2*exp(-z^2)     erf(z)
1561  *       ------------ - --------
1562  *       sqrt(Pi)*z^2      z^3
1563  *
1564  * 5. Multiply the entire expression by beta^3. This will get you
1565  *
1566  *       beta^3*2*exp(-z^2)     beta^3*erf(z)
1567  *       ------------------  - ---------------
1568  *          sqrt(Pi)*z^2            z^3
1569  *
1570  *    or, switching back to r (z=r*beta):
1571  *
1572  *       2*beta*exp(-r^2*beta^2)   erf(r*beta)
1573  *       ----------------------- - -----------
1574  *            sqrt(Pi)*r^2            r^3
1575  *
1576  *
1577  *    With a bit of math exercise you should be able to confirm that
1578  *    this is exactly D[Erf[beta*r]/r,r] divided by r another time.
1579  *
1580  * 6. Add the result to 1/r^3, multiply by the product of the charges,
1581  *    and you have your force (divided by r). A final multiplication
1582  *    with the vector connecting the two particles and you have your
1583  *    vectorial force to add to the particles.
1584  *
1585  */
1586 static __m256d
1587 gmx_mm256_pmecorrF_pd(__m256d z2)
1588 {
1589     const __m256d  FN10     = _mm256_set1_pd(-8.0072854618360083154e-14);
1590     const __m256d  FN9      = _mm256_set1_pd(1.1859116242260148027e-11);
1591     const __m256d  FN8      = _mm256_set1_pd(-8.1490406329798423616e-10);
1592     const __m256d  FN7      = _mm256_set1_pd(3.4404793543907847655e-8);
1593     const __m256d  FN6      = _mm256_set1_pd(-9.9471420832602741006e-7);
1594     const __m256d  FN5      = _mm256_set1_pd(0.000020740315999115847456);
1595     const __m256d  FN4      = _mm256_set1_pd(-0.00031991745139313364005);
1596     const __m256d  FN3      = _mm256_set1_pd(0.0035074449373659008203);
1597     const __m256d  FN2      = _mm256_set1_pd(-0.031750380176100813405);
1598     const __m256d  FN1      = _mm256_set1_pd(0.13884101728898463426);
1599     const __m256d  FN0      = _mm256_set1_pd(-0.75225277815249618847);
1600
1601     const __m256d  FD5      = _mm256_set1_pd(0.000016009278224355026701);
1602     const __m256d  FD4      = _mm256_set1_pd(0.00051055686934806966046);
1603     const __m256d  FD3      = _mm256_set1_pd(0.0081803507497974289008);
1604     const __m256d  FD2      = _mm256_set1_pd(0.077181146026670287235);
1605     const __m256d  FD1      = _mm256_set1_pd(0.41543303143712535988);
1606     const __m256d  FD0      = _mm256_set1_pd(1.0);
1607
1608     __m256d        z4;
1609     __m256d        polyFN0, polyFN1, polyFD0, polyFD1;
1610
1611     z4             = _mm256_mul_pd(z2, z2);
1612
1613     polyFD1        = _mm256_mul_pd(FD5, z4);
1614     polyFD0        = _mm256_mul_pd(FD4, z4);
1615     polyFD1        = _mm256_add_pd(polyFD1, FD3);
1616     polyFD0        = _mm256_add_pd(polyFD0, FD2);
1617     polyFD1        = _mm256_mul_pd(polyFD1, z4);
1618     polyFD0        = _mm256_mul_pd(polyFD0, z4);
1619     polyFD1        = _mm256_add_pd(polyFD1, FD1);
1620     polyFD0        = _mm256_add_pd(polyFD0, FD0);
1621     polyFD1        = _mm256_mul_pd(polyFD1, z2);
1622     polyFD0        = _mm256_add_pd(polyFD0, polyFD1);
1623
1624     polyFD0        = gmx_mm256_inv_pd(polyFD0);
1625
1626     polyFN0        = _mm256_mul_pd(FN10, z4);
1627     polyFN1        = _mm256_mul_pd(FN9, z4);
1628     polyFN0        = _mm256_add_pd(polyFN0, FN8);
1629     polyFN1        = _mm256_add_pd(polyFN1, FN7);
1630     polyFN0        = _mm256_mul_pd(polyFN0, z4);
1631     polyFN1        = _mm256_mul_pd(polyFN1, z4);
1632     polyFN0        = _mm256_add_pd(polyFN0, FN6);
1633     polyFN1        = _mm256_add_pd(polyFN1, FN5);
1634     polyFN0        = _mm256_mul_pd(polyFN0, z4);
1635     polyFN1        = _mm256_mul_pd(polyFN1, z4);
1636     polyFN0        = _mm256_add_pd(polyFN0, FN4);
1637     polyFN1        = _mm256_add_pd(polyFN1, FN3);
1638     polyFN0        = _mm256_mul_pd(polyFN0, z4);
1639     polyFN1        = _mm256_mul_pd(polyFN1, z4);
1640     polyFN0        = _mm256_add_pd(polyFN0, FN2);
1641     polyFN1        = _mm256_add_pd(polyFN1, FN1);
1642     polyFN0        = _mm256_mul_pd(polyFN0, z4);
1643     polyFN1        = _mm256_mul_pd(polyFN1, z2);
1644     polyFN0        = _mm256_add_pd(polyFN0, FN0);
1645     polyFN0        = _mm256_add_pd(polyFN0, polyFN1);
1646
1647     return _mm256_mul_pd(polyFN0, polyFD0);
1648 }
1649
1650
1651 static __m128d
1652 gmx_mm_pmecorrF_pd(__m128d z2)
1653 {
1654     const __m128d  FN10     = _mm_set1_pd(-8.0072854618360083154e-14);
1655     const __m128d  FN9      = _mm_set1_pd(1.1859116242260148027e-11);
1656     const __m128d  FN8      = _mm_set1_pd(-8.1490406329798423616e-10);
1657     const __m128d  FN7      = _mm_set1_pd(3.4404793543907847655e-8);
1658     const __m128d  FN6      = _mm_set1_pd(-9.9471420832602741006e-7);
1659     const __m128d  FN5      = _mm_set1_pd(0.000020740315999115847456);
1660     const __m128d  FN4      = _mm_set1_pd(-0.00031991745139313364005);
1661     const __m128d  FN3      = _mm_set1_pd(0.0035074449373659008203);
1662     const __m128d  FN2      = _mm_set1_pd(-0.031750380176100813405);
1663     const __m128d  FN1      = _mm_set1_pd(0.13884101728898463426);
1664     const __m128d  FN0      = _mm_set1_pd(-0.75225277815249618847);
1665
1666     const __m128d  FD5      = _mm_set1_pd(0.000016009278224355026701);
1667     const __m128d  FD4      = _mm_set1_pd(0.00051055686934806966046);
1668     const __m128d  FD3      = _mm_set1_pd(0.0081803507497974289008);
1669     const __m128d  FD2      = _mm_set1_pd(0.077181146026670287235);
1670     const __m128d  FD1      = _mm_set1_pd(0.41543303143712535988);
1671     const __m128d  FD0      = _mm_set1_pd(1.0);
1672
1673     __m128d        z4;
1674     __m128d        polyFN0, polyFN1, polyFD0, polyFD1;
1675
1676     z4             = _mm_mul_pd(z2, z2);
1677
1678     polyFD1        = _mm_mul_pd(FD5, z4);
1679     polyFD0        = _mm_mul_pd(FD4, z4);
1680     polyFD1        = _mm_add_pd(polyFD1, FD3);
1681     polyFD0        = _mm_add_pd(polyFD0, FD2);
1682     polyFD1        = _mm_mul_pd(polyFD1, z4);
1683     polyFD0        = _mm_mul_pd(polyFD0, z4);
1684     polyFD1        = _mm_add_pd(polyFD1, FD1);
1685     polyFD0        = _mm_add_pd(polyFD0, FD0);
1686     polyFD1        = _mm_mul_pd(polyFD1, z2);
1687     polyFD0        = _mm_add_pd(polyFD0, polyFD1);
1688
1689     polyFD0        = gmx_mm_inv_pd(polyFD0);
1690
1691     polyFN0        = _mm_mul_pd(FN10, z4);
1692     polyFN1        = _mm_mul_pd(FN9, z4);
1693     polyFN0        = _mm_add_pd(polyFN0, FN8);
1694     polyFN1        = _mm_add_pd(polyFN1, FN7);
1695     polyFN0        = _mm_mul_pd(polyFN0, z4);
1696     polyFN1        = _mm_mul_pd(polyFN1, z4);
1697     polyFN0        = _mm_add_pd(polyFN0, FN6);
1698     polyFN1        = _mm_add_pd(polyFN1, FN5);
1699     polyFN0        = _mm_mul_pd(polyFN0, z4);
1700     polyFN1        = _mm_mul_pd(polyFN1, z4);
1701     polyFN0        = _mm_add_pd(polyFN0, FN4);
1702     polyFN1        = _mm_add_pd(polyFN1, FN3);
1703     polyFN0        = _mm_mul_pd(polyFN0, z4);
1704     polyFN1        = _mm_mul_pd(polyFN1, z4);
1705     polyFN0        = _mm_add_pd(polyFN0, FN2);
1706     polyFN1        = _mm_add_pd(polyFN1, FN1);
1707     polyFN0        = _mm_mul_pd(polyFN0, z4);
1708     polyFN1        = _mm_mul_pd(polyFN1, z2);
1709     polyFN0        = _mm_add_pd(polyFN0, FN0);
1710     polyFN0        = _mm_add_pd(polyFN0, polyFN1);
1711
1712     return _mm_mul_pd(polyFN0, polyFD0);
1713 }
1714
1715
1716
1717
1718 /* Calculate the potential correction due to PME analytically.
1719  *
1720  * See gmx_mm256_pmecorrF_ps() for details about the approximation.
1721  *
1722  * This routine calculates Erf(z)/z, although you should provide z^2
1723  * as the input argument.
1724  *
1725  * Here's how it should be used:
1726  *
1727  * 1. Calculate r^2.
1728  * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
1729  * 3. Evaluate this routine with z^2 as the argument.
1730  * 4. The return value is the expression:
1731  *
1732  *
1733  *        erf(z)
1734  *       --------
1735  *          z
1736  *
1737  * 5. Multiply the entire expression by beta and switching back to r (z=r*beta):
1738  *
1739  *       erf(r*beta)
1740  *       -----------
1741  *           r
1742  *
1743  * 6. Subtract the result from 1/r, multiply by the product of the charges,
1744  *    and you have your potential.
1745  *
1746  */
1747 static __m256d
1748 gmx_mm256_pmecorrV_pd(__m256d z2)
1749 {
1750     const __m256d  VN9      = _mm256_set1_pd(-9.3723776169321855475e-13);
1751     const __m256d  VN8      = _mm256_set1_pd(1.2280156762674215741e-10);
1752     const __m256d  VN7      = _mm256_set1_pd(-7.3562157912251309487e-9);
1753     const __m256d  VN6      = _mm256_set1_pd(2.6215886208032517509e-7);
1754     const __m256d  VN5      = _mm256_set1_pd(-4.9532491651265819499e-6);
1755     const __m256d  VN4      = _mm256_set1_pd(0.00025907400778966060389);
1756     const __m256d  VN3      = _mm256_set1_pd(0.0010585044856156469792);
1757     const __m256d  VN2      = _mm256_set1_pd(0.045247661136833092885);
1758     const __m256d  VN1      = _mm256_set1_pd(0.11643931522926034421);
1759     const __m256d  VN0      = _mm256_set1_pd(1.1283791671726767970);
1760
1761     const __m256d  VD5      = _mm256_set1_pd(0.000021784709867336150342);
1762     const __m256d  VD4      = _mm256_set1_pd(0.00064293662010911388448);
1763     const __m256d  VD3      = _mm256_set1_pd(0.0096311444822588683504);
1764     const __m256d  VD2      = _mm256_set1_pd(0.085608012351550627051);
1765     const __m256d  VD1      = _mm256_set1_pd(0.43652499166614811084);
1766     const __m256d  VD0      = _mm256_set1_pd(1.0);
1767
1768     __m256d        z4;
1769     __m256d        polyVN0, polyVN1, polyVD0, polyVD1;
1770
1771     z4             = _mm256_mul_pd(z2, z2);
1772
1773     polyVD1        = _mm256_mul_pd(VD5, z4);
1774     polyVD0        = _mm256_mul_pd(VD4, z4);
1775     polyVD1        = _mm256_add_pd(polyVD1, VD3);
1776     polyVD0        = _mm256_add_pd(polyVD0, VD2);
1777     polyVD1        = _mm256_mul_pd(polyVD1, z4);
1778     polyVD0        = _mm256_mul_pd(polyVD0, z4);
1779     polyVD1        = _mm256_add_pd(polyVD1, VD1);
1780     polyVD0        = _mm256_add_pd(polyVD0, VD0);
1781     polyVD1        = _mm256_mul_pd(polyVD1, z2);
1782     polyVD0        = _mm256_add_pd(polyVD0, polyVD1);
1783
1784     polyVD0        = gmx_mm256_inv_pd(polyVD0);
1785
1786     polyVN1        = _mm256_mul_pd(VN9, z4);
1787     polyVN0        = _mm256_mul_pd(VN8, z4);
1788     polyVN1        = _mm256_add_pd(polyVN1, VN7);
1789     polyVN0        = _mm256_add_pd(polyVN0, VN6);
1790     polyVN1        = _mm256_mul_pd(polyVN1, z4);
1791     polyVN0        = _mm256_mul_pd(polyVN0, z4);
1792     polyVN1        = _mm256_add_pd(polyVN1, VN5);
1793     polyVN0        = _mm256_add_pd(polyVN0, VN4);
1794     polyVN1        = _mm256_mul_pd(polyVN1, z4);
1795     polyVN0        = _mm256_mul_pd(polyVN0, z4);
1796     polyVN1        = _mm256_add_pd(polyVN1, VN3);
1797     polyVN0        = _mm256_add_pd(polyVN0, VN2);
1798     polyVN1        = _mm256_mul_pd(polyVN1, z4);
1799     polyVN0        = _mm256_mul_pd(polyVN0, z4);
1800     polyVN1        = _mm256_add_pd(polyVN1, VN1);
1801     polyVN0        = _mm256_add_pd(polyVN0, VN0);
1802     polyVN1        = _mm256_mul_pd(polyVN1, z2);
1803     polyVN0        = _mm256_add_pd(polyVN0, polyVN1);
1804
1805     return _mm256_mul_pd(polyVN0, polyVD0);
1806 }
1807
1808
1809 static __m128d
1810 gmx_mm_pmecorrV_pd(__m128d z2)
1811 {
1812     const __m128d  VN9      = _mm_set1_pd(-9.3723776169321855475e-13);
1813     const __m128d  VN8      = _mm_set1_pd(1.2280156762674215741e-10);
1814     const __m128d  VN7      = _mm_set1_pd(-7.3562157912251309487e-9);
1815     const __m128d  VN6      = _mm_set1_pd(2.6215886208032517509e-7);
1816     const __m128d  VN5      = _mm_set1_pd(-4.9532491651265819499e-6);
1817     const __m128d  VN4      = _mm_set1_pd(0.00025907400778966060389);
1818     const __m128d  VN3      = _mm_set1_pd(0.0010585044856156469792);
1819     const __m128d  VN2      = _mm_set1_pd(0.045247661136833092885);
1820     const __m128d  VN1      = _mm_set1_pd(0.11643931522926034421);
1821     const __m128d  VN0      = _mm_set1_pd(1.1283791671726767970);
1822
1823     const __m128d  VD5      = _mm_set1_pd(0.000021784709867336150342);
1824     const __m128d  VD4      = _mm_set1_pd(0.00064293662010911388448);
1825     const __m128d  VD3      = _mm_set1_pd(0.0096311444822588683504);
1826     const __m128d  VD2      = _mm_set1_pd(0.085608012351550627051);
1827     const __m128d  VD1      = _mm_set1_pd(0.43652499166614811084);
1828     const __m128d  VD0      = _mm_set1_pd(1.0);
1829
1830     __m128d        z4;
1831     __m128d        polyVN0, polyVN1, polyVD0, polyVD1;
1832
1833     z4             = _mm_mul_pd(z2, z2);
1834
1835     polyVD1        = _mm_mul_pd(VD5, z4);
1836     polyVD0        = _mm_mul_pd(VD4, z4);
1837     polyVD1        = _mm_add_pd(polyVD1, VD3);
1838     polyVD0        = _mm_add_pd(polyVD0, VD2);
1839     polyVD1        = _mm_mul_pd(polyVD1, z4);
1840     polyVD0        = _mm_mul_pd(polyVD0, z4);
1841     polyVD1        = _mm_add_pd(polyVD1, VD1);
1842     polyVD0        = _mm_add_pd(polyVD0, VD0);
1843     polyVD1        = _mm_mul_pd(polyVD1, z2);
1844     polyVD0        = _mm_add_pd(polyVD0, polyVD1);
1845
1846     polyVD0        = gmx_mm_inv_pd(polyVD0);
1847
1848     polyVN1        = _mm_mul_pd(VN9, z4);
1849     polyVN0        = _mm_mul_pd(VN8, z4);
1850     polyVN1        = _mm_add_pd(polyVN1, VN7);
1851     polyVN0        = _mm_add_pd(polyVN0, VN6);
1852     polyVN1        = _mm_mul_pd(polyVN1, z4);
1853     polyVN0        = _mm_mul_pd(polyVN0, z4);
1854     polyVN1        = _mm_add_pd(polyVN1, VN5);
1855     polyVN0        = _mm_add_pd(polyVN0, VN4);
1856     polyVN1        = _mm_mul_pd(polyVN1, z4);
1857     polyVN0        = _mm_mul_pd(polyVN0, z4);
1858     polyVN1        = _mm_add_pd(polyVN1, VN3);
1859     polyVN0        = _mm_add_pd(polyVN0, VN2);
1860     polyVN1        = _mm_mul_pd(polyVN1, z4);
1861     polyVN0        = _mm_mul_pd(polyVN0, z4);
1862     polyVN1        = _mm_add_pd(polyVN1, VN1);
1863     polyVN0        = _mm_add_pd(polyVN0, VN0);
1864     polyVN1        = _mm_mul_pd(polyVN1, z2);
1865     polyVN0        = _mm_add_pd(polyVN0, polyVN1);
1866
1867     return _mm_mul_pd(polyVN0, polyVD0);
1868 }
1869
1870
1871
1872 static int
1873 gmx_mm256_sincos_pd(__m256d  x,
1874                     __m256d *sinval,
1875                     __m256d *cosval)
1876 {
1877 #ifdef _MSC_VER
1878     __declspec(align(16))
1879     const double sintable[34] =
1880     {
1881         1.00000000000000000e+00, 0.00000000000000000e+00,
1882         9.95184726672196929e-01, 9.80171403295606036e-02,
1883         9.80785280403230431e-01, 1.95090322016128248e-01,
1884         9.56940335732208824e-01, 2.90284677254462331e-01,
1885         9.23879532511286738e-01, 3.82683432365089782e-01,
1886         8.81921264348355050e-01, 4.71396736825997642e-01,
1887         8.31469612302545236e-01, 5.55570233019602178e-01,
1888         7.73010453362736993e-01, 6.34393284163645488e-01,
1889         7.07106781186547573e-01, 7.07106781186547462e-01,
1890         6.34393284163645599e-01, 7.73010453362736882e-01,
1891         5.55570233019602289e-01, 8.31469612302545125e-01,
1892         4.71396736825997809e-01, 8.81921264348354939e-01,
1893         3.82683432365089837e-01, 9.23879532511286738e-01,
1894         2.90284677254462276e-01, 9.56940335732208935e-01,
1895         1.95090322016128304e-01, 9.80785280403230431e-01,
1896         9.80171403295607702e-02, 9.95184726672196818e-01,
1897         0.0, 1.00000000000000000e+00
1898     };
1899 #else
1900     const __m128d sintable[17] =
1901     {
1902         _mm_set_pd( 0.0, 1.0 ),
1903         _mm_set_pd( sin(  1.0 * (M_PI/2.0) / 16.0), cos(  1.0 * (M_PI/2.0) / 16.0) ),
1904         _mm_set_pd( sin(  2.0 * (M_PI/2.0) / 16.0), cos(  2.0 * (M_PI/2.0) / 16.0) ),
1905         _mm_set_pd( sin(  3.0 * (M_PI/2.0) / 16.0), cos(  3.0 * (M_PI/2.0) / 16.0) ),
1906         _mm_set_pd( sin(  4.0 * (M_PI/2.0) / 16.0), cos(  4.0 * (M_PI/2.0) / 16.0) ),
1907         _mm_set_pd( sin(  5.0 * (M_PI/2.0) / 16.0), cos(  5.0 * (M_PI/2.0) / 16.0) ),
1908         _mm_set_pd( sin(  6.0 * (M_PI/2.0) / 16.0), cos(  6.0 * (M_PI/2.0) / 16.0) ),
1909         _mm_set_pd( sin(  7.0 * (M_PI/2.0) / 16.0), cos(  7.0 * (M_PI/2.0) / 16.0) ),
1910         _mm_set_pd( sin(  8.0 * (M_PI/2.0) / 16.0), cos(  8.0 * (M_PI/2.0) / 16.0) ),
1911         _mm_set_pd( sin(  9.0 * (M_PI/2.0) / 16.0), cos(  9.0 * (M_PI/2.0) / 16.0) ),
1912         _mm_set_pd( sin( 10.0 * (M_PI/2.0) / 16.0), cos( 10.0 * (M_PI/2.0) / 16.0) ),
1913         _mm_set_pd( sin( 11.0 * (M_PI/2.0) / 16.0), cos( 11.0 * (M_PI/2.0) / 16.0) ),
1914         _mm_set_pd( sin( 12.0 * (M_PI/2.0) / 16.0), cos( 12.0 * (M_PI/2.0) / 16.0) ),
1915         _mm_set_pd( sin( 13.0 * (M_PI/2.0) / 16.0), cos( 13.0 * (M_PI/2.0) / 16.0) ),
1916         _mm_set_pd( sin( 14.0 * (M_PI/2.0) / 16.0), cos( 14.0 * (M_PI/2.0) / 16.0) ),
1917         _mm_set_pd( sin( 15.0 * (M_PI/2.0) / 16.0), cos( 15.0 * (M_PI/2.0) / 16.0) ),
1918         _mm_set_pd(  1.0, 0.0 )
1919     };
1920 #endif
1921
1922     const __m256d signmask    = _mm256_castsi256_pd( _mm256_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF,
1923                                                                       0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
1924
1925     const __m256d tabscale      = _mm256_set1_pd(32.0/M_PI);
1926     const __m256d invtabscale0  = _mm256_set1_pd(9.81747508049011230469e-02);
1927     const __m256d invtabscale1  = _mm256_set1_pd(1.96197799156550576057e-08);
1928     const __m128i ione          = _mm_set1_epi32(1);
1929     const __m128i i32           = _mm_set1_epi32(32);
1930     const __m128i i16           = _mm_set1_epi32(16);
1931     const __m128i tabmask       = _mm_set1_epi32(0x3F);
1932     const __m256d sinP7         = _mm256_set1_pd(-1.0/5040.0);
1933     const __m256d sinP5         = _mm256_set1_pd(1.0/120.0);
1934     const __m256d sinP3         = _mm256_set1_pd(-1.0/6.0);
1935     const __m256d sinP1         = _mm256_set1_pd(1.0);
1936
1937     const __m256d cosP6         = _mm256_set1_pd(-1.0/720.0);
1938     const __m256d cosP4         = _mm256_set1_pd(1.0/24.0);
1939     const __m256d cosP2         = _mm256_set1_pd(-1.0/2.0);
1940     const __m256d cosP0         = _mm256_set1_pd(1.0);
1941
1942     __m256d       scalex;
1943     __m128i       tabidx, corridx;
1944     __m256d       xabs, z, z2, polySin, polyCos;
1945     __m256d       xpoint;
1946     __m256d       t1, t2;
1947
1948     __m256d       sinpoint, cospoint;
1949     __m256d       xsign, ssign, csign;
1950     __m128i       imask, sswapsign, cswapsign;
1951
1952     xsign    = _mm256_andnot_pd(signmask, x);
1953     xabs     = _mm256_and_pd(x, signmask);
1954
1955     scalex   = _mm256_mul_pd(tabscale, xabs);
1956     tabidx   = _mm256_cvtpd_epi32(scalex);
1957
1958     xpoint   = _mm256_round_pd(scalex, _MM_FROUND_TO_NEAREST_INT);
1959
1960     /* Extended precision arithmetics */
1961     z        = _mm256_sub_pd(xabs, _mm256_mul_pd(invtabscale0, xpoint));
1962     z        = _mm256_sub_pd(z, _mm256_mul_pd(invtabscale1, xpoint));
1963
1964     /* Range reduction to 0..2*Pi */
1965     tabidx   = _mm_and_si128(tabidx, tabmask);
1966
1967     /* tabidx is now in range [0,..,64] */
1968     imask     = _mm_cmpgt_epi32(tabidx, i32);
1969     sswapsign = imask;
1970     cswapsign = imask;
1971     corridx   = _mm_and_si128(imask, i32);
1972     tabidx    = _mm_sub_epi32(tabidx, corridx);
1973
1974     /* tabidx is now in range [0..32] */
1975     imask     = _mm_cmpgt_epi32(tabidx, i16);
1976     cswapsign = _mm_xor_si128(cswapsign, imask);
1977     corridx   = _mm_sub_epi32(i32, tabidx);
1978     tabidx    = _mm_blendv_epi8(tabidx, corridx, imask);
1979     /* tabidx is now in range [0..16] */
1980     ssign     = _mm256_cvtepi32_pd( _mm_or_si128( sswapsign, ione ) );
1981     csign     = _mm256_cvtepi32_pd( _mm_or_si128( cswapsign, ione ) );
1982
1983     /* First lookup into table */
1984 #ifdef _MSC_VER
1985     t1       = _mm256_insertf128_pd(_mm256_castpd128_pd256(_mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx, 0))),
1986                                     _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx, 2)), 0x1);
1987     t2       = _mm256_insertf128_pd(_mm256_castpd128_pd256(_mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx, 1))),
1988                                     _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx, 3)), 0x1);
1989 #else
1990     t1       = _mm256_insertf128_pd(_mm256_castpd128_pd256(sintable[_mm_extract_epi32(tabidx, 0)]),
1991                                     sintable[_mm_extract_epi32(tabidx, 2)], 0x1);
1992     t2       = _mm256_insertf128_pd(_mm256_castpd128_pd256(sintable[_mm_extract_epi32(tabidx, 1)]),
1993                                     sintable[_mm_extract_epi32(tabidx, 3)], 0x1);
1994 #endif
1995
1996     sinpoint  = _mm256_unpackhi_pd(t1, t2);
1997     cospoint  = _mm256_unpacklo_pd(t1, t2);
1998
1999     sinpoint = _mm256_mul_pd(sinpoint, ssign);
2000     cospoint = _mm256_mul_pd(cospoint, csign);
2001
2002     z2       = _mm256_mul_pd(z, z);
2003
2004     polySin  = _mm256_mul_pd(sinP7, z2);
2005     polySin  = _mm256_add_pd(polySin, sinP5);
2006     polySin  = _mm256_mul_pd(polySin, z2);
2007     polySin  = _mm256_add_pd(polySin, sinP3);
2008     polySin  = _mm256_mul_pd(polySin, z2);
2009     polySin  = _mm256_add_pd(polySin, sinP1);
2010     polySin  = _mm256_mul_pd(polySin, z);
2011
2012     polyCos  = _mm256_mul_pd(cosP6, z2);
2013     polyCos  = _mm256_add_pd(polyCos, cosP4);
2014     polyCos  = _mm256_mul_pd(polyCos, z2);
2015     polyCos  = _mm256_add_pd(polyCos, cosP2);
2016     polyCos  = _mm256_mul_pd(polyCos, z2);
2017     polyCos  = _mm256_add_pd(polyCos, cosP0);
2018
2019     *sinval  = _mm256_xor_pd(_mm256_add_pd( _mm256_mul_pd(sinpoint, polyCos), _mm256_mul_pd(cospoint, polySin) ), xsign);
2020     *cosval  = _mm256_sub_pd( _mm256_mul_pd(cospoint, polyCos), _mm256_mul_pd(sinpoint, polySin) );
2021
2022     return 0;
2023 }
2024
2025 static int
2026 gmx_mm_sincos_pd(__m128d  x,
2027                  __m128d *sinval,
2028                  __m128d *cosval)
2029 {
2030 #ifdef _MSC_VER
2031     __declspec(align(16))
2032     const double sintable[34] =
2033     {
2034         1.00000000000000000e+00, 0.00000000000000000e+00,
2035         9.95184726672196929e-01, 9.80171403295606036e-02,
2036         9.80785280403230431e-01, 1.95090322016128248e-01,
2037         9.56940335732208824e-01, 2.90284677254462331e-01,
2038         9.23879532511286738e-01, 3.82683432365089782e-01,
2039         8.81921264348355050e-01, 4.71396736825997642e-01,
2040         8.31469612302545236e-01, 5.55570233019602178e-01,
2041         7.73010453362736993e-01, 6.34393284163645488e-01,
2042         7.07106781186547573e-01, 7.07106781186547462e-01,
2043         6.34393284163645599e-01, 7.73010453362736882e-01,
2044         5.55570233019602289e-01, 8.31469612302545125e-01,
2045         4.71396736825997809e-01, 8.81921264348354939e-01,
2046         3.82683432365089837e-01, 9.23879532511286738e-01,
2047         2.90284677254462276e-01, 9.56940335732208935e-01,
2048         1.95090322016128304e-01, 9.80785280403230431e-01,
2049         9.80171403295607702e-02, 9.95184726672196818e-01,
2050         0.0, 1.00000000000000000e+00
2051     };
2052 #else
2053     const __m128d sintable[17] =
2054     {
2055         _mm_set_pd( 0.0, 1.0 ),
2056         _mm_set_pd( sin(  1.0 * (M_PI/2.0) / 16.0), cos(  1.0 * (M_PI/2.0) / 16.0) ),
2057         _mm_set_pd( sin(  2.0 * (M_PI/2.0) / 16.0), cos(  2.0 * (M_PI/2.0) / 16.0) ),
2058         _mm_set_pd( sin(  3.0 * (M_PI/2.0) / 16.0), cos(  3.0 * (M_PI/2.0) / 16.0) ),
2059         _mm_set_pd( sin(  4.0 * (M_PI/2.0) / 16.0), cos(  4.0 * (M_PI/2.0) / 16.0) ),
2060         _mm_set_pd( sin(  5.0 * (M_PI/2.0) / 16.0), cos(  5.0 * (M_PI/2.0) / 16.0) ),
2061         _mm_set_pd( sin(  6.0 * (M_PI/2.0) / 16.0), cos(  6.0 * (M_PI/2.0) / 16.0) ),
2062         _mm_set_pd( sin(  7.0 * (M_PI/2.0) / 16.0), cos(  7.0 * (M_PI/2.0) / 16.0) ),
2063         _mm_set_pd( sin(  8.0 * (M_PI/2.0) / 16.0), cos(  8.0 * (M_PI/2.0) / 16.0) ),
2064         _mm_set_pd( sin(  9.0 * (M_PI/2.0) / 16.0), cos(  9.0 * (M_PI/2.0) / 16.0) ),
2065         _mm_set_pd( sin( 10.0 * (M_PI/2.0) / 16.0), cos( 10.0 * (M_PI/2.0) / 16.0) ),
2066         _mm_set_pd( sin( 11.0 * (M_PI/2.0) / 16.0), cos( 11.0 * (M_PI/2.0) / 16.0) ),
2067         _mm_set_pd( sin( 12.0 * (M_PI/2.0) / 16.0), cos( 12.0 * (M_PI/2.0) / 16.0) ),
2068         _mm_set_pd( sin( 13.0 * (M_PI/2.0) / 16.0), cos( 13.0 * (M_PI/2.0) / 16.0) ),
2069         _mm_set_pd( sin( 14.0 * (M_PI/2.0) / 16.0), cos( 14.0 * (M_PI/2.0) / 16.0) ),
2070         _mm_set_pd( sin( 15.0 * (M_PI/2.0) / 16.0), cos( 15.0 * (M_PI/2.0) / 16.0) ),
2071         _mm_set_pd(  1.0, 0.0 )
2072     };
2073 #endif
2074
2075     const __m128d signmask       = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
2076
2077     const __m128d tabscale      = _mm_set1_pd(32.0/M_PI);
2078     const __m128d invtabscale0  = _mm_set1_pd(9.81747508049011230469e-02);
2079     const __m128d invtabscale1  = _mm_set1_pd(1.96197799156550576057e-08);
2080     const __m128i ione          = _mm_set1_epi32(1);
2081     const __m128i i32           = _mm_set1_epi32(32);
2082     const __m128i i16           = _mm_set1_epi32(16);
2083     const __m128i tabmask       = _mm_set1_epi32(0x3F);
2084     const __m128d sinP7         = _mm_set1_pd(-1.0/5040.0);
2085     const __m128d sinP5         = _mm_set1_pd(1.0/120.0);
2086     const __m128d sinP3         = _mm_set1_pd(-1.0/6.0);
2087     const __m128d sinP1         = _mm_set1_pd(1.0);
2088
2089     const __m128d cosP6         = _mm_set1_pd(-1.0/720.0);
2090     const __m128d cosP4         = _mm_set1_pd(1.0/24.0);
2091     const __m128d cosP2         = _mm_set1_pd(-1.0/2.0);
2092     const __m128d cosP0         = _mm_set1_pd(1.0);
2093
2094     __m128d       scalex;
2095     __m128i       tabidx, corridx;
2096     __m128d       xabs, z, z2, polySin, polyCos;
2097     __m128d       xpoint;
2098     __m128d       ypoint0, ypoint1;
2099
2100     __m128d       sinpoint, cospoint;
2101     __m128d       xsign, ssign, csign;
2102     __m128i       imask, sswapsign, cswapsign;
2103
2104     xsign    = _mm_andnot_pd(signmask, x);
2105     xabs     = _mm_and_pd(x, signmask);
2106
2107     scalex   = _mm_mul_pd(tabscale, xabs);
2108     tabidx   = _mm_cvtpd_epi32(scalex);
2109
2110     xpoint   = _mm_round_pd(scalex, _MM_FROUND_TO_NEAREST_INT);
2111
2112     /* Extended precision arithmetics */
2113     z        = _mm_sub_pd(xabs, _mm_mul_pd(invtabscale0, xpoint));
2114     z        = _mm_sub_pd(z, _mm_mul_pd(invtabscale1, xpoint));
2115
2116     /* Range reduction to 0..2*Pi */
2117     tabidx   = _mm_and_si128(tabidx, tabmask);
2118
2119     /* tabidx is now in range [0,..,64] */
2120     imask     = _mm_cmpgt_epi32(tabidx, i32);
2121     sswapsign = imask;
2122     cswapsign = imask;
2123     corridx   = _mm_and_si128(imask, i32);
2124     tabidx    = _mm_sub_epi32(tabidx, corridx);
2125
2126     /* tabidx is now in range [0..32] */
2127     imask     = _mm_cmpgt_epi32(tabidx, i16);
2128     cswapsign = _mm_xor_si128(cswapsign, imask);
2129     corridx   = _mm_sub_epi32(i32, tabidx);
2130     tabidx    = _mm_blendv_epi8(tabidx, corridx, imask);
2131     /* tabidx is now in range [0..16] */
2132     ssign     = _mm_cvtepi32_pd( _mm_or_si128( sswapsign, ione ) );
2133     csign     = _mm_cvtepi32_pd( _mm_or_si128( cswapsign, ione ) );
2134
2135 #ifdef _MSC_VER
2136     ypoint0  = _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx, 0));
2137     ypoint1  = _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx, 1));
2138 #else
2139     ypoint0  = sintable[_mm_extract_epi32(tabidx, 0)];
2140     ypoint1  = sintable[_mm_extract_epi32(tabidx, 1)];
2141 #endif
2142     sinpoint = _mm_unpackhi_pd(ypoint0, ypoint1);
2143     cospoint = _mm_unpacklo_pd(ypoint0, ypoint1);
2144
2145     sinpoint = _mm_mul_pd(sinpoint, ssign);
2146     cospoint = _mm_mul_pd(cospoint, csign);
2147
2148     z2       = _mm_mul_pd(z, z);
2149
2150     polySin  = _mm_mul_pd(sinP7, z2);
2151     polySin  = _mm_add_pd(polySin, sinP5);
2152     polySin  = _mm_mul_pd(polySin, z2);
2153     polySin  = _mm_add_pd(polySin, sinP3);
2154     polySin  = _mm_mul_pd(polySin, z2);
2155     polySin  = _mm_add_pd(polySin, sinP1);
2156     polySin  = _mm_mul_pd(polySin, z);
2157
2158     polyCos  = _mm_mul_pd(cosP6, z2);
2159     polyCos  = _mm_add_pd(polyCos, cosP4);
2160     polyCos  = _mm_mul_pd(polyCos, z2);
2161     polyCos  = _mm_add_pd(polyCos, cosP2);
2162     polyCos  = _mm_mul_pd(polyCos, z2);
2163     polyCos  = _mm_add_pd(polyCos, cosP0);
2164
2165     *sinval  = _mm_xor_pd(_mm_add_pd( _mm_mul_pd(sinpoint, polyCos), _mm_mul_pd(cospoint, polySin) ), xsign);
2166     *cosval  = _mm_sub_pd( _mm_mul_pd(cospoint, polyCos), _mm_mul_pd(sinpoint, polySin) );
2167
2168     return 0;
2169 }
2170
2171
2172 /*
2173  * IMPORTANT: Do NOT call both sin & cos if you need both results, since each of them
2174  * will then call the sincos() routine and waste a factor 2 in performance!
2175  */
2176 static __m256d
2177 gmx_mm256_sin_pd(__m256d x)
2178 {
2179     __m256d s, c;
2180     gmx_mm256_sincos_pd(x, &s, &c);
2181     return s;
2182 }
2183
2184 static __m128d
2185 gmx_mm_sin_pd(__m128d x)
2186 {
2187     __m128d s, c;
2188     gmx_mm_sincos_pd(x, &s, &c);
2189     return s;
2190 }
2191
2192 /*
2193  * IMPORTANT: Do NOT call both sin & cos if you need both results, since each of them
2194  * will then call the sincos() routine and waste a factor 2 in performance!
2195  */
2196 static __m256d
2197 gmx_mm256_cos_pd(__m256d x)
2198 {
2199     __m256d s, c;
2200     gmx_mm256_sincos_pd(x, &s, &c);
2201     return c;
2202 }
2203
2204 static __m128d
2205 gmx_mm_cos_pd(__m128d x)
2206 {
2207     __m128d s, c;
2208     gmx_mm_sincos_pd(x, &s, &c);
2209     return c;
2210 }
2211
2212
2213 static __m256d
2214 gmx_mm256_tan_pd(__m256d x)
2215 {
2216     __m256d sinval, cosval;
2217     __m256d tanval;
2218
2219     gmx_mm256_sincos_pd(x, &sinval, &cosval);
2220
2221     tanval = _mm256_mul_pd(sinval, gmx_mm256_inv_pd(cosval));
2222
2223     return tanval;
2224 }
2225
2226 static __m128d
2227 gmx_mm_tan_pd(__m128d x)
2228 {
2229     __m128d sinval, cosval;
2230     __m128d tanval;
2231
2232     gmx_mm_sincos_pd(x, &sinval, &cosval);
2233
2234     tanval = _mm_mul_pd(sinval, gmx_mm_inv_pd(cosval));
2235
2236     return tanval;
2237 }
2238
2239
2240 static __m256d
2241 gmx_mm256_asin_pd(__m256d x)
2242 {
2243     /* Same algorithm as cephes library */
2244     const __m256d signmask  = _mm256_castsi256_pd( _mm256_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF,
2245                                                                     0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
2246     const __m256d limit1    = _mm256_set1_pd(0.625);
2247     const __m256d limit2    = _mm256_set1_pd(1e-8);
2248     const __m256d one       = _mm256_set1_pd(1.0);
2249     const __m256d quarterpi = _mm256_set1_pd(M_PI/4.0);
2250     const __m256d morebits  = _mm256_set1_pd(6.123233995736765886130e-17);
2251
2252     const __m256d P5        = _mm256_set1_pd(4.253011369004428248960e-3);
2253     const __m256d P4        = _mm256_set1_pd(-6.019598008014123785661e-1);
2254     const __m256d P3        = _mm256_set1_pd(5.444622390564711410273e0);
2255     const __m256d P2        = _mm256_set1_pd(-1.626247967210700244449e1);
2256     const __m256d P1        = _mm256_set1_pd(1.956261983317594739197e1);
2257     const __m256d P0        = _mm256_set1_pd(-8.198089802484824371615e0);
2258
2259     const __m256d Q4        = _mm256_set1_pd(-1.474091372988853791896e1);
2260     const __m256d Q3        = _mm256_set1_pd(7.049610280856842141659e1);
2261     const __m256d Q2        = _mm256_set1_pd(-1.471791292232726029859e2);
2262     const __m256d Q1        = _mm256_set1_pd(1.395105614657485689735e2);
2263     const __m256d Q0        = _mm256_set1_pd(-4.918853881490881290097e1);
2264
2265     const __m256d R4        = _mm256_set1_pd(2.967721961301243206100e-3);
2266     const __m256d R3        = _mm256_set1_pd(-5.634242780008963776856e-1);
2267     const __m256d R2        = _mm256_set1_pd(6.968710824104713396794e0);
2268     const __m256d R1        = _mm256_set1_pd(-2.556901049652824852289e1);
2269     const __m256d R0        = _mm256_set1_pd(2.853665548261061424989e1);
2270
2271     const __m256d S3        = _mm256_set1_pd(-2.194779531642920639778e1);
2272     const __m256d S2        = _mm256_set1_pd(1.470656354026814941758e2);
2273     const __m256d S1        = _mm256_set1_pd(-3.838770957603691357202e2);
2274     const __m256d S0        = _mm256_set1_pd(3.424398657913078477438e2);
2275
2276     __m256d       sign;
2277     __m256d       mask;
2278     __m256d       xabs;
2279     __m256d       zz, ww, z, q, w, zz2, ww2;
2280     __m256d       PA, PB;
2281     __m256d       QA, QB;
2282     __m256d       RA, RB;
2283     __m256d       SA, SB;
2284     __m256d       nom, denom;
2285
2286     sign  = _mm256_andnot_pd(signmask, x);
2287     xabs  = _mm256_and_pd(x, signmask);
2288
2289     mask  = _mm256_cmp_pd(xabs, limit1, _CMP_GT_OQ);
2290
2291     zz    = _mm256_sub_pd(one, xabs);
2292     ww    = _mm256_mul_pd(xabs, xabs);
2293     zz2   = _mm256_mul_pd(zz, zz);
2294     ww2   = _mm256_mul_pd(ww, ww);
2295
2296     /* R */
2297     RA    = _mm256_mul_pd(R4, zz2);
2298     RB    = _mm256_mul_pd(R3, zz2);
2299     RA    = _mm256_add_pd(RA, R2);
2300     RB    = _mm256_add_pd(RB, R1);
2301     RA    = _mm256_mul_pd(RA, zz2);
2302     RB    = _mm256_mul_pd(RB, zz);
2303     RA    = _mm256_add_pd(RA, R0);
2304     RA    = _mm256_add_pd(RA, RB);
2305
2306     /* S, SA = zz2 */
2307     SB    = _mm256_mul_pd(S3, zz2);
2308     SA    = _mm256_add_pd(zz2, S2);
2309     SB    = _mm256_add_pd(SB, S1);
2310     SA    = _mm256_mul_pd(SA, zz2);
2311     SB    = _mm256_mul_pd(SB, zz);
2312     SA    = _mm256_add_pd(SA, S0);
2313     SA    = _mm256_add_pd(SA, SB);
2314
2315     /* P */
2316     PA    = _mm256_mul_pd(P5, ww2);
2317     PB    = _mm256_mul_pd(P4, ww2);
2318     PA    = _mm256_add_pd(PA, P3);
2319     PB    = _mm256_add_pd(PB, P2);
2320     PA    = _mm256_mul_pd(PA, ww2);
2321     PB    = _mm256_mul_pd(PB, ww2);
2322     PA    = _mm256_add_pd(PA, P1);
2323     PB    = _mm256_add_pd(PB, P0);
2324     PA    = _mm256_mul_pd(PA, ww);
2325     PA    = _mm256_add_pd(PA, PB);
2326
2327     /* Q, QA = ww2 */
2328     QB    = _mm256_mul_pd(Q4, ww2);
2329     QA    = _mm256_add_pd(ww2, Q3);
2330     QB    = _mm256_add_pd(QB, Q2);
2331     QA    = _mm256_mul_pd(QA, ww2);
2332     QB    = _mm256_mul_pd(QB, ww2);
2333     QA    = _mm256_add_pd(QA, Q1);
2334     QB    = _mm256_add_pd(QB, Q0);
2335     QA    = _mm256_mul_pd(QA, ww);
2336     QA    = _mm256_add_pd(QA, QB);
2337
2338     RA    = _mm256_mul_pd(RA, zz);
2339     PA    = _mm256_mul_pd(PA, ww);
2340
2341     nom   = _mm256_blendv_pd( PA, RA, mask );
2342     denom = _mm256_blendv_pd( QA, SA, mask );
2343
2344     q     = _mm256_mul_pd( nom, gmx_mm256_inv_pd(denom) );
2345
2346     zz    = _mm256_add_pd(zz, zz);
2347     zz    = gmx_mm256_sqrt_pd(zz);
2348     z     = _mm256_sub_pd(quarterpi, zz);
2349     zz    = _mm256_mul_pd(zz, q);
2350     zz    = _mm256_sub_pd(zz, morebits);
2351     z     = _mm256_sub_pd(z, zz);
2352     z     = _mm256_add_pd(z, quarterpi);
2353
2354     w     = _mm256_mul_pd(xabs, q);
2355     w     = _mm256_add_pd(w, xabs);
2356
2357     z     = _mm256_blendv_pd( w, z, mask );
2358
2359     mask  = _mm256_cmp_pd(xabs, limit2, _CMP_GT_OQ);
2360     z     = _mm256_blendv_pd( xabs, z, mask );
2361
2362     z = _mm256_xor_pd(z, sign);
2363
2364     return z;
2365 }
2366
2367 static __m128d
2368 gmx_mm_asin_pd(__m128d x)
2369 {
2370     /* Same algorithm as cephes library */
2371     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
2372     const __m128d limit1    = _mm_set1_pd(0.625);
2373     const __m128d limit2    = _mm_set1_pd(1e-8);
2374     const __m128d one       = _mm_set1_pd(1.0);
2375     const __m128d quarterpi = _mm_set1_pd(M_PI/4.0);
2376     const __m128d morebits  = _mm_set1_pd(6.123233995736765886130e-17);
2377
2378     const __m128d P5        = _mm_set1_pd(4.253011369004428248960e-3);
2379     const __m128d P4        = _mm_set1_pd(-6.019598008014123785661e-1);
2380     const __m128d P3        = _mm_set1_pd(5.444622390564711410273e0);
2381     const __m128d P2        = _mm_set1_pd(-1.626247967210700244449e1);
2382     const __m128d P1        = _mm_set1_pd(1.956261983317594739197e1);
2383     const __m128d P0        = _mm_set1_pd(-8.198089802484824371615e0);
2384
2385     const __m128d Q4        = _mm_set1_pd(-1.474091372988853791896e1);
2386     const __m128d Q3        = _mm_set1_pd(7.049610280856842141659e1);
2387     const __m128d Q2        = _mm_set1_pd(-1.471791292232726029859e2);
2388     const __m128d Q1        = _mm_set1_pd(1.395105614657485689735e2);
2389     const __m128d Q0        = _mm_set1_pd(-4.918853881490881290097e1);
2390
2391     const __m128d R4        = _mm_set1_pd(2.967721961301243206100e-3);
2392     const __m128d R3        = _mm_set1_pd(-5.634242780008963776856e-1);
2393     const __m128d R2        = _mm_set1_pd(6.968710824104713396794e0);
2394     const __m128d R1        = _mm_set1_pd(-2.556901049652824852289e1);
2395     const __m128d R0        = _mm_set1_pd(2.853665548261061424989e1);
2396
2397     const __m128d S3        = _mm_set1_pd(-2.194779531642920639778e1);
2398     const __m128d S2        = _mm_set1_pd(1.470656354026814941758e2);
2399     const __m128d S1        = _mm_set1_pd(-3.838770957603691357202e2);
2400     const __m128d S0        = _mm_set1_pd(3.424398657913078477438e2);
2401
2402     __m128d       sign;
2403     __m128d       mask;
2404     __m128d       xabs;
2405     __m128d       zz, ww, z, q, w, zz2, ww2;
2406     __m128d       PA, PB;
2407     __m128d       QA, QB;
2408     __m128d       RA, RB;
2409     __m128d       SA, SB;
2410     __m128d       nom, denom;
2411
2412     sign  = _mm_andnot_pd(signmask, x);
2413     xabs  = _mm_and_pd(x, signmask);
2414
2415     mask  = _mm_cmpgt_pd(xabs, limit1);
2416
2417     zz    = _mm_sub_pd(one, xabs);
2418     ww    = _mm_mul_pd(xabs, xabs);
2419     zz2   = _mm_mul_pd(zz, zz);
2420     ww2   = _mm_mul_pd(ww, ww);
2421
2422     /* R */
2423     RA    = _mm_mul_pd(R4, zz2);
2424     RB    = _mm_mul_pd(R3, zz2);
2425     RA    = _mm_add_pd(RA, R2);
2426     RB    = _mm_add_pd(RB, R1);
2427     RA    = _mm_mul_pd(RA, zz2);
2428     RB    = _mm_mul_pd(RB, zz);
2429     RA    = _mm_add_pd(RA, R0);
2430     RA    = _mm_add_pd(RA, RB);
2431
2432     /* S, SA = zz2 */
2433     SB    = _mm_mul_pd(S3, zz2);
2434     SA    = _mm_add_pd(zz2, S2);
2435     SB    = _mm_add_pd(SB, S1);
2436     SA    = _mm_mul_pd(SA, zz2);
2437     SB    = _mm_mul_pd(SB, zz);
2438     SA    = _mm_add_pd(SA, S0);
2439     SA    = _mm_add_pd(SA, SB);
2440
2441     /* P */
2442     PA    = _mm_mul_pd(P5, ww2);
2443     PB    = _mm_mul_pd(P4, ww2);
2444     PA    = _mm_add_pd(PA, P3);
2445     PB    = _mm_add_pd(PB, P2);
2446     PA    = _mm_mul_pd(PA, ww2);
2447     PB    = _mm_mul_pd(PB, ww2);
2448     PA    = _mm_add_pd(PA, P1);
2449     PB    = _mm_add_pd(PB, P0);
2450     PA    = _mm_mul_pd(PA, ww);
2451     PA    = _mm_add_pd(PA, PB);
2452
2453     /* Q, QA = ww2 */
2454     QB    = _mm_mul_pd(Q4, ww2);
2455     QA    = _mm_add_pd(ww2, Q3);
2456     QB    = _mm_add_pd(QB, Q2);
2457     QA    = _mm_mul_pd(QA, ww2);
2458     QB    = _mm_mul_pd(QB, ww2);
2459     QA    = _mm_add_pd(QA, Q1);
2460     QB    = _mm_add_pd(QB, Q0);
2461     QA    = _mm_mul_pd(QA, ww);
2462     QA    = _mm_add_pd(QA, QB);
2463
2464     RA    = _mm_mul_pd(RA, zz);
2465     PA    = _mm_mul_pd(PA, ww);
2466
2467     nom   = _mm_blendv_pd( PA, RA, mask );
2468     denom = _mm_blendv_pd( QA, SA, mask );
2469
2470     q     = _mm_mul_pd( nom, gmx_mm_inv_pd(denom) );
2471
2472     zz    = _mm_add_pd(zz, zz);
2473     zz    = gmx_mm_sqrt_pd(zz);
2474     z     = _mm_sub_pd(quarterpi, zz);
2475     zz    = _mm_mul_pd(zz, q);
2476     zz    = _mm_sub_pd(zz, morebits);
2477     z     = _mm_sub_pd(z, zz);
2478     z     = _mm_add_pd(z, quarterpi);
2479
2480     w     = _mm_mul_pd(xabs, q);
2481     w     = _mm_add_pd(w, xabs);
2482
2483     z     = _mm_blendv_pd( w, z, mask );
2484
2485     mask  = _mm_cmpgt_pd(xabs, limit2);
2486     z     = _mm_blendv_pd( xabs, z, mask );
2487
2488     z = _mm_xor_pd(z, sign);
2489
2490     return z;
2491 }
2492
2493
2494 static __m256d
2495 gmx_mm256_acos_pd(__m256d x)
2496 {
2497     const __m256d one        = _mm256_set1_pd(1.0);
2498     const __m256d half       = _mm256_set1_pd(0.5);
2499     const __m256d quarterpi0 = _mm256_set1_pd(7.85398163397448309616e-1);
2500     const __m256d quarterpi1 = _mm256_set1_pd(6.123233995736765886130e-17);
2501
2502
2503     __m256d mask1;
2504
2505     __m256d z, z1, z2;
2506
2507     mask1 = _mm256_cmp_pd(x, half, _CMP_GT_OQ);
2508     z1    = _mm256_mul_pd(half, _mm256_sub_pd(one, x));
2509     z1    = gmx_mm256_sqrt_pd(z1);
2510     z     = _mm256_blendv_pd( x, z1, mask1 );
2511
2512     z     = gmx_mm256_asin_pd(z);
2513
2514     z1    = _mm256_add_pd(z, z);
2515
2516     z2    = _mm256_sub_pd(quarterpi0, z);
2517     z2    = _mm256_add_pd(z2, quarterpi1);
2518     z2    = _mm256_add_pd(z2, quarterpi0);
2519
2520     z     = _mm256_blendv_pd(z2, z1, mask1);
2521
2522     return z;
2523 }
2524
2525 static __m128d
2526 gmx_mm_acos_pd(__m128d x)
2527 {
2528     const __m128d one        = _mm_set1_pd(1.0);
2529     const __m128d half       = _mm_set1_pd(0.5);
2530     const __m128d quarterpi0 = _mm_set1_pd(7.85398163397448309616e-1);
2531     const __m128d quarterpi1 = _mm_set1_pd(6.123233995736765886130e-17);
2532
2533
2534     __m128d mask1;
2535
2536     __m128d z, z1, z2;
2537
2538     mask1 = _mm_cmpgt_pd(x, half);
2539     z1    = _mm_mul_pd(half, _mm_sub_pd(one, x));
2540     z1    = gmx_mm_sqrt_pd(z1);
2541     z     = _mm_blendv_pd( x, z1, mask1 );
2542
2543     z     = gmx_mm_asin_pd(z);
2544
2545     z1    = _mm_add_pd(z, z);
2546
2547     z2    = _mm_sub_pd(quarterpi0, z);
2548     z2    = _mm_add_pd(z2, quarterpi1);
2549     z2    = _mm_add_pd(z2, quarterpi0);
2550
2551     z     = _mm_blendv_pd(z2, z1, mask1);
2552
2553     return z;
2554 }
2555
2556
2557 static __m256d
2558 gmx_mm256_atan_pd(__m256d x)
2559 {
2560     /* Same algorithm as cephes library */
2561     const __m256d signmask  = _mm256_castsi256_pd( _mm256_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF,
2562                                                                     0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
2563     const __m256d limit1    = _mm256_set1_pd(0.66);
2564     const __m256d limit2    = _mm256_set1_pd(2.41421356237309504880);
2565     const __m256d quarterpi = _mm256_set1_pd(M_PI/4.0);
2566     const __m256d halfpi    = _mm256_set1_pd(M_PI/2.0);
2567     const __m256d mone      = _mm256_set1_pd(-1.0);
2568     const __m256d morebits1 = _mm256_set1_pd(0.5*6.123233995736765886130E-17);
2569     const __m256d morebits2 = _mm256_set1_pd(6.123233995736765886130E-17);
2570
2571     const __m256d P4        = _mm256_set1_pd(-8.750608600031904122785E-1);
2572     const __m256d P3        = _mm256_set1_pd(-1.615753718733365076637E1);
2573     const __m256d P2        = _mm256_set1_pd(-7.500855792314704667340E1);
2574     const __m256d P1        = _mm256_set1_pd(-1.228866684490136173410E2);
2575     const __m256d P0        = _mm256_set1_pd(-6.485021904942025371773E1);
2576
2577     const __m256d Q4        = _mm256_set1_pd(2.485846490142306297962E1);
2578     const __m256d Q3        = _mm256_set1_pd(1.650270098316988542046E2);
2579     const __m256d Q2        = _mm256_set1_pd(4.328810604912902668951E2);
2580     const __m256d Q1        = _mm256_set1_pd(4.853903996359136964868E2);
2581     const __m256d Q0        = _mm256_set1_pd(1.945506571482613964425E2);
2582
2583     __m256d       sign;
2584     __m256d       mask1, mask2;
2585     __m256d       y, t1, t2;
2586     __m256d       z, z2;
2587     __m256d       P_A, P_B, Q_A, Q_B;
2588
2589     sign   = _mm256_andnot_pd(signmask, x);
2590     x      = _mm256_and_pd(x, signmask);
2591
2592     mask1  = _mm256_cmp_pd(x, limit1, _CMP_GT_OQ);
2593     mask2  = _mm256_cmp_pd(x, limit2, _CMP_GT_OQ);
2594
2595     t1     = _mm256_mul_pd(_mm256_add_pd(x, mone), gmx_mm256_inv_pd(_mm256_sub_pd(x, mone)));
2596     t2     = _mm256_mul_pd(mone, gmx_mm256_inv_pd(x));
2597
2598     y      = _mm256_and_pd(mask1, quarterpi);
2599     y      = _mm256_or_pd( _mm256_and_pd(mask2, halfpi), _mm256_andnot_pd(mask2, y) );
2600
2601     x      = _mm256_or_pd( _mm256_and_pd(mask1, t1), _mm256_andnot_pd(mask1, x) );
2602     x      = _mm256_or_pd( _mm256_and_pd(mask2, t2), _mm256_andnot_pd(mask2, x) );
2603
2604     z      = _mm256_mul_pd(x, x);
2605     z2     = _mm256_mul_pd(z, z);
2606
2607     P_A    = _mm256_mul_pd(P4, z2);
2608     P_B    = _mm256_mul_pd(P3, z2);
2609     P_A    = _mm256_add_pd(P_A, P2);
2610     P_B    = _mm256_add_pd(P_B, P1);
2611     P_A    = _mm256_mul_pd(P_A, z2);
2612     P_B    = _mm256_mul_pd(P_B, z);
2613     P_A    = _mm256_add_pd(P_A, P0);
2614     P_A    = _mm256_add_pd(P_A, P_B);
2615
2616     /* Q_A = z2 */
2617     Q_B    = _mm256_mul_pd(Q4, z2);
2618     Q_A    = _mm256_add_pd(z2, Q3);
2619     Q_B    = _mm256_add_pd(Q_B, Q2);
2620     Q_A    = _mm256_mul_pd(Q_A, z2);
2621     Q_B    = _mm256_mul_pd(Q_B, z2);
2622     Q_A    = _mm256_add_pd(Q_A, Q1);
2623     Q_B    = _mm256_add_pd(Q_B, Q0);
2624     Q_A    = _mm256_mul_pd(Q_A, z);
2625     Q_A    = _mm256_add_pd(Q_A, Q_B);
2626
2627     z      = _mm256_mul_pd(z, P_A);
2628     z      = _mm256_mul_pd(z, gmx_mm256_inv_pd(Q_A));
2629     z      = _mm256_mul_pd(z, x);
2630     z      = _mm256_add_pd(z, x);
2631
2632     t1     = _mm256_and_pd(mask1, morebits1);
2633     t1     = _mm256_or_pd( _mm256_and_pd(mask2, morebits2), _mm256_andnot_pd(mask2, t1) );
2634
2635     z      = _mm256_add_pd(z, t1);
2636     y      = _mm256_add_pd(y, z);
2637
2638     y      = _mm256_xor_pd(y, sign);
2639
2640     return y;
2641 }
2642
2643 static __m128d
2644 gmx_mm_atan_pd(__m128d x)
2645 {
2646     /* Same algorithm as cephes library */
2647     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
2648     const __m128d limit1    = _mm_set1_pd(0.66);
2649     const __m128d limit2    = _mm_set1_pd(2.41421356237309504880);
2650     const __m128d quarterpi = _mm_set1_pd(M_PI/4.0);
2651     const __m128d halfpi    = _mm_set1_pd(M_PI/2.0);
2652     const __m128d mone      = _mm_set1_pd(-1.0);
2653     const __m128d morebits1 = _mm_set1_pd(0.5*6.123233995736765886130E-17);
2654     const __m128d morebits2 = _mm_set1_pd(6.123233995736765886130E-17);
2655
2656     const __m128d P4        = _mm_set1_pd(-8.750608600031904122785E-1);
2657     const __m128d P3        = _mm_set1_pd(-1.615753718733365076637E1);
2658     const __m128d P2        = _mm_set1_pd(-7.500855792314704667340E1);
2659     const __m128d P1        = _mm_set1_pd(-1.228866684490136173410E2);
2660     const __m128d P0        = _mm_set1_pd(-6.485021904942025371773E1);
2661
2662     const __m128d Q4        = _mm_set1_pd(2.485846490142306297962E1);
2663     const __m128d Q3        = _mm_set1_pd(1.650270098316988542046E2);
2664     const __m128d Q2        = _mm_set1_pd(4.328810604912902668951E2);
2665     const __m128d Q1        = _mm_set1_pd(4.853903996359136964868E2);
2666     const __m128d Q0        = _mm_set1_pd(1.945506571482613964425E2);
2667
2668     __m128d       sign;
2669     __m128d       mask1, mask2;
2670     __m128d       y, t1, t2;
2671     __m128d       z, z2;
2672     __m128d       P_A, P_B, Q_A, Q_B;
2673
2674     sign   = _mm_andnot_pd(signmask, x);
2675     x      = _mm_and_pd(x, signmask);
2676
2677     mask1  = _mm_cmpgt_pd(x, limit1);
2678     mask2  = _mm_cmpgt_pd(x, limit2);
2679
2680     t1     = _mm_mul_pd(_mm_add_pd(x, mone), gmx_mm_inv_pd(_mm_sub_pd(x, mone)));
2681     t2     = _mm_mul_pd(mone, gmx_mm_inv_pd(x));
2682
2683     y      = _mm_and_pd(mask1, quarterpi);
2684     y      = _mm_or_pd( _mm_and_pd(mask2, halfpi), _mm_andnot_pd(mask2, y) );
2685
2686     x      = _mm_or_pd( _mm_and_pd(mask1, t1), _mm_andnot_pd(mask1, x) );
2687     x      = _mm_or_pd( _mm_and_pd(mask2, t2), _mm_andnot_pd(mask2, x) );
2688
2689     z      = _mm_mul_pd(x, x);
2690     z2     = _mm_mul_pd(z, z);
2691
2692     P_A    = _mm_mul_pd(P4, z2);
2693     P_B    = _mm_mul_pd(P3, z2);
2694     P_A    = _mm_add_pd(P_A, P2);
2695     P_B    = _mm_add_pd(P_B, P1);
2696     P_A    = _mm_mul_pd(P_A, z2);
2697     P_B    = _mm_mul_pd(P_B, z);
2698     P_A    = _mm_add_pd(P_A, P0);
2699     P_A    = _mm_add_pd(P_A, P_B);
2700
2701     /* Q_A = z2 */
2702     Q_B    = _mm_mul_pd(Q4, z2);
2703     Q_A    = _mm_add_pd(z2, Q3);
2704     Q_B    = _mm_add_pd(Q_B, Q2);
2705     Q_A    = _mm_mul_pd(Q_A, z2);
2706     Q_B    = _mm_mul_pd(Q_B, z2);
2707     Q_A    = _mm_add_pd(Q_A, Q1);
2708     Q_B    = _mm_add_pd(Q_B, Q0);
2709     Q_A    = _mm_mul_pd(Q_A, z);
2710     Q_A    = _mm_add_pd(Q_A, Q_B);
2711
2712     z      = _mm_mul_pd(z, P_A);
2713     z      = _mm_mul_pd(z, gmx_mm_inv_pd(Q_A));
2714     z      = _mm_mul_pd(z, x);
2715     z      = _mm_add_pd(z, x);
2716
2717     t1     = _mm_and_pd(mask1, morebits1);
2718     t1     = _mm_or_pd( _mm_and_pd(mask2, morebits2), _mm_andnot_pd(mask2, t1) );
2719
2720     z      = _mm_add_pd(z, t1);
2721     y      = _mm_add_pd(y, z);
2722
2723     y      = _mm_xor_pd(y, sign);
2724
2725     return y;
2726 }
2727
2728
2729
2730 static __m256d
2731 gmx_mm256_atan2_pd(__m256d y, __m256d x)
2732 {
2733     const __m256d pi          = _mm256_set1_pd(M_PI);
2734     const __m256d minuspi     = _mm256_set1_pd(-M_PI);
2735     const __m256d halfpi      = _mm256_set1_pd(M_PI/2.0);
2736     const __m256d minushalfpi = _mm256_set1_pd(-M_PI/2.0);
2737
2738     __m256d       z, z1, z3, z4;
2739     __m256d       w;
2740     __m256d       maskx_lt, maskx_eq;
2741     __m256d       masky_lt, masky_eq;
2742     __m256d       mask1, mask2, mask3, mask4, maskall;
2743
2744     maskx_lt  = _mm256_cmp_pd(x, _mm256_setzero_pd(), _CMP_LT_OQ);
2745     masky_lt  = _mm256_cmp_pd(y, _mm256_setzero_pd(), _CMP_LT_OQ);
2746     maskx_eq  = _mm256_cmp_pd(x, _mm256_setzero_pd(), _CMP_EQ_OQ);
2747     masky_eq  = _mm256_cmp_pd(y, _mm256_setzero_pd(), _CMP_EQ_OQ);
2748
2749     z         = _mm256_mul_pd(y, gmx_mm256_inv_pd(x));
2750     z         = gmx_mm256_atan_pd(z);
2751
2752     mask1     = _mm256_and_pd(maskx_eq, masky_lt);
2753     mask2     = _mm256_andnot_pd(maskx_lt, masky_eq);
2754     mask3     = _mm256_andnot_pd( _mm256_or_pd(masky_lt, masky_eq), maskx_eq);
2755     mask4     = _mm256_and_pd(masky_eq, maskx_lt);
2756
2757     maskall   = _mm256_or_pd( _mm256_or_pd(mask1, mask2), _mm256_or_pd(mask3, mask4) );
2758
2759     z         = _mm256_andnot_pd(maskall, z);
2760     z1        = _mm256_and_pd(mask1, minushalfpi);
2761     z3        = _mm256_and_pd(mask3, halfpi);
2762     z4        = _mm256_and_pd(mask4, pi);
2763
2764     z         = _mm256_or_pd( _mm256_or_pd(z, z1), _mm256_or_pd(z3, z4) );
2765
2766     w         = _mm256_blendv_pd(pi, minuspi, masky_lt);
2767     w         = _mm256_and_pd(w, maskx_lt);
2768
2769     w         = _mm256_andnot_pd(maskall, w);
2770
2771     z         = _mm256_add_pd(z, w);
2772
2773     return z;
2774 }
2775
2776 static __m128d
2777 gmx_mm_atan2_pd(__m128d y, __m128d x)
2778 {
2779     const __m128d pi          = _mm_set1_pd(M_PI);
2780     const __m128d minuspi     = _mm_set1_pd(-M_PI);
2781     const __m128d halfpi      = _mm_set1_pd(M_PI/2.0);
2782     const __m128d minushalfpi = _mm_set1_pd(-M_PI/2.0);
2783
2784     __m128d       z, z1, z3, z4;
2785     __m128d       w;
2786     __m128d       maskx_lt, maskx_eq;
2787     __m128d       masky_lt, masky_eq;
2788     __m128d       mask1, mask2, mask3, mask4, maskall;
2789
2790     maskx_lt  = _mm_cmplt_pd(x, _mm_setzero_pd());
2791     masky_lt  = _mm_cmplt_pd(y, _mm_setzero_pd());
2792     maskx_eq  = _mm_cmpeq_pd(x, _mm_setzero_pd());
2793     masky_eq  = _mm_cmpeq_pd(y, _mm_setzero_pd());
2794
2795     z         = _mm_mul_pd(y, gmx_mm_inv_pd(x));
2796     z         = gmx_mm_atan_pd(z);
2797
2798     mask1     = _mm_and_pd(maskx_eq, masky_lt);
2799     mask2     = _mm_andnot_pd(maskx_lt, masky_eq);
2800     mask3     = _mm_andnot_pd( _mm_or_pd(masky_lt, masky_eq), maskx_eq);
2801     mask4     = _mm_and_pd(masky_eq, maskx_lt);
2802
2803     maskall   = _mm_or_pd( _mm_or_pd(mask1, mask2), _mm_or_pd(mask3, mask4) );
2804
2805     z         = _mm_andnot_pd(maskall, z);
2806     z1        = _mm_and_pd(mask1, minushalfpi);
2807     z3        = _mm_and_pd(mask3, halfpi);
2808     z4        = _mm_and_pd(mask4, pi);
2809
2810     z         = _mm_or_pd( _mm_or_pd(z, z1), _mm_or_pd(z3, z4) );
2811
2812     w         = _mm_blendv_pd(pi, minuspi, masky_lt);
2813     w         = _mm_and_pd(w, maskx_lt);
2814
2815     w         = _mm_andnot_pd(maskall, w);
2816
2817     z         = _mm_add_pd(z, w);
2818
2819     return z;
2820 }
2821
2822 #endif /* _gmx_math_x86_avx_256_double_h_ */