src/gromacs/legacyheaders/gmx_math_x86_avx_128_fma_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_128_fma_double_h_
  22 #define _gmx_math_x86_avx_128_fma_double_h_
  23
  24 #include <immintrin.h> /* AVX */
  25 #ifdef HAVE_X86INTRIN_H
  26 #include <x86intrin.h> /* FMA */
  27 #endif
  28 #ifdef HAVE_INTRIN_H
  29 #include <intrin.h> /* FMA MSVC */
  30 #endif
  31
  32 #include <math.h>
  33
  34 #include "gmx_x86_avx_128_fma.h"
  35
  36
  37 #ifndef M_PI
  38 #  define M_PI 3.14159265358979323846264338327950288
  39 #endif
  40
  41
  42 /************************
  43  *                      *
  44  * Simple math routines *
  45  *                      *
  46  ************************/
  47
  48 /* 1.0/sqrt(x) */
  49 static gmx_inline __m128d
  50 gmx_mm_invsqrt_pd(__m128d x)
  51 {
  52     const __m128d half  = _mm_set1_pd(0.5);
  53     const __m128d three = _mm_set1_pd(3.0);
  54
  55     /* Lookup instruction only exists in single precision, convert back and forth... */
  56     __m128d lu = _mm_cvtps_pd(_mm_rsqrt_ps( _mm_cvtpd_ps(x)));
  57
  58     lu = _mm_mul_pd(_mm_mul_pd(half,lu),_mm_nmacc_pd(_mm_mul_pd(lu,lu),x,three));
  59     return _mm_mul_pd(_mm_mul_pd(half,lu),_mm_nmacc_pd(_mm_mul_pd(lu,lu),x,three));
  60 }
  61
  62 /* 1.0/sqrt(x), done for a pair of arguments to improve throughput */
  63 static void
  64 gmx_mm_invsqrt_pair_pd(__m128d x1, __m128d x2, __m128d *invsqrt1, __m128d *invsqrt2)
  65 {
  66     const __m128d half   = _mm_set1_pd(0.5);
  67     const __m128d three  = _mm_set1_pd(3.0);
  68     const __m128  halff  = _mm_set1_ps(0.5f);
  69     const __m128  threef = _mm_set1_ps(3.0f);
  70
  71     __m128 xf,luf;
  72     __m128d lu1,lu2;
  73
  74     /* Do first N-R step in float for 2x throughput */
  75     xf  = _mm_shuffle_ps(_mm_cvtpd_ps(x1),_mm_cvtpd_ps(x2),_MM_SHUFFLE(1,0,1,0));
  76     luf = _mm_rsqrt_ps(xf);
  77
  78     luf = _mm_mul_ps(_mm_mul_ps(halff,luf),_mm_nmacc_ps(_mm_mul_ps(luf,luf),xf,threef));
  79
  80
  81     lu2 = _mm_cvtps_pd(_mm_shuffle_ps(luf,luf,_MM_SHUFFLE(3,2,3,2)));
  82     lu1 = _mm_cvtps_pd(luf);
  83
  84     *invsqrt1 = _mm_mul_pd(_mm_mul_pd(half,lu1),_mm_nmacc_pd(_mm_mul_pd(lu1,lu1),x1,three));
  85     *invsqrt2 = _mm_mul_pd(_mm_mul_pd(half,lu2),_mm_nmacc_pd(_mm_mul_pd(lu2,lu2),x2,three));
  86 }
  87
  88 /* sqrt(x) - Do NOT use this (but rather invsqrt) if you actually need 1.0/sqrt(x) */
  89 static gmx_inline __m128d
  90 gmx_mm_sqrt_pd(__m128d x)
  91 {
  92     __m128d mask;
  93     __m128d res;
  94
  95     mask = _mm_cmpeq_pd(x,_mm_setzero_pd());
  96     res  = _mm_andnot_pd(mask,gmx_mm_invsqrt_pd(x));
  97
  98     res  = _mm_mul_pd(x,res);
  99
 100     return res;
 101 }
 102
 103 /* 1.0/x */
 104 static gmx_inline __m128d
 105 gmx_mm_inv_pd(__m128d x)
 106 {
 107     const __m128d two  = _mm_set1_pd(2.0);
 108
 109     /* Lookup instruction only exists in single precision, convert back and forth... */
 110     __m128d lu = _mm_cvtps_pd(_mm_rcp_ps( _mm_cvtpd_ps(x)));
 111
 112     /* Perform two N-R steps for double precision */
 113     lu         = _mm_mul_pd(lu,_mm_nmacc_pd(lu,x,two));
 114     return _mm_mul_pd(lu,_mm_nmacc_pd(lu,x,two));
 115 }
 116
 117 static gmx_inline __m128d
 118 gmx_mm_abs_pd(__m128d x)
 119 {
 120     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF, 0xFFFFFFFF, 0x7FFFFFFF, 0xFFFFFFFF) );
 121
 122     return _mm_and_pd(x,signmask);
 123 }
 124
 125
 126 /*
 127  * 2^x function.
 128  *
 129  * The 2^w term is calculated from a (6,0)-th order (no denominator) Minimax polynomia on the interval
 130  * [-0.5,0.5].
 131  *
 132  * The approximation on [-0.5,0.5] is a rational Padé approximation, 1+2*P(x^2)/(Q(x^2)-P(x^2)),
 133  * according to the same algorithm as used in the Cephes/netlib math routines.
 134  */
 135 static __m128d
 136 gmx_mm_exp2_pd(__m128d x)
 137 {
 138     /* Lower bound: We do not allow numbers that would lead to an IEEE fp representation exponent smaller than -126. */
 139     const __m128d arglimit = _mm_set1_pd(1022.0);
 140     const __m128i expbase  = _mm_set1_epi32(1023);
 141
 142     const __m128d P2       = _mm_set1_pd(2.30933477057345225087e-2);
 143     const __m128d P1       = _mm_set1_pd(2.02020656693165307700e1);
 144     const __m128d P0       = _mm_set1_pd(1.51390680115615096133e3);
 145     /* Q2 == 1.0 */
 146     const __m128d Q1       = _mm_set1_pd(2.33184211722314911771e2);
 147     const __m128d Q0       = _mm_set1_pd(4.36821166879210612817e3);
 148     const __m128d one      = _mm_set1_pd(1.0);
 149     const __m128d two      = _mm_set1_pd(2.0);
 150
 151     __m128d valuemask;
 152     __m128i iexppart;
 153     __m128d fexppart;
 154     __m128d intpart;
 155     __m128d z,z2;
 156     __m128d PolyP,PolyQ;
 157
 158     iexppart  = _mm_cvtpd_epi32(x);
 159     intpart   = _mm_round_pd(x,_MM_FROUND_TO_NEAREST_INT);
 160
 161     /* The two lowest elements of iexppart now contains 32-bit numbers with a correctly biased exponent.
 162      * To be able to shift it into the exponent for a double precision number we first need to
 163      * shuffle so that the lower half contains the first element, and the upper half the second.
 164      * This should really be done as a zero-extension, but since the next instructions will shift
 165      * the registers left by 52 bits it doesn't matter what we put there - it will be shifted out.
 166      * (thus we just use element 2 from iexppart).
 167      */
 168     iexppart  = _mm_shuffle_epi32(iexppart,_MM_SHUFFLE(2,1,2,0));
 169
 170     /* Do the shift operation on the 64-bit registers */
 171     iexppart  = _mm_add_epi32(iexppart,expbase);
 172     iexppart  = _mm_slli_epi64(iexppart,52);
 173
 174     valuemask = _mm_cmpge_pd(arglimit,gmx_mm_abs_pd(x));
 175     fexppart  = _mm_and_pd(valuemask,gmx_mm_castsi128_pd(iexppart));
 176
 177     z         = _mm_sub_pd(x,intpart);
 178     z2        = _mm_mul_pd(z,z);
 179
 180     PolyP     = _mm_macc_pd(P2,z2,P1);
 181     PolyQ     = _mm_add_pd(z2,Q1);
 182     PolyP     = _mm_macc_pd(PolyP,z2,P0);
 183     PolyQ     = _mm_macc_pd(PolyQ,z2,Q0);
 184     PolyP     = _mm_mul_pd(PolyP,z);
 185
 186     z         = _mm_mul_pd(PolyP,gmx_mm_inv_pd(_mm_sub_pd(PolyQ,PolyP)));
 187     z         = _mm_macc_pd(two,z,one);
 188
 189     z         = _mm_mul_pd(z,fexppart);
 190
 191     return  z;
 192 }
 193
 194 /* Exponential function. This could be calculated from 2^x as Exp(x)=2^(y), where y=log2(e)*x,
 195  * but there will then be a small rounding error since we lose some precision due to the
 196  * multiplication. This will then be magnified a lot by the exponential.
 197  *
 198  * Instead, we calculate the fractional part directly as a Padé approximation of
 199  * Exp(z) on [-0.5,0.5]. We use extended precision arithmetics to calculate the fraction
 200  * remaining after 2^y, which avoids the precision-loss.
 201  */
 202 static __m128d
 203 gmx_mm_exp_pd(__m128d exparg)
 204 {
 205     const __m128d argscale = _mm_set1_pd(1.4426950408889634073599);
 206     /* Lower bound: We do not allow numbers that would lead to an IEEE fp representation exponent smaller than -126. */
 207     const __m128d arglimit = _mm_set1_pd(1022.0);
 208     const __m128i expbase  = _mm_set1_epi32(1023);
 209
 210     const __m128d invargscale0  = _mm_set1_pd(6.93145751953125e-1);
 211     const __m128d invargscale1  = _mm_set1_pd(1.42860682030941723212e-6);
 212
 213     const __m128d P2       = _mm_set1_pd(1.26177193074810590878e-4);
 214     const __m128d P1       = _mm_set1_pd(3.02994407707441961300e-2);
 215     /* P0 == 1.0 */
 216     const __m128d Q3       = _mm_set1_pd(3.00198505138664455042E-6);
 217     const __m128d Q2       = _mm_set1_pd(2.52448340349684104192E-3);
 218     const __m128d Q1       = _mm_set1_pd(2.27265548208155028766E-1);
 219     /* Q0 == 2.0 */
 220     const __m128d one      = _mm_set1_pd(1.0);
 221     const __m128d two      = _mm_set1_pd(2.0);
 222
 223     __m128d valuemask;
 224     __m128i iexppart;
 225     __m128d fexppart;
 226     __m128d intpart;
 227     __m128d x,z,z2;
 228     __m128d PolyP,PolyQ;
 229
 230     x             = _mm_mul_pd(exparg,argscale);
 231
 232     iexppart  = _mm_cvtpd_epi32(x);
 233     intpart   = _mm_round_pd(x,_MM_FROUND_TO_NEAREST_INT);
 234
 235     /* The two lowest elements of iexppart now contains 32-bit numbers with a correctly biased exponent.
 236      * To be able to shift it into the exponent for a double precision number we first need to
 237      * shuffle so that the lower half contains the first element, and the upper half the second.
 238      * This should really be done as a zero-extension, but since the next instructions will shift
 239      * the registers left by 52 bits it doesn't matter what we put there - it will be shifted out.
 240      * (thus we just use element 2 from iexppart).
 241      */
 242     iexppart  = _mm_shuffle_epi32(iexppart,_MM_SHUFFLE(2,1,2,0));
 243
 244     /* Do the shift operation on the 64-bit registers */
 245     iexppart  = _mm_add_epi32(iexppart,expbase);
 246     iexppart  = _mm_slli_epi64(iexppart,52);
 247
 248     valuemask = _mm_cmpge_pd(arglimit,gmx_mm_abs_pd(x));
 249     fexppart  = _mm_and_pd(valuemask,gmx_mm_castsi128_pd(iexppart));
 250
 251     z         = _mm_sub_pd(exparg,_mm_mul_pd(invargscale0,intpart));
 252     z         = _mm_sub_pd(z,_mm_mul_pd(invargscale1,intpart));
 253
 254     z2        = _mm_mul_pd(z,z);
 255
 256     PolyQ     = _mm_macc_pd(Q3,z2,Q2);
 257     PolyP     = _mm_macc_pd(P2,z2,P1);
 258     PolyQ     = _mm_macc_pd(PolyQ,z2,Q1);
 259
 260     PolyP     = _mm_macc_pd(PolyP,z2,one);
 261     PolyQ     = _mm_macc_pd(PolyQ,z2,two);
 262
 263     PolyP     = _mm_mul_pd(PolyP,z);
 264
 265     z         = _mm_mul_pd(PolyP,gmx_mm_inv_pd(_mm_sub_pd(PolyQ,PolyP)));
 266     z         = _mm_macc_pd(two,z,one);
 267
 268     z         = _mm_mul_pd(z,fexppart);
 269
 270     return  z;
 271 }
 272
 273
 274
 275 static __m128d
 276 gmx_mm_log_pd(__m128d x)
 277 {
 278     /* Same algorithm as cephes library */
 279     const __m128d expmask    = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FF00000, 0x00000000, 0x7FF00000, 0x00000000) );
 280
 281     const __m128i expbase_m1 = _mm_set1_epi32(1023-1); /* We want non-IEEE format */
 282
 283     const __m128d half       = _mm_set1_pd(0.5);
 284     const __m128d one        = _mm_set1_pd(1.0);
 285     const __m128d two        = _mm_set1_pd(2.0);
 286     const __m128d invsq2     = _mm_set1_pd(1.0/sqrt(2.0));
 287
 288     const __m128d corr1      = _mm_set1_pd(-2.121944400546905827679e-4);
 289     const __m128d corr2      = _mm_set1_pd(0.693359375);
 290
 291     const __m128d P5         = _mm_set1_pd(1.01875663804580931796e-4);
 292     const __m128d P4         = _mm_set1_pd(4.97494994976747001425e-1);
 293     const __m128d P3         = _mm_set1_pd(4.70579119878881725854e0);
 294     const __m128d P2         = _mm_set1_pd(1.44989225341610930846e1);
 295     const __m128d P1         = _mm_set1_pd(1.79368678507819816313e1);
 296     const __m128d P0         = _mm_set1_pd(7.70838733755885391666e0);
 297
 298     const __m128d Q4         = _mm_set1_pd(1.12873587189167450590e1);
 299     const __m128d Q3         = _mm_set1_pd(4.52279145837532221105e1);
 300     const __m128d Q2         = _mm_set1_pd(8.29875266912776603211e1);
 301     const __m128d Q1         = _mm_set1_pd(7.11544750618563894466e1);
 302     const __m128d Q0         = _mm_set1_pd(2.31251620126765340583e1);
 303
 304     const __m128d R2         = _mm_set1_pd(-7.89580278884799154124e-1);
 305     const __m128d R1         = _mm_set1_pd(1.63866645699558079767e1);
 306     const __m128d R0         = _mm_set1_pd(-6.41409952958715622951e1);
 307
 308     const __m128d S2         = _mm_set1_pd(-3.56722798256324312549E1);
 309     const __m128d S1         = _mm_set1_pd(3.12093766372244180303E2);
 310     const __m128d S0         = _mm_set1_pd(-7.69691943550460008604E2);
 311
 312     __m128d fexp;
 313     __m128i iexp;
 314
 315     __m128d mask1,mask2;
 316     __m128d corr,t1,t2,q;
 317     __m128d zA,yA,xA,zB,yB,xB,z;
 318     __m128d polyR,polyS;
 319     __m128d polyP1,polyP2,polyQ1,polyQ2;
 320
 321     /* Separate x into exponent and mantissa, with a mantissa in the range [0.5..1[ (not IEEE754 standard!) */
 322     fexp   = _mm_and_pd(x,expmask);
 323     iexp   = gmx_mm_castpd_si128(fexp);
 324     iexp   = _mm_srli_epi64(iexp,52);
 325     iexp   = _mm_sub_epi32(iexp,expbase_m1);
 326     iexp   = _mm_shuffle_epi32(iexp, _MM_SHUFFLE(1,1,2,0) );
 327     fexp   = _mm_cvtepi32_pd(iexp);
 328
 329     x      = _mm_andnot_pd(expmask,x);
 330     x      = _mm_or_pd(x,one);
 331     x      = _mm_mul_pd(x,half);
 332
 333     mask1     = _mm_cmpgt_pd(gmx_mm_abs_pd(fexp),two);
 334     mask2     = _mm_cmplt_pd(x,invsq2);
 335
 336     fexp   = _mm_sub_pd(fexp,_mm_and_pd(mask2,one));
 337
 338     /* If mask1 is set ('A') */
 339     zA     = _mm_sub_pd(x,half);
 340     t1     = _mm_blendv_pd( zA,x,mask2 );
 341     zA     = _mm_sub_pd(t1,half);
 342     t2     = _mm_blendv_pd( x,zA,mask2 );
 343     yA     = _mm_mul_pd(half,_mm_add_pd(t2,one));
 344
 345     xA     = _mm_mul_pd(zA,gmx_mm_inv_pd(yA));
 346     zA     = _mm_mul_pd(xA,xA);
 347
 348     /* EVALUATE POLY */
 349     polyR  = _mm_macc_pd(R2,zA,R1);
 350     polyR  = _mm_macc_pd(polyR,zA,R0);
 351
 352     polyS  = _mm_add_pd(zA,S2);
 353     polyS  = _mm_macc_pd(polyS,zA,S1);
 354     polyS  = _mm_macc_pd(polyS,zA,S0);
 355
 356     q      = _mm_mul_pd(polyR,gmx_mm_inv_pd(polyS));
 357     zA     = _mm_mul_pd(_mm_mul_pd(xA,zA),q);
 358
 359     zA     = _mm_macc_pd(corr1,fexp,zA);
 360     zA     = _mm_add_pd(zA,xA);
 361     zA     = _mm_macc_pd(corr2,fexp,zA);
 362
 363     /* If mask1 is not set ('B') */
 364     corr   = _mm_and_pd(mask2,x);
 365     xB     = _mm_add_pd(x,corr);
 366     xB     = _mm_sub_pd(xB,one);
 367     zB     = _mm_mul_pd(xB,xB);
 368
 369     polyP1 = _mm_macc_pd(P5,zB,P3);
 370     polyP2 = _mm_macc_pd(P4,zB,P2);
 371     polyP1 = _mm_macc_pd(polyP1,zB,P1);
 372     polyP2 = _mm_macc_pd(polyP2,zB,P0);
 373     polyP1 = _mm_macc_pd(polyP1,xB,polyP2);
 374
 375     polyQ2 = _mm_macc_pd(Q4,zB,Q2);
 376     polyQ1 = _mm_add_pd(zB,Q3);
 377     polyQ1 = _mm_macc_pd(polyQ1,zB,Q1);
 378     polyQ2 = _mm_macc_pd(polyQ2,zB,Q0);
 379     polyQ1 = _mm_macc_pd(polyQ1,xB,polyQ2);
 380
 381     fexp   = _mm_and_pd(fexp,_mm_cmpneq_pd(fexp,_mm_setzero_pd()));
 382
 383     q      = _mm_mul_pd(polyP1,gmx_mm_inv_pd(polyQ1));
 384     yB     = _mm_macc_pd(_mm_mul_pd(xB,zB),q,_mm_mul_pd(corr1,fexp));
 385
 386     yB     = _mm_nmacc_pd(half,zB,yB);
 387     zB     = _mm_add_pd(xB,yB);
 388     zB     = _mm_macc_pd(corr2,fexp,zB);
 389
 390     z      = _mm_blendv_pd( zB,zA,mask1 );
 391
 392     return z;
 393 }
 394
 395
 396
 397 static __m128d
 398 gmx_mm_erf_pd(__m128d x)
 399 {
 400     /* Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75] */
 401     const __m128d CAP4      = _mm_set1_pd(-0.431780540597889301512e-4);
 402     const __m128d CAP3      = _mm_set1_pd(-0.00578562306260059236059);
 403     const __m128d CAP2      = _mm_set1_pd(-0.028593586920219752446);
 404     const __m128d CAP1      = _mm_set1_pd(-0.315924962948621698209);
 405     const __m128d CAP0      = _mm_set1_pd(0.14952975608477029151);
 406
 407     const __m128d CAQ5      = _mm_set1_pd(-0.374089300177174709737e-5);
 408     const __m128d CAQ4      = _mm_set1_pd(0.00015126584532155383535);
 409     const __m128d CAQ3      = _mm_set1_pd(0.00536692680669480725423);
 410     const __m128d CAQ2      = _mm_set1_pd(0.0668686825594046122636);
 411     const __m128d CAQ1      = _mm_set1_pd(0.402604990869284362773);
 412     /* CAQ0 == 1.0 */
 413     const __m128d CAoffset  = _mm_set1_pd(0.9788494110107421875);
 414
 415     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5] */
 416     const __m128d CBP6      = _mm_set1_pd(2.49650423685462752497647637088e-10);
 417     const __m128d CBP5      = _mm_set1_pd(0.00119770193298159629350136085658);
 418     const __m128d CBP4      = _mm_set1_pd(0.0164944422378370965881008942733);
 419     const __m128d CBP3      = _mm_set1_pd(0.0984581468691775932063932439252);
 420     const __m128d CBP2      = _mm_set1_pd(0.317364595806937763843589437418);
 421     const __m128d CBP1      = _mm_set1_pd(0.554167062641455850932670067075);
 422     const __m128d CBP0      = _mm_set1_pd(0.427583576155807163756925301060);
 423     const __m128d CBQ7      = _mm_set1_pd(0.00212288829699830145976198384930);
 424     const __m128d CBQ6      = _mm_set1_pd(0.0334810979522685300554606393425);
 425     const __m128d CBQ5      = _mm_set1_pd(0.2361713785181450957579508850717);
 426     const __m128d CBQ4      = _mm_set1_pd(0.955364736493055670530981883072);
 427     const __m128d CBQ3      = _mm_set1_pd(2.36815675631420037315349279199);
 428     const __m128d CBQ2      = _mm_set1_pd(3.55261649184083035537184223542);
 429     const __m128d CBQ1      = _mm_set1_pd(2.93501136050160872574376997993);
 430     /* CBQ0 == 1.0 */
 431
 432     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf] */
 433     const __m128d CCP6      = _mm_set1_pd(-2.8175401114513378771);
 434     const __m128d CCP5      = _mm_set1_pd(-3.22729451764143718517);
 435     const __m128d CCP4      = _mm_set1_pd(-2.5518551727311523996);
 436     const __m128d CCP3      = _mm_set1_pd(-0.687717681153649930619);
 437     const __m128d CCP2      = _mm_set1_pd(-0.212652252872804219852);
 438     const __m128d CCP1      = _mm_set1_pd(0.0175389834052493308818);
 439     const __m128d CCP0      = _mm_set1_pd(0.00628057170626964891937);
 440
 441     const __m128d CCQ6      = _mm_set1_pd(5.48409182238641741584);
 442     const __m128d CCQ5      = _mm_set1_pd(13.5064170191802889145);
 443     const __m128d CCQ4      = _mm_set1_pd(22.9367376522880577224);
 444     const __m128d CCQ3      = _mm_set1_pd(15.930646027911794143);
 445     const __m128d CCQ2      = _mm_set1_pd(11.0567237927800161565);
 446     const __m128d CCQ1      = _mm_set1_pd(2.79257750980575282228);
 447     /* CCQ0 == 1.0 */
 448     const __m128d CCoffset  = _mm_set1_pd(0.5579090118408203125);
 449
 450     const __m128d one       = _mm_set1_pd(1.0);
 451     const __m128d two       = _mm_set1_pd(2.0);
 452
 453     const __m128d signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 454
 455     __m128d xabs,x2,x4,t,t2,w,w2;
 456     __m128d PolyAP0,PolyAP1,PolyAQ0,PolyAQ1;
 457     __m128d PolyBP0,PolyBP1,PolyBQ0,PolyBQ1;
 458     __m128d PolyCP0,PolyCP1,PolyCQ0,PolyCQ1;
 459     __m128d res_erf,res_erfcB,res_erfcC,res_erfc,res;
 460     __m128d mask,expmx2;
 461
 462     /* Calculate erf() */
 463     xabs     = gmx_mm_abs_pd(x);
 464     x2       = _mm_mul_pd(x,x);
 465     x4       = _mm_mul_pd(x2,x2);
 466
 467     PolyAP0  = _mm_macc_pd(CAP4,x4,CAP2);
 468     PolyAP1  = _mm_macc_pd(CAP3,x4,CAP1);
 469     PolyAP0  = _mm_macc_pd(PolyAP0,x4,CAP0);
 470     PolyAP0  = _mm_macc_pd(PolyAP1,x2,PolyAP0);
 471
 472     PolyAQ1  = _mm_macc_pd(CAQ5,x4,CAQ3);
 473     PolyAQ0  = _mm_macc_pd(CAQ4,x4,CAQ2);
 474     PolyAQ1  = _mm_macc_pd(PolyAQ1,x4,CAQ1);
 475     PolyAQ0  = _mm_macc_pd(PolyAQ0,x4,one);
 476     PolyAQ0  = _mm_macc_pd(PolyAQ1,x2,PolyAQ0);
 477
 478     res_erf  = _mm_macc_pd(PolyAP0,gmx_mm_inv_pd(PolyAQ0),CAoffset);
 479     res_erf  = _mm_mul_pd(x,res_erf);
 480
 481     /* Calculate erfc() in range [1,4.5] */
 482     t       = _mm_sub_pd(xabs,one);
 483     t2      = _mm_mul_pd(t,t);
 484
 485     PolyBP0  = _mm_macc_pd(CBP6,t2,CBP4);
 486     PolyBP1  = _mm_macc_pd(CBP5,t2,CBP3);
 487     PolyBP0  = _mm_macc_pd(PolyBP0,t2,CBP2);
 488     PolyBP1  = _mm_macc_pd(PolyBP1,t2,CBP1);
 489     PolyBP0  = _mm_macc_pd(PolyBP0,t2,CBP0);
 490     PolyBP0  = _mm_macc_pd(PolyBP1,t,PolyBP0);
 491
 492     PolyBQ1 = _mm_macc_pd(CBQ7,t2,CBQ5);
 493     PolyBQ0 = _mm_macc_pd(CBQ6,t2,CBQ4);
 494     PolyBQ1 = _mm_macc_pd(PolyBQ1,t2,CBQ3);
 495     PolyBQ0 = _mm_macc_pd(PolyBQ0,t2,CBQ2);
 496     PolyBQ1 = _mm_macc_pd(PolyBQ1,t2,CBQ1);
 497     PolyBQ0 = _mm_macc_pd(PolyBQ0,t2,one);
 498     PolyBQ0 = _mm_macc_pd(PolyBQ1,t,PolyBQ0);
 499
 500     res_erfcB = _mm_mul_pd(PolyBP0,gmx_mm_inv_pd(PolyBQ0));
 501
 502     res_erfcB = _mm_mul_pd(res_erfcB,xabs);
 503
 504     /* Calculate erfc() in range [4.5,inf] */
 505     w       = gmx_mm_inv_pd(xabs);
 506     w2      = _mm_mul_pd(w,w);
 507
 508     PolyCP0  = _mm_macc_pd(CCP6,w2,CCP4);
 509     PolyCP1  = _mm_macc_pd(CCP5,w2,CCP3);
 510     PolyCP0  = _mm_macc_pd(PolyCP0,w2,CCP2);
 511     PolyCP1  = _mm_macc_pd(PolyCP1,w2,CCP1);
 512     PolyCP0  = _mm_macc_pd(PolyCP0,w2,CCP0);
 513     PolyCP0  = _mm_macc_pd(PolyCP1,w,PolyCP0);
 514
 515     PolyCQ0  = _mm_macc_pd(CCQ6,w2,CCQ4);
 516     PolyCQ1  = _mm_macc_pd(CCQ5,w2,CCQ3);
 517     PolyCQ0  = _mm_macc_pd(PolyCQ0,w2,CCQ2);
 518     PolyCQ1  = _mm_macc_pd(PolyCQ1,w2,CCQ1);
 519     PolyCQ0  = _mm_macc_pd(PolyCQ0,w2,one);
 520     PolyCQ0  = _mm_macc_pd(PolyCQ1,w,PolyCQ0);
 521
 522     expmx2   = gmx_mm_exp_pd( _mm_or_pd(signbit, x2) );
 523
 524     res_erfcC = _mm_macc_pd(PolyCP0,gmx_mm_inv_pd(PolyCQ0),CCoffset);
 525     res_erfcC = _mm_mul_pd(res_erfcC,w);
 526
 527     mask = _mm_cmpgt_pd(xabs,_mm_set1_pd(4.5));
 528     res_erfc = _mm_blendv_pd(res_erfcB,res_erfcC,mask);
 529
 530     res_erfc = _mm_mul_pd(res_erfc,expmx2);
 531
 532     /* erfc(x<0) = 2-erfc(|x|) */
 533     mask = _mm_cmplt_pd(x,_mm_setzero_pd());
 534     res_erfc = _mm_blendv_pd(res_erfc,_mm_sub_pd(two,res_erfc),mask);
 535
 536     /* Select erf() or erfc() */
 537     mask = _mm_cmplt_pd(xabs,one);
 538     res  = _mm_blendv_pd(_mm_sub_pd(one,res_erfc),res_erf,mask);
 539
 540     return res;
 541 }
 542
 543
 544 static __m128d
 545 gmx_mm_erfc_pd(__m128d x)
 546 {
 547     /* Coefficients for minimax approximation of erf(x)=x*(CAoffset + P(x^2)/Q(x^2)) in range [-0.75,0.75] */
 548     const __m128d CAP4      = _mm_set1_pd(-0.431780540597889301512e-4);
 549     const __m128d CAP3      = _mm_set1_pd(-0.00578562306260059236059);
 550     const __m128d CAP2      = _mm_set1_pd(-0.028593586920219752446);
 551     const __m128d CAP1      = _mm_set1_pd(-0.315924962948621698209);
 552     const __m128d CAP0      = _mm_set1_pd(0.14952975608477029151);
 553
 554     const __m128d CAQ5      = _mm_set1_pd(-0.374089300177174709737e-5);
 555     const __m128d CAQ4      = _mm_set1_pd(0.00015126584532155383535);
 556     const __m128d CAQ3      = _mm_set1_pd(0.00536692680669480725423);
 557     const __m128d CAQ2      = _mm_set1_pd(0.0668686825594046122636);
 558     const __m128d CAQ1      = _mm_set1_pd(0.402604990869284362773);
 559     /* CAQ0 == 1.0 */
 560     const __m128d CAoffset  = _mm_set1_pd(0.9788494110107421875);
 561
 562     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)*x*(P(x-1)/Q(x-1)) in range [1.0,4.5] */
 563     const __m128d CBP6      = _mm_set1_pd(2.49650423685462752497647637088e-10);
 564     const __m128d CBP5      = _mm_set1_pd(0.00119770193298159629350136085658);
 565     const __m128d CBP4      = _mm_set1_pd(0.0164944422378370965881008942733);
 566     const __m128d CBP3      = _mm_set1_pd(0.0984581468691775932063932439252);
 567     const __m128d CBP2      = _mm_set1_pd(0.317364595806937763843589437418);
 568     const __m128d CBP1      = _mm_set1_pd(0.554167062641455850932670067075);
 569     const __m128d CBP0      = _mm_set1_pd(0.427583576155807163756925301060);
 570     const __m128d CBQ7      = _mm_set1_pd(0.00212288829699830145976198384930);
 571     const __m128d CBQ6      = _mm_set1_pd(0.0334810979522685300554606393425);
 572     const __m128d CBQ5      = _mm_set1_pd(0.2361713785181450957579508850717);
 573     const __m128d CBQ4      = _mm_set1_pd(0.955364736493055670530981883072);
 574     const __m128d CBQ3      = _mm_set1_pd(2.36815675631420037315349279199);
 575     const __m128d CBQ2      = _mm_set1_pd(3.55261649184083035537184223542);
 576     const __m128d CBQ1      = _mm_set1_pd(2.93501136050160872574376997993);
 577     /* CBQ0 == 1.0 */
 578
 579     /* Coefficients for minimax approximation of erfc(x)=exp(-x^2)/x*(P(1/x)/Q(1/x)) in range [4.5,inf] */
 580     const __m128d CCP6      = _mm_set1_pd(-2.8175401114513378771);
 581     const __m128d CCP5      = _mm_set1_pd(-3.22729451764143718517);
 582     const __m128d CCP4      = _mm_set1_pd(-2.5518551727311523996);
 583     const __m128d CCP3      = _mm_set1_pd(-0.687717681153649930619);
 584     const __m128d CCP2      = _mm_set1_pd(-0.212652252872804219852);
 585     const __m128d CCP1      = _mm_set1_pd(0.0175389834052493308818);
 586     const __m128d CCP0      = _mm_set1_pd(0.00628057170626964891937);
 587
 588     const __m128d CCQ6      = _mm_set1_pd(5.48409182238641741584);
 589     const __m128d CCQ5      = _mm_set1_pd(13.5064170191802889145);
 590     const __m128d CCQ4      = _mm_set1_pd(22.9367376522880577224);
 591     const __m128d CCQ3      = _mm_set1_pd(15.930646027911794143);
 592     const __m128d CCQ2      = _mm_set1_pd(11.0567237927800161565);
 593     const __m128d CCQ1      = _mm_set1_pd(2.79257750980575282228);
 594     /* CCQ0 == 1.0 */
 595     const __m128d CCoffset  = _mm_set1_pd(0.5579090118408203125);
 596
 597     const __m128d one       = _mm_set1_pd(1.0);
 598     const __m128d two       = _mm_set1_pd(2.0);
 599
 600     const __m128d signbit   = gmx_mm_castsi128_pd( _mm_set_epi32(0x80000000,0x00000000,0x80000000,0x00000000) );
 601
 602     __m128d xabs,x2,x4,t,t2,w,w2;
 603     __m128d PolyAP0,PolyAP1,PolyAQ0,PolyAQ1;
 604     __m128d PolyBP0,PolyBP1,PolyBQ0,PolyBQ1;
 605     __m128d PolyCP0,PolyCP1,PolyCQ0,PolyCQ1;
 606     __m128d res_erf,res_erfcB,res_erfcC,res_erfc,res;
 607     __m128d mask,expmx2;
 608
 609     /* Calculate erf() */
 610     xabs     = gmx_mm_abs_pd(x);
 611     x2       = _mm_mul_pd(x,x);
 612     x4       = _mm_mul_pd(x2,x2);
 613
 614     PolyAP0  = _mm_macc_pd(CAP4,x4,CAP2);
 615     PolyAP1  = _mm_macc_pd(CAP3,x4,CAP1);
 616     PolyAP0  = _mm_macc_pd(PolyAP0,x4,CAP0);
 617     PolyAP0  = _mm_macc_pd(PolyAP1,x2,PolyAP0);
 618
 619     PolyAQ1  = _mm_macc_pd(CAQ5,x4,CAQ3);
 620     PolyAQ0  = _mm_macc_pd(CAQ4,x4,CAQ2);
 621     PolyAQ1  = _mm_macc_pd(PolyAQ1,x4,CAQ1);
 622     PolyAQ0  = _mm_macc_pd(PolyAQ0,x4,one);
 623     PolyAQ0  = _mm_macc_pd(PolyAQ1,x2,PolyAQ0);
 624
 625     res_erf  = _mm_macc_pd(PolyAP0,gmx_mm_inv_pd(PolyAQ0),CAoffset);
 626     res_erf  = _mm_mul_pd(x,res_erf);
 627
 628     /* Calculate erfc() in range [1,4.5] */
 629     t       = _mm_sub_pd(xabs,one);
 630     t2      = _mm_mul_pd(t,t);
 631
 632     PolyBP0  = _mm_macc_pd(CBP6,t2,CBP4);
 633     PolyBP1  = _mm_macc_pd(CBP5,t2,CBP3);
 634     PolyBP0  = _mm_macc_pd(PolyBP0,t2,CBP2);
 635     PolyBP1  = _mm_macc_pd(PolyBP1,t2,CBP1);
 636     PolyBP0  = _mm_macc_pd(PolyBP0,t2,CBP0);
 637     PolyBP0  = _mm_macc_pd(PolyBP1,t,PolyBP0);
 638
 639     PolyBQ1 = _mm_macc_pd(CBQ7,t2,CBQ5);
 640     PolyBQ0 = _mm_macc_pd(CBQ6,t2,CBQ4);
 641     PolyBQ1 = _mm_macc_pd(PolyBQ1,t2,CBQ3);
 642     PolyBQ0 = _mm_macc_pd(PolyBQ0,t2,CBQ2);
 643     PolyBQ1 = _mm_macc_pd(PolyBQ1,t2,CBQ1);
 644     PolyBQ0 = _mm_macc_pd(PolyBQ0,t2,one);
 645     PolyBQ0 = _mm_macc_pd(PolyBQ1,t,PolyBQ0);
 646
 647     res_erfcB = _mm_mul_pd(PolyBP0,gmx_mm_inv_pd(PolyBQ0));
 648
 649     res_erfcB = _mm_mul_pd(res_erfcB,xabs);
 650
 651     /* Calculate erfc() in range [4.5,inf] */
 652     w       = gmx_mm_inv_pd(xabs);
 653     w2      = _mm_mul_pd(w,w);
 654
 655     PolyCP0  = _mm_macc_pd(CCP6,w2,CCP4);
 656     PolyCP1  = _mm_macc_pd(CCP5,w2,CCP3);
 657     PolyCP0  = _mm_macc_pd(PolyCP0,w2,CCP2);
 658     PolyCP1  = _mm_macc_pd(PolyCP1,w2,CCP1);
 659     PolyCP0  = _mm_macc_pd(PolyCP0,w2,CCP0);
 660     PolyCP0  = _mm_macc_pd(PolyCP1,w,PolyCP0);
 661
 662     PolyCQ0  = _mm_macc_pd(CCQ6,w2,CCQ4);
 663     PolyCQ1  = _mm_macc_pd(CCQ5,w2,CCQ3);
 664     PolyCQ0  = _mm_macc_pd(PolyCQ0,w2,CCQ2);
 665     PolyCQ1  = _mm_macc_pd(PolyCQ1,w2,CCQ1);
 666     PolyCQ0  = _mm_macc_pd(PolyCQ0,w2,one);
 667     PolyCQ0  = _mm_macc_pd(PolyCQ1,w,PolyCQ0);
 668
 669     expmx2   = gmx_mm_exp_pd( _mm_or_pd(signbit, x2) );
 670
 671     res_erfcC = _mm_macc_pd(PolyCP0,gmx_mm_inv_pd(PolyCQ0),CCoffset);
 672     res_erfcC = _mm_mul_pd(res_erfcC,w);
 673
 674     mask = _mm_cmpgt_pd(xabs,_mm_set1_pd(4.5));
 675     res_erfc = _mm_blendv_pd(res_erfcB,res_erfcC,mask);
 676
 677     res_erfc = _mm_mul_pd(res_erfc,expmx2);
 678
 679     /* erfc(x<0) = 2-erfc(|x|) */
 680     mask = _mm_cmplt_pd(x,_mm_setzero_pd());
 681     res_erfc = _mm_blendv_pd(res_erfc,_mm_sub_pd(two,res_erfc),mask);
 682
 683     /* Select erf() or erfc() */
 684     mask = _mm_cmplt_pd(xabs,one);
 685     res  = _mm_blendv_pd(res_erfc,_mm_sub_pd(one,res_erf),mask);
 686
 687     return res;
 688 }
 689
 690
 691
 692 /* Calculate the force correction due to PME analytically.
 693  *
 694  * This routine is meant to enable analytical evaluation of the
 695  * direct-space PME electrostatic force to avoid tables.
 696  *
 697  * The direct-space potential should be Erfc(beta*r)/r, but there
 698  * are some problems evaluating that:
 699  *
 700  * First, the error function is difficult (read: expensive) to
 701  * approxmiate accurately for intermediate to large arguments, and
 702  * this happens already in ranges of beta*r that occur in simulations.
 703  * Second, we now try to avoid calculating potentials in Gromacs but
 704  * use forces directly.
 705  *
 706  * We can simply things slight by noting that the PME part is really
 707  * a correction to the normal Coulomb force since Erfc(z)=1-Erf(z), i.e.
 708  *
 709  * V= 1/r - Erf(beta*r)/r
 710  *
 711  * The first term we already have from the inverse square root, so
 712  * that we can leave out of this routine.
 713  *
 714  * For pme tolerances of 1e-3 to 1e-8 and cutoffs of 0.5nm to 1.8nm,
 715  * the argument beta*r will be in the range 0.15 to ~4. Use your
 716  * favorite plotting program to realize how well-behaved Erf(z)/z is
 717  * in this range!
 718  *
 719  * We approximate f(z)=erf(z)/z with a rational minimax polynomial.
 720  * However, it turns out it is more efficient to approximate f(z)/z and
 721  * then only use even powers. This is another minor optimization, since
 722  * we actually WANT f(z)/z, because it is going to be multiplied by
 723  * the vector between the two atoms to get the vectorial force. The
 724  * fastest flops are the ones we can avoid calculating!
 725  *
 726  * So, here's how it should be used:
 727  *
 728  * 1. Calculate r^2.
 729  * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
 730  * 3. Evaluate this routine with z^2 as the argument.
 731  * 4. The return value is the expression:
 732  *
 733  *
 734  *       2*exp(-z^2)     erf(z)
 735  *       ------------ - --------
 736  *       sqrt(Pi)*z^2      z^3
 737  *
 738  * 5. Multiply the entire expression by beta^3. This will get you
 739  *
 740  *       beta^3*2*exp(-z^2)     beta^3*erf(z)
 741  *       ------------------  - ---------------
 742  *          sqrt(Pi)*z^2            z^3
 743  *
 744  *    or, switching back to r (z=r*beta):
 745  *
 746  *       2*beta*exp(-r^2*beta^2)   erf(r*beta)
 747  *       ----------------------- - -----------
 748  *            sqrt(Pi)*r^2            r^3
 749  *
 750  *
 751  *    With a bit of math exercise you should be able to confirm that
 752  *    this is exactly D[Erf[beta*r]/r,r] divided by r another time.
 753  *
 754  * 6. Add the result to 1/r^3, multiply by the product of the charges,
 755  *    and you have your force (divided by r). A final multiplication
 756  *    with the vector connecting the two particles and you have your
 757  *    vectorial force to add to the particles.
 758  *
 759  */
 760 static __m128d
 761 gmx_mm_pmecorrF_pd(__m128d z2)
 762 {
 763     const __m128d  FN10     = _mm_set1_pd(-8.0072854618360083154e-14);
 764     const __m128d  FN9      = _mm_set1_pd(1.1859116242260148027e-11);
 765     const __m128d  FN8      = _mm_set1_pd(-8.1490406329798423616e-10);
 766     const __m128d  FN7      = _mm_set1_pd(3.4404793543907847655e-8);
 767     const __m128d  FN6      = _mm_set1_pd(-9.9471420832602741006e-7);
 768     const __m128d  FN5      = _mm_set1_pd(0.000020740315999115847456);
 769     const __m128d  FN4      = _mm_set1_pd(-0.00031991745139313364005);
 770     const __m128d  FN3      = _mm_set1_pd(0.0035074449373659008203);
 771     const __m128d  FN2      = _mm_set1_pd(-0.031750380176100813405);
 772     const __m128d  FN1      = _mm_set1_pd(0.13884101728898463426);
 773     const __m128d  FN0      = _mm_set1_pd(-0.75225277815249618847);
 774
 775     const __m128d  FD5      = _mm_set1_pd(0.000016009278224355026701);
 776     const __m128d  FD4      = _mm_set1_pd(0.00051055686934806966046);
 777     const __m128d  FD3      = _mm_set1_pd(0.0081803507497974289008);
 778     const __m128d  FD2      = _mm_set1_pd(0.077181146026670287235);
 779     const __m128d  FD1      = _mm_set1_pd(0.41543303143712535988);
 780     const __m128d  FD0      = _mm_set1_pd(1.0);
 781
 782     __m128d z4;
 783     __m128d polyFN0,polyFN1,polyFD0,polyFD1;
 784
 785     z4             = _mm_mul_pd(z2,z2);
 786
 787     polyFD1        = _mm_macc_pd(FD5,z4,FD3);
 788     polyFD1        = _mm_macc_pd(polyFD1,z4,FD1);
 789     polyFD1        = _mm_mul_pd(polyFD1,z2);
 790     polyFD0        = _mm_macc_pd(FD4,z4,FD2);
 791     polyFD0        = _mm_macc_pd(polyFD0,z4,FD0);
 792     polyFD0        = _mm_add_pd(polyFD0,polyFD1);
 793
 794     polyFD0        = gmx_mm_inv_pd(polyFD0);
 795
 796     polyFN0        = _mm_macc_pd(FN10,z4,FN8);
 797     polyFN0        = _mm_macc_pd(polyFN0,z4,FN6);
 798     polyFN0        = _mm_macc_pd(polyFN0,z4,FN4);
 799     polyFN0        = _mm_macc_pd(polyFN0,z4,FN2);
 800     polyFN0        = _mm_macc_pd(polyFN0,z4,FN0);
 801     polyFN1        = _mm_macc_pd(FN9,z4,FN7);
 802     polyFN1        = _mm_macc_pd(polyFN1,z4,FN5);
 803     polyFN1        = _mm_macc_pd(polyFN1,z4,FN3);
 804     polyFN1        = _mm_macc_pd(polyFN1,z4,FN1);
 805     polyFN0        = _mm_macc_pd(polyFN1,z2,polyFN0);
 806
 807     return   _mm_mul_pd(polyFN0,polyFD0);
 808 }
 809
 810
 811 /* Calculate the potential correction due to PME analytically.
 812  *
 813  * This routine calculates Erf(z)/z, although you should provide z^2
 814  * as the input argument.
 815  *
 816  * Here's how it should be used:
 817  *
 818  * 1. Calculate r^2.
 819  * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
 820  * 3. Evaluate this routine with z^2 as the argument.
 821  * 4. The return value is the expression:
 822  *
 823  *
 824  *        erf(z)
 825  *       --------
 826  *          z
 827  *
 828  * 5. Multiply the entire expression by beta and switching back to r (z=r*beta):
 829  *
 830  *       erf(r*beta)
 831  *       -----------
 832  *           r
 833  *
 834  * 6. Subtract the result from 1/r, multiply by the product of the charges,
 835  *    and you have your potential.
 836  *
 837  */
 838 static __m128d
 839 gmx_mm_pmecorrV_pd(__m128d z2)
 840 {
 841     const __m128d  VN9      = _mm_set1_pd(-9.3723776169321855475e-13);
 842     const __m128d  VN8      = _mm_set1_pd(1.2280156762674215741e-10);
 843     const __m128d  VN7      = _mm_set1_pd(-7.3562157912251309487e-9);
 844     const __m128d  VN6      = _mm_set1_pd(2.6215886208032517509e-7);
 845     const __m128d  VN5      = _mm_set1_pd(-4.9532491651265819499e-6);
 846     const __m128d  VN4      = _mm_set1_pd(0.00025907400778966060389);
 847     const __m128d  VN3      = _mm_set1_pd(0.0010585044856156469792);
 848     const __m128d  VN2      = _mm_set1_pd(0.045247661136833092885);
 849     const __m128d  VN1      = _mm_set1_pd(0.11643931522926034421);
 850     const __m128d  VN0      = _mm_set1_pd(1.1283791671726767970);
 851
 852     const __m128d  VD5      = _mm_set1_pd(0.000021784709867336150342);
 853     const __m128d  VD4      = _mm_set1_pd(0.00064293662010911388448);
 854     const __m128d  VD3      = _mm_set1_pd(0.0096311444822588683504);
 855     const __m128d  VD2      = _mm_set1_pd(0.085608012351550627051);
 856     const __m128d  VD1      = _mm_set1_pd(0.43652499166614811084);
 857     const __m128d  VD0      = _mm_set1_pd(1.0);
 858
 859     __m128d z4;
 860     __m128d polyVN0,polyVN1,polyVD0,polyVD1;
 861
 862     z4             = _mm_mul_pd(z2,z2);
 863
 864     polyVD1        = _mm_macc_pd(VD5,z4,VD3);
 865     polyVD0        = _mm_macc_pd(VD4,z4,VD2);
 866     polyVD1        = _mm_macc_pd(polyVD1,z4,VD1);
 867     polyVD0        = _mm_macc_pd(polyVD0,z4,VD0);
 868     polyVD0        = _mm_macc_pd(polyVD1,z2,polyVD0);
 869
 870     polyVD0        = gmx_mm_inv_pd(polyVD0);
 871
 872     polyVN1        = _mm_macc_pd(VN9,z4,VN7);
 873     polyVN0        = _mm_macc_pd(VN8,z4,VN6);
 874     polyVN1        = _mm_macc_pd(polyVN1,z4,VN5);
 875     polyVN0        = _mm_macc_pd(polyVN0,z4,VN4);
 876     polyVN1        = _mm_macc_pd(polyVN1,z4,VN3);
 877     polyVN0        = _mm_macc_pd(polyVN0,z4,VN2);
 878     polyVN1        = _mm_macc_pd(polyVN1,z4,VN1);
 879     polyVN0        = _mm_macc_pd(polyVN0,z4,VN0);
 880     polyVN0        = _mm_macc_pd(polyVN1,z2,polyVN0);
 881
 882     return   _mm_mul_pd(polyVN0,polyVD0);
 883 }
 884
 885
 886 static int
 887 gmx_mm_sincos_pd(__m128d x,
 888                  __m128d *sinval,
 889                  __m128d *cosval)
 890 {
 891 #ifdef _MSC_VER
 892     __declspec(align(16))
 893     const double sintable[34] =
 894     {
 895         1.00000000000000000e+00 , 0.00000000000000000e+00 ,
 896         9.95184726672196929e-01 , 9.80171403295606036e-02 ,
 897         9.80785280403230431e-01 , 1.95090322016128248e-01 ,
 898         9.56940335732208824e-01 , 2.90284677254462331e-01 ,
 899         9.23879532511286738e-01 , 3.82683432365089782e-01 ,
 900         8.81921264348355050e-01 , 4.71396736825997642e-01 ,
 901         8.31469612302545236e-01 , 5.55570233019602178e-01 ,
 902         7.73010453362736993e-01 , 6.34393284163645488e-01 ,
 903         7.07106781186547573e-01 , 7.07106781186547462e-01 ,
 904         6.34393284163645599e-01 , 7.73010453362736882e-01 ,
 905         5.55570233019602289e-01 , 8.31469612302545125e-01 ,
 906         4.71396736825997809e-01 , 8.81921264348354939e-01 ,
 907         3.82683432365089837e-01 , 9.23879532511286738e-01 ,
 908         2.90284677254462276e-01 , 9.56940335732208935e-01 ,
 909         1.95090322016128304e-01 , 9.80785280403230431e-01 ,
 910         9.80171403295607702e-02 , 9.95184726672196818e-01 ,
 911         0.0 , 1.00000000000000000e+00
 912     };
 913 #else
 914     const __m128d sintable[17] =
 915     {
 916         _mm_set_pd( 0.0 , 1.0 ),
 917         _mm_set_pd( sin(  1.0 * (M_PI/2.0) / 16.0) , cos(  1.0 * (M_PI/2.0) / 16.0) ),
 918         _mm_set_pd( sin(  2.0 * (M_PI/2.0) / 16.0) , cos(  2.0 * (M_PI/2.0) / 16.0) ),
 919         _mm_set_pd( sin(  3.0 * (M_PI/2.0) / 16.0) , cos(  3.0 * (M_PI/2.0) / 16.0) ),
 920         _mm_set_pd( sin(  4.0 * (M_PI/2.0) / 16.0) , cos(  4.0 * (M_PI/2.0) / 16.0) ),
 921         _mm_set_pd( sin(  5.0 * (M_PI/2.0) / 16.0) , cos(  5.0 * (M_PI/2.0) / 16.0) ),
 922         _mm_set_pd( sin(  6.0 * (M_PI/2.0) / 16.0) , cos(  6.0 * (M_PI/2.0) / 16.0) ),
 923         _mm_set_pd( sin(  7.0 * (M_PI/2.0) / 16.0) , cos(  7.0 * (M_PI/2.0) / 16.0) ),
 924         _mm_set_pd( sin(  8.0 * (M_PI/2.0) / 16.0) , cos(  8.0 * (M_PI/2.0) / 16.0) ),
 925         _mm_set_pd( sin(  9.0 * (M_PI/2.0) / 16.0) , cos(  9.0 * (M_PI/2.0) / 16.0) ),
 926         _mm_set_pd( sin( 10.0 * (M_PI/2.0) / 16.0) , cos( 10.0 * (M_PI/2.0) / 16.0) ),
 927         _mm_set_pd( sin( 11.0 * (M_PI/2.0) / 16.0) , cos( 11.0 * (M_PI/2.0) / 16.0) ),
 928         _mm_set_pd( sin( 12.0 * (M_PI/2.0) / 16.0) , cos( 12.0 * (M_PI/2.0) / 16.0) ),
 929         _mm_set_pd( sin( 13.0 * (M_PI/2.0) / 16.0) , cos( 13.0 * (M_PI/2.0) / 16.0) ),
 930         _mm_set_pd( sin( 14.0 * (M_PI/2.0) / 16.0) , cos( 14.0 * (M_PI/2.0) / 16.0) ),
 931         _mm_set_pd( sin( 15.0 * (M_PI/2.0) / 16.0) , cos( 15.0 * (M_PI/2.0) / 16.0) ),
 932         _mm_set_pd(  1.0 , 0.0 )
 933     };
 934 #endif
 935
 936     const __m128d signmask    = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
 937     const __m128i signbit_epi32  = _mm_set1_epi32(0x80000000);
 938
 939     const __m128d tabscale      = _mm_set1_pd(32.0/M_PI);
 940     const __m128d invtabscale0  = _mm_set1_pd(9.81747508049011230469e-02);
 941     const __m128d invtabscale1  = _mm_set1_pd(1.96197799156550576057e-08);
 942     const __m128i ione          = _mm_set1_epi32(1);
 943     const __m128i i32           = _mm_set1_epi32(32);
 944     const __m128i i16           = _mm_set1_epi32(16);
 945     const __m128i tabmask       = _mm_set1_epi32(0x3F);
 946     const __m128d sinP7         = _mm_set1_pd(-1.0/5040.0);
 947     const __m128d sinP5         = _mm_set1_pd(1.0/120.0);
 948     const __m128d sinP3         = _mm_set1_pd(-1.0/6.0);
 949     const __m128d sinP1         = _mm_set1_pd(1.0);
 950
 951     const __m128d cosP6         = _mm_set1_pd(-1.0/720.0);
 952     const __m128d cosP4         = _mm_set1_pd(1.0/24.0);
 953     const __m128d cosP2         = _mm_set1_pd(-1.0/2.0);
 954     const __m128d cosP0         = _mm_set1_pd(1.0);
 955
 956     __m128d scalex;
 957     __m128i tabidx,corridx;
 958     __m128d xabs,z,z2,polySin,polyCos;
 959     __m128d xpoint;
 960     __m128d ypoint0,ypoint1;
 961
 962     __m128d sinpoint,cospoint;
 963     __m128d xsign,ssign,csign;
 964     __m128i imask,sswapsign,cswapsign;
 965     __m128d minusone;
 966
 967     xsign    = _mm_andnot_pd(signmask,x);
 968     xabs     = _mm_and_pd(x,signmask);
 969
 970     scalex   = _mm_mul_pd(tabscale,xabs);
 971     tabidx   = _mm_cvtpd_epi32(scalex);
 972
 973     xpoint   = _mm_round_pd(scalex,_MM_FROUND_TO_NEAREST_INT);
 974
 975     /* Extended precision arithmetics */
 976     z        = _mm_nmacc_pd(invtabscale0,xpoint,xabs);
 977     z        = _mm_nmacc_pd(invtabscale1,xpoint,z);
 978
 979     /* Range reduction to 0..2*Pi */
 980     tabidx   = _mm_and_si128(tabidx,tabmask);
 981
 982     /* tabidx is now in range [0,..,64] */
 983     imask     = _mm_cmpgt_epi32(tabidx,i32);
 984     sswapsign = imask;
 985     cswapsign = imask;
 986     corridx   = _mm_and_si128(imask,i32);
 987     tabidx    = _mm_sub_epi32(tabidx,corridx);
 988
 989     /* tabidx is now in range [0..32] */
 990     imask     = _mm_cmpgt_epi32(tabidx,i16);
 991     cswapsign = _mm_xor_si128(cswapsign,imask);
 992     corridx   = _mm_sub_epi32(i32,tabidx);
 993     tabidx    = _mm_blendv_epi8(tabidx,corridx,imask);
 994     /* tabidx is now in range [0..16] */
 995     ssign     = _mm_cvtepi32_pd( _mm_or_si128( sswapsign , ione ) );
 996     csign     = _mm_cvtepi32_pd( _mm_or_si128( cswapsign , ione ) );
 997
 998 #ifdef _MSC_VER
 999     ypoint0  = _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx,0));
1000     ypoint1  = _mm_load_pd(sintable + 2*_mm_extract_epi32(tabidx,1));
1001 #else
1002     ypoint0  = sintable[_mm_extract_epi32(tabidx,0)];
1003     ypoint1  = sintable[_mm_extract_epi32(tabidx,1)];
1004 #endif
1005     sinpoint = _mm_unpackhi_pd(ypoint0,ypoint1);
1006     cospoint = _mm_unpacklo_pd(ypoint0,ypoint1);
1007
1008     sinpoint = _mm_mul_pd(sinpoint,ssign);
1009     cospoint = _mm_mul_pd(cospoint,csign);
1010
1011     z2       = _mm_mul_pd(z,z);
1012
1013     polySin  = _mm_macc_pd(sinP7,z2,sinP5);
1014     polySin  = _mm_macc_pd(polySin,z2,sinP3);
1015     polySin  = _mm_macc_pd(polySin,z2,sinP1);
1016     polySin  = _mm_mul_pd(polySin,z);
1017
1018     polyCos  = _mm_macc_pd(cosP6,z2,cosP4);
1019     polyCos  = _mm_macc_pd(polyCos,z2,cosP2);
1020     polyCos  = _mm_macc_pd(polyCos,z2,cosP0);
1021
1022     *sinval  = _mm_xor_pd(_mm_add_pd( _mm_mul_pd(sinpoint,polyCos) , _mm_mul_pd(cospoint,polySin) ),xsign);
1023     *cosval  = _mm_sub_pd( _mm_mul_pd(cospoint,polyCos) , _mm_mul_pd(sinpoint,polySin) );
1024
1025     return 0;
1026 }
1027
1028 /*
1029  * IMPORTANT: Do NOT call both sin & cos if you need both results, since each of them
1030  * will then call the sincos() routine and waste a factor 2 in performance!
1031  */
1032 static __m128d
1033 gmx_mm_sin_pd(__m128d x)
1034 {
1035     __m128d s,c;
1036     gmx_mm_sincos_pd(x,&s,&c);
1037     return s;
1038 }
1039
1040 /*
1041  * IMPORTANT: Do NOT call both sin & cos if you need both results, since each of them
1042  * will then call the sincos() routine and waste a factor 2 in performance!
1043  */
1044 static __m128d
1045 gmx_mm_cos_pd(__m128d x)
1046 {
1047     __m128d s,c;
1048     gmx_mm_sincos_pd(x,&s,&c);
1049     return c;
1050 }
1051
1052
1053
1054 static __m128d
1055 gmx_mm_tan_pd(__m128d x)
1056 {
1057     __m128d sinval,cosval;
1058     __m128d tanval;
1059
1060     gmx_mm_sincos_pd(x,&sinval,&cosval);
1061
1062     tanval = _mm_mul_pd(sinval,gmx_mm_inv_pd(cosval));
1063
1064     return tanval;
1065 }
1066
1067
1068
1069 static __m128d
1070 gmx_mm_asin_pd(__m128d x)
1071 {
1072     /* Same algorithm as cephes library */
1073     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
1074     const __m128d limit1    = _mm_set1_pd(0.625);
1075     const __m128d limit2    = _mm_set1_pd(1e-8);
1076     const __m128d one       = _mm_set1_pd(1.0);
1077     const __m128d halfpi    = _mm_set1_pd(M_PI/2.0);
1078     const __m128d quarterpi = _mm_set1_pd(M_PI/4.0);
1079     const __m128d morebits  = _mm_set1_pd(6.123233995736765886130e-17);
1080
1081     const __m128d P5        = _mm_set1_pd(4.253011369004428248960e-3);
1082     const __m128d P4        = _mm_set1_pd(-6.019598008014123785661e-1);
1083     const __m128d P3        = _mm_set1_pd(5.444622390564711410273e0);
1084     const __m128d P2        = _mm_set1_pd(-1.626247967210700244449e1);
1085     const __m128d P1        = _mm_set1_pd(1.956261983317594739197e1);
1086     const __m128d P0        = _mm_set1_pd(-8.198089802484824371615e0);
1087
1088     const __m128d Q4        = _mm_set1_pd(-1.474091372988853791896e1);
1089     const __m128d Q3        = _mm_set1_pd(7.049610280856842141659e1);
1090     const __m128d Q2        = _mm_set1_pd(-1.471791292232726029859e2);
1091     const __m128d Q1        = _mm_set1_pd(1.395105614657485689735e2);
1092     const __m128d Q0        = _mm_set1_pd(-4.918853881490881290097e1);
1093
1094     const __m128d R4        = _mm_set1_pd(2.967721961301243206100e-3);
1095     const __m128d R3        = _mm_set1_pd(-5.634242780008963776856e-1);
1096     const __m128d R2        = _mm_set1_pd(6.968710824104713396794e0);
1097     const __m128d R1        = _mm_set1_pd(-2.556901049652824852289e1);
1098     const __m128d R0        = _mm_set1_pd(2.853665548261061424989e1);
1099
1100     const __m128d S3        = _mm_set1_pd(-2.194779531642920639778e1);
1101     const __m128d S2        = _mm_set1_pd(1.470656354026814941758e2);
1102     const __m128d S1        = _mm_set1_pd(-3.838770957603691357202e2);
1103     const __m128d S0        = _mm_set1_pd(3.424398657913078477438e2);
1104
1105     __m128d sign;
1106     __m128d mask;
1107     __m128d xabs;
1108     __m128d zz,ww,z,q,w,y,zz2,ww2;
1109     __m128d PA,PB;
1110     __m128d QA,QB;
1111     __m128d RA,RB;
1112     __m128d SA,SB;
1113     __m128d nom,denom;
1114
1115     sign  = _mm_andnot_pd(signmask,x);
1116     xabs  = _mm_and_pd(x,signmask);
1117
1118     mask  = _mm_cmpgt_pd(xabs,limit1);
1119
1120     zz    = _mm_sub_pd(one,xabs);
1121     ww    = _mm_mul_pd(xabs,xabs);
1122     zz2   = _mm_mul_pd(zz,zz);
1123     ww2   = _mm_mul_pd(ww,ww);
1124
1125     /* R */
1126     RA    = _mm_macc_pd(R4,zz2,R2);
1127     RB    = _mm_macc_pd(R3,zz2,R1);
1128     RA    = _mm_macc_pd(RA,zz2,R0);
1129     RA    = _mm_macc_pd(RB,zz,RA);
1130
1131     /* S, SA = zz2 */
1132     SB    = _mm_macc_pd(S3,zz2,S1);
1133     SA    = _mm_add_pd(zz2,S2);
1134     SA    = _mm_macc_pd(SA,zz2,S0);
1135     SA    = _mm_macc_pd(SB,zz,SA);
1136
1137     /* P */
1138     PA    = _mm_macc_pd(P5,ww2,P3);
1139     PB    = _mm_macc_pd(P4,ww2,P2);
1140     PA    = _mm_macc_pd(PA,ww2,P1);
1141     PB    = _mm_macc_pd(PB,ww2,P0);
1142     PA    = _mm_macc_pd(PA,ww,PB);
1143
1144     /* Q, QA = ww2 */
1145     QB    = _mm_macc_pd(Q4,ww2,Q2);
1146     QA    = _mm_add_pd(ww2,Q3);
1147     QA    = _mm_macc_pd(QA,ww2,Q1);
1148     QB    = _mm_macc_pd(QB,ww2,Q0);
1149     QA    = _mm_macc_pd(QA,ww,QB);
1150
1151     RA    = _mm_mul_pd(RA,zz);
1152     PA    = _mm_mul_pd(PA,ww);
1153
1154     nom   = _mm_blendv_pd( PA,RA,mask );
1155     denom = _mm_blendv_pd( QA,SA,mask );
1156
1157     q     = _mm_mul_pd( nom, gmx_mm_inv_pd(denom) );
1158
1159     zz    = _mm_add_pd(zz,zz);
1160     zz    = gmx_mm_sqrt_pd(zz);
1161     z     = _mm_sub_pd(quarterpi,zz);
1162     zz    = _mm_mul_pd(zz,q);
1163     zz    = _mm_sub_pd(zz,morebits);
1164     z     = _mm_sub_pd(z,zz);
1165     z     = _mm_add_pd(z,quarterpi);
1166
1167     w     = _mm_macc_pd(xabs,q,xabs);
1168
1169     z     = _mm_blendv_pd( w,z,mask );
1170
1171     mask  = _mm_cmpgt_pd(xabs,limit2);
1172     z     = _mm_blendv_pd( xabs,z,mask );
1173
1174     z = _mm_xor_pd(z,sign);
1175
1176     return z;
1177 }
1178
1179
1180 static __m128d
1181 gmx_mm_acos_pd(__m128d x)
1182 {
1183     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
1184     const __m128d one        = _mm_set1_pd(1.0);
1185     const __m128d half       = _mm_set1_pd(0.5);
1186     const __m128d pi         = _mm_set1_pd(M_PI);
1187     const __m128d quarterpi0 = _mm_set1_pd(7.85398163397448309616e-1);
1188     const __m128d quarterpi1 = _mm_set1_pd(6.123233995736765886130e-17);
1189
1190
1191     __m128d mask1;
1192
1193     __m128d z,z1,z2;
1194
1195     mask1 = _mm_cmpgt_pd(x,half);
1196     z1    = _mm_mul_pd(half,_mm_sub_pd(one,x));
1197     z1    = gmx_mm_sqrt_pd(z1);
1198     z     = _mm_blendv_pd( x,z1,mask1 );
1199
1200     z     = gmx_mm_asin_pd(z);
1201
1202     z1    = _mm_add_pd(z,z);
1203
1204     z2    = _mm_sub_pd(quarterpi0,z);
1205     z2    = _mm_add_pd(z2,quarterpi1);
1206     z2    = _mm_add_pd(z2,quarterpi0);
1207
1208     z     = _mm_blendv_pd(z2,z1,mask1);
1209
1210     return z;
1211 }
1212
1213 static __m128d
1214 gmx_mm_atan_pd(__m128d x)
1215 {
1216     /* Same algorithm as cephes library */
1217     const __m128d signmask  = gmx_mm_castsi128_pd( _mm_set_epi32(0x7FFFFFFF,0xFFFFFFFF,0x7FFFFFFF,0xFFFFFFFF) );
1218     const __m128d limit1    = _mm_set1_pd(0.66);
1219     const __m128d limit2    = _mm_set1_pd(2.41421356237309504880);
1220     const __m128d quarterpi = _mm_set1_pd(M_PI/4.0);
1221     const __m128d halfpi    = _mm_set1_pd(M_PI/2.0);
1222     const __m128d mone      = _mm_set1_pd(-1.0);
1223     const __m128d morebits1 = _mm_set1_pd(0.5*6.123233995736765886130E-17);
1224     const __m128d morebits2 = _mm_set1_pd(6.123233995736765886130E-17);
1225
1226     const __m128d P4        = _mm_set1_pd(-8.750608600031904122785E-1);
1227     const __m128d P3        = _mm_set1_pd(-1.615753718733365076637E1);
1228     const __m128d P2        = _mm_set1_pd(-7.500855792314704667340E1);
1229     const __m128d P1        = _mm_set1_pd(-1.228866684490136173410E2);
1230     const __m128d P0        = _mm_set1_pd(-6.485021904942025371773E1);
1231
1232     const __m128d Q4        = _mm_set1_pd(2.485846490142306297962E1);
1233     const __m128d Q3        = _mm_set1_pd(1.650270098316988542046E2);
1234     const __m128d Q2        = _mm_set1_pd(4.328810604912902668951E2);
1235     const __m128d Q1        = _mm_set1_pd(4.853903996359136964868E2);
1236     const __m128d Q0        = _mm_set1_pd(1.945506571482613964425E2);
1237
1238     __m128d sign;
1239     __m128d mask1,mask2;
1240     __m128d y,t1,t2;
1241     __m128d z,z2;
1242     __m128d P_A,P_B,Q_A,Q_B;
1243
1244     sign   = _mm_andnot_pd(signmask,x);
1245     x      = _mm_and_pd(x,signmask);
1246
1247     mask1  = _mm_cmpgt_pd(x,limit1);
1248     mask2  = _mm_cmpgt_pd(x,limit2);
1249
1250     t1     = _mm_mul_pd(_mm_add_pd(x,mone),gmx_mm_inv_pd(_mm_sub_pd(x,mone)));
1251     t2     = _mm_mul_pd(mone,gmx_mm_inv_pd(x));
1252
1253     y      = _mm_and_pd(mask1,quarterpi);
1254     y      = _mm_or_pd( _mm_and_pd(mask2,halfpi) , _mm_andnot_pd(mask2,y) );
1255
1256     x      = _mm_or_pd( _mm_and_pd(mask1,t1) , _mm_andnot_pd(mask1,x) );
1257     x      = _mm_or_pd( _mm_and_pd(mask2,t2) , _mm_andnot_pd(mask2,x) );
1258
1259     z      = _mm_mul_pd(x,x);
1260     z2     = _mm_mul_pd(z,z);
1261
1262     P_A    = _mm_macc_pd(P4,z2,P2);
1263     P_B    = _mm_macc_pd(P3,z2,P1);
1264     P_A    = _mm_macc_pd(P_A,z2,P0);
1265     P_A    = _mm_macc_pd(P_B,z,P_A);
1266
1267     /* Q_A = z2 */
1268     Q_B    = _mm_macc_pd(Q4,z2,Q2);
1269     Q_A    = _mm_add_pd(z2,Q3);
1270     Q_A    = _mm_macc_pd(Q_A,z2,Q1);
1271     Q_B    = _mm_macc_pd(Q_B,z2,Q0);
1272     Q_A    = _mm_macc_pd(Q_A,z,Q_B);
1273
1274     z      = _mm_mul_pd(z,P_A);
1275     z      = _mm_mul_pd(z,gmx_mm_inv_pd(Q_A));
1276     z      = _mm_macc_pd(z,x,x);
1277
1278     t1     = _mm_and_pd(mask1,morebits1);
1279     t1     = _mm_or_pd( _mm_and_pd(mask2,morebits2) , _mm_andnot_pd(mask2,t1) );
1280
1281     z      = _mm_add_pd(z,t1);
1282     y      = _mm_add_pd(y,z);
1283
1284     y      = _mm_xor_pd(y,sign);
1285
1286     return y;
1287 }
1288
1289
1290 static __m128d
1291 gmx_mm_atan2_pd(__m128d y, __m128d x)
1292 {
1293     const __m128d pi          = _mm_set1_pd(M_PI);
1294     const __m128d minuspi     = _mm_set1_pd(-M_PI);
1295     const __m128d halfpi      = _mm_set1_pd(M_PI/2.0);
1296     const __m128d minushalfpi = _mm_set1_pd(-M_PI/2.0);
1297
1298     __m128d z,z1,z3,z4;
1299     __m128d w;
1300     __m128d maskx_lt,maskx_eq;
1301     __m128d masky_lt,masky_eq;
1302     __m128d mask1,mask2,mask3,mask4,maskall;
1303
1304     maskx_lt  = _mm_cmplt_pd(x,_mm_setzero_pd());
1305     masky_lt  = _mm_cmplt_pd(y,_mm_setzero_pd());
1306     maskx_eq  = _mm_cmpeq_pd(x,_mm_setzero_pd());
1307     masky_eq  = _mm_cmpeq_pd(y,_mm_setzero_pd());
1308
1309     z         = _mm_mul_pd(y,gmx_mm_inv_pd(x));
1310     z         = gmx_mm_atan_pd(z);
1311
1312     mask1     = _mm_and_pd(maskx_eq,masky_lt);
1313     mask2     = _mm_andnot_pd(maskx_lt,masky_eq);
1314     mask3     = _mm_andnot_pd( _mm_or_pd(masky_lt,masky_eq) , maskx_eq);
1315     mask4     = _mm_and_pd(masky_eq,maskx_lt);
1316
1317     maskall   = _mm_or_pd( _mm_or_pd(mask1,mask2), _mm_or_pd(mask3,mask4) );
1318
1319     z         = _mm_andnot_pd(maskall,z);
1320     z1        = _mm_and_pd(mask1,minushalfpi);
1321     z3        = _mm_and_pd(mask3,halfpi);
1322     z4        = _mm_and_pd(mask4,pi);
1323
1324     z         = _mm_or_pd( _mm_or_pd(z,z1), _mm_or_pd(z3,z4) );
1325
1326     w         = _mm_blendv_pd(pi,minuspi,masky_lt);
1327     w         = _mm_and_pd(w,maskx_lt);
1328
1329     w         = _mm_andnot_pd(maskall,w);
1330
1331     z         = _mm_add_pd(z,w);
1332
1333     return z;
1334 }
1335
1336 #endif /*_gmx_math_x86_avx_128_fma_double_h_ */