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