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