?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Added .gitignore for kernel generation
[gromacs.git] / include / gmx_math_x86_avx_256_single.h
blobb25b644ae868057f309d4ba3cc5d3c5f190008c9
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This file is part of GROMACS.
5 * Copyright (c) 2012-
6 *
7 * Written by the Gromacs development team under coordination of
8 * David van der Spoel, Berk Hess, and Erik Lindahl.
9 *
10 * This library is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2
13 * of the License, or (at your option) any later version.
14 *
15 * To help us fund GROMACS development, we humbly ask that you cite
16 * the research papers on the package. Check out http://www.gromacs.org
17 *
18 * And Hey:
19 * Gnomes, ROck Monsters And Chili Sauce
20 */
21 #ifndef _gmx_math_x86_avx_256_single_h_
22 #define _gmx_math_x86_avx_256_single_h_
23
24 #include "gmx_x86_avx_256.h"
25
26
27 #ifndef M_PI
28 # define M_PI 3.14159265358979323846264338327950288
29 #endif
30
31
32
33 /************************
34 * *
35 * Simple math routines *
36 * *
37 ************************/
38
39 /* 1.0/sqrt(x), 256-bit wide version */
40 static gmx_inline __m256
41 gmx_mm256_invsqrt_ps(__m256 x)
42 {
43 const __m256 half = _mm256_set1_ps(0.5f);
44 const __m256 three = _mm256_set1_ps(3.0f);
45
46 __m256 lu = _mm256_rsqrt_ps(x);
47
48 return _mm256_mul_ps(half,_mm256_mul_ps(_mm256_sub_ps(three,_mm256_mul_ps(_mm256_mul_ps(lu,lu),x)),lu));
49 }
50
51 /* 1.0/sqrt(x), 128-bit wide version */
52 static gmx_inline __m128
53 gmx_mm_invsqrt_ps(__m128 x)
54 {
55 const __m128 half = _mm_set_ps(0.5,0.5,0.5,0.5);
56 const __m128 three = _mm_set_ps(3.0,3.0,3.0,3.0);
57
58 __m128 lu = _mm_rsqrt_ps(x);
59
60 return _mm_mul_ps(half,_mm_mul_ps(_mm_sub_ps(three,_mm_mul_ps(_mm_mul_ps(lu,lu),x)),lu));
61 }
62
63
64 /* sqrt(x) (256 bit) - Do NOT use this (but rather invsqrt) if you actually need 1.0/sqrt(x) */
65 static gmx_inline __m256
66 gmx_mm256_sqrt_ps(__m256 x)
67 {
68 __m256 mask;
69 __m256 res;
70
71 mask = _mm256_cmp_ps(x,_mm256_setzero_ps(),_CMP_EQ_OQ);
72 res = _mm256_andnot_ps(mask,gmx_mm256_invsqrt_ps(x));
73
74 res = _mm256_mul_ps(x,res);
75
76 return res;
77 }
78
79 /* sqrt(x) (128 bit) - Do NOT use this (but rather invsqrt) if you actually need 1.0/sqrt(x) */
80 static gmx_inline __m128
81 gmx_mm_sqrt_ps(__m128 x)
82 {
83 __m128 mask;
84 __m128 res;
85
86 mask = _mm_cmpeq_ps(x,_mm_setzero_ps());
87 res = _mm_andnot_ps(mask,gmx_mm_invsqrt_ps(x));
88
89 res = _mm_mul_ps(x,res);
90
91 return res;
92 }
93
94
95 /* 1.0/x, 256-bit wide */
96 static gmx_inline __m256
97 gmx_mm256_inv_ps(__m256 x)
98 {
99 const __m256 two = _mm256_set1_ps(2.0f);
100
101 __m256 lu = _mm256_rcp_ps(x);
102
103 return _mm256_mul_ps(lu,_mm256_sub_ps(two,_mm256_mul_ps(lu,x)));
104 }
105
106 /* 1.0/x, 128-bit wide */
107 static gmx_inline __m128
108 gmx_mm_inv_ps(__m128 x)
109 {
110 const __m128 two = _mm_set_ps(2.0f,2.0f,2.0f,2.0f);
111
112 __m128 lu = _mm_rcp_ps(x);
113
114 return _mm_mul_ps(lu,_mm_sub_ps(two,_mm_mul_ps(lu,x)));
115 }
116
117
118 static gmx_inline __m256
119 gmx_mm256_abs_ps(__m256 x)
120 {
121 const __m256 signmask = _mm256_castsi256_ps( _mm256_set1_epi32(0x7FFFFFFF) );
122
123 return _mm256_and_ps(x,signmask);
124 }
125
126 static gmx_inline __m128
127 gmx_mm_abs_ps(__m128 x)
128 {
129 const __m128 signmask = gmx_mm_castsi128_ps( _mm_set1_epi32(0x7FFFFFFF) );
130
131 return _mm_and_ps(x,signmask);
132 }
133
134
135 static __m256
136 gmx_mm256_log_ps(__m256 x)
137 {
138 const __m256 expmask = _mm256_castsi256_ps( _mm256_set1_epi32(0x7F800000) );
139 const __m128i expbase_m1 = _mm_set1_epi32(127-1); /* We want non-IEEE format */
140 const __m256 half = _mm256_set1_ps(0.5f);
141 const __m256 one = _mm256_set1_ps(1.0f);
142 const __m256 invsq2 = _mm256_set1_ps(1.0f/sqrt(2.0f));
143 const __m256 corr1 = _mm256_set1_ps(-2.12194440e-4f);
144 const __m256 corr2 = _mm256_set1_ps(0.693359375f);
145
146 const __m256 CA_1 = _mm256_set1_ps(0.070376836292f);
147 const __m256 CB_0 = _mm256_set1_ps(1.6714950086782716f);
148 const __m256 CB_1 = _mm256_set1_ps(-2.452088066061482f);
149 const __m256 CC_0 = _mm256_set1_ps(1.5220770854701728f);
150 const __m256 CC_1 = _mm256_set1_ps(-1.3422238433233642f);
151 const __m256 CD_0 = _mm256_set1_ps(1.386218787509749f);
152 const __m256 CD_1 = _mm256_set1_ps(0.35075468953796346f);
153 const __m256 CE_0 = _mm256_set1_ps(1.3429983063133937f);
154 const __m256 CE_1 = _mm256_set1_ps(1.807420826584643f);
155
156 __m256 fexp,fexp1;
157 __m256i iexp;
158 __m128i iexp128a,iexp128b;
159 __m256 mask;
160 __m256i imask;
161 __m128i imask128a,imask128b;
162 __m256 x1,x2;
163 __m256 y;
164 __m256 pA,pB,pC,pD,pE,tB,tC,tD,tE;
165
166 /* Separate x into exponent and mantissa, with a mantissa in the range [0.5..1[ (not IEEE754 standard!) */
167 fexp = _mm256_and_ps(x,expmask);
168 iexp = _mm256_castps_si256(fexp);
169
170 iexp128b = _mm256_extractf128_si256(iexp,0x1);
171 iexp128a = _mm256_castsi256_si128(iexp);
172
173 iexp128a = _mm_srli_epi32(iexp128a,23);
174 iexp128b = _mm_srli_epi32(iexp128b,23);
175 iexp128a = _mm_sub_epi32(iexp128a,expbase_m1);
176 iexp128b = _mm_sub_epi32(iexp128b,expbase_m1);
177
178 x = _mm256_andnot_ps(expmask,x);
179 x = _mm256_or_ps(x,one);
180 x = _mm256_mul_ps(x,half);
181
182 mask = _mm256_cmp_ps(x,invsq2,_CMP_LT_OQ);
183
184 x = _mm256_add_ps(x,_mm256_and_ps(mask,x));
185 x = _mm256_sub_ps(x,one);
186
187 imask = _mm256_castps_si256(mask);
188
189 imask128b = _mm256_extractf128_si256(imask,0x1);
190 imask128a = _mm256_castsi256_si128(imask);
191
192 iexp128a = _mm_add_epi32(iexp128a,imask128a);
193 iexp128b = _mm_add_epi32(iexp128b,imask128b);
194
195 iexp = _mm256_castsi128_si256(iexp128a);
196 iexp = _mm256_insertf128_si256(iexp,iexp128b,0x1);
197
198 x2 = _mm256_mul_ps(x,x);
199
200 pA = _mm256_mul_ps(CA_1,x);
201 pB = _mm256_mul_ps(CB_1,x);
202 pC = _mm256_mul_ps(CC_1,x);
203 pD = _mm256_mul_ps(CD_1,x);
204 pE = _mm256_mul_ps(CE_1,x);
205 tB = _mm256_add_ps(CB_0,x2);
206 tC = _mm256_add_ps(CC_0,x2);
207 tD = _mm256_add_ps(CD_0,x2);
208 tE = _mm256_add_ps(CE_0,x2);
209 pB = _mm256_add_ps(pB,tB);
210 pC = _mm256_add_ps(pC,tC);
211 pD = _mm256_add_ps(pD,tD);
212 pE = _mm256_add_ps(pE,tE);
213
214 pA = _mm256_mul_ps(pA,pB);
215 pC = _mm256_mul_ps(pC,pD);
216 pE = _mm256_mul_ps(pE,x2);
217 pA = _mm256_mul_ps(pA,pC);
218 y = _mm256_mul_ps(pA,pE);
219
220 fexp = _mm256_cvtepi32_ps(iexp);
221 y = _mm256_add_ps(y,_mm256_mul_ps(fexp,corr1));
222
223 y = _mm256_sub_ps(y, _mm256_mul_ps(half,x2));
224 x2 = _mm256_add_ps(x,y);
225
226 x2 = _mm256_add_ps(x2,_mm256_mul_ps(fexp,corr2));
227
228 return x2;
229 }
230
231
232 static __m128
233 gmx_mm_log_ps(__m128 x)
234 {
235 /* Same algorithm as cephes library */
236 const __m128 expmask = gmx_mm_castsi128_ps( _mm_set_epi32(0x7F800000, 0x7F800000, 0x7F800000, 0x7F800000) );
237 const __m128i expbase_m1 = _mm_set1_epi32(127-1); /* We want non-IEEE format */
238 const __m128 half = _mm_set1_ps(0.5f);
239 const __m128 one = _mm_set1_ps(1.0f);
240 const __m128 invsq2 = _mm_set1_ps(1.0f/sqrt(2.0f));
241 const __m128 corr1 = _mm_set1_ps(-2.12194440e-4f);
242 const __m128 corr2 = _mm_set1_ps(0.693359375f);
243
244 const __m128 CA_1 = _mm_set1_ps(0.070376836292f);
245 const __m128 CB_0 = _mm_set1_ps(1.6714950086782716f);
246 const __m128 CB_1 = _mm_set1_ps(-2.452088066061482f);
247 const __m128 CC_0 = _mm_set1_ps(1.5220770854701728f);
248 const __m128 CC_1 = _mm_set1_ps(-1.3422238433233642f);
249 const __m128 CD_0 = _mm_set1_ps(1.386218787509749f);
250 const __m128 CD_1 = _mm_set1_ps(0.35075468953796346f);
251 const __m128 CE_0 = _mm_set1_ps(1.3429983063133937f);
252 const __m128 CE_1 = _mm_set1_ps(1.807420826584643f);
253
254 __m128 fexp,fexp1;
255 __m128i iexp;
256 __m128 mask;
257 __m128 x1,x2;
258 __m128 y;
259 __m128 pA,pB,pC,pD,pE,tB,tC,tD,tE;
260
261 /* Separate x into exponent and mantissa, with a mantissa in the range [0.5..1[ (not IEEE754 standard!) */
262 fexp = _mm_and_ps(x,expmask);
263 iexp = gmx_mm_castps_si128(fexp);
264 iexp = _mm_srli_epi32(iexp,23);
265 iexp = _mm_sub_epi32(iexp,expbase_m1);
266
267 x = _mm_andnot_ps(expmask,x);
268 x = _mm_or_ps(x,one);
269 x = _mm_mul_ps(x,half);
270
271 mask = _mm_cmplt_ps(x,invsq2);
272
273 x = _mm_add_ps(x,_mm_and_ps(mask,x));
274 x = _mm_sub_ps(x,one);
275 iexp = _mm_add_epi32(iexp,gmx_mm_castps_si128(mask)); /* 0xFFFFFFFF = -1 as int */
276
277 x2 = _mm_mul_ps(x,x);
278
279 pA = _mm_mul_ps(CA_1,x);
280 pB = _mm_mul_ps(CB_1,x);
281 pC = _mm_mul_ps(CC_1,x);
282 pD = _mm_mul_ps(CD_1,x);
283 pE = _mm_mul_ps(CE_1,x);
284 tB = _mm_add_ps(CB_0,x2);
285 tC = _mm_add_ps(CC_0,x2);
286 tD = _mm_add_ps(CD_0,x2);
287 tE = _mm_add_ps(CE_0,x2);
288 pB = _mm_add_ps(pB,tB);
289 pC = _mm_add_ps(pC,tC);
290 pD = _mm_add_ps(pD,tD);
291 pE = _mm_add_ps(pE,tE);
292
293 pA = _mm_mul_ps(pA,pB);
294 pC = _mm_mul_ps(pC,pD);
295 pE = _mm_mul_ps(pE,x2);
296 pA = _mm_mul_ps(pA,pC);
297 y = _mm_mul_ps(pA,pE);
298
299 fexp = _mm_cvtepi32_ps(iexp);
300 y = _mm_add_ps(y,_mm_mul_ps(fexp,corr1));
301
302 y = _mm_sub_ps(y, _mm_mul_ps(half,x2));
303 x2 = _mm_add_ps(x,y);
304
305 x2 = _mm_add_ps(x2,_mm_mul_ps(fexp,corr2));
306
307 return x2;
308 }
309
310
311 /*
312 * 2^x function, 256-bit wide
313 *
314 * The 2^w term is calculated from a (6,0)-th order (no denominator) Minimax polynomia on the interval
315 * [-0.5,0.5]. The coefficiencts of this was derived in Mathematica using the command:
316 *
317 * MiniMaxApproximation[(2^x), {x, {-0.5, 0.5}, 6, 0}, WorkingPrecision -> 15]
318 *
319 * The largest-magnitude exponent we can represent in IEEE single-precision binary format
320 * is 2^-126 for small numbers and 2^127 for large ones. To avoid wrap-around problems, we set the
321 * result to zero if the argument falls outside this range. For small numbers this is just fine, but
322 * for large numbers you could be fancy and return the smallest/largest IEEE single-precision
323 * number instead. That would take a few extra cycles and not really help, since something is
324 * wrong if you are using single precision to work with numbers that cannot really be represented
325 * in single precision.
326 *
327 * The accuracy is at least 23 bits.
328 */
329 static __m256
330 gmx_mm256_exp2_ps(__m256 x)
331 {
332 /* Lower bound: Disallow numbers that would lead to an IEEE fp exponent reaching +-127. */
333 const __m256 arglimit = _mm256_set1_ps(126.0f);
334
335 const __m128i expbase = _mm_set1_epi32(127);
336 const __m256 CC6 = _mm256_set1_ps(1.535336188319500E-004);
337 const __m256 CC5 = _mm256_set1_ps(1.339887440266574E-003);
338 const __m256 CC4 = _mm256_set1_ps(9.618437357674640E-003);
339 const __m256 CC3 = _mm256_set1_ps(5.550332471162809E-002);
340 const __m256 CC2 = _mm256_set1_ps(2.402264791363012E-001);
341 const __m256 CC1 = _mm256_set1_ps(6.931472028550421E-001);
342 const __m256 CC0 = _mm256_set1_ps(1.0f);
343
344 __m256 p0,p1;
345 __m256 valuemask;
346 __m256i iexppart;
347 __m128i iexppart128a,iexppart128b;
348 __m256 fexppart;
349 __m256 intpart;
350 __m256 x2;
351
352
353 iexppart = _mm256_cvtps_epi32(x);
354 intpart = _mm256_round_ps(x,_MM_FROUND_TO_NEAREST_INT);
355
356 iexppart128b = _mm256_extractf128_si256(iexppart,0x1);
357 iexppart128a = _mm256_castsi256_si128(iexppart);
358
359 iexppart128a = _mm_slli_epi32(_mm_add_epi32(iexppart128a,expbase),23);
360 iexppart128b = _mm_slli_epi32(_mm_add_epi32(iexppart128b,expbase),23);
361
362 iexppart = _mm256_castsi128_si256(iexppart128a);
363 iexppart = _mm256_insertf128_si256(iexppart,iexppart128b,0x1);
364 valuemask = _mm256_cmp_ps(arglimit,gmx_mm256_abs_ps(x),_CMP_GE_OQ);
365 fexppart = _mm256_and_ps(valuemask,_mm256_castsi256_ps(iexppart));
366
367 x = _mm256_sub_ps(x,intpart);
368 x2 = _mm256_mul_ps(x,x);
369
370 p0 = _mm256_mul_ps(CC6,x2);
371 p1 = _mm256_mul_ps(CC5,x2);
372 p0 = _mm256_add_ps(p0,CC4);
373 p1 = _mm256_add_ps(p1,CC3);
374 p0 = _mm256_mul_ps(p0,x2);
375 p1 = _mm256_mul_ps(p1,x2);
376 p0 = _mm256_add_ps(p0,CC2);
377 p1 = _mm256_add_ps(p1,CC1);
378 p0 = _mm256_mul_ps(p0,x2);
379 p1 = _mm256_mul_ps(p1,x);
380 p0 = _mm256_add_ps(p0,CC0);
381 p0 = _mm256_add_ps(p0,p1);
382 x = _mm256_mul_ps(p0,fexppart);
383
384 return x;
385 }
386
387
388 /* 2^x, 128 bit wide */
389 static __m128
390 gmx_mm_exp2_ps(__m128 x)
391 {
392 /* Lower bound: We do not allow numbers that would lead to an IEEE fp representation exponent smaller than -126. */
393 const __m128 arglimit = _mm_set1_ps(126.0f);
394
395 const __m128i expbase = _mm_set1_epi32(127);
396 const __m128 CA6 = _mm_set1_ps(1.535336188319500E-004);
397 const __m128 CA5 = _mm_set1_ps(1.339887440266574E-003);
398 const __m128 CA4 = _mm_set1_ps(9.618437357674640E-003);
399 const __m128 CA3 = _mm_set1_ps(5.550332471162809E-002);
400 const __m128 CA2 = _mm_set1_ps(2.402264791363012E-001);
401 const __m128 CA1 = _mm_set1_ps(6.931472028550421E-001);
402 const __m128 CA0 = _mm_set1_ps(1.0f);
403
404 __m128 valuemask;
405 __m128i iexppart;
406 __m128 fexppart;
407 __m128 intpart;
408 __m128 x2;
409 __m128 p0,p1;
410
411 iexppart = _mm_cvtps_epi32(x);
412 intpart = _mm_round_ps(x,_MM_FROUND_TO_NEAREST_INT);
413 iexppart = _mm_slli_epi32(_mm_add_epi32(iexppart,expbase),23);
414 valuemask = _mm_cmpge_ps(arglimit,gmx_mm_abs_ps(x));
415 fexppart = _mm_and_ps(valuemask,gmx_mm_castsi128_ps(iexppart));
416
417 x = _mm_sub_ps(x,intpart);
418 x2 = _mm_mul_ps(x,x);
419
420 p0 = _mm_mul_ps(CA6,x2);
421 p1 = _mm_mul_ps(CA5,x2);
422 p0 = _mm_add_ps(p0,CA4);
423 p1 = _mm_add_ps(p1,CA3);
424 p0 = _mm_mul_ps(p0,x2);
425 p1 = _mm_mul_ps(p1,x2);
426 p0 = _mm_add_ps(p0,CA2);
427 p1 = _mm_add_ps(p1,CA1);
428 p0 = _mm_mul_ps(p0,x2);
429 p1 = _mm_mul_ps(p1,x);
430 p0 = _mm_add_ps(p0,CA0);
431 p0 = _mm_add_ps(p0,p1);
432 x = _mm_mul_ps(p0,fexppart);
433
434 return x;
435 }
436
437
438 /* Exponential function, 256 bit wide. This could be calculated from 2^x as Exp(x)=2^(y),
439 * where y=log2(e)*x, but there will then be a small rounding error since we lose some
440 * precision due to the multiplication. This will then be magnified a lot by the exponential.
441 *
442 * Instead, we calculate the fractional part directly as a minimax approximation of
443 * Exp(z) on [-0.5,0.5]. We use extended precision arithmetics to calculate the fraction
444 * remaining after 2^y, which avoids the precision-loss.
445 * The final result is correct to within 1 LSB over the entire argument range.
446 */
447 static __m256
448 gmx_mm256_exp_ps(__m256 exparg)
449 {
450 const __m256 argscale = _mm256_set1_ps(1.44269504088896341f);
451 /* Lower bound: Disallow numbers that would lead to an IEEE fp exponent reaching +-127. */
452 const __m256 arglimit = _mm256_set1_ps(126.0f);
453 const __m128i expbase = _mm_set1_epi32(127);
454
455 const __m256 invargscale0 = _mm256_set1_ps(0.693359375f);
456 const __m256 invargscale1 = _mm256_set1_ps(-2.12194440e-4f);
457
458 const __m256 CE5 = _mm256_set1_ps(1.9875691500e-4f);
459 const __m256 CE4 = _mm256_set1_ps(1.3981999507e-3f);
460 const __m256 CE3 = _mm256_set1_ps(8.3334519073e-3f);
461 const __m256 CE2 = _mm256_set1_ps(4.1665795894e-2f);
462 const __m256 CE1 = _mm256_set1_ps(1.6666665459e-1f);
463 const __m256 CE0 = _mm256_set1_ps(5.0000001201e-1f);
464 const __m256 one = _mm256_set1_ps(1.0f);
465
466 __m256 exparg2,exp2arg;
467 __m256 pE0,pE1;
468 __m256 valuemask;
469 __m256i iexppart;
470 __m128i iexppart128a,iexppart128b;
471 __m256 fexppart;
472 __m256 intpart;
473
474 exp2arg = _mm256_mul_ps(exparg,argscale);
475
476 iexppart = _mm256_cvtps_epi32(exp2arg);
477 intpart = _mm256_round_ps(exp2arg,_MM_FROUND_TO_NEAREST_INT);
478
479 iexppart128b = _mm256_extractf128_si256(iexppart,0x1);
480 iexppart128a = _mm256_castsi256_si128(iexppart);
481
482 iexppart128a = _mm_slli_epi32(_mm_add_epi32(iexppart128a,expbase),23);
483 iexppart128b = _mm_slli_epi32(_mm_add_epi32(iexppart128b,expbase),23);
484
485 iexppart = _mm256_castsi128_si256(iexppart128a);
486 iexppart = _mm256_insertf128_si256(iexppart,iexppart128b,0x1);
487 valuemask = _mm256_cmp_ps(arglimit,gmx_mm256_abs_ps(exp2arg),_CMP_GE_OQ);
488 fexppart = _mm256_and_ps(valuemask,_mm256_castsi256_ps(iexppart));
489
490 /* Extended precision arithmetics */
491 exparg = _mm256_sub_ps(exparg,_mm256_mul_ps(invargscale0,intpart));
492 exparg = _mm256_sub_ps(exparg,_mm256_mul_ps(invargscale1,intpart));
493
494 exparg2 = _mm256_mul_ps(exparg,exparg);
495
496 pE1 = _mm256_mul_ps(CE5,exparg2);
497 pE0 = _mm256_mul_ps(CE4,exparg2);
498 pE1 = _mm256_add_ps(pE1,CE3);
499 pE0 = _mm256_add_ps(pE0,CE2);
500 pE1 = _mm256_mul_ps(pE1,exparg2);
501 pE0 = _mm256_mul_ps(pE0,exparg2);
502 pE1 = _mm256_add_ps(pE1,CE1);
503 pE0 = _mm256_add_ps(pE0,CE0);
504 pE1 = _mm256_mul_ps(pE1,exparg);
505 pE0 = _mm256_add_ps(pE0,pE1);
506 pE0 = _mm256_mul_ps(pE0,exparg2);
507 exparg = _mm256_add_ps(exparg,one);
508 exparg = _mm256_add_ps(exparg,pE0);
509
510 exparg = _mm256_mul_ps(exparg,fexppart);
511
512 return exparg;
513 }
514
515
516 /* exp(), 128 bit wide. */
517 static __m128
518 gmx_mm_exp_ps(__m128 x)
519 {
520 const __m128 argscale = _mm_set1_ps(1.44269504088896341f);
521 /* Lower bound: Disallow numbers that would lead to an IEEE fp exponent reaching +-127. */
522 const __m128 arglimit = _mm_set1_ps(126.0f);
523 const __m128i expbase = _mm_set1_epi32(127);
524
525 const __m128 invargscale0 = _mm_set1_ps(0.693359375f);
526 const __m128 invargscale1 = _mm_set1_ps(-2.12194440e-4f);
527
528 const __m128 CC5 = _mm_set1_ps(1.9875691500e-4f);
529 const __m128 CC4 = _mm_set1_ps(1.3981999507e-3f);
530 const __m128 CC3 = _mm_set1_ps(8.3334519073e-3f);
531 const __m128 CC2 = _mm_set1_ps(4.1665795894e-2f);
532 const __m128 CC1 = _mm_set1_ps(1.6666665459e-1f);
533 const __m128 CC0 = _mm_set1_ps(5.0000001201e-1f);
534 const __m128 one = _mm_set1_ps(1.0f);
535
536 __m128 y,x2;
537 __m128 p0,p1;
538 __m128 valuemask;
539 __m128i iexppart;
540 __m128 fexppart;
541 __m128 intpart;
542
543 y = _mm_mul_ps(x,argscale);
544
545 iexppart = _mm_cvtps_epi32(y);
546 intpart = _mm_round_ps(y,_MM_FROUND_TO_NEAREST_INT);
547
548 iexppart = _mm_slli_epi32(_mm_add_epi32(iexppart,expbase),23);
549 valuemask = _mm_cmpge_ps(arglimit,gmx_mm_abs_ps(y));
550 fexppart = _mm_and_ps(valuemask,gmx_mm_castsi128_ps(iexppart));
551
552 /* Extended precision arithmetics */
553 x = _mm_sub_ps(x,_mm_mul_ps(invargscale0,intpart));
554 x = _mm_sub_ps(x,_mm_mul_ps(invargscale1,intpart));
555
556 x2 = _mm_mul_ps(x,x);
557
558 p1 = _mm_mul_ps(CC5,x2);
559 p0 = _mm_mul_ps(CC4,x2);
560 p1 = _mm_add_ps(p1,CC3);
561 p0 = _mm_add_ps(p0,CC2);
562 p1 = _mm_mul_ps(p1,x2);
563 p0 = _mm_mul_ps(p0,x2);
564 p1 = _mm_add_ps(p1,CC1);
565 p0 = _mm_add_ps(p0,CC0);
566 p1 = _mm_mul_ps(p1,x);
567 p0 = _mm_add_ps(p0,p1);
568 p0 = _mm_mul_ps(p0,x2);
569 x = _mm_add_ps(x,one);
570 x = _mm_add_ps(x,p0);
571
572 x = _mm_mul_ps(x,fexppart);
573
574 return x;
575 }
576
577
578
579 /* FULL precision erf(), 256-bit wide. Only errors in LSB */
580 static __m256
581 gmx_mm256_erf_ps(__m256 x)
582 {
583 /* Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1] */
584 const __m256 CA6 = _mm256_set1_ps(7.853861353153693e-5f);
585 const __m256 CA5 = _mm256_set1_ps(-8.010193625184903e-4f);
586 const __m256 CA4 = _mm256_set1_ps(5.188327685732524e-3f);
587 const __m256 CA3 = _mm256_set1_ps(-2.685381193529856e-2f);
588 const __m256 CA2 = _mm256_set1_ps(1.128358514861418e-1f);
589 const __m256 CA1 = _mm256_set1_ps(-3.761262582423300e-1f);
590 const __m256 CA0 = _mm256_set1_ps(1.128379165726710f);
591 /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2] */
592 const __m256 CB9 = _mm256_set1_ps(-0.0018629930017603923f);
593 const __m256 CB8 = _mm256_set1_ps(0.003909821287598495f);
594 const __m256 CB7 = _mm256_set1_ps(-0.0052094582210355615f);
595 const __m256 CB6 = _mm256_set1_ps(0.005685614362160572f);
596 const __m256 CB5 = _mm256_set1_ps(-0.0025367682853477272f);
597 const __m256 CB4 = _mm256_set1_ps(-0.010199799682318782f);
598 const __m256 CB3 = _mm256_set1_ps(0.04369575504816542f);
599 const __m256 CB2 = _mm256_set1_ps(-0.11884063474674492f);
600 const __m256 CB1 = _mm256_set1_ps(0.2732120154030589f);
601 const __m256 CB0 = _mm256_set1_ps(0.42758357702025784f);
602 /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19] */
603 const __m256 CC10 = _mm256_set1_ps(-0.0445555913112064f);
604 const __m256 CC9 = _mm256_set1_ps(0.21376355144663348f);
605 const __m256 CC8 = _mm256_set1_ps(-0.3473187200259257f);
606 const __m256 CC7 = _mm256_set1_ps(0.016690861551248114f);
607 const __m256 CC6 = _mm256_set1_ps(0.7560973182491192f);
608 const __m256 CC5 = _mm256_set1_ps(-1.2137903600145787f);
609 const __m256 CC4 = _mm256_set1_ps(0.8411872321232948f);
610 const __m256 CC3 = _mm256_set1_ps(-0.08670413896296343f);
611 const __m256 CC2 = _mm256_set1_ps(-0.27124782687240334f);
612 const __m256 CC1 = _mm256_set1_ps(-0.0007502488047806069f);
613 const __m256 CC0 = _mm256_set1_ps(0.5642114853803148f);
614
615 /* Coefficients for expansion of exp(x) in [0,0.1] */
616 /* CD0 and CD1 are both 1.0, so no need to declare them separately */
617 const __m256 CD2 = _mm256_set1_ps(0.5000066608081202f);
618 const __m256 CD3 = _mm256_set1_ps(0.1664795422874624f);
619 const __m256 CD4 = _mm256_set1_ps(0.04379839977652482f);
620
621 const __m256 sieve = _mm256_castsi256_ps( _mm256_set1_epi32(0xfffff000) );
622 const __m256 signbit = _mm256_castsi256_ps( _mm256_set1_epi32(0x80000000) );
623 const __m256 one = _mm256_set1_ps(1.0f);
624 const __m256 two = _mm256_set1_ps(2.0f);
625
626 __m256 x2,x4,y;
627 __m256 z,q,t,t2,w,w2;
628 __m256 pA0,pA1,pB0,pB1,pC0,pC1;
629 __m256 expmx2,corr;
630 __m256 res_erf,res_erfc,res;
631 __m256 mask;
632
633 /* Calculate erf() */
634 x2 = _mm256_mul_ps(x,x);
635 x4 = _mm256_mul_ps(x2,x2);
636
637 pA0 = _mm256_mul_ps(CA6,x4);
638 pA1 = _mm256_mul_ps(CA5,x4);
639 pA0 = _mm256_add_ps(pA0,CA4);
640 pA1 = _mm256_add_ps(pA1,CA3);
641 pA0 = _mm256_mul_ps(pA0,x4);
642 pA1 = _mm256_mul_ps(pA1,x4);
643 pA0 = _mm256_add_ps(pA0,CA2);
644 pA1 = _mm256_add_ps(pA1,CA1);
645 pA0 = _mm256_mul_ps(pA0,x4);
646 pA1 = _mm256_mul_ps(pA1,x2);
647 pA0 = _mm256_add_ps(pA0,pA1);
648 pA0 = _mm256_add_ps(pA0,CA0);
649
650 res_erf = _mm256_mul_ps(x,pA0);
651
652 /* Calculate erfc */
653
654 y = gmx_mm256_abs_ps(x);
655 t = gmx_mm256_inv_ps(y);
656 w = _mm256_sub_ps(t,one);
657 t2 = _mm256_mul_ps(t,t);
658 w2 = _mm256_mul_ps(w,w);
659 /*
660 * We cannot simply calculate exp(-x2) directly in single precision, since
661 * that will lose a couple of bits of precision due to the multiplication.
662 * Instead, we introduce x=z+w, where the last 12 bits of precision are in w.
663 * Then we get exp(-x2) = exp(-z2)*exp((z-x)*(z+x)).
664 *
665 * The only drawback with this is that it requires TWO separate exponential
666 * evaluations, which would be horrible performance-wise. However, the argument
667 * for the second exp() call is always small, so there we simply use a
668 * low-order minimax expansion on [0,0.1].
669 */
670
671 z = _mm256_and_ps(y,sieve);
672 q = _mm256_mul_ps( _mm256_sub_ps(z,y) , _mm256_add_ps(z,y) );
673
674 corr = _mm256_mul_ps(CD4,q);
675 corr = _mm256_add_ps(corr,CD3);
676 corr = _mm256_mul_ps(corr,q);
677 corr = _mm256_add_ps(corr,CD2);
678 corr = _mm256_mul_ps(corr,q);
679 corr = _mm256_add_ps(corr,one);
680 corr = _mm256_mul_ps(corr,q);
681 corr = _mm256_add_ps(corr,one);
682
683 expmx2 = gmx_mm256_exp_ps( _mm256_or_ps( signbit , _mm256_mul_ps(z,z) ) );
684 expmx2 = _mm256_mul_ps(expmx2,corr);
685
686 pB1 = _mm256_mul_ps(CB9,w2);
687 pB0 = _mm256_mul_ps(CB8,w2);
688 pB1 = _mm256_add_ps(pB1,CB7);
689 pB0 = _mm256_add_ps(pB0,CB6);
690 pB1 = _mm256_mul_ps(pB1,w2);
691 pB0 = _mm256_mul_ps(pB0,w2);
692 pB1 = _mm256_add_ps(pB1,CB5);
693 pB0 = _mm256_add_ps(pB0,CB4);
694 pB1 = _mm256_mul_ps(pB1,w2);
695 pB0 = _mm256_mul_ps(pB0,w2);
696 pB1 = _mm256_add_ps(pB1,CB3);
697 pB0 = _mm256_add_ps(pB0,CB2);
698 pB1 = _mm256_mul_ps(pB1,w2);
699 pB0 = _mm256_mul_ps(pB0,w2);
700 pB1 = _mm256_add_ps(pB1,CB1);
701 pB1 = _mm256_mul_ps(pB1,w);
702 pB0 = _mm256_add_ps(pB0,pB1);
703 pB0 = _mm256_add_ps(pB0,CB0);
704
705 pC0 = _mm256_mul_ps(CC10,t2);
706 pC1 = _mm256_mul_ps(CC9,t2);
707 pC0 = _mm256_add_ps(pC0,CC8);
708 pC1 = _mm256_add_ps(pC1,CC7);
709 pC0 = _mm256_mul_ps(pC0,t2);
710 pC1 = _mm256_mul_ps(pC1,t2);
711 pC0 = _mm256_add_ps(pC0,CC6);
712 pC1 = _mm256_add_ps(pC1,CC5);
713 pC0 = _mm256_mul_ps(pC0,t2);
714 pC1 = _mm256_mul_ps(pC1,t2);
715 pC0 = _mm256_add_ps(pC0,CC4);
716 pC1 = _mm256_add_ps(pC1,CC3);
717 pC0 = _mm256_mul_ps(pC0,t2);
718 pC1 = _mm256_mul_ps(pC1,t2);
719 pC0 = _mm256_add_ps(pC0,CC2);
720 pC1 = _mm256_add_ps(pC1,CC1);
721 pC0 = _mm256_mul_ps(pC0,t2);
722 pC1 = _mm256_mul_ps(pC1,t);
723 pC0 = _mm256_add_ps(pC0,pC1);
724 pC0 = _mm256_add_ps(pC0,CC0);
725 pC0 = _mm256_mul_ps(pC0,t);
726
727 /* SELECT pB0 or pC0 for erfc() */
728 mask = _mm256_cmp_ps(two,y,_CMP_LT_OQ);
729 res_erfc = _mm256_blendv_ps(pB0,pC0,mask);
730 res_erfc = _mm256_mul_ps(res_erfc,expmx2);
731
732 /* erfc(x<0) = 2-erfc(|x|) */
733 mask = _mm256_cmp_ps(x,_mm256_setzero_ps(),_CMP_LT_OQ);
734 res_erfc = _mm256_blendv_ps(res_erfc,_mm256_sub_ps(two,res_erfc),mask);
735
736 /* Select erf() or erfc() */
737 mask = _mm256_cmp_ps(y,_mm256_set1_ps(0.75f),_CMP_LT_OQ);
738 res = _mm256_blendv_ps(_mm256_sub_ps(one,res_erfc),res_erf,mask);
739
740 return res;
741 }
742
743
744 /* erf(), 128 bit wide */
745 static __m128
746 gmx_mm_erf_ps(__m128 x)
747 {
748 /* Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1] */
749 const __m128 CA6 = _mm_set1_ps(7.853861353153693e-5f);
750 const __m128 CA5 = _mm_set1_ps(-8.010193625184903e-4f);
751 const __m128 CA4 = _mm_set1_ps(5.188327685732524e-3f);
752 const __m128 CA3 = _mm_set1_ps(-2.685381193529856e-2f);
753 const __m128 CA2 = _mm_set1_ps(1.128358514861418e-1f);
754 const __m128 CA1 = _mm_set1_ps(-3.761262582423300e-1f);
755 const __m128 CA0 = _mm_set1_ps(1.128379165726710f);
756 /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2] */
757 const __m128 CB9 = _mm_set1_ps(-0.0018629930017603923f);
758 const __m128 CB8 = _mm_set1_ps(0.003909821287598495f);
759 const __m128 CB7 = _mm_set1_ps(-0.0052094582210355615f);
760 const __m128 CB6 = _mm_set1_ps(0.005685614362160572f);
761 const __m128 CB5 = _mm_set1_ps(-0.0025367682853477272f);
762 const __m128 CB4 = _mm_set1_ps(-0.010199799682318782f);
763 const __m128 CB3 = _mm_set1_ps(0.04369575504816542f);
764 const __m128 CB2 = _mm_set1_ps(-0.11884063474674492f);
765 const __m128 CB1 = _mm_set1_ps(0.2732120154030589f);
766 const __m128 CB0 = _mm_set1_ps(0.42758357702025784f);
767 /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19] */
768 const __m128 CC10 = _mm_set1_ps(-0.0445555913112064f);
769 const __m128 CC9 = _mm_set1_ps(0.21376355144663348f);
770 const __m128 CC8 = _mm_set1_ps(-0.3473187200259257f);
771 const __m128 CC7 = _mm_set1_ps(0.016690861551248114f);
772 const __m128 CC6 = _mm_set1_ps(0.7560973182491192f);
773 const __m128 CC5 = _mm_set1_ps(-1.2137903600145787f);
774 const __m128 CC4 = _mm_set1_ps(0.8411872321232948f);
775 const __m128 CC3 = _mm_set1_ps(-0.08670413896296343f);
776 const __m128 CC2 = _mm_set1_ps(-0.27124782687240334f);
777 const __m128 CC1 = _mm_set1_ps(-0.0007502488047806069f);
778 const __m128 CC0 = _mm_set1_ps(0.5642114853803148f);
779
780 /* Coefficients for expansion of exp(x) in [0,0.1] */
781 /* CD0 and CD1 are both 1.0, so no need to declare them separately */
782 const __m128 CD2 = _mm_set1_ps(0.5000066608081202f);
783 const __m128 CD3 = _mm_set1_ps(0.1664795422874624f);
784 const __m128 CD4 = _mm_set1_ps(0.04379839977652482f);
785
786 const __m128 sieve = gmx_mm_castsi128_ps( _mm_set1_epi32(0xfffff000) );
787 const __m128 signbit = gmx_mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
788 const __m128 one = _mm_set1_ps(1.0f);
789 const __m128 two = _mm_set1_ps(2.0f);
790
791 __m128 x2,x4,y;
792 __m128 z,q,t,t2,w,w2;
793 __m128 pA0,pA1,pB0,pB1,pC0,pC1;
794 __m128 expmx2,corr;
795 __m128 res_erf,res_erfc,res;
796 __m128 mask;
797
798 /* Calculate erf() */
799 x2 = _mm_mul_ps(x,x);
800 x4 = _mm_mul_ps(x2,x2);
801
802 pA0 = _mm_mul_ps(CA6,x4);
803 pA1 = _mm_mul_ps(CA5,x4);
804 pA0 = _mm_add_ps(pA0,CA4);
805 pA1 = _mm_add_ps(pA1,CA3);
806 pA0 = _mm_mul_ps(pA0,x4);
807 pA1 = _mm_mul_ps(pA1,x4);
808 pA0 = _mm_add_ps(pA0,CA2);
809 pA1 = _mm_add_ps(pA1,CA1);
810 pA0 = _mm_mul_ps(pA0,x4);
811 pA1 = _mm_mul_ps(pA1,x2);
812 pA0 = _mm_add_ps(pA0,pA1);
813 pA0 = _mm_add_ps(pA0,CA0);
814
815 res_erf = _mm_mul_ps(x,pA0);
816
817 /* Calculate erfc */
818
819 y = gmx_mm_abs_ps(x);
820 t = gmx_mm_inv_ps(y);
821 w = _mm_sub_ps(t,one);
822 t2 = _mm_mul_ps(t,t);
823 w2 = _mm_mul_ps(w,w);
824 /*
825 * We cannot simply calculate exp(-x2) directly in single precision, since
826 * that will lose a couple of bits of precision due to the multiplication.
827 * Instead, we introduce x=z+w, where the last 12 bits of precision are in w.
828 * Then we get exp(-x2) = exp(-z2)*exp((z-x)*(z+x)).
829 *
830 * The only drawback with this is that it requires TWO separate exponential
831 * evaluations, which would be horrible performance-wise. However, the argument
832 * for the second exp() call is always small, so there we simply use a
833 * low-order minimax expansion on [0,0.1].
834 */
835
836 z = _mm_and_ps(y,sieve);
837 q = _mm_mul_ps( _mm_sub_ps(z,y) , _mm_add_ps(z,y) );
838
839 corr = _mm_mul_ps(CD4,q);
840 corr = _mm_add_ps(corr,CD3);
841 corr = _mm_mul_ps(corr,q);
842 corr = _mm_add_ps(corr,CD2);
843 corr = _mm_mul_ps(corr,q);
844 corr = _mm_add_ps(corr,one);
845 corr = _mm_mul_ps(corr,q);
846 corr = _mm_add_ps(corr,one);
847
848 expmx2 = gmx_mm_exp_ps( _mm_or_ps( signbit , _mm_mul_ps(z,z) ) );
849 expmx2 = _mm_mul_ps(expmx2,corr);
850
851 pB1 = _mm_mul_ps(CB9,w2);
852 pB0 = _mm_mul_ps(CB8,w2);
853 pB1 = _mm_add_ps(pB1,CB7);
854 pB0 = _mm_add_ps(pB0,CB6);
855 pB1 = _mm_mul_ps(pB1,w2);
856 pB0 = _mm_mul_ps(pB0,w2);
857 pB1 = _mm_add_ps(pB1,CB5);
858 pB0 = _mm_add_ps(pB0,CB4);
859 pB1 = _mm_mul_ps(pB1,w2);
860 pB0 = _mm_mul_ps(pB0,w2);
861 pB1 = _mm_add_ps(pB1,CB3);
862 pB0 = _mm_add_ps(pB0,CB2);
863 pB1 = _mm_mul_ps(pB1,w2);
864 pB0 = _mm_mul_ps(pB0,w2);
865 pB1 = _mm_add_ps(pB1,CB1);
866 pB1 = _mm_mul_ps(pB1,w);
867 pB0 = _mm_add_ps(pB0,pB1);
868 pB0 = _mm_add_ps(pB0,CB0);
869
870 pC0 = _mm_mul_ps(CC10,t2);
871 pC1 = _mm_mul_ps(CC9,t2);
872 pC0 = _mm_add_ps(pC0,CC8);
873 pC1 = _mm_add_ps(pC1,CC7);
874 pC0 = _mm_mul_ps(pC0,t2);
875 pC1 = _mm_mul_ps(pC1,t2);
876 pC0 = _mm_add_ps(pC0,CC6);
877 pC1 = _mm_add_ps(pC1,CC5);
878 pC0 = _mm_mul_ps(pC0,t2);
879 pC1 = _mm_mul_ps(pC1,t2);
880 pC0 = _mm_add_ps(pC0,CC4);
881 pC1 = _mm_add_ps(pC1,CC3);
882 pC0 = _mm_mul_ps(pC0,t2);
883 pC1 = _mm_mul_ps(pC1,t2);
884 pC0 = _mm_add_ps(pC0,CC2);
885 pC1 = _mm_add_ps(pC1,CC1);
886 pC0 = _mm_mul_ps(pC0,t2);
887 pC1 = _mm_mul_ps(pC1,t);
888 pC0 = _mm_add_ps(pC0,pC1);
889 pC0 = _mm_add_ps(pC0,CC0);
890 pC0 = _mm_mul_ps(pC0,t);
891
892 /* SELECT pB0 or pC0 for erfc() */
893 mask = _mm_cmplt_ps(two,y);
894 res_erfc = _mm_blendv_ps(pB0,pC0,mask);
895 res_erfc = _mm_mul_ps(res_erfc,expmx2);
896
897 /* erfc(x<0) = 2-erfc(|x|) */
898 mask = _mm_cmplt_ps(x,_mm_setzero_ps());
899 res_erfc = _mm_blendv_ps(res_erfc,_mm_sub_ps(two,res_erfc),mask);
900
901 /* Select erf() or erfc() */
902 mask = _mm_cmplt_ps(y,_mm_set1_ps(0.75f));
903 res = _mm_blendv_ps(_mm_sub_ps(one,res_erfc),res_erf,mask);
904
905 return res;
906 }
907
908
909
910
911 /* FULL precision erfc(), 256 bit wide. Only errors in LSB */
912 static __m256
913 gmx_mm256_erfc_ps(__m256 x)
914 {
915 /* Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1] */
916 const __m256 CA6 = _mm256_set1_ps(7.853861353153693e-5f);
917 const __m256 CA5 = _mm256_set1_ps(-8.010193625184903e-4f);
918 const __m256 CA4 = _mm256_set1_ps(5.188327685732524e-3f);
919 const __m256 CA3 = _mm256_set1_ps(-2.685381193529856e-2f);
920 const __m256 CA2 = _mm256_set1_ps(1.128358514861418e-1f);
921 const __m256 CA1 = _mm256_set1_ps(-3.761262582423300e-1f);
922 const __m256 CA0 = _mm256_set1_ps(1.128379165726710f);
923 /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2] */
924 const __m256 CB9 = _mm256_set1_ps(-0.0018629930017603923f);
925 const __m256 CB8 = _mm256_set1_ps(0.003909821287598495f);
926 const __m256 CB7 = _mm256_set1_ps(-0.0052094582210355615f);
927 const __m256 CB6 = _mm256_set1_ps(0.005685614362160572f);
928 const __m256 CB5 = _mm256_set1_ps(-0.0025367682853477272f);
929 const __m256 CB4 = _mm256_set1_ps(-0.010199799682318782f);
930 const __m256 CB3 = _mm256_set1_ps(0.04369575504816542f);
931 const __m256 CB2 = _mm256_set1_ps(-0.11884063474674492f);
932 const __m256 CB1 = _mm256_set1_ps(0.2732120154030589f);
933 const __m256 CB0 = _mm256_set1_ps(0.42758357702025784f);
934 /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19] */
935 const __m256 CC10 = _mm256_set1_ps(-0.0445555913112064f);
936 const __m256 CC9 = _mm256_set1_ps(0.21376355144663348f);
937 const __m256 CC8 = _mm256_set1_ps(-0.3473187200259257f);
938 const __m256 CC7 = _mm256_set1_ps(0.016690861551248114f);
939 const __m256 CC6 = _mm256_set1_ps(0.7560973182491192f);
940 const __m256 CC5 = _mm256_set1_ps(-1.2137903600145787f);
941 const __m256 CC4 = _mm256_set1_ps(0.8411872321232948f);
942 const __m256 CC3 = _mm256_set1_ps(-0.08670413896296343f);
943 const __m256 CC2 = _mm256_set1_ps(-0.27124782687240334f);
944 const __m256 CC1 = _mm256_set1_ps(-0.0007502488047806069f);
945 const __m256 CC0 = _mm256_set1_ps(0.5642114853803148f);
946
947 /* Coefficients for expansion of exp(x) in [0,0.1] */
948 /* CD0 and CD1 are both 1.0, so no need to declare them separately */
949 const __m256 CD2 = _mm256_set1_ps(0.5000066608081202f);
950 const __m256 CD3 = _mm256_set1_ps(0.1664795422874624f);
951 const __m256 CD4 = _mm256_set1_ps(0.04379839977652482f);
952
953 const __m256 sieve = _mm256_castsi256_ps( _mm256_set1_epi32(0xfffff000) );
954 const __m256 signbit = _mm256_castsi256_ps( _mm256_set1_epi32(0x80000000) );
955 const __m256 one = _mm256_set1_ps(1.0f);
956 const __m256 two = _mm256_set1_ps(2.0f);
957
958 __m256 x2,x4,y;
959 __m256 z,q,t,t2,w,w2;
960 __m256 pA0,pA1,pB0,pB1,pC0,pC1;
961 __m256 expmx2,corr;
962 __m256 res_erf,res_erfc,res;
963 __m256 mask;
964
965 /* Calculate erf() */
966 x2 = _mm256_mul_ps(x,x);
967 x4 = _mm256_mul_ps(x2,x2);
968
969 pA0 = _mm256_mul_ps(CA6,x4);
970 pA1 = _mm256_mul_ps(CA5,x4);
971 pA0 = _mm256_add_ps(pA0,CA4);
972 pA1 = _mm256_add_ps(pA1,CA3);
973 pA0 = _mm256_mul_ps(pA0,x4);
974 pA1 = _mm256_mul_ps(pA1,x4);
975 pA0 = _mm256_add_ps(pA0,CA2);
976 pA1 = _mm256_add_ps(pA1,CA1);
977 pA0 = _mm256_mul_ps(pA0,x4);
978 pA1 = _mm256_mul_ps(pA1,x2);
979 pA0 = _mm256_add_ps(pA0,pA1);
980 pA0 = _mm256_add_ps(pA0,CA0);
981
982 res_erf = _mm256_mul_ps(x,pA0);
983
984 /* Calculate erfc */
985 y = gmx_mm256_abs_ps(x);
986 t = gmx_mm256_inv_ps(y);
987 w = _mm256_sub_ps(t,one);
988 t2 = _mm256_mul_ps(t,t);
989 w2 = _mm256_mul_ps(w,w);
990 /*
991 * We cannot simply calculate exp(-x2) directly in single precision, since
992 * that will lose a couple of bits of precision due to the multiplication.
993 * Instead, we introduce x=z+w, where the last 12 bits of precision are in w.
994 * Then we get exp(-x2) = exp(-z2)*exp((z-x)*(z+x)).
995 *
996 * The only drawback with this is that it requires TWO separate exponential
997 * evaluations, which would be horrible performance-wise. However, the argument
998 * for the second exp() call is always small, so there we simply use a
999 * low-order minimax expansion on [0,0.1].
1000 */
1001
1002 z = _mm256_and_ps(y,sieve);
1003 q = _mm256_mul_ps( _mm256_sub_ps(z,y) , _mm256_add_ps(z,y) );
1004
1005 corr = _mm256_mul_ps(CD4,q);
1006 corr = _mm256_add_ps(corr,CD3);
1007 corr = _mm256_mul_ps(corr,q);
1008 corr = _mm256_add_ps(corr,CD2);
1009 corr = _mm256_mul_ps(corr,q);
1010 corr = _mm256_add_ps(corr,one);
1011 corr = _mm256_mul_ps(corr,q);
1012 corr = _mm256_add_ps(corr,one);
1013
1014 expmx2 = gmx_mm256_exp_ps( _mm256_or_ps( signbit , _mm256_mul_ps(z,z) ) );
1015 expmx2 = _mm256_mul_ps(expmx2,corr);
1016
1017 pB1 = _mm256_mul_ps(CB9,w2);
1018 pB0 = _mm256_mul_ps(CB8,w2);
1019 pB1 = _mm256_add_ps(pB1,CB7);
1020 pB0 = _mm256_add_ps(pB0,CB6);
1021 pB1 = _mm256_mul_ps(pB1,w2);
1022 pB0 = _mm256_mul_ps(pB0,w2);
1023 pB1 = _mm256_add_ps(pB1,CB5);
1024 pB0 = _mm256_add_ps(pB0,CB4);
1025 pB1 = _mm256_mul_ps(pB1,w2);
1026 pB0 = _mm256_mul_ps(pB0,w2);
1027 pB1 = _mm256_add_ps(pB1,CB3);
1028 pB0 = _mm256_add_ps(pB0,CB2);
1029 pB1 = _mm256_mul_ps(pB1,w2);
1030 pB0 = _mm256_mul_ps(pB0,w2);
1031 pB1 = _mm256_add_ps(pB1,CB1);
1032 pB1 = _mm256_mul_ps(pB1,w);
1033 pB0 = _mm256_add_ps(pB0,pB1);
1034 pB0 = _mm256_add_ps(pB0,CB0);
1035
1036 pC0 = _mm256_mul_ps(CC10,t2);
1037 pC1 = _mm256_mul_ps(CC9,t2);
1038 pC0 = _mm256_add_ps(pC0,CC8);
1039 pC1 = _mm256_add_ps(pC1,CC7);
1040 pC0 = _mm256_mul_ps(pC0,t2);
1041 pC1 = _mm256_mul_ps(pC1,t2);
1042 pC0 = _mm256_add_ps(pC0,CC6);
1043 pC1 = _mm256_add_ps(pC1,CC5);
1044 pC0 = _mm256_mul_ps(pC0,t2);
1045 pC1 = _mm256_mul_ps(pC1,t2);
1046 pC0 = _mm256_add_ps(pC0,CC4);
1047 pC1 = _mm256_add_ps(pC1,CC3);
1048 pC0 = _mm256_mul_ps(pC0,t2);
1049 pC1 = _mm256_mul_ps(pC1,t2);
1050 pC0 = _mm256_add_ps(pC0,CC2);
1051 pC1 = _mm256_add_ps(pC1,CC1);
1052 pC0 = _mm256_mul_ps(pC0,t2);
1053 pC1 = _mm256_mul_ps(pC1,t);
1054 pC0 = _mm256_add_ps(pC0,pC1);
1055 pC0 = _mm256_add_ps(pC0,CC0);
1056 pC0 = _mm256_mul_ps(pC0,t);
1057
1058 /* SELECT pB0 or pC0 for erfc() */
1059 mask = _mm256_cmp_ps(two,y,_CMP_LT_OQ);
1060 res_erfc = _mm256_blendv_ps(pB0,pC0,mask);
1061 res_erfc = _mm256_mul_ps(res_erfc,expmx2);
1062
1063 /* erfc(x<0) = 2-erfc(|x|) */
1064 mask = _mm256_cmp_ps(x,_mm256_setzero_ps(),_CMP_LT_OQ);
1065 res_erfc = _mm256_blendv_ps(res_erfc,_mm256_sub_ps(two,res_erfc),mask);
1066
1067 /* Select erf() or erfc() */
1068 mask = _mm256_cmp_ps(y,_mm256_set1_ps(0.75f),_CMP_LT_OQ);
1069 res = _mm256_blendv_ps(res_erfc,_mm256_sub_ps(one,res_erf),mask);
1070
1071 return res;
1072 }
1073
1074
1075 /* erfc(), 128 bit wide */
1076 static __m128
1077 gmx_mm_erfc_ps(__m128 x)
1078 {
1079 /* Coefficients for minimax approximation of erf(x)=x*P(x^2) in range [-1,1] */
1080 const __m128 CA6 = _mm_set1_ps(7.853861353153693e-5f);
1081 const __m128 CA5 = _mm_set1_ps(-8.010193625184903e-4f);
1082 const __m128 CA4 = _mm_set1_ps(5.188327685732524e-3f);
1083 const __m128 CA3 = _mm_set1_ps(-2.685381193529856e-2f);
1084 const __m128 CA2 = _mm_set1_ps(1.128358514861418e-1f);
1085 const __m128 CA1 = _mm_set1_ps(-3.761262582423300e-1f);
1086 const __m128 CA0 = _mm_set1_ps(1.128379165726710f);
1087 /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*P((1/(x-1))^2) in range [0.67,2] */
1088 const __m128 CB9 = _mm_set1_ps(-0.0018629930017603923f);
1089 const __m128 CB8 = _mm_set1_ps(0.003909821287598495f);
1090 const __m128 CB7 = _mm_set1_ps(-0.0052094582210355615f);
1091 const __m128 CB6 = _mm_set1_ps(0.005685614362160572f);
1092 const __m128 CB5 = _mm_set1_ps(-0.0025367682853477272f);
1093 const __m128 CB4 = _mm_set1_ps(-0.010199799682318782f);
1094 const __m128 CB3 = _mm_set1_ps(0.04369575504816542f);
1095 const __m128 CB2 = _mm_set1_ps(-0.11884063474674492f);
1096 const __m128 CB1 = _mm_set1_ps(0.2732120154030589f);
1097 const __m128 CB0 = _mm_set1_ps(0.42758357702025784f);
1098 /* Coefficients for minimax approximation of erfc(x)=Exp(-x^2)*(1/x)*P((1/x)^2) in range [2,9.19] */
1099 const __m128 CC10 = _mm_set1_ps(-0.0445555913112064f);
1100 const __m128 CC9 = _mm_set1_ps(0.21376355144663348f);
1101 const __m128 CC8 = _mm_set1_ps(-0.3473187200259257f);
1102 const __m128 CC7 = _mm_set1_ps(0.016690861551248114f);
1103 const __m128 CC6 = _mm_set1_ps(0.7560973182491192f);
1104 const __m128 CC5 = _mm_set1_ps(-1.2137903600145787f);
1105 const __m128 CC4 = _mm_set1_ps(0.8411872321232948f);
1106 const __m128 CC3 = _mm_set1_ps(-0.08670413896296343f);
1107 const __m128 CC2 = _mm_set1_ps(-0.27124782687240334f);
1108 const __m128 CC1 = _mm_set1_ps(-0.0007502488047806069f);
1109 const __m128 CC0 = _mm_set1_ps(0.5642114853803148f);
1110
1111 /* Coefficients for expansion of exp(x) in [0,0.1] */
1112 /* CD0 and CD1 are both 1.0, so no need to declare them separately */
1113 const __m128 CD2 = _mm_set1_ps(0.5000066608081202f);
1114 const __m128 CD3 = _mm_set1_ps(0.1664795422874624f);
1115 const __m128 CD4 = _mm_set1_ps(0.04379839977652482f);
1116
1117 const __m128 sieve = gmx_mm_castsi128_ps( _mm_set1_epi32(0xfffff000) );
1118 const __m128 signbit = gmx_mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
1119 const __m128 one = _mm_set1_ps(1.0f);
1120 const __m128 two = _mm_set1_ps(2.0f);
1121
1122 __m128 x2,x4,y;
1123 __m128 z,q,t,t2,w,w2;
1124 __m128 pA0,pA1,pB0,pB1,pC0,pC1;
1125 __m128 expmx2,corr;
1126 __m128 res_erf,res_erfc,res;
1127 __m128 mask;
1128
1129 /* Calculate erf() */
1130 x2 = _mm_mul_ps(x,x);
1131 x4 = _mm_mul_ps(x2,x2);
1132
1133 pA0 = _mm_mul_ps(CA6,x4);
1134 pA1 = _mm_mul_ps(CA5,x4);
1135 pA0 = _mm_add_ps(pA0,CA4);
1136 pA1 = _mm_add_ps(pA1,CA3);
1137 pA0 = _mm_mul_ps(pA0,x4);
1138 pA1 = _mm_mul_ps(pA1,x4);
1139 pA0 = _mm_add_ps(pA0,CA2);
1140 pA1 = _mm_add_ps(pA1,CA1);
1141 pA0 = _mm_mul_ps(pA0,x4);
1142 pA1 = _mm_mul_ps(pA1,x2);
1143 pA0 = _mm_add_ps(pA0,pA1);
1144 pA0 = _mm_add_ps(pA0,CA0);
1145
1146 res_erf = _mm_mul_ps(x,pA0);
1147
1148 /* Calculate erfc */
1149 y = gmx_mm_abs_ps(x);
1150 t = gmx_mm_inv_ps(y);
1151 w = _mm_sub_ps(t,one);
1152 t2 = _mm_mul_ps(t,t);
1153 w2 = _mm_mul_ps(w,w);
1154 /*
1155 * We cannot simply calculate exp(-x2) directly in single precision, since
1156 * that will lose a couple of bits of precision due to the multiplication.
1157 * Instead, we introduce x=z+w, where the last 12 bits of precision are in w.
1158 * Then we get exp(-x2) = exp(-z2)*exp((z-x)*(z+x)).
1159 *
1160 * The only drawback with this is that it requires TWO separate exponential
1161 * evaluations, which would be horrible performance-wise. However, the argument
1162 * for the second exp() call is always small, so there we simply use a
1163 * low-order minimax expansion on [0,0.1].
1164 */
1165
1166 z = _mm_and_ps(y,sieve);
1167 q = _mm_mul_ps( _mm_sub_ps(z,y) , _mm_add_ps(z,y) );
1168
1169 corr = _mm_mul_ps(CD4,q);
1170 corr = _mm_add_ps(corr,CD3);
1171 corr = _mm_mul_ps(corr,q);
1172 corr = _mm_add_ps(corr,CD2);
1173 corr = _mm_mul_ps(corr,q);
1174 corr = _mm_add_ps(corr,one);
1175 corr = _mm_mul_ps(corr,q);
1176 corr = _mm_add_ps(corr,one);
1177
1178 expmx2 = gmx_mm_exp_ps( _mm_or_ps( signbit , _mm_mul_ps(z,z) ) );
1179 expmx2 = _mm_mul_ps(expmx2,corr);
1180
1181 pB1 = _mm_mul_ps(CB9,w2);
1182 pB0 = _mm_mul_ps(CB8,w2);
1183 pB1 = _mm_add_ps(pB1,CB7);
1184 pB0 = _mm_add_ps(pB0,CB6);
1185 pB1 = _mm_mul_ps(pB1,w2);
1186 pB0 = _mm_mul_ps(pB0,w2);
1187 pB1 = _mm_add_ps(pB1,CB5);
1188 pB0 = _mm_add_ps(pB0,CB4);
1189 pB1 = _mm_mul_ps(pB1,w2);
1190 pB0 = _mm_mul_ps(pB0,w2);
1191 pB1 = _mm_add_ps(pB1,CB3);
1192 pB0 = _mm_add_ps(pB0,CB2);
1193 pB1 = _mm_mul_ps(pB1,w2);
1194 pB0 = _mm_mul_ps(pB0,w2);
1195 pB1 = _mm_add_ps(pB1,CB1);
1196 pB1 = _mm_mul_ps(pB1,w);
1197 pB0 = _mm_add_ps(pB0,pB1);
1198 pB0 = _mm_add_ps(pB0,CB0);
1199
1200 pC0 = _mm_mul_ps(CC10,t2);
1201 pC1 = _mm_mul_ps(CC9,t2);
1202 pC0 = _mm_add_ps(pC0,CC8);
1203 pC1 = _mm_add_ps(pC1,CC7);
1204 pC0 = _mm_mul_ps(pC0,t2);
1205 pC1 = _mm_mul_ps(pC1,t2);
1206 pC0 = _mm_add_ps(pC0,CC6);
1207 pC1 = _mm_add_ps(pC1,CC5);
1208 pC0 = _mm_mul_ps(pC0,t2);
1209 pC1 = _mm_mul_ps(pC1,t2);
1210 pC0 = _mm_add_ps(pC0,CC4);
1211 pC1 = _mm_add_ps(pC1,CC3);
1212 pC0 = _mm_mul_ps(pC0,t2);
1213 pC1 = _mm_mul_ps(pC1,t2);
1214 pC0 = _mm_add_ps(pC0,CC2);
1215 pC1 = _mm_add_ps(pC1,CC1);
1216 pC0 = _mm_mul_ps(pC0,t2);
1217 pC1 = _mm_mul_ps(pC1,t);
1218 pC0 = _mm_add_ps(pC0,pC1);
1219 pC0 = _mm_add_ps(pC0,CC0);
1220 pC0 = _mm_mul_ps(pC0,t);
1221
1222 /* SELECT pB0 or pC0 for erfc() */
1223 mask = _mm_cmplt_ps(two,y);
1224 res_erfc = _mm_blendv_ps(pB0,pC0,mask);
1225 res_erfc = _mm_mul_ps(res_erfc,expmx2);
1226
1227 /* erfc(x<0) = 2-erfc(|x|) */
1228 mask = _mm_cmplt_ps(x,_mm_setzero_ps());
1229 res_erfc = _mm_blendv_ps(res_erfc,_mm_sub_ps(two,res_erfc),mask);
1230
1231 /* Select erf() or erfc() */
1232 mask = _mm_cmplt_ps(y,_mm_set1_ps(0.75f));
1233 res = _mm_blendv_ps(res_erfc,_mm_sub_ps(one,res_erf),mask);
1234
1235 return res;
1236 }
1237
1238
1239
1240 /* Calculate the force correction due to PME analytically.
1241 *
1242 * This routine is meant to enable analytical evaluation of the
1243 * direct-space PME electrostatic force to avoid tables.
1244 *
1245 * The direct-space potential should be Erfc(beta*r)/r, but there
1246 * are some problems evaluating that:
1247 *
1248 * First, the error function is difficult (read: expensive) to
1249 * approxmiate accurately for intermediate to large arguments, and
1250 * this happens already in ranges of beta*r that occur in simulations.
1251 * Second, we now try to avoid calculating potentials in Gromacs but
1252 * use forces directly.
1253 *
1254 * We can simply things slight by noting that the PME part is really
1255 * a correction to the normal Coulomb force since Erfc(z)=1-Erf(z), i.e.
1256 *
1257 * V= 1/r - Erf(beta*r)/r
1258 *
1259 * The first term we already have from the inverse square root, so
1260 * that we can leave out of this routine.
1261 *
1262 * For pme tolerances of 1e-3 to 1e-8 and cutoffs of 0.5nm to 1.8nm,
1263 * the argument beta*r will be in the range 0.15 to ~4. Use your
1264 * favorite plotting program to realize how well-behaved Erf(z)/z is
1265 * in this range!
1266 *
1267 * We approximate f(z)=erf(z)/z with a rational minimax polynomial.
1268 * However, it turns out it is more efficient to approximate f(z)/z and
1269 * then only use even powers. This is another minor optimization, since
1270 * we actually WANT f(z)/z, because it is going to be multiplied by
1271 * the vector between the two atoms to get the vectorial force. The
1272 * fastest flops are the ones we can avoid calculating!
1273 *
1274 * So, here's how it should be used:
1275 *
1276 * 1. Calculate r^2.
1277 * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
1278 * 3. Evaluate this routine with z^2 as the argument.
1279 * 4. The return value is the expression:
1280 *
1281 *
1282 * 2*exp(-z^2) erf(z)
1283 * ------------ - --------
1284 * sqrt(Pi)*z^2 z^3
1285 *
1286 * 5. Multiply the entire expression by beta^3. This will get you
1287 *
1288 * beta^3*2*exp(-z^2) beta^3*erf(z)
1289 * ------------------ - ---------------
1290 * sqrt(Pi)*z^2 z^3
1291 *
1292 * or, switching back to r (z=r*beta):
1293 *
1294 * 2*beta*exp(-r^2*beta^2) erf(r*beta)
1295 * ----------------------- - -----------
1296 * sqrt(Pi)*r^2 r^3
1297 *
1298 *
1299 * With a bit of math exercise you should be able to confirm that
1300 * this is exactly D[Erf[beta*r]/r,r] divided by r another time.
1301 *
1302 * 6. Add the result to 1/r^3, multiply by the product of the charges,
1303 * and you have your force (divided by r). A final multiplication
1304 * with the vector connecting the two particles and you have your
1305 * vectorial force to add to the particles.
1306 *
1307 */
1308 static __m256
1309 gmx_mm256_pmecorrF_ps(__m256 z2)
1310 {
1311 const __m256 FN6 = _mm256_set1_ps(-1.7357322914161492954e-8f);
1312 const __m256 FN5 = _mm256_set1_ps(1.4703624142580877519e-6f);
1313 const __m256 FN4 = _mm256_set1_ps(-0.000053401640219807709149f);
1314 const __m256 FN3 = _mm256_set1_ps(0.0010054721316683106153f);
1315 const __m256 FN2 = _mm256_set1_ps(-0.019278317264888380590f);
1316 const __m256 FN1 = _mm256_set1_ps(0.069670166153766424023f);
1317 const __m256 FN0 = _mm256_set1_ps(-0.75225204789749321333f);
1318
1319 const __m256 FD4 = _mm256_set1_ps(0.0011193462567257629232f);
1320 const __m256 FD3 = _mm256_set1_ps(0.014866955030185295499f);
1321 const __m256 FD2 = _mm256_set1_ps(0.11583842382862377919f);
1322 const __m256 FD1 = _mm256_set1_ps(0.50736591960530292870f);
1323 const __m256 FD0 = _mm256_set1_ps(1.0f);
1324
1325 __m256 z4;
1326 __m256 polyFN0,polyFN1,polyFD0,polyFD1;
1327
1328 z4 = _mm256_mul_ps(z2,z2);
1329
1330 polyFD0 = _mm256_mul_ps(FD4,z4);
1331 polyFD1 = _mm256_mul_ps(FD3,z4);
1332 polyFD0 = _mm256_add_ps(polyFD0,FD2);
1333 polyFD1 = _mm256_add_ps(polyFD1,FD1);
1334 polyFD0 = _mm256_mul_ps(polyFD0,z4);
1335 polyFD1 = _mm256_mul_ps(polyFD1,z2);
1336 polyFD0 = _mm256_add_ps(polyFD0,FD0);
1337 polyFD0 = _mm256_add_ps(polyFD0,polyFD1);
1338
1339 polyFD0 = gmx_mm256_inv_ps(polyFD0);
1340
1341 polyFN0 = _mm256_mul_ps(FN6,z4);
1342 polyFN1 = _mm256_mul_ps(FN5,z4);
1343 polyFN0 = _mm256_add_ps(polyFN0,FN4);
1344 polyFN1 = _mm256_add_ps(polyFN1,FN3);
1345 polyFN0 = _mm256_mul_ps(polyFN0,z4);
1346 polyFN1 = _mm256_mul_ps(polyFN1,z4);
1347 polyFN0 = _mm256_add_ps(polyFN0,FN2);
1348 polyFN1 = _mm256_add_ps(polyFN1,FN1);
1349 polyFN0 = _mm256_mul_ps(polyFN0,z4);
1350 polyFN1 = _mm256_mul_ps(polyFN1,z2);
1351 polyFN0 = _mm256_add_ps(polyFN0,FN0);
1352 polyFN0 = _mm256_add_ps(polyFN0,polyFN1);
1353
1354 return _mm256_mul_ps(polyFN0,polyFD0);
1355 }
1356
1357
1358 static __m128
1359 gmx_mm_pmecorrF_ps(__m128 z2)
1360 {
1361 const __m128 FN6 = _mm_set1_ps(-1.7357322914161492954e-8f);
1362 const __m128 FN5 = _mm_set1_ps(1.4703624142580877519e-6f);
1363 const __m128 FN4 = _mm_set1_ps(-0.000053401640219807709149f);
1364 const __m128 FN3 = _mm_set1_ps(0.0010054721316683106153f);
1365 const __m128 FN2 = _mm_set1_ps(-0.019278317264888380590f);
1366 const __m128 FN1 = _mm_set1_ps(0.069670166153766424023f);
1367 const __m128 FN0 = _mm_set1_ps(-0.75225204789749321333f);
1368
1369 const __m128 FD4 = _mm_set1_ps(0.0011193462567257629232f);
1370 const __m128 FD3 = _mm_set1_ps(0.014866955030185295499f);
1371 const __m128 FD2 = _mm_set1_ps(0.11583842382862377919f);
1372 const __m128 FD1 = _mm_set1_ps(0.50736591960530292870f);
1373 const __m128 FD0 = _mm_set1_ps(1.0f);
1374
1375 __m128 z4;
1376 __m128 polyFN0,polyFN1,polyFD0,polyFD1;
1377
1378 z4 = _mm_mul_ps(z2,z2);
1379
1380 polyFD0 = _mm_mul_ps(FD4,z4);
1381 polyFD1 = _mm_mul_ps(FD3,z4);
1382 polyFD0 = _mm_add_ps(polyFD0,FD2);
1383 polyFD1 = _mm_add_ps(polyFD1,FD1);
1384 polyFD0 = _mm_mul_ps(polyFD0,z4);
1385 polyFD1 = _mm_mul_ps(polyFD1,z2);
1386 polyFD0 = _mm_add_ps(polyFD0,FD0);
1387 polyFD0 = _mm_add_ps(polyFD0,polyFD1);
1388
1389 polyFD0 = gmx_mm_inv_ps(polyFD0);
1390
1391 polyFN0 = _mm_mul_ps(FN6,z4);
1392 polyFN1 = _mm_mul_ps(FN5,z4);
1393 polyFN0 = _mm_add_ps(polyFN0,FN4);
1394 polyFN1 = _mm_add_ps(polyFN1,FN3);
1395 polyFN0 = _mm_mul_ps(polyFN0,z4);
1396 polyFN1 = _mm_mul_ps(polyFN1,z4);
1397 polyFN0 = _mm_add_ps(polyFN0,FN2);
1398 polyFN1 = _mm_add_ps(polyFN1,FN1);
1399 polyFN0 = _mm_mul_ps(polyFN0,z4);
1400 polyFN1 = _mm_mul_ps(polyFN1,z2);
1401 polyFN0 = _mm_add_ps(polyFN0,FN0);
1402 polyFN0 = _mm_add_ps(polyFN0,polyFN1);
1403
1404 return _mm_mul_ps(polyFN0,polyFD0);
1405 }
1406
1407
1408
1409 /* Calculate the potential correction due to PME analytically.
1410 *
1411 * See gmx_mm256_pmecorrF_ps() for details about the approximation.
1412 *
1413 * This routine calculates Erf(z)/z, although you should provide z^2
1414 * as the input argument.
1415 *
1416 * Here's how it should be used:
1417 *
1418 * 1. Calculate r^2.
1419 * 2. Multiply by beta^2, so you get z^2=beta^2*r^2.
1420 * 3. Evaluate this routine with z^2 as the argument.
1421 * 4. The return value is the expression:
1422 *
1423 *
1424 * erf(z)
1425 * --------
1426 * z
1427 *
1428 * 5. Multiply the entire expression by beta and switching back to r (z=r*beta):
1429 *
1430 * erf(r*beta)
1431 * -----------
1432 * r
1433 *
1434 * 6. Subtract the result from 1/r, multiply by the product of the charges,
1435 * and you have your potential.
1436 */
1437 static __m256
1438 gmx_mm256_pmecorrV_ps(__m256 z2)
1439 {
1440 const __m256 VN6 = _mm256_set1_ps(1.9296833005951166339e-8f);
1441 const __m256 VN5 = _mm256_set1_ps(-1.4213390571557850962e-6f);
1442 const __m256 VN4 = _mm256_set1_ps(0.000041603292906656984871f);
1443 const __m256 VN3 = _mm256_set1_ps(-0.00013134036773265025626f);
1444 const __m256 VN2 = _mm256_set1_ps(0.038657983986041781264f);
1445 const __m256 VN1 = _mm256_set1_ps(0.11285044772717598220f);
1446 const __m256 VN0 = _mm256_set1_ps(1.1283802385263030286f);
1447
1448 const __m256 VD3 = _mm256_set1_ps(0.0066752224023576045451f);
1449 const __m256 VD2 = _mm256_set1_ps(0.078647795836373922256f);
1450 const __m256 VD1 = _mm256_set1_ps(0.43336185284710920150f);
1451 const __m256 VD0 = _mm256_set1_ps(1.0f);
1452
1453 __m256 z4;
1454 __m256 polyVN0,polyVN1,polyVD0,polyVD1;
1455
1456 z4 = _mm256_mul_ps(z2,z2);
1457
1458 polyVD1 = _mm256_mul_ps(VD3,z4);
1459 polyVD0 = _mm256_mul_ps(VD2,z4);
1460 polyVD1 = _mm256_add_ps(polyVD1,VD1);
1461 polyVD0 = _mm256_add_ps(polyVD0,VD0);
1462 polyVD1 = _mm256_mul_ps(polyVD1,z2);
1463 polyVD0 = _mm256_add_ps(polyVD0,polyVD1);
1464
1465 polyVD0 = gmx_mm256_inv_ps(polyVD0);
1466
1467 polyVN0 = _mm256_mul_ps(VN6,z4);
1468 polyVN1 = _mm256_mul_ps(VN5,z4);
1469 polyVN0 = _mm256_add_ps(polyVN0,VN4);
1470 polyVN1 = _mm256_add_ps(polyVN1,VN3);
1471 polyVN0 = _mm256_mul_ps(polyVN0,z4);
1472 polyVN1 = _mm256_mul_ps(polyVN1,z4);
1473 polyVN0 = _mm256_add_ps(polyVN0,VN2);
1474 polyVN1 = _mm256_add_ps(polyVN1,VN1);
1475 polyVN0 = _mm256_mul_ps(polyVN0,z4);
1476 polyVN1 = _mm256_mul_ps(polyVN1,z2);
1477 polyVN0 = _mm256_add_ps(polyVN0,VN0);
1478 polyVN0 = _mm256_add_ps(polyVN0,polyVN1);
1479
1480 return _mm256_mul_ps(polyVN0,polyVD0);
1481 }
1482
1483
1484 static __m128
1485 gmx_mm_pmecorrV_ps(__m128 z2)
1486 {
1487 const __m128 VN6 = _mm_set1_ps(1.9296833005951166339e-8f);
1488 const __m128 VN5 = _mm_set1_ps(-1.4213390571557850962e-6f);
1489 const __m128 VN4 = _mm_set1_ps(0.000041603292906656984871f);
1490 const __m128 VN3 = _mm_set1_ps(-0.00013134036773265025626f);
1491 const __m128 VN2 = _mm_set1_ps(0.038657983986041781264f);
1492 const __m128 VN1 = _mm_set1_ps(0.11285044772717598220f);
1493 const __m128 VN0 = _mm_set1_ps(1.1283802385263030286f);
1494
1495 const __m128 VD3 = _mm_set1_ps(0.0066752224023576045451f);
1496 const __m128 VD2 = _mm_set1_ps(0.078647795836373922256f);
1497 const __m128 VD1 = _mm_set1_ps(0.43336185284710920150f);
1498 const __m128 VD0 = _mm_set1_ps(1.0f);
1499
1500 __m128 z4;
1501 __m128 polyVN0,polyVN1,polyVD0,polyVD1;
1502
1503 z4 = _mm_mul_ps(z2,z2);
1504
1505 polyVD1 = _mm_mul_ps(VD3,z4);
1506 polyVD0 = _mm_mul_ps(VD2,z4);
1507 polyVD1 = _mm_add_ps(polyVD1,VD1);
1508 polyVD0 = _mm_add_ps(polyVD0,VD0);
1509 polyVD1 = _mm_mul_ps(polyVD1,z2);
1510 polyVD0 = _mm_add_ps(polyVD0,polyVD1);
1511
1512 polyVD0 = gmx_mm_inv_ps(polyVD0);
1513
1514 polyVN0 = _mm_mul_ps(VN6,z4);
1515 polyVN1 = _mm_mul_ps(VN5,z4);
1516 polyVN0 = _mm_add_ps(polyVN0,VN4);
1517 polyVN1 = _mm_add_ps(polyVN1,VN3);
1518 polyVN0 = _mm_mul_ps(polyVN0,z4);
1519 polyVN1 = _mm_mul_ps(polyVN1,z4);
1520 polyVN0 = _mm_add_ps(polyVN0,VN2);
1521 polyVN1 = _mm_add_ps(polyVN1,VN1);
1522 polyVN0 = _mm_mul_ps(polyVN0,z4);
1523 polyVN1 = _mm_mul_ps(polyVN1,z2);
1524 polyVN0 = _mm_add_ps(polyVN0,VN0);
1525 polyVN0 = _mm_add_ps(polyVN0,polyVN1);
1526
1527 return _mm_mul_ps(polyVN0,polyVD0);
1528 }
1529
1530
1531 static int
1532 gmx_mm256_sincos_ps(__m256 x,
1533 __m256 *sinval,
1534 __m256 *cosval)
1535 {
1536 const __m256 two_over_pi = _mm256_set1_ps(2.0f/(float)M_PI);
1537 const __m256 half = _mm256_set1_ps(0.5f);
1538 const __m256 one = _mm256_set1_ps(1.0f);
1539 const __m256 zero = _mm256_setzero_ps();
1540
1541 const __m128i ione = _mm_set1_epi32(1);
1542
1543 const __m256 mask_one = _mm256_castsi256_ps(_mm256_set1_epi32(1));
1544 const __m256 mask_two = _mm256_castsi256_ps(_mm256_set1_epi32(2));
1545 const __m256 mask_three = _mm256_castsi256_ps(_mm256_set1_epi32(3));
1546
1547 const __m256 CA1 = _mm256_set1_ps(1.5703125f);
1548 const __m256 CA2 = _mm256_set1_ps(4.837512969970703125e-4f);
1549 const __m256 CA3 = _mm256_set1_ps(7.54978995489188216e-8f);
1550
1551 const __m256 CC0 = _mm256_set1_ps(-0.0013602249f);
1552 const __m256 CC1 = _mm256_set1_ps(0.0416566950f);
1553 const __m256 CC2 = _mm256_set1_ps(-0.4999990225f);
1554 const __m256 CS0 = _mm256_set1_ps(-0.0001950727f);
1555 const __m256 CS1 = _mm256_set1_ps(0.0083320758f);
1556 const __m256 CS2 = _mm256_set1_ps(-0.1666665247f);
1557
1558 const __m256 signbit = _mm256_castsi256_ps( _mm256_set1_epi32(0x80000000) );
1559
1560 __m256 y,y2;
1561 __m256 z;
1562 __m256i iz;
1563 __m128i iz_high,iz_low;
1564 __m256 offset_sin,offset_cos;
1565 __m256 mask_sin,mask_cos;
1566 __m256 tmp1,tmp2;
1567 __m256 tmp_sin,tmp_cos;
1568
1569 y = _mm256_mul_ps(x,two_over_pi);
1570 y = _mm256_add_ps(y,_mm256_or_ps(_mm256_and_ps(y,signbit),half));
1571
1572 iz = _mm256_cvttps_epi32(y);
1573 z = _mm256_round_ps(y,_MM_FROUND_TO_ZERO);
1574
1575 offset_sin = _mm256_and_ps(_mm256_castsi256_ps(iz),mask_three);
1576
1577 iz_high = _mm256_extractf128_si256(iz,0x1);
1578 iz_low = _mm256_castsi256_si128(iz);
1579 iz_low = _mm_add_epi32(iz_low,ione);
1580 iz_high = _mm_add_epi32(iz_high,ione);
1581 iz = _mm256_castsi128_si256(iz_low);
1582 iz = _mm256_insertf128_si256(iz,iz_high,0x1);
1583 offset_cos = _mm256_castsi256_ps(iz);
1584
1585 /* Extended precision arithmethic to achieve full precision */
1586 y = _mm256_mul_ps(z,CA1);
1587 tmp1 = _mm256_mul_ps(z,CA2);
1588 tmp2 = _mm256_mul_ps(z,CA3);
1589 y = _mm256_sub_ps(x,y);
1590 y = _mm256_sub_ps(y,tmp1);
1591 y = _mm256_sub_ps(y,tmp2);
1592
1593 y2 = _mm256_mul_ps(y,y);
1594
1595 tmp1 = _mm256_mul_ps(CC0,y2);
1596 tmp1 = _mm256_add_ps(tmp1,CC1);
1597 tmp2 = _mm256_mul_ps(CS0,y2);
1598 tmp2 = _mm256_add_ps(tmp2,CS1);
1599 tmp1 = _mm256_mul_ps(tmp1,y2);
1600 tmp1 = _mm256_add_ps(tmp1,CC2);
1601 tmp2 = _mm256_mul_ps(tmp2,y2);
1602 tmp2 = _mm256_add_ps(tmp2,CS2);
1603
1604 tmp1 = _mm256_mul_ps(tmp1,y2);
1605 tmp1 = _mm256_add_ps(tmp1,one);
1606
1607 tmp2 = _mm256_mul_ps(tmp2,_mm256_mul_ps(y,y2));
1608 tmp2 = _mm256_add_ps(tmp2,y);
1609
1610 #ifdef __INTEL_COMPILER
1611 /* Intel Compiler version 12.1.3 20120130 is buggy if optimization is enabled unless we cast explicitly! */
1612 mask_sin = _mm256_cmp_ps(_mm256_cvtepi32_ps(_mm256_castps_si256(_mm256_and_ps(offset_sin,mask_one))),zero,_CMP_EQ_OQ);
1613 mask_cos = _mm256_cmp_ps(_mm256_cvtepi32_ps(_mm256_castps_si256(_mm256_and_ps(offset_cos,mask_one))),zero,_CMP_EQ_OQ);
1614 #else
1615 mask_sin = _mm256_cmp_ps( _mm256_and_ps(offset_sin,mask_one), zero, _CMP_EQ_OQ);
1616 mask_cos = _mm256_cmp_ps( _mm256_and_ps(offset_cos,mask_one), zero, _CMP_EQ_OQ);
1617 #endif
1618 tmp_sin = _mm256_blendv_ps(tmp1,tmp2,mask_sin);
1619 tmp_cos = _mm256_blendv_ps(tmp1,tmp2,mask_cos);
1620
1621 tmp1 = _mm256_xor_ps(signbit,tmp_sin);
1622 tmp2 = _mm256_xor_ps(signbit,tmp_cos);
1623
1624 #ifdef __INTEL_COMPILER
1625 /* Intel Compiler version 12.1.3 20120130 is buggy if optimization is enabled unless we cast explicitly! */
1626 mask_sin = _mm256_cmp_ps(_mm256_cvtepi32_ps(_mm256_castps_si256(_mm256_and_ps(offset_sin,mask_two))),zero,_CMP_EQ_OQ);
1627 mask_cos = _mm256_cmp_ps(_mm256_cvtepi32_ps(_mm256_castps_si256(_mm256_and_ps(offset_cos,mask_two))),zero,_CMP_EQ_OQ);
1628 #else
1629 mask_sin = _mm256_cmp_ps( _mm256_and_ps(offset_sin,mask_two), zero, _CMP_EQ_OQ);
1630 mask_cos = _mm256_cmp_ps( _mm256_and_ps(offset_cos,mask_two), zero, _CMP_EQ_OQ);
1631
1632 #endif
1633 *sinval = _mm256_blendv_ps(tmp1,tmp_sin,mask_sin);
1634 *cosval = _mm256_blendv_ps(tmp2,tmp_cos,mask_cos);
1635
1636 return 0;
1637 }
1638
1639 static int
1640 gmx_mm_sincos_ps(__m128 x,
1641 __m128 *sinval,
1642 __m128 *cosval)
1643 {
1644 const __m128 two_over_pi = _mm_set1_ps(2.0/M_PI);
1645 const __m128 half = _mm_set1_ps(0.5);
1646 const __m128 one = _mm_set1_ps(1.0);
1647
1648 const __m128i izero = _mm_set1_epi32(0);
1649 const __m128i ione = _mm_set1_epi32(1);
1650 const __m128i itwo = _mm_set1_epi32(2);
1651 const __m128i ithree = _mm_set1_epi32(3);
1652 const __m128 signbit = gmx_mm_castsi128_ps( _mm_set1_epi32(0x80000000) );
1653
1654 const __m128 CA1 = _mm_set1_ps(1.5703125f);
1655 const __m128 CA2 = _mm_set1_ps(4.837512969970703125e-4f);
1656 const __m128 CA3 = _mm_set1_ps(7.54978995489188216e-8f);
1657
1658 const __m128 CC0 = _mm_set1_ps(-0.0013602249f);
1659 const __m128 CC1 = _mm_set1_ps(0.0416566950f);
1660 const __m128 CC2 = _mm_set1_ps(-0.4999990225f);
1661 const __m128 CS0 = _mm_set1_ps(-0.0001950727f);
1662 const __m128 CS1 = _mm_set1_ps(0.0083320758f);
1663 const __m128 CS2 = _mm_set1_ps(-0.1666665247f);
1664
1665 __m128 y,y2;
1666 __m128 z;
1667 __m128i iz;
1668 __m128i offset_sin,offset_cos;
1669 __m128 tmp1,tmp2;
1670 __m128 mask_sin,mask_cos;
1671 __m128 tmp_sin,tmp_cos;
1672
1673 y = _mm_mul_ps(x,two_over_pi);
1674 y = _mm_add_ps(y,_mm_or_ps(_mm_and_ps(y,signbit),half));
1675
1676 iz = _mm_cvttps_epi32(y);
1677 z = _mm_round_ps(y,_MM_FROUND_TO_ZERO);
1678
1679 offset_sin = _mm_and_si128(iz,ithree);
1680 offset_cos = _mm_add_epi32(iz,ione);
1681
1682 /* Extended precision arithmethic to achieve full precision */
1683 y = _mm_mul_ps(z,CA1);
1684 tmp1 = _mm_mul_ps(z,CA2);
1685 tmp2 = _mm_mul_ps(z,CA3);
1686 y = _mm_sub_ps(x,y);
1687 y = _mm_sub_ps(y,tmp1);
1688 y = _mm_sub_ps(y,tmp2);
1689
1690 y2 = _mm_mul_ps(y,y);
1691
1692 tmp1 = _mm_mul_ps(CC0,y2);
1693 tmp1 = _mm_add_ps(tmp1,CC1);
1694 tmp2 = _mm_mul_ps(CS0,y2);
1695 tmp2 = _mm_add_ps(tmp2,CS1);
1696 tmp1 = _mm_mul_ps(tmp1,y2);
1697 tmp1 = _mm_add_ps(tmp1,CC2);
1698 tmp2 = _mm_mul_ps(tmp2,y2);
1699 tmp2 = _mm_add_ps(tmp2,CS2);
1700
1701 tmp1 = _mm_mul_ps(tmp1,y2);
1702 tmp1 = _mm_add_ps(tmp1,one);
1703
1704 tmp2 = _mm_mul_ps(tmp2,_mm_mul_ps(y,y2));
1705 tmp2 = _mm_add_ps(tmp2,y);
1706
1707 mask_sin = gmx_mm_castsi128_ps(_mm_cmpeq_epi32( _mm_and_si128(offset_sin,ione), izero));
1708 mask_cos = gmx_mm_castsi128_ps(_mm_cmpeq_epi32( _mm_and_si128(offset_cos,ione), izero));
1709
1710 tmp_sin = _mm_blendv_ps(tmp1,tmp2,mask_sin);
1711 tmp_cos = _mm_blendv_ps(tmp1,tmp2,mask_cos);
1712
1713 mask_sin = gmx_mm_castsi128_ps(_mm_cmpeq_epi32( _mm_and_si128(offset_sin,itwo), izero));
1714 mask_cos = gmx_mm_castsi128_ps(_mm_cmpeq_epi32( _mm_and_si128(offset_cos,itwo), izero));
1715
1716 tmp1 = _mm_xor_ps(signbit,tmp_sin);
1717 tmp2 = _mm_xor_ps(signbit,tmp_cos);
1718
1719 *sinval = _mm_blendv_ps(tmp1,tmp_sin,mask_sin);
1720 *cosval = _mm_blendv_ps(tmp2,tmp_cos,mask_cos);
1721
1722 return 0;
1723 }
1724
1725
1726
1727
1728 /*
1729 * IMPORTANT: Do NOT call both sin & cos if you need both results, since each of them
1730 * will then call the sincos() routine and waste a factor 2 in performance!
1731 */
1732 static __m256
1733 gmx_mm256_sin_ps(__m256 x)
1734 {
1735 __m256 s,c;
1736 gmx_mm256_sincos_ps(x,&s,&c);
1737 return s;
1738 }
1739
1740 static __m128
1741 gmx_mm_sin_ps(__m128 x)
1742 {
1743 __m128 s,c;
1744 gmx_mm_sincos_ps(x,&s,&c);
1745 return s;
1746 }
1747
1748
1749 /*
1750 * IMPORTANT: Do NOT call both sin & cos if you need both results, since each of them
1751 * will then call the sincos() routine and waste a factor 2 in performance!
1752 */
1753 static __m256
1754 gmx_mm256_cos_ps(__m256 x)
1755 {
1756 __m256 s,c;
1757 gmx_mm256_sincos_ps(x,&s,&c);
1758 return c;
1759 }
1760
1761 static __m128
1762 gmx_mm_cos_ps(__m128 x)
1763 {
1764 __m128 s,c;
1765 gmx_mm_sincos_ps(x,&s,&c);
1766 return c;
1767 }
1768
1769
1770 static __m256
1771 gmx_mm256_tan_ps(__m256 x)
1772 {
1773 __m256 sinval,cosval;
1774 __m256 tanval;
1775
1776 gmx_mm256_sincos_ps(x,&sinval,&cosval);
1777
1778 tanval = _mm256_mul_ps(sinval,gmx_mm256_inv_ps(cosval));
1779
1780 return tanval;
1781 }
1782
1783 static __m128
1784 gmx_mm_tan_ps(__m128 x)
1785 {
1786 __m128 sinval,cosval;
1787 __m128 tanval;
1788
1789 gmx_mm_sincos_ps(x,&sinval,&cosval);
1790
1791 tanval = _mm_mul_ps(sinval,gmx_mm_inv_ps(cosval));
1792
1793 return tanval;
1794 }
1795
1796
1797 static __m256
1798 gmx_mm256_asin_ps(__m256 x)
1799 {
1800 const __m256 signmask = _mm256_castsi256_ps( _mm256_set1_epi32(0x7FFFFFFF) );
1801 const __m256 limitlow = _mm256_set1_ps(1e-4f);
1802 const __m256 half = _mm256_set1_ps(0.5f);
1803 const __m256 one = _mm256_set1_ps(1.0f);
1804 const __m256 halfpi = _mm256_set1_ps((float)M_PI/2.0f);
1805
1806 const __m256 CC5 = _mm256_set1_ps(4.2163199048E-2f);
1807 const __m256 CC4 = _mm256_set1_ps(2.4181311049E-2f);
1808 const __m256 CC3 = _mm256_set1_ps(4.5470025998E-2f);
1809 const __m256 CC2 = _mm256_set1_ps(7.4953002686E-2f);
1810 const __m256 CC1 = _mm256_set1_ps(1.6666752422E-1f);
1811
1812 __m256 sign;
1813 __m256 mask;
1814 __m256 xabs;
1815 __m256 z,z1,z2,q,q1,q2;
1816 __m256 pA,pB;
1817
1818 sign = _mm256_andnot_ps(signmask,x);
1819 xabs = _mm256_and_ps(x,signmask);
1820
1821 mask = _mm256_cmp_ps(xabs,half,_CMP_GT_OQ);
1822
1823 z1 = _mm256_mul_ps(half, _mm256_sub_ps(one,xabs));
1824 q1 = _mm256_mul_ps(z1,gmx_mm256_invsqrt_ps(z1));
1825 q1 = _mm256_andnot_ps(_mm256_cmp_ps(xabs,one,_CMP_EQ_OQ),q1);
1826
1827 q2 = xabs;
1828 z2 = _mm256_mul_ps(q2,q2);
1829
1830 z = _mm256_blendv_ps(z2,z1,mask);
1831 q = _mm256_blendv_ps(q2,q1,mask);
1832
1833 z2 = _mm256_mul_ps(z,z);
1834
1835 pA = _mm256_mul_ps(CC5,z2);
1836 pB = _mm256_mul_ps(CC4,z2);
1837
1838 pA = _mm256_add_ps(pA,CC3);
1839 pB = _mm256_add_ps(pB,CC2);
1840
1841 pA = _mm256_mul_ps(pA,z2);
1842 pB = _mm256_mul_ps(pB,z2);
1843
1844 pA = _mm256_add_ps(pA,CC1);
1845 pA = _mm256_mul_ps(pA,z);
1846
1847 z = _mm256_add_ps(pA,pB);
1848 z = _mm256_mul_ps(z,q);
1849 z = _mm256_add_ps(z,q);
1850
1851 q2 = _mm256_sub_ps(halfpi,z);
1852 q2 = _mm256_sub_ps(q2,z);
1853
1854 z = _mm256_blendv_ps(z,q2,mask);
1855
1856 mask = _mm256_cmp_ps(xabs,limitlow,_CMP_GT_OQ);
1857 z = _mm256_blendv_ps(xabs,z,mask);
1858
1859 z = _mm256_xor_ps(z,sign);
1860
1861 return z;
1862 }
1863
1864 static __m128
1865 gmx_mm_asin_ps(__m128 x)
1866 {
1867 /* Same algorithm as cephes library */
1868 const __m128 signmask = gmx_mm_castsi128_ps( _mm_set1_epi32(0x7FFFFFFF) );
1869 const __m128 limitlow = _mm_set1_ps(1e-4f);
1870 const __m128 half = _mm_set1_ps(0.5f);
1871 const __m128 one = _mm_set1_ps(1.0f);
1872 const __m128 halfpi = _mm_set1_ps(M_PI/2.0f);
1873
1874 const __m128 CC5 = _mm_set1_ps(4.2163199048E-2f);
1875 const __m128 CC4 = _mm_set1_ps(2.4181311049E-2f);
1876 const __m128 CC3 = _mm_set1_ps(4.5470025998E-2f);
1877 const __m128 CC2 = _mm_set1_ps(7.4953002686E-2f);
1878 const __m128 CC1 = _mm_set1_ps(1.6666752422E-1f);
1879
1880 __m128 sign;
1881 __m128 mask;
1882 __m128 xabs;
1883 __m128 z,z1,z2,q,q1,q2;
1884 __m128 pA,pB;
1885
1886 sign = _mm_andnot_ps(signmask,x);
1887 xabs = _mm_and_ps(x,signmask);
1888
1889 mask = _mm_cmpgt_ps(xabs,half);
1890
1891 z1 = _mm_mul_ps(half, _mm_sub_ps(one,xabs));
1892 q1 = _mm_mul_ps(z1,gmx_mm_invsqrt_ps(z1));
1893 q1 = _mm_andnot_ps(_mm_cmpeq_ps(xabs,one),q1);
1894
1895 q2 = xabs;
1896 z2 = _mm_mul_ps(q2,q2);
1897
1898 z = _mm_or_ps( _mm_and_ps(mask,z1) , _mm_andnot_ps(mask,z2) );
1899 q = _mm_or_ps( _mm_and_ps(mask,q1) , _mm_andnot_ps(mask,q2) );
1900
1901 z2 = _mm_mul_ps(z,z);
1902
1903 pA = _mm_mul_ps(CC5,z2);
1904 pB = _mm_mul_ps(CC4,z2);
1905
1906 pA = _mm_add_ps(pA,CC3);
1907 pB = _mm_add_ps(pB,CC2);
1908
1909 pA = _mm_mul_ps(pA,z2);
1910 pB = _mm_mul_ps(pB,z2);
1911
1912 pA = _mm_add_ps(pA,CC1);
1913 pA = _mm_mul_ps(pA,z);
1914
1915 z = _mm_add_ps(pA,pB);
1916 z = _mm_mul_ps(z,q);
1917 z = _mm_add_ps(z,q);
1918
1919 q2 = _mm_sub_ps(halfpi,z);
1920 q2 = _mm_sub_ps(q2,z);
1921
1922 z = _mm_or_ps( _mm_and_ps(mask,q2) , _mm_andnot_ps(mask,z) );
1923
1924 mask = _mm_cmpgt_ps(xabs,limitlow);
1925 z = _mm_or_ps( _mm_and_ps(mask,z) , _mm_andnot_ps(mask,xabs) );
1926
1927 z = _mm_xor_ps(z,sign);
1928
1929 return z;
1930 }
1931
1932
1933 static __m256
1934 gmx_mm256_acos_ps(__m256 x)
1935 {
1936 const __m256 signmask = _mm256_castsi256_ps( _mm256_set1_epi32(0x7FFFFFFF) );
1937 const __m256 one_ps = _mm256_set1_ps(1.0f);
1938 const __m256 half_ps = _mm256_set1_ps(0.5f);
1939 const __m256 pi_ps = _mm256_set1_ps((float)M_PI);
1940 const __m256 halfpi_ps = _mm256_set1_ps((float)M_PI/2.0f);
1941
1942 __m256 mask1;
1943 __m256 mask2;
1944 __m256 xabs;
1945 __m256 z,z1,z2,z3;
1946
1947 xabs = _mm256_and_ps(x,signmask);
1948 mask1 = _mm256_cmp_ps(xabs,half_ps,_CMP_GT_OQ);
1949 mask2 = _mm256_cmp_ps(x,_mm256_setzero_ps(),_CMP_GT_OQ);
1950
1951 z = _mm256_mul_ps(half_ps,_mm256_sub_ps(one_ps,xabs));
1952 z = _mm256_mul_ps(z,gmx_mm256_invsqrt_ps(z));
1953 z = _mm256_andnot_ps(_mm256_cmp_ps(xabs,one_ps,_CMP_EQ_OQ),z);
1954
1955 z = _mm256_blendv_ps(x,z,mask1);
1956 z = gmx_mm256_asin_ps(z);
1957
1958 z2 = _mm256_add_ps(z,z);
1959 z1 = _mm256_sub_ps(pi_ps,z2);
1960 z3 = _mm256_sub_ps(halfpi_ps,z);
1961
1962 z = _mm256_blendv_ps(z1,z2,mask2);
1963 z = _mm256_blendv_ps(z3,z,mask1);
1964
1965 return z;
1966 }
1967
1968 static __m128
1969 gmx_mm_acos_ps(__m128 x)
1970 {
1971 const __m128 signmask = gmx_mm_castsi128_ps( _mm_set1_epi32(0x7FFFFFFF) );
1972 const __m128 one_ps = _mm_set1_ps(1.0f);
1973 const __m128 half_ps = _mm_set1_ps(0.5f);
1974 const __m128 pi_ps = _mm_set1_ps(M_PI);
1975 const __m128 halfpi_ps = _mm_set1_ps(M_PI/2.0f);
1976
1977 __m128 mask1;
1978 __m128 mask2;
1979 __m128 xabs;
1980 __m128 z,z1,z2,z3;
1981
1982 xabs = _mm_and_ps(x,signmask);
1983 mask1 = _mm_cmpgt_ps(xabs,half_ps);
1984 mask2 = _mm_cmpgt_ps(x,_mm_setzero_ps());
1985
1986 z = _mm_mul_ps(half_ps,_mm_sub_ps(one_ps,xabs));
1987 z = _mm_mul_ps(z,gmx_mm_invsqrt_ps(z));
1988 z = _mm_andnot_ps(_mm_cmpeq_ps(xabs,one_ps),z);
1989
1990 z = _mm_blendv_ps(x,z,mask1);
1991 z = gmx_mm_asin_ps(z);
1992
1993 z2 = _mm_add_ps(z,z);
1994 z1 = _mm_sub_ps(pi_ps,z2);
1995 z3 = _mm_sub_ps(halfpi_ps,z);
1996
1997 z = _mm_blendv_ps(z1,z2,mask2);
1998 z = _mm_blendv_ps(z3,z,mask1);
1999
2000 return z;
2001 }
2002
2003
2004 static __m256
2005 gmx_mm256_atan_ps(__m256 x)
2006 {
2007 const __m256 signmask = _mm256_castsi256_ps( _mm256_set1_epi32(0x7FFFFFFF) );
2008 const __m256 limit1 = _mm256_set1_ps(0.414213562373095f);
2009 const __m256 limit2 = _mm256_set1_ps(2.414213562373095f);
2010 const __m256 quarterpi = _mm256_set1_ps(0.785398163397448f);
2011 const __m256 halfpi = _mm256_set1_ps(1.570796326794896f);
2012 const __m256 mone = _mm256_set1_ps(-1.0f);
2013 const __m256 CC3 = _mm256_set1_ps(-3.33329491539E-1f);
2014 const __m256 CC5 = _mm256_set1_ps(1.99777106478E-1f);
2015 const __m256 CC7 = _mm256_set1_ps(-1.38776856032E-1);
2016 const __m256 CC9 = _mm256_set1_ps(8.05374449538e-2f);
2017
2018 __m256 sign;
2019 __m256 mask1,mask2;
2020 __m256 y,z1,z2;
2021 __m256 x2,x4;
2022 __m256 sum1,sum2;
2023
2024 sign = _mm256_andnot_ps(signmask,x);
2025 x = _mm256_and_ps(x,signmask);
2026
2027 mask1 = _mm256_cmp_ps(x,limit1,_CMP_GT_OQ);
2028 mask2 = _mm256_cmp_ps(x,limit2,_CMP_GT_OQ);
2029
2030 z1 = _mm256_mul_ps(_mm256_add_ps(x,mone),gmx_mm256_inv_ps(_mm256_sub_ps(x,mone)));
2031 z2 = _mm256_mul_ps(mone,gmx_mm256_inv_ps(x));
2032
2033 y = _mm256_and_ps(mask1,quarterpi);
2034 y = _mm256_blendv_ps(y,halfpi,mask2);
2035
2036 x = _mm256_blendv_ps(x,z1,mask1);
2037 x = _mm256_blendv_ps(x,z2,mask2);
2038
2039 x2 = _mm256_mul_ps(x,x);
2040 x4 = _mm256_mul_ps(x2,x2);
2041
2042 sum1 = _mm256_mul_ps(CC9,x4);
2043 sum2 = _mm256_mul_ps(CC7,x4);
2044 sum1 = _mm256_add_ps(sum1,CC5);
2045 sum2 = _mm256_add_ps(sum2,CC3);
2046 sum1 = _mm256_mul_ps(sum1,x4);
2047 sum2 = _mm256_mul_ps(sum2,x2);
2048
2049 sum1 = _mm256_add_ps(sum1,sum2);
2050 sum1 = _mm256_sub_ps(sum1,mone);
2051 sum1 = _mm256_mul_ps(sum1,x);
2052 y = _mm256_add_ps(y,sum1);
2053
2054 y = _mm256_xor_ps(y,sign);
2055
2056 return y;
2057 }
2058
2059 static __m128
2060 gmx_mm_atan_ps(__m128 x)
2061 {
2062 /* Same algorithm as cephes library */
2063 const __m128 signmask = gmx_mm_castsi128_ps( _mm_set1_epi32(0x7FFFFFFF) );
2064 const __m128 limit1 = _mm_set1_ps(0.414213562373095f);
2065 const __m128 limit2 = _mm_set1_ps(2.414213562373095f);
2066 const __m128 quarterpi = _mm_set1_ps(0.785398163397448f);
2067 const __m128 halfpi = _mm_set1_ps(1.570796326794896f);
2068 const __m128 mone = _mm_set1_ps(-1.0f);
2069 const __m128 CC3 = _mm_set1_ps(-3.33329491539E-1f);
2070 const __m128 CC5 = _mm_set1_ps(1.99777106478E-1f);
2071 const __m128 CC7 = _mm_set1_ps(-1.38776856032E-1);
2072 const __m128 CC9 = _mm_set1_ps(8.05374449538e-2f);
2073
2074 __m128 sign;
2075 __m128 mask1,mask2;
2076 __m128 y,z1,z2;
2077 __m128 x2,x4;
2078 __m128 sum1,sum2;
2079
2080 sign = _mm_andnot_ps(signmask,x);
2081 x = _mm_and_ps(x,signmask);
2082
2083 mask1 = _mm_cmpgt_ps(x,limit1);
2084 mask2 = _mm_cmpgt_ps(x,limit2);
2085
2086 z1 = _mm_mul_ps(_mm_add_ps(x,mone),gmx_mm_inv_ps(_mm_sub_ps(x,mone)));
2087 z2 = _mm_mul_ps(mone,gmx_mm_inv_ps(x));
2088
2089 y = _mm_and_ps(mask1,quarterpi);
2090 y = _mm_blendv_ps(y,halfpi,mask2);
2091
2092 x = _mm_blendv_ps(x,z1,mask1);
2093 x = _mm_blendv_ps(x,z2,mask2);
2094
2095 x2 = _mm_mul_ps(x,x);
2096 x4 = _mm_mul_ps(x2,x2);
2097
2098 sum1 = _mm_mul_ps(CC9,x4);
2099 sum2 = _mm_mul_ps(CC7,x4);
2100 sum1 = _mm_add_ps(sum1,CC5);
2101 sum2 = _mm_add_ps(sum2,CC3);
2102 sum1 = _mm_mul_ps(sum1,x4);
2103 sum2 = _mm_mul_ps(sum2,x2);
2104
2105 sum1 = _mm_add_ps(sum1,sum2);
2106 sum1 = _mm_sub_ps(sum1,mone);
2107 sum1 = _mm_mul_ps(sum1,x);
2108 y = _mm_add_ps(y,sum1);
2109
2110 y = _mm_xor_ps(y,sign);
2111
2112 return y;
2113 }
2114
2115
2116 static __m256
2117 gmx_mm256_atan2_ps(__m256 y, __m256 x)
2118 {
2119 const __m256 pi = _mm256_set1_ps( (float) M_PI);
2120 const __m256 minuspi = _mm256_set1_ps( (float) -M_PI);
2121 const __m256 halfpi = _mm256_set1_ps( (float) M_PI/2.0f);
2122 const __m256 minushalfpi = _mm256_set1_ps( (float) -M_PI/2.0f);
2123
2124 __m256 z,z1,z3,z4;
2125 __m256 w;
2126 __m256 maskx_lt,maskx_eq;
2127 __m256 masky_lt,masky_eq;
2128 __m256 mask1,mask2,mask3,mask4,maskall;
2129
2130 maskx_lt = _mm256_cmp_ps(x,_mm256_setzero_ps(),_CMP_LT_OQ);
2131 masky_lt = _mm256_cmp_ps(y,_mm256_setzero_ps(),_CMP_LT_OQ);
2132 maskx_eq = _mm256_cmp_ps(x,_mm256_setzero_ps(),_CMP_EQ_OQ);
2133 masky_eq = _mm256_cmp_ps(y,_mm256_setzero_ps(),_CMP_EQ_OQ);
2134
2135 z = _mm256_mul_ps(y,gmx_mm256_inv_ps(x));
2136 z = gmx_mm256_atan_ps(z);
2137
2138 mask1 = _mm256_and_ps(maskx_eq,masky_lt);
2139 mask2 = _mm256_andnot_ps(maskx_lt,masky_eq);
2140 mask3 = _mm256_andnot_ps( _mm256_or_ps(masky_lt,masky_eq) , maskx_eq);
2141 mask4 = _mm256_and_ps(maskx_lt,masky_eq);
2142 maskall = _mm256_or_ps( _mm256_or_ps(mask1,mask2), _mm256_or_ps(mask3,mask4) );
2143
2144 z = _mm256_andnot_ps(maskall,z);
2145 z1 = _mm256_and_ps(mask1,minushalfpi);
2146 z3 = _mm256_and_ps(mask3,halfpi);
2147 z4 = _mm256_and_ps(mask4,pi);
2148
2149 z = _mm256_or_ps( _mm256_or_ps(z,z1), _mm256_or_ps(z3,z4) );
2150
2151 w = _mm256_blendv_ps(pi,minuspi,masky_lt);
2152 w = _mm256_and_ps(w,maskx_lt);
2153
2154 w = _mm256_andnot_ps(maskall,w);
2155
2156 z = _mm256_add_ps(z,w);
2157
2158 return z;
2159 }
2160
2161 static __m128
2162 gmx_mm_atan2_ps(__m128 y, __m128 x)
2163 {
2164 const __m128 pi = _mm_set1_ps(M_PI);
2165 const __m128 minuspi = _mm_set1_ps(-M_PI);
2166 const __m128 halfpi = _mm_set1_ps(M_PI/2.0);
2167 const __m128 minushalfpi = _mm_set1_ps(-M_PI/2.0);
2168
2169 __m128 z,z1,z3,z4;
2170 __m128 w;
2171 __m128 maskx_lt,maskx_eq;
2172 __m128 masky_lt,masky_eq;
2173 __m128 mask1,mask2,mask3,mask4,maskall;
2174
2175 maskx_lt = _mm_cmplt_ps(x,_mm_setzero_ps());
2176 masky_lt = _mm_cmplt_ps(y,_mm_setzero_ps());
2177 maskx_eq = _mm_cmpeq_ps(x,_mm_setzero_ps());
2178 masky_eq = _mm_cmpeq_ps(y,_mm_setzero_ps());
2179
2180 z = _mm_mul_ps(y,gmx_mm_inv_ps(x));
2181 z = gmx_mm_atan_ps(z);
2182
2183 mask1 = _mm_and_ps(maskx_eq,masky_lt);
2184 mask2 = _mm_andnot_ps(maskx_lt,masky_eq);
2185 mask3 = _mm_andnot_ps( _mm_or_ps(masky_lt,masky_eq) , maskx_eq);
2186 mask4 = _mm_and_ps(masky_eq,maskx_lt);
2187
2188 maskall = _mm_or_ps( _mm_or_ps(mask1,mask2), _mm_or_ps(mask3,mask4) );
2189
2190 z = _mm_andnot_ps(maskall,z);
2191 z1 = _mm_and_ps(mask1,minushalfpi);
2192 z3 = _mm_and_ps(mask3,halfpi);
2193 z4 = _mm_and_ps(mask4,pi);
2194
2195 z = _mm_or_ps( _mm_or_ps(z,z1), _mm_or_ps(z3,z4) );
2196
2197 mask1 = _mm_andnot_ps(masky_lt,maskx_lt);
2198 mask2 = _mm_and_ps(maskx_lt,masky_lt);
2199
2200 w = _mm_or_ps( _mm_and_ps(mask1,pi), _mm_and_ps(mask2,minuspi) );
2201 w = _mm_andnot_ps(maskall,w);
2202
2203 z = _mm_add_ps(z,w);
2204
2205 return z;
2206 }
2207
2208 #endif /* _gmx_math_x86_avx_256_single_h_ */
The ultimate molecular dynamics simulation package