include/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 "gmx_x86_avx_256.h"
  25
  26
  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. Add the result to 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     const __m128i signbit_epi32  = _mm_set1_epi32(0x80000000);
1925
1926     const __m256d tabscale      = _mm256_set1_pd(32.0/M_PI);
1927     const __m256d invtabscale0  = _mm256_set1_pd(9.81747508049011230469e-02);
1928     const __m256d invtabscale1  = _mm256_set1_pd(1.96197799156550576057e-08);
1929     const __m128i ione          = _mm_set1_epi32(1);
1930     const __m128i i32           = _mm_set1_epi32(32);
1931     const __m128i i16           = _mm_set1_epi32(16);
1932     const __m128i tabmask       = _mm_set1_epi32(0x3F);
1933     const __m256d sinP7         = _mm256_set1_pd(-1.0/5040.0);
1934     const __m256d sinP5         = _mm256_set1_pd(1.0/120.0);
1935     const __m256d sinP3         = _mm256_set1_pd(-1.0/6.0);
1936     const __m256d sinP1         = _mm256_set1_pd(1.0);
1937
1938     const __m256d cosP6         = _mm256_set1_pd(-1.0/720.0);
1939     const __m256d cosP4         = _mm256_set1_pd(1.0/24.0);
1940     const __m256d cosP2         = _mm256_set1_pd(-1.0/2.0);
1941     const __m256d cosP0         = _mm256_set1_pd(1.0);
1942
1943     __m256d scalex;
1944     __m128i tabidx,corridx;
1945     __m256d xabs,z,z2,polySin,polyCos;
1946     __m256d xpoint;
1947     __m256d t1,t2,t3,t4;
1948
1949     __m256d sinpoint,cospoint;
1950     __m256d xsign,ssign,csign;
1951     __m128i imask,sswapsign,cswapsign;
1952     __m256d minusone;
1953
1954     xsign    = _mm256_andnot_pd(signmask,x);
1955     xabs     = _mm256_and_pd(x,signmask);
1956
1957     scalex   = _mm256_mul_pd(tabscale,xabs);
1958     tabidx   = _mm256_cvtpd_epi32(scalex);
1959
1960     xpoint   = _mm256_round_pd(scalex,_MM_FROUND_TO_NEAREST_INT);
1961
1962     /* Extended precision arithmetics */
1963     z        = _mm256_sub_pd(xabs,_mm256_mul_pd(invtabscale0,xpoint));
1964     z        = _mm256_sub_pd(z,_mm256_mul_pd(invtabscale1,xpoint));
1965
1966     /* Range reduction to 0..2*Pi */
1967     tabidx   = _mm_and_si128(tabidx,tabmask);
1968
1969     /* tabidx is now in range [0,..,64] */
1970     imask     = _mm_cmpgt_epi32(tabidx,i32);
1971     sswapsign = imask;
1972     cswapsign = imask;
1973     corridx   = _mm_and_si128(imask,i32);
1974     tabidx    = _mm_sub_epi32(tabidx,corridx);
1975
1976     /* tabidx is now in range [0..32] */
1977     imask     = _mm_cmpgt_epi32(tabidx,i16);
1978     cswapsign = _mm_xor_si128(cswapsign,imask);
1979     corridx   = _mm_sub_epi32(i32,tabidx);
1980     tabidx    = _mm_blendv_epi8(tabidx,corridx,imask);
1981     /* tabidx is now in range [0..16] */
1982     ssign     = _mm256_cvtepi32_pd( _mm_or_si128( sswapsign , ione ) );
1983     csign     = _mm256_cvtepi32_pd( _mm_or_si128( cswapsign , ione ) );
1984
1985     /* First lookup into table */
1986 #ifdef _MSC_VER
1987     t1       = _mm256_insertf128_pd(_mm256_castpd128_pd256(_mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx,0))),
1988                                     _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx,2)), 0x1);
1989     t2       = _mm256_insertf128_pd(_mm256_castpd128_pd256(_mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx,1))),
1990                                     _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx,3)), 0x1);
1991 #else
1992     t1       = _mm256_insertf128_pd(_mm256_castpd128_pd256(sintable[_mm_extract_epi32(tabidx,0)]),
1993                                     sintable[_mm_extract_epi32(tabidx,2)], 0x1);
1994     t2       = _mm256_insertf128_pd(_mm256_castpd128_pd256(sintable[_mm_extract_epi32(tabidx,1)]),
1995                                     sintable[_mm_extract_epi32(tabidx,3)], 0x1);
1996 #endif
1997
1998     sinpoint  = _mm256_unpackhi_pd(t1,t2);
1999     cospoint  = _mm256_unpacklo_pd(t1,t2);
2000
2001     sinpoint = _mm256_mul_pd(sinpoint,ssign);
2002     cospoint = _mm256_mul_pd(cospoint,csign);
2003
2004     z2       = _mm256_mul_pd(z,z);
2005
2006     polySin  = _mm256_mul_pd(sinP7,z2);
2007     polySin  = _mm256_add_pd(polySin,sinP5);
2008     polySin  = _mm256_mul_pd(polySin,z2);
2009     polySin  = _mm256_add_pd(polySin,sinP3);
2010     polySin  = _mm256_mul_pd(polySin,z2);
2011     polySin  = _mm256_add_pd(polySin,sinP1);
2012     polySin  = _mm256_mul_pd(polySin,z);
2013
2014     polyCos  = _mm256_mul_pd(cosP6,z2);
2015     polyCos  = _mm256_add_pd(polyCos,cosP4);
2016     polyCos  = _mm256_mul_pd(polyCos,z2);
2017     polyCos  = _mm256_add_pd(polyCos,cosP2);
2018     polyCos  = _mm256_mul_pd(polyCos,z2);
2019     polyCos  = _mm256_add_pd(polyCos,cosP0);
2020
2021     *sinval  = _mm256_xor_pd(_mm256_add_pd( _mm256_mul_pd(sinpoint,polyCos) , _mm256_mul_pd(cospoint,polySin) ),xsign);
2022     *cosval  = _mm256_sub_pd( _mm256_mul_pd(cospoint,polyCos) , _mm256_mul_pd(sinpoint,polySin) );
2023
2024     return 0;
2025 }
2026
2027 static int
2028 gmx_mm_sincos_pd(__m128d x,
2029                  __m128d *sinval,
2030                  __m128d *cosval)
2031 {
2032 #ifdef _MSC_VER
2033     __declspec(align(16))
2034     const double sintable[34] =
2035     {
2036         1.00000000000000000e+00 , 0.00000000000000000e+00 ,
2037         9.95184726672196929e-01 , 9.80171403295606036e-02 ,
2038         9.80785280403230431e-01 , 1.95090322016128248e-01 ,
2039         9.56940335732208824e-01 , 2.90284677254462331e-01 ,
2040         9.23879532511286738e-01 , 3.82683432365089782e-01 ,
2041         8.81921264348355050e-01 , 4.71396736825997642e-01 ,
2042         8.31469612302545236e-01 , 5.55570233019602178e-01 ,
2043         7.73010453362736993e-01 , 6.34393284163645488e-01 ,
2044         7.07106781186547573e-01 , 7.07106781186547462e-01 ,
2045         6.34393284163645599e-01 , 7.73010453362736882e-01 ,
2046         5.55570233019602289e-01 , 8.31469612302545125e-01 ,
2047         4.71396736825997809e-01 , 8.81921264348354939e-01 ,
2048         3.82683432365089837e-01 , 9.23879532511286738e-01 ,
2049         2.90284677254462276e-01 , 9.56940335732208935e-01 ,
2050         1.95090322016128304e-01 , 9.80785280403230431e-01 ,
2051         9.80171403295607702e-02 , 9.95184726672196818e-01 ,
2052         0.0 , 1.00000000000000000e+00
2053     };
2054 #else
2055     const __m128d sintable[17] =
2056     {
2057         _mm_set_pd( 0.0 , 1.0 ),
2058         _mm_set_pd( sin(  1.0 * (M_PI/2.0) / 16.0) , cos(  1.0 * (M_PI/2.0) / 16.0) ),
2059         _mm_set_pd( sin(  2.0 * (M_PI/2.0) / 16.0) , cos(  2.0 * (M_PI/2.0) / 16.0) ),
2060         _mm_set_pd( sin(  3.0 * (M_PI/2.0) / 16.0) , cos(  3.0 * (M_PI/2.0) / 16.0) ),
2061         _mm_set_pd( sin(  4.0 * (M_PI/2.0) / 16.0) , cos(  4.0 * (M_PI/2.0) / 16.0) ),
2062         _mm_set_pd( sin(  5.0 * (M_PI/2.0) / 16.0) , cos(  5.0 * (M_PI/2.0) / 16.0) ),
2063         _mm_set_pd( sin(  6.0 * (M_PI/2.0) / 16.0) , cos(  6.0 * (M_PI/2.0) / 16.0) ),
2064         _mm_set_pd( sin(  7.0 * (M_PI/2.0) / 16.0) , cos(  7.0 * (M_PI/2.0) / 16.0) ),
2065         _mm_set_pd( sin(  8.0 * (M_PI/2.0) / 16.0) , cos(  8.0 * (M_PI/2.0) / 16.0) ),
2066         _mm_set_pd( sin(  9.0 * (M_PI/2.0) / 16.0) , cos(  9.0 * (M_PI/2.0) / 16.0) ),
2067         _mm_set_pd( sin( 10.0 * (M_PI/2.0) / 16.0) , cos( 10.0 * (M_PI/2.0) / 16.0) ),
2068         _mm_set_pd( sin( 11.0 * (M_PI/2.0) / 16.0) , cos( 11.0 * (M_PI/2.0) / 16.0) ),
2069         _mm_set_pd( sin( 12.0 * (M_PI/2.0) / 16.0) , cos( 12.0 * (M_PI/2.0) / 16.0) ),
2070         _mm_set_pd( sin( 13.0 * (M_PI/2.0) / 16.0) , cos( 13.0 * (M_PI/2.0) / 16.0) ),
2071         _mm_set_pd( sin( 14.0 * (M_PI/2.0) / 16.0) , cos( 14.0 * (M_PI/2.0) / 16.0) ),
2072         _mm_set_pd( sin( 15.0 * (M_PI/2.0) / 16.0) , cos( 15.0 * (M_PI/2.0) / 16.0) ),
2073         _mm_set_pd(  1.0 , 0.0 )
2074     };
2075 #endif
2076
2077     const __m128d signmask    = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
2078     const __m128i signbit_epi32  = _mm_set1_epi32(0x80000000);
2079
2080     const __m128d tabscale      = _mm_set1_pd(32.0/M_PI);
2081     const __m128d invtabscale0  = _mm_set1_pd(9.81747508049011230469e-02);
2082     const __m128d invtabscale1  = _mm_set1_pd(1.96197799156550576057e-08);
2083     const __m128i ione          = _mm_set1_epi32(1);
2084     const __m128i i32           = _mm_set1_epi32(32);
2085     const __m128i i16           = _mm_set1_epi32(16);
2086     const __m128i tabmask       = _mm_set1_epi32(0x3F);
2087     const __m128d sinP7         = _mm_set1_pd(-1.0/5040.0);
2088     const __m128d sinP5         = _mm_set1_pd(1.0/120.0);
2089     const __m128d sinP3         = _mm_set1_pd(-1.0/6.0);
2090     const __m128d sinP1         = _mm_set1_pd(1.0);
2091
2092     const __m128d cosP6         = _mm_set1_pd(-1.0/720.0);
2093     const __m128d cosP4         = _mm_set1_pd(1.0/24.0);
2094     const __m128d cosP2         = _mm_set1_pd(-1.0/2.0);
2095     const __m128d cosP0         = _mm_set1_pd(1.0);
2096
2097     __m128d scalex;
2098     __m128i tabidx,corridx;
2099     __m128d xabs,z,z2,polySin,polyCos;
2100     __m128d xpoint;
2101     __m128d ypoint0,ypoint1;
2102
2103     __m128d sinpoint,cospoint;
2104     __m128d xsign,ssign,csign;
2105     __m128i imask,sswapsign,cswapsign;
2106     __m128d minusone;
2107
2108     xsign    = _mm_andnot_pd(signmask,x);
2109     xabs     = _mm_and_pd(x,signmask);
2110
2111     scalex   = _mm_mul_pd(tabscale,xabs);
2112     tabidx   = _mm_cvtpd_epi32(scalex);
2113
2114     xpoint   = _mm_round_pd(scalex,_MM_FROUND_TO_NEAREST_INT);
2115
2116     /* Extended precision arithmetics */
2117     z        = _mm_sub_pd(xabs,_mm_mul_pd(invtabscale0,xpoint));
2118     z        = _mm_sub_pd(z,_mm_mul_pd(invtabscale1,xpoint));
2119
2120     /* Range reduction to 0..2*Pi */
2121     tabidx   = _mm_and_si128(tabidx,tabmask);
2122
2123     /* tabidx is now in range [0,..,64] */
2124     imask     = _mm_cmpgt_epi32(tabidx,i32);
2125     sswapsign = imask;
2126     cswapsign = imask;
2127     corridx   = _mm_and_si128(imask,i32);
2128     tabidx    = _mm_sub_epi32(tabidx,corridx);
2129
2130     /* tabidx is now in range [0..32] */
2131     imask     = _mm_cmpgt_epi32(tabidx,i16);
2132     cswapsign = _mm_xor_si128(cswapsign,imask);
2133     corridx   = _mm_sub_epi32(i32,tabidx);
2134     tabidx    = _mm_blendv_epi8(tabidx,corridx,imask);
2135     /* tabidx is now in range [0..16] */
2136     ssign     = _mm_cvtepi32_pd( _mm_or_si128( sswapsign , ione ) );
2137     csign     = _mm_cvtepi32_pd( _mm_or_si128( cswapsign , ione ) );
2138
2139 #ifdef _MSC_VER
2140     ypoint0  = _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx,0));
2141     ypoint1  = _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx,1));
2142 #else
2143     ypoint0  = sintable[_mm_extract_epi32(tabidx,0)];
2144     ypoint1  = sintable[_mm_extract_epi32(tabidx,1)];
2145 #endif
2146     sinpoint = _mm_unpackhi_pd(ypoint0,ypoint1);
2147     cospoint = _mm_unpacklo_pd(ypoint0,ypoint1);
2148
2149     sinpoint = _mm_mul_pd(sinpoint,ssign);
2150     cospoint = _mm_mul_pd(cospoint,csign);
2151
2152     z2       = _mm_mul_pd(z,z);
2153
2154     polySin  = _mm_mul_pd(sinP7,z2);
2155     polySin  = _mm_add_pd(polySin,sinP5);
2156     polySin  = _mm_mul_pd(polySin,z2);
2157     polySin  = _mm_add_pd(polySin,sinP3);
2158     polySin  = _mm_mul_pd(polySin,z2);
2159     polySin  = _mm_add_pd(polySin,sinP1);
2160     polySin  = _mm_mul_pd(polySin,z);
2161
2162     polyCos  = _mm_mul_pd(cosP6,z2);
2163     polyCos  = _mm_add_pd(polyCos,cosP4);
2164     polyCos  = _mm_mul_pd(polyCos,z2);
2165     polyCos  = _mm_add_pd(polyCos,cosP2);
2166     polyCos  = _mm_mul_pd(polyCos,z2);
2167     polyCos  = _mm_add_pd(polyCos,cosP0);
2168
2169     *sinval  = _mm_xor_pd(_mm_add_pd( _mm_mul_pd(sinpoint,polyCos) , _mm_mul_pd(cospoint,polySin) ),xsign);
2170     *cosval  = _mm_sub_pd( _mm_mul_pd(cospoint,polyCos) , _mm_mul_pd(sinpoint,polySin) );
2171
2172     return 0;
2173 }
2174
2175
2176 /*
2177  * IMPORTANT: Do NOT call both sin & cos if you need both results, since each of them
2178  * will then call the sincos() routine and waste a factor 2 in performance!
2179  */
2180 static __m256d
2181 gmx_mm256_sin_pd(__m256d x)
2182 {
2183     __m256d s,c;
2184     gmx_mm256_sincos_pd(x,&s,&c);
2185     return s;
2186 }
2187
2188 static __m128d
2189 gmx_mm_sin_pd(__m128d x)
2190 {
2191     __m128d s,c;
2192     gmx_mm_sincos_pd(x,&s,&c);
2193     return s;
2194 }
2195
2196 /*
2197  * IMPORTANT: Do NOT call both sin & cos if you need both results, since each of them
2198  * will then call the sincos() routine and waste a factor 2 in performance!
2199  */
2200 static __m256d
2201 gmx_mm256_cos_pd(__m256d x)
2202 {
2203     __m256d s,c;
2204     gmx_mm256_sincos_pd(x,&s,&c);
2205     return c;
2206 }
2207
2208 static __m128d
2209 gmx_mm_cos_pd(__m128d x)
2210 {
2211     __m128d s,c;
2212     gmx_mm_sincos_pd(x,&s,&c);
2213     return c;
2214 }
2215
2216
2217 static __m256d
2218 gmx_mm256_tan_pd(__m256d x)
2219 {
2220     __m256d sinval,cosval;
2221     __m256d tanval;
2222
2223     gmx_mm256_sincos_pd(x,&sinval,&cosval);
2224
2225     tanval = _mm256_mul_pd(sinval,gmx_mm256_inv_pd(cosval));
2226
2227     return tanval;
2228 }
2229
2230 static __m128d
2231 gmx_mm_tan_pd(__m128d x)
2232 {
2233     __m128d sinval,cosval;
2234     __m128d tanval;
2235
2236     gmx_mm_sincos_pd(x,&sinval,&cosval);
2237
2238     tanval = _mm_mul_pd(sinval,gmx_mm_inv_pd(cosval));
2239
2240     return tanval;
2241 }
2242
2243
2244 static __m256d
2245 gmx_mm256_asin_pd(__m256d x)
2246 {
2247     /* Same algorithm as cephes library */
2248     const __m256d signmask  = _mm256_castsi256_pd( _mm256_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF,
2249                                                                     0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
2250     const __m256d limit1    = _mm256_set1_pd(0.625);
2251     const __m256d limit2    = _mm256_set1_pd(1e-8);
2252     const __m256d one       = _mm256_set1_pd(1.0);
2253     const __m256d halfpi    = _mm256_set1_pd(M_PI/2.0);
2254     const __m256d quarterpi = _mm256_set1_pd(M_PI/4.0);
2255     const __m256d morebits  = _mm256_set1_pd(6.123233995736765886130e-17);
2256
2257     const __m256d P5        = _mm256_set1_pd(4.253011369004428248960e-3);
2258     const __m256d P4        = _mm256_set1_pd(-6.019598008014123785661e-1);
2259     const __m256d P3        = _mm256_set1_pd(5.444622390564711410273e0);
2260     const __m256d P2        = _mm256_set1_pd(-1.626247967210700244449e1);
2261     const __m256d P1        = _mm256_set1_pd(1.956261983317594739197e1);
2262     const __m256d P0        = _mm256_set1_pd(-8.198089802484824371615e0);
2263
2264     const __m256d Q4        = _mm256_set1_pd(-1.474091372988853791896e1);
2265     const __m256d Q3        = _mm256_set1_pd(7.049610280856842141659e1);
2266     const __m256d Q2        = _mm256_set1_pd(-1.471791292232726029859e2);
2267     const __m256d Q1        = _mm256_set1_pd(1.395105614657485689735e2);
2268     const __m256d Q0        = _mm256_set1_pd(-4.918853881490881290097e1);
2269
2270     const __m256d R4        = _mm256_set1_pd(2.967721961301243206100e-3);
2271     const __m256d R3        = _mm256_set1_pd(-5.634242780008963776856e-1);
2272     const __m256d R2        = _mm256_set1_pd(6.968710824104713396794e0);
2273     const __m256d R1        = _mm256_set1_pd(-2.556901049652824852289e1);
2274     const __m256d R0        = _mm256_set1_pd(2.853665548261061424989e1);
2275
2276     const __m256d S3        = _mm256_set1_pd(-2.194779531642920639778e1);
2277     const __m256d S2        = _mm256_set1_pd(1.470656354026814941758e2);
2278     const __m256d S1        = _mm256_set1_pd(-3.838770957603691357202e2);
2279     const __m256d S0        = _mm256_set1_pd(3.424398657913078477438e2);
2280
2281     __m256d sign;
2282     __m256d mask;
2283     __m256d xabs;
2284     __m256d zz,ww,z,q,w,y,zz2,ww2;
2285     __m256d PA,PB;
2286     __m256d QA,QB;
2287     __m256d RA,RB;
2288     __m256d SA,SB;
2289     __m256d nom,denom;
2290
2291     sign  = _mm256_andnot_pd(signmask,x);
2292     xabs  = _mm256_and_pd(x,signmask);
2293
2294     mask  = _mm256_cmp_pd(xabs,limit1,_CMP_GT_OQ);
2295
2296     zz    = _mm256_sub_pd(one,xabs);
2297     ww    = _mm256_mul_pd(xabs,xabs);
2298     zz2   = _mm256_mul_pd(zz,zz);
2299     ww2   = _mm256_mul_pd(ww,ww);
2300
2301     /* R */
2302     RA    = _mm256_mul_pd(R4,zz2);
2303     RB    = _mm256_mul_pd(R3,zz2);
2304     RA    = _mm256_add_pd(RA,R2);
2305     RB    = _mm256_add_pd(RB,R1);
2306     RA    = _mm256_mul_pd(RA,zz2);
2307     RB    = _mm256_mul_pd(RB,zz);
2308     RA    = _mm256_add_pd(RA,R0);
2309     RA    = _mm256_add_pd(RA,RB);
2310
2311     /* S, SA = zz2 */
2312     SB    = _mm256_mul_pd(S3,zz2);
2313     SA    = _mm256_add_pd(zz2,S2);
2314     SB    = _mm256_add_pd(SB,S1);
2315     SA    = _mm256_mul_pd(SA,zz2);
2316     SB    = _mm256_mul_pd(SB,zz);
2317     SA    = _mm256_add_pd(SA,S0);
2318     SA    = _mm256_add_pd(SA,SB);
2319
2320     /* P */
2321     PA    = _mm256_mul_pd(P5,ww2);
2322     PB    = _mm256_mul_pd(P4,ww2);
2323     PA    = _mm256_add_pd(PA,P3);
2324     PB    = _mm256_add_pd(PB,P2);
2325     PA    = _mm256_mul_pd(PA,ww2);
2326     PB    = _mm256_mul_pd(PB,ww2);
2327     PA    = _mm256_add_pd(PA,P1);
2328     PB    = _mm256_add_pd(PB,P0);
2329     PA    = _mm256_mul_pd(PA,ww);
2330     PA    = _mm256_add_pd(PA,PB);
2331
2332     /* Q, QA = ww2 */
2333     QB    = _mm256_mul_pd(Q4,ww2);
2334     QA    = _mm256_add_pd(ww2,Q3);
2335     QB    = _mm256_add_pd(QB,Q2);
2336     QA    = _mm256_mul_pd(QA,ww2);
2337     QB    = _mm256_mul_pd(QB,ww2);
2338     QA    = _mm256_add_pd(QA,Q1);
2339     QB    = _mm256_add_pd(QB,Q0);
2340     QA    = _mm256_mul_pd(QA,ww);
2341     QA    = _mm256_add_pd(QA,QB);
2342
2343     RA    = _mm256_mul_pd(RA,zz);
2344     PA    = _mm256_mul_pd(PA,ww);
2345
2346     nom   = _mm256_blendv_pd( PA,RA,mask );
2347     denom = _mm256_blendv_pd( QA,SA,mask );
2348
2349     q     = _mm256_mul_pd( nom, gmx_mm256_inv_pd(denom) );
2350
2351     zz    = _mm256_add_pd(zz,zz);
2352     zz    = gmx_mm256_sqrt_pd(zz);
2353     z     = _mm256_sub_pd(quarterpi,zz);
2354     zz    = _mm256_mul_pd(zz,q);
2355     zz    = _mm256_sub_pd(zz,morebits);
2356     z     = _mm256_sub_pd(z,zz);
2357     z     = _mm256_add_pd(z,quarterpi);
2358
2359     w     = _mm256_mul_pd(xabs,q);
2360     w     = _mm256_add_pd(w,xabs);
2361
2362     z     = _mm256_blendv_pd( w,z,mask );
2363
2364     mask  = _mm256_cmp_pd(xabs,limit2,_CMP_GT_OQ);
2365     z     = _mm256_blendv_pd( xabs,z,mask );
2366
2367     z = _mm256_xor_pd(z,sign);
2368
2369     return z;
2370 }
2371
2372 static __m128d
2373 gmx_mm_asin_pd(__m128d x)
2374 {
2375     /* Same algorithm as cephes library */
2376     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
2377     const __m128d limit1    = _mm_set1_pd(0.625);
2378     const __m128d limit2    = _mm_set1_pd(1e-8);
2379     const __m128d one       = _mm_set1_pd(1.0);
2380     const __m128d halfpi    = _mm_set1_pd(M_PI/2.0);
2381     const __m128d quarterpi = _mm_set1_pd(M_PI/4.0);
2382     const __m128d morebits  = _mm_set1_pd(6.123233995736765886130e-17);
2383
2384     const __m128d P5        = _mm_set1_pd(4.253011369004428248960e-3);
2385     const __m128d P4        = _mm_set1_pd(-6.019598008014123785661e-1);
2386     const __m128d P3        = _mm_set1_pd(5.444622390564711410273e0);
2387     const __m128d P2        = _mm_set1_pd(-1.626247967210700244449e1);
2388     const __m128d P1        = _mm_set1_pd(1.956261983317594739197e1);
2389     const __m128d P0        = _mm_set1_pd(-8.198089802484824371615e0);
2390
2391     const __m128d Q4        = _mm_set1_pd(-1.474091372988853791896e1);
2392     const __m128d Q3        = _mm_set1_pd(7.049610280856842141659e1);
2393     const __m128d Q2        = _mm_set1_pd(-1.471791292232726029859e2);
2394     const __m128d Q1        = _mm_set1_pd(1.395105614657485689735e2);
2395     const __m128d Q0        = _mm_set1_pd(-4.918853881490881290097e1);
2396
2397     const __m128d R4        = _mm_set1_pd(2.967721961301243206100e-3);
2398     const __m128d R3        = _mm_set1_pd(-5.634242780008963776856e-1);
2399     const __m128d R2        = _mm_set1_pd(6.968710824104713396794e0);
2400     const __m128d R1        = _mm_set1_pd(-2.556901049652824852289e1);
2401     const __m128d R0        = _mm_set1_pd(2.853665548261061424989e1);
2402
2403     const __m128d S3        = _mm_set1_pd(-2.194779531642920639778e1);
2404     const __m128d S2        = _mm_set1_pd(1.470656354026814941758e2);
2405     const __m128d S1        = _mm_set1_pd(-3.838770957603691357202e2);
2406     const __m128d S0        = _mm_set1_pd(3.424398657913078477438e2);
2407
2408     __m128d sign;
2409     __m128d mask;
2410     __m128d xabs;
2411     __m128d zz,ww,z,q,w,y,zz2,ww2;
2412     __m128d PA,PB;
2413     __m128d QA,QB;
2414     __m128d RA,RB;
2415     __m128d SA,SB;
2416     __m128d nom,denom;
2417
2418     sign  = _mm_andnot_pd(signmask,x);
2419     xabs  = _mm_and_pd(x,signmask);
2420
2421     mask  = _mm_cmpgt_pd(xabs,limit1);
2422
2423     zz    = _mm_sub_pd(one,xabs);
2424     ww    = _mm_mul_pd(xabs,xabs);
2425     zz2   = _mm_mul_pd(zz,zz);
2426     ww2   = _mm_mul_pd(ww,ww);
2427
2428     /* R */
2429     RA    = _mm_mul_pd(R4,zz2);
2430     RB    = _mm_mul_pd(R3,zz2);
2431     RA    = _mm_add_pd(RA,R2);
2432     RB    = _mm_add_pd(RB,R1);
2433     RA    = _mm_mul_pd(RA,zz2);
2434     RB    = _mm_mul_pd(RB,zz);
2435     RA    = _mm_add_pd(RA,R0);
2436     RA    = _mm_add_pd(RA,RB);
2437
2438     /* S, SA = zz2 */
2439     SB    = _mm_mul_pd(S3,zz2);
2440     SA    = _mm_add_pd(zz2,S2);
2441     SB    = _mm_add_pd(SB,S1);
2442     SA    = _mm_mul_pd(SA,zz2);
2443     SB    = _mm_mul_pd(SB,zz);
2444     SA    = _mm_add_pd(SA,S0);
2445     SA    = _mm_add_pd(SA,SB);
2446
2447     /* P */
2448     PA    = _mm_mul_pd(P5,ww2);
2449     PB    = _mm_mul_pd(P4,ww2);
2450     PA    = _mm_add_pd(PA,P3);
2451     PB    = _mm_add_pd(PB,P2);
2452     PA    = _mm_mul_pd(PA,ww2);
2453     PB    = _mm_mul_pd(PB,ww2);
2454     PA    = _mm_add_pd(PA,P1);
2455     PB    = _mm_add_pd(PB,P0);
2456     PA    = _mm_mul_pd(PA,ww);
2457     PA    = _mm_add_pd(PA,PB);
2458
2459     /* Q, QA = ww2 */
2460     QB    = _mm_mul_pd(Q4,ww2);
2461     QA    = _mm_add_pd(ww2,Q3);
2462     QB    = _mm_add_pd(QB,Q2);
2463     QA    = _mm_mul_pd(QA,ww2);
2464     QB    = _mm_mul_pd(QB,ww2);
2465     QA    = _mm_add_pd(QA,Q1);
2466     QB    = _mm_add_pd(QB,Q0);
2467     QA    = _mm_mul_pd(QA,ww);
2468     QA    = _mm_add_pd(QA,QB);
2469
2470     RA    = _mm_mul_pd(RA,zz);
2471     PA    = _mm_mul_pd(PA,ww);
2472
2473     nom   = _mm_blendv_pd( PA,RA,mask );
2474     denom = _mm_blendv_pd( QA,SA,mask );
2475
2476     q     = _mm_mul_pd( nom, gmx_mm_inv_pd(denom) );
2477
2478     zz    = _mm_add_pd(zz,zz);
2479     zz    = gmx_mm_sqrt_pd(zz);
2480     z     = _mm_sub_pd(quarterpi,zz);
2481     zz    = _mm_mul_pd(zz,q);
2482     zz    = _mm_sub_pd(zz,morebits);
2483     z     = _mm_sub_pd(z,zz);
2484     z     = _mm_add_pd(z,quarterpi);
2485
2486     w     = _mm_mul_pd(xabs,q);
2487     w     = _mm_add_pd(w,xabs);
2488
2489     z     = _mm_blendv_pd( w,z,mask );
2490
2491     mask  = _mm_cmpgt_pd(xabs,limit2);
2492     z     = _mm_blendv_pd( xabs,z,mask );
2493
2494     z = _mm_xor_pd(z,sign);
2495
2496     return z;
2497 }
2498
2499
2500 static __m256d
2501 gmx_mm256_acos_pd(__m256d x)
2502 {
2503     const __m256d signmask  = _mm256_castsi256_pd( _mm256_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF,
2504                                                                     0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
2505     const __m256d one        = _mm256_set1_pd(1.0);
2506     const __m256d half       = _mm256_set1_pd(0.5);
2507     const __m256d pi         = _mm256_set1_pd(M_PI);
2508     const __m256d quarterpi0 = _mm256_set1_pd(7.85398163397448309616e-1);
2509     const __m256d quarterpi1 = _mm256_set1_pd(6.123233995736765886130e-17);
2510
2511
2512     __m256d mask1;
2513
2514     __m256d z,z1,z2;
2515
2516     mask1 = _mm256_cmp_pd(x,half,_CMP_GT_OQ);
2517     z1    = _mm256_mul_pd(half,_mm256_sub_pd(one,x));
2518     z1    = gmx_mm256_sqrt_pd(z1);
2519     z     = _mm256_blendv_pd( x,z1,mask1 );
2520
2521     z     = gmx_mm256_asin_pd(z);
2522
2523     z1    = _mm256_add_pd(z,z);
2524
2525     z2    = _mm256_sub_pd(quarterpi0,z);
2526     z2    = _mm256_add_pd(z2,quarterpi1);
2527     z2    = _mm256_add_pd(z2,quarterpi0);
2528
2529     z     = _mm256_blendv_pd(z2,z1,mask1);
2530
2531     return z;
2532 }
2533
2534 static __m128d
2535 gmx_mm_acos_pd(__m128d x)
2536 {
2537     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
2538     const __m128d one        = _mm_set1_pd(1.0);
2539     const __m128d half       = _mm_set1_pd(0.5);
2540     const __m128d pi         = _mm_set1_pd(M_PI);
2541     const __m128d quarterpi0 = _mm_set1_pd(7.85398163397448309616e-1);
2542     const __m128d quarterpi1 = _mm_set1_pd(6.123233995736765886130e-17);
2543
2544
2545     __m128d mask1;
2546
2547     __m128d z,z1,z2;
2548
2549     mask1 = _mm_cmpgt_pd(x,half);
2550     z1    = _mm_mul_pd(half,_mm_sub_pd(one,x));
2551     z1    = gmx_mm_sqrt_pd(z1);
2552     z     = _mm_blendv_pd( x,z1,mask1 );
2553
2554     z     = gmx_mm_asin_pd(z);
2555
2556     z1    = _mm_add_pd(z,z);
2557
2558     z2    = _mm_sub_pd(quarterpi0,z);
2559     z2    = _mm_add_pd(z2,quarterpi1);
2560     z2    = _mm_add_pd(z2,quarterpi0);
2561
2562     z     = _mm_blendv_pd(z2,z1,mask1);
2563
2564     return z;
2565 }
2566
2567
2568 static __m256d
2569 gmx_mm256_atan_pd(__m256d x)
2570 {
2571     /* Same algorithm as cephes library */
2572     const __m256d signmask  = _mm256_castsi256_pd( _mm256_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF,
2573                                                                     0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
2574     const __m256d limit1    = _mm256_set1_pd(0.66);
2575     const __m256d limit2    = _mm256_set1_pd(2.41421356237309504880);
2576     const __m256d quarterpi = _mm256_set1_pd(M_PI/4.0);
2577     const __m256d halfpi    = _mm256_set1_pd(M_PI/2.0);
2578     const __m256d mone      = _mm256_set1_pd(-1.0);
2579     const __m256d morebits1 = _mm256_set1_pd(0.5*6.123233995736765886130E-17);
2580     const __m256d morebits2 = _mm256_set1_pd(6.123233995736765886130E-17);
2581
2582     const __m256d P4        = _mm256_set1_pd(-8.750608600031904122785E-1);
2583     const __m256d P3        = _mm256_set1_pd(-1.615753718733365076637E1);
2584     const __m256d P2        = _mm256_set1_pd(-7.500855792314704667340E1);
2585     const __m256d P1        = _mm256_set1_pd(-1.228866684490136173410E2);
2586     const __m256d P0        = _mm256_set1_pd(-6.485021904942025371773E1);
2587
2588     const __m256d Q4        = _mm256_set1_pd(2.485846490142306297962E1);
2589     const __m256d Q3        = _mm256_set1_pd(1.650270098316988542046E2);
2590     const __m256d Q2        = _mm256_set1_pd(4.328810604912902668951E2);
2591     const __m256d Q1        = _mm256_set1_pd(4.853903996359136964868E2);
2592     const __m256d Q0        = _mm256_set1_pd(1.945506571482613964425E2);
2593
2594     __m256d sign;
2595     __m256d mask1,mask2;
2596     __m256d y,t1,t2;
2597     __m256d z,z2;
2598     __m256d P_A,P_B,Q_A,Q_B;
2599
2600     sign   = _mm256_andnot_pd(signmask,x);
2601     x      = _mm256_and_pd(x,signmask);
2602
2603     mask1  = _mm256_cmp_pd(x,limit1,_CMP_GT_OQ);
2604     mask2  = _mm256_cmp_pd(x,limit2,_CMP_GT_OQ);
2605
2606     t1     = _mm256_mul_pd(_mm256_add_pd(x,mone),gmx_mm256_inv_pd(_mm256_sub_pd(x,mone)));
2607     t2     = _mm256_mul_pd(mone,gmx_mm256_inv_pd(x));
2608
2609     y      = _mm256_and_pd(mask1,quarterpi);
2610     y      = _mm256_or_pd( _mm256_and_pd(mask2,halfpi) , _mm256_andnot_pd(mask2,y) );
2611
2612     x      = _mm256_or_pd( _mm256_and_pd(mask1,t1) , _mm256_andnot_pd(mask1,x) );
2613     x      = _mm256_or_pd( _mm256_and_pd(mask2,t2) , _mm256_andnot_pd(mask2,x) );
2614
2615     z      = _mm256_mul_pd(x,x);
2616     z2     = _mm256_mul_pd(z,z);
2617
2618     P_A    = _mm256_mul_pd(P4,z2);
2619     P_B    = _mm256_mul_pd(P3,z2);
2620     P_A    = _mm256_add_pd(P_A,P2);
2621     P_B    = _mm256_add_pd(P_B,P1);
2622     P_A    = _mm256_mul_pd(P_A,z2);
2623     P_B    = _mm256_mul_pd(P_B,z);
2624     P_A    = _mm256_add_pd(P_A,P0);
2625     P_A    = _mm256_add_pd(P_A,P_B);
2626
2627     /* Q_A = z2 */
2628     Q_B    = _mm256_mul_pd(Q4,z2);
2629     Q_A    = _mm256_add_pd(z2,Q3);
2630     Q_B    = _mm256_add_pd(Q_B,Q2);
2631     Q_A    = _mm256_mul_pd(Q_A,z2);
2632     Q_B    = _mm256_mul_pd(Q_B,z2);
2633     Q_A    = _mm256_add_pd(Q_A,Q1);
2634     Q_B    = _mm256_add_pd(Q_B,Q0);
2635     Q_A    = _mm256_mul_pd(Q_A,z);
2636     Q_A    = _mm256_add_pd(Q_A,Q_B);
2637
2638     z      = _mm256_mul_pd(z,P_A);
2639     z      = _mm256_mul_pd(z,gmx_mm256_inv_pd(Q_A));
2640     z      = _mm256_mul_pd(z,x);
2641     z      = _mm256_add_pd(z,x);
2642
2643     t1     = _mm256_and_pd(mask1,morebits1);
2644     t1     = _mm256_or_pd( _mm256_and_pd(mask2,morebits2) , _mm256_andnot_pd(mask2,t1) );
2645
2646     z      = _mm256_add_pd(z,t1);
2647     y      = _mm256_add_pd(y,z);
2648
2649     y      = _mm256_xor_pd(y,sign);
2650
2651     return y;
2652 }
2653
2654 static __m128d
2655 gmx_mm_atan_pd(__m128d x)
2656 {
2657     /* Same algorithm as cephes library */
2658     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
2659     const __m128d limit1    = _mm_set1_pd(0.66);
2660     const __m128d limit2    = _mm_set1_pd(2.41421356237309504880);
2661     const __m128d quarterpi = _mm_set1_pd(M_PI/4.0);
2662     const __m128d halfpi    = _mm_set1_pd(M_PI/2.0);
2663     const __m128d mone      = _mm_set1_pd(-1.0);
2664     const __m128d morebits1 = _mm_set1_pd(0.5*6.123233995736765886130E-17);
2665     const __m128d morebits2 = _mm_set1_pd(6.123233995736765886130E-17);
2666
2667     const __m128d P4        = _mm_set1_pd(-8.750608600031904122785E-1);
2668     const __m128d P3        = _mm_set1_pd(-1.615753718733365076637E1);
2669     const __m128d P2        = _mm_set1_pd(-7.500855792314704667340E1);
2670     const __m128d P1        = _mm_set1_pd(-1.228866684490136173410E2);
2671     const __m128d P0        = _mm_set1_pd(-6.485021904942025371773E1);
2672
2673     const __m128d Q4        = _mm_set1_pd(2.485846490142306297962E1);
2674     const __m128d Q3        = _mm_set1_pd(1.650270098316988542046E2);
2675     const __m128d Q2        = _mm_set1_pd(4.328810604912902668951E2);
2676     const __m128d Q1        = _mm_set1_pd(4.853903996359136964868E2);
2677     const __m128d Q0        = _mm_set1_pd(1.945506571482613964425E2);
2678
2679     __m128d sign;
2680     __m128d mask1,mask2;
2681     __m128d y,t1,t2;
2682     __m128d z,z2;
2683     __m128d P_A,P_B,Q_A,Q_B;
2684
2685     sign   = _mm_andnot_pd(signmask,x);
2686     x      = _mm_and_pd(x,signmask);
2687
2688     mask1  = _mm_cmpgt_pd(x,limit1);
2689     mask2  = _mm_cmpgt_pd(x,limit2);
2690
2691     t1     = _mm_mul_pd(_mm_add_pd(x,mone),gmx_mm_inv_pd(_mm_sub_pd(x,mone)));
2692     t2     = _mm_mul_pd(mone,gmx_mm_inv_pd(x));
2693
2694     y      = _mm_and_pd(mask1,quarterpi);
2695     y      = _mm_or_pd( _mm_and_pd(mask2,halfpi) , _mm_andnot_pd(mask2,y) );
2696
2697     x      = _mm_or_pd( _mm_and_pd(mask1,t1) , _mm_andnot_pd(mask1,x) );
2698     x      = _mm_or_pd( _mm_and_pd(mask2,t2) , _mm_andnot_pd(mask2,x) );
2699
2700     z      = _mm_mul_pd(x,x);
2701     z2     = _mm_mul_pd(z,z);
2702
2703     P_A    = _mm_mul_pd(P4,z2);
2704     P_B    = _mm_mul_pd(P3,z2);
2705     P_A    = _mm_add_pd(P_A,P2);
2706     P_B    = _mm_add_pd(P_B,P1);
2707     P_A    = _mm_mul_pd(P_A,z2);
2708     P_B    = _mm_mul_pd(P_B,z);
2709     P_A    = _mm_add_pd(P_A,P0);
2710     P_A    = _mm_add_pd(P_A,P_B);
2711
2712     /* Q_A = z2 */
2713     Q_B    = _mm_mul_pd(Q4,z2);
2714     Q_A    = _mm_add_pd(z2,Q3);
2715     Q_B    = _mm_add_pd(Q_B,Q2);
2716     Q_A    = _mm_mul_pd(Q_A,z2);
2717     Q_B    = _mm_mul_pd(Q_B,z2);
2718     Q_A    = _mm_add_pd(Q_A,Q1);
2719     Q_B    = _mm_add_pd(Q_B,Q0);
2720     Q_A    = _mm_mul_pd(Q_A,z);
2721     Q_A    = _mm_add_pd(Q_A,Q_B);
2722
2723     z      = _mm_mul_pd(z,P_A);
2724     z      = _mm_mul_pd(z,gmx_mm_inv_pd(Q_A));
2725     z      = _mm_mul_pd(z,x);
2726     z      = _mm_add_pd(z,x);
2727
2728     t1     = _mm_and_pd(mask1,morebits1);
2729     t1     = _mm_or_pd( _mm_and_pd(mask2,morebits2) , _mm_andnot_pd(mask2,t1) );
2730
2731     z      = _mm_add_pd(z,t1);
2732     y      = _mm_add_pd(y,z);
2733
2734     y      = _mm_xor_pd(y,sign);
2735
2736     return y;
2737 }
2738
2739
2740
2741 static __m256d
2742 gmx_mm256_atan2_pd(__m256d y, __m256d x)
2743 {
2744     const __m256d pi          = _mm256_set1_pd(M_PI);
2745     const __m256d minuspi     = _mm256_set1_pd(-M_PI);
2746     const __m256d halfpi      = _mm256_set1_pd(M_PI/2.0);
2747     const __m256d minushalfpi = _mm256_set1_pd(-M_PI/2.0);
2748
2749     __m256d z,z1,z3,z4;
2750     __m256d w;
2751     __m256d maskx_lt,maskx_eq;
2752     __m256d masky_lt,masky_eq;
2753     __m256d mask1,mask2,mask3,mask4,maskall;
2754
2755     maskx_lt  = _mm256_cmp_pd(x,_mm256_setzero_pd(),_CMP_LT_OQ);
2756     masky_lt  = _mm256_cmp_pd(y,_mm256_setzero_pd(),_CMP_LT_OQ);
2757     maskx_eq  = _mm256_cmp_pd(x,_mm256_setzero_pd(),_CMP_EQ_OQ);
2758     masky_eq  = _mm256_cmp_pd(y,_mm256_setzero_pd(),_CMP_EQ_OQ);
2759
2760     z         = _mm256_mul_pd(y,gmx_mm256_inv_pd(x));
2761     z         = gmx_mm256_atan_pd(z);
2762
2763     mask1     = _mm256_and_pd(maskx_eq,masky_lt);
2764     mask2     = _mm256_andnot_pd(maskx_lt,masky_eq);
2765     mask3     = _mm256_andnot_pd( _mm256_or_pd(masky_lt,masky_eq) , maskx_eq);
2766     mask4     = _mm256_and_pd(masky_eq,maskx_lt);
2767
2768     maskall   = _mm256_or_pd( _mm256_or_pd(mask1,mask2), _mm256_or_pd(mask3,mask4) );
2769
2770     z         = _mm256_andnot_pd(maskall,z);
2771     z1        = _mm256_and_pd(mask1,minushalfpi);
2772     z3        = _mm256_and_pd(mask3,halfpi);
2773     z4        = _mm256_and_pd(mask4,pi);
2774
2775     z         = _mm256_or_pd( _mm256_or_pd(z,z1), _mm256_or_pd(z3,z4) );
2776
2777     w         = _mm256_blendv_pd(pi,minuspi,masky_lt);
2778     w         = _mm256_and_pd(w,maskx_lt);
2779
2780     w         = _mm256_andnot_pd(maskall,w);
2781
2782     z         = _mm256_add_pd(z,w);
2783
2784     return z;
2785 }
2786
2787 static __m128d
2788 gmx_mm_atan2_pd(__m128d y, __m128d x)
2789 {
2790     const __m128d pi          = _mm_set1_pd(M_PI);
2791     const __m128d minuspi     = _mm_set1_pd(-M_PI);
2792     const __m128d halfpi      = _mm_set1_pd(M_PI/2.0);
2793     const __m128d minushalfpi = _mm_set1_pd(-M_PI/2.0);
2794
2795     __m128d z,z1,z3,z4;
2796     __m128d w;
2797     __m128d maskx_lt,maskx_eq;
2798     __m128d masky_lt,masky_eq;
2799     __m128d mask1,mask2,mask3,mask4,maskall;
2800
2801     maskx_lt  = _mm_cmplt_pd(x,_mm_setzero_pd());
2802     masky_lt  = _mm_cmplt_pd(y,_mm_setzero_pd());
2803     maskx_eq  = _mm_cmpeq_pd(x,_mm_setzero_pd());
2804     masky_eq  = _mm_cmpeq_pd(y,_mm_setzero_pd());
2805
2806     z         = _mm_mul_pd(y,gmx_mm_inv_pd(x));
2807     z         = gmx_mm_atan_pd(z);
2808
2809     mask1     = _mm_and_pd(maskx_eq,masky_lt);
2810     mask2     = _mm_andnot_pd(maskx_lt,masky_eq);
2811     mask3     = _mm_andnot_pd( _mm_or_pd(masky_lt,masky_eq) , maskx_eq);
2812     mask4     = _mm_and_pd(masky_eq,maskx_lt);
2813
2814     maskall   = _mm_or_pd( _mm_or_pd(mask1,mask2), _mm_or_pd(mask3,mask4) );
2815
2816     z         = _mm_andnot_pd(maskall,z);
2817     z1        = _mm_and_pd(mask1,minushalfpi);
2818     z3        = _mm_and_pd(mask3,halfpi);
2819     z4        = _mm_and_pd(mask4,pi);
2820
2821     z         = _mm_or_pd( _mm_or_pd(z,z1), _mm_or_pd(z3,z4) );
2822
2823     w         = _mm_blendv_pd(pi,minuspi,masky_lt);
2824     w         = _mm_and_pd(w,maskx_lt);
2825
2826     w         = _mm_andnot_pd(maskall,w);
2827
2828     z         = _mm_add_pd(z,w);
2829
2830     return z;
2831 }
2832
2833 #endif /* _gmx_math_x86_avx_256_double_h_ */