EigenRand  0.2.0
MorePacketMath.h
Go to the documentation of this file.
1 
12 #ifndef EIGENRAND_MORE_PACKET_MATH_H
13 #define EIGENRAND_MORE_PACKET_MATH_H
14 
15 #include <Eigen/Dense>
16 
17 namespace Eigen
18 {
19  namespace internal
20  {
21  template<typename Packet>
22  struct reinterpreter
23  {
24  };
25 
26  template<typename Packet>
27  inline auto reinterpret_to_float(const Packet& x)
28  -> decltype(reinterpreter<Packet>{}.to_float(x))
29  {
30  return reinterpreter<Packet>{}.to_float(x);
31  }
32 
33  template<typename Packet>
34  inline auto reinterpret_to_double(const Packet& x)
35  -> decltype(reinterpreter<Packet>{}.to_double(x))
36  {
37  return reinterpreter<Packet>{}.to_double(x);
38  }
39 
40  template<typename Packet>
41  inline auto reinterpret_to_int(const Packet& x)
42  -> decltype(reinterpreter<Packet>{}.to_int(x))
43  {
44  return reinterpreter<Packet>{}.to_int(x);
45  }
46 
47  template<typename Packet>
48  EIGEN_STRONG_INLINE Packet pseti64(uint64_t a);
49 
50  template<typename Packet>
51  EIGEN_STRONG_INLINE Packet pcmpeq(const Packet& a, const Packet& b);
52 
53  template<typename Packet>
54  EIGEN_STRONG_INLINE Packet psll(const Packet& a, int b);
55 
56  template<typename Packet>
57  EIGEN_STRONG_INLINE Packet psrl(const Packet& a, int b);
58 
59  template<typename Packet>
60  EIGEN_STRONG_INLINE Packet psll64(const Packet& a, int b);
61 
62  template<typename Packet>
63  EIGEN_STRONG_INLINE Packet psrl64(const Packet& a, int b);
64 
65  template<typename Packet>
66  EIGEN_STRONG_INLINE int pmovemask(const Packet& a);
67 
68  template<typename Packet>
69  EIGEN_STRONG_INLINE void psincos(Packet x, Packet &s, Packet &c)
70  {
71  Packet xmm1, xmm2, xmm3 = pset1<Packet>(0), sign_bit_sin, y;
72  using IntPacket = decltype(reinterpret_to_int(x));
73  IntPacket emm0, emm2, emm4;
74 
75  sign_bit_sin = x;
76  /* take the absolute value */
77  x = pabs(x);
78  /* extract the sign bit (upper one) */
79  sign_bit_sin = reinterpret_to_float(
80  pand(reinterpret_to_int(sign_bit_sin), pset1<IntPacket>(0x80000000))
81  );
82 
83  /* scale by 4/Pi */
84  y = pmul(x, pset1<Packet>(1.27323954473516));
85 
86  /* store the integer part of y in emm2 */
87  emm2 = pcast<Packet, IntPacket>(y);
88 
89  /* j=(j+1) & (~1) (see the cephes sources) */
90  emm2 = padd(emm2, pset1<IntPacket>(1));
91  emm2 = pand(emm2, pset1<IntPacket>(~1));
92  y = pcast<IntPacket, Packet>(emm2);
93 
94  emm4 = emm2;
95 
96  /* get the swap sign flag for the sine */
97  emm0 = pand(emm2, pset1<IntPacket>(4));
98  emm0 = psll(emm0, 29);
99  Packet swap_sign_bit_sin = reinterpret_to_float(emm0);
100 
101  /* get the polynom selection mask for the sine*/
102  emm2 = pand(emm2, pset1<IntPacket>(2));
103 
104  emm2 = pcmpeq(emm2, pset1<IntPacket>(0));
105  Packet poly_mask = reinterpret_to_float(emm2);
106 
107  /* The magic pass: "Extended precision modular arithmetic"
108  x = ((x - y * DP1) - y * DP2) - y * DP3; */
109  xmm1 = pset1<Packet>(-0.78515625);
110  xmm2 = pset1<Packet>(-2.4187564849853515625e-4);
111  xmm3 = pset1<Packet>(-3.77489497744594108e-8);
112  xmm1 = pmul(y, xmm1);
113  xmm2 = pmul(y, xmm2);
114  xmm3 = pmul(y, xmm3);
115  x = padd(x, xmm1);
116  x = padd(x, xmm2);
117  x = padd(x, xmm3);
118 
119  emm4 = psub(emm4, pset1<IntPacket>(2));
120  emm4 = pandnot(emm4, pset1<IntPacket>(4));
121  emm4 = psll(emm4, 29);
122  Packet sign_bit_cos = reinterpret_to_float(emm4);
123  sign_bit_sin = pxor(sign_bit_sin, swap_sign_bit_sin);
124 
125 
126  /* Evaluate the first polynom (0 <= x <= Pi/4) */
127  Packet z = pmul(x, x);
128  y = pset1<Packet>(2.443315711809948E-005);
129 
130  y = pmul(y, z);
131  y = padd(y, pset1<Packet>(-1.388731625493765E-003));
132  y = pmul(y, z);
133  y = padd(y, pset1<Packet>(4.166664568298827E-002));
134  y = pmul(y, z);
135  y = pmul(y, z);
136  Packet tmp = pmul(z, pset1<Packet>(0.5));
137  y = psub(y, tmp);
138  y = padd(y, pset1<Packet>(1));
139 
140  /* Evaluate the second polynom (Pi/4 <= x <= 0) */
141 
142  Packet y2 = pset1<Packet>(-1.9515295891E-4);
143  y2 = pmul(y2, z);
144  y2 = padd(y2, pset1<Packet>(8.3321608736E-3));
145  y2 = pmul(y2, z);
146  y2 = padd(y2, pset1<Packet>(-1.6666654611E-1));
147  y2 = pmul(y2, z);
148  y2 = pmul(y2, x);
149  y2 = padd(y2, x);
150 
151  /* select the correct result from the two polynoms */
152  xmm3 = poly_mask;
153  Packet ysin2 = pand(xmm3, y2);
154  Packet ysin1 = pandnot(xmm3, y);
155  y2 = psub(y2, ysin2);
156  y = psub(y, ysin1);
157 
158  xmm1 = padd(ysin1, ysin2);
159  xmm2 = padd(y, y2);
160 
161  /* update the sign */
162  s = pxor(xmm1, sign_bit_sin);
163  c = pxor(xmm2, sign_bit_cos);
164  }
165 
166  // approximation : lgamma(z) ~= (z+2.5)ln(z+3) - z - 3 + 0.5 ln (2pi) + 1/12/(z + 3) - ln (z(z+1)(z+2))
167  template<typename Packet>
168  EIGEN_STRONG_INLINE Packet plgamma(const Packet& x)
169  {
170  auto x_3 = padd(x, pset1<Packet>(3));
171  auto ret = pmul(padd(x_3, pset1<Packet>(-0.5)), plog(x_3));
172  ret = psub(ret, x_3);
173  ret = padd(ret, pset1<Packet>(0.9189385332046727));
174  ret = padd(ret, pdiv(pset1<Packet>(1 / 12.), x_3));
175  ret = psub(ret, plog(pmul(
176  pmul(psub(x_3, pset1<Packet>(1)), psub(x_3, pset1<Packet>(2))), x)));
177  return ret;
178  }
179 
180  template<typename Packet>
181  EIGEN_STRONG_INLINE Packet pcmplt(const Packet& a, const Packet& b);
182 
183  template<typename Packet>
184  EIGEN_STRONG_INLINE Packet pcmple(const Packet& a, const Packet& b);
185 
186  template<typename Packet>
187  EIGEN_STRONG_INLINE Packet pblendv(const Packet& ifPacket, const Packet& thenPacket, const Packet& elsePacket);
188 
189  template<typename Packet>
190  EIGEN_STRONG_INLINE Packet pgather(const int* addr, const Packet& index);
191 
192  template<typename Packet>
193  EIGEN_STRONG_INLINE auto pgather(const float* addr, const Packet& index) -> decltype(reinterpret_to_float(std::declval<Packet>()));
194 
195  template<typename Packet>
196  EIGEN_STRONG_INLINE auto pgather(const double* addr, const Packet& index, bool upperhalf = false) -> decltype(reinterpret_to_double(std::declval<Packet>()));
197 
198  template<typename Packet>
199  EIGEN_STRONG_INLINE Packet ptruncate(const Packet& a);
200 
201  template<typename IntPacket>
202  EIGEN_STRONG_INLINE auto bit_to_ur_float(const IntPacket& x) -> decltype(reinterpret_to_float(x))
203  {
204  using FloatPacket = decltype(reinterpret_to_float(x));
205 
206  const IntPacket lower = pset1<IntPacket>(0x7FFFFF),
207  upper = pset1<IntPacket>(127 << 23);
208  const FloatPacket one = pset1<FloatPacket>(1);
209 
210  return psub(reinterpret_to_float(por(pand(x, lower), upper)), one);
211  }
212 
213  template<typename IntPacket>
214  EIGEN_STRONG_INLINE auto bit_to_ur_double(const IntPacket& x) -> decltype(reinterpret_to_double(x))
215  {
216  using DoublePacket = decltype(reinterpret_to_double(x));
217 
218  const IntPacket lower = pseti64<IntPacket>(0xFFFFFFFFFFFFFull),
219  upper = pseti64<IntPacket>(1023ull << 52);
220  const DoublePacket one = pset1<DoublePacket>(1);
221 
222  return psub(reinterpret_to_double(por(pand(x, lower), upper)), one);
223  }
224 
225  template<typename Scalar>
226  struct bit_scalar;
227 
228  template<>
229  struct bit_scalar<float>
230  {
231  float to_ur(uint32_t x)
232  {
233  union
234  {
235  uint32_t u;
236  float f;
237  };
238  u = (x & 0x7FFFFF) | (127 << 23);
239  return f - 1.f;
240  }
241 
242  float to_nzur(uint32_t x)
243  {
244  return to_ur(x) + std::numeric_limits<float>::epsilon() / 8;
245  }
246  };
247 
248  template<>
249  struct bit_scalar<double>
250  {
251  double to_ur(uint64_t x)
252  {
253  union
254  {
255  uint64_t u;
256  double f;
257  };
258  u = (x & 0xFFFFFFFFFFFFFull) | (1023ull << 52);
259  return f - 1.;
260  }
261 
262  double to_nzur(uint64_t x)
263  {
264  return to_ur(x) + std::numeric_limits<double>::epsilon() / 8;
265  }
266  };
267  }
268 }
269 
270 #ifdef EIGEN_VECTORIZE_AVX
271 #include <immintrin.h>
272 
273 namespace Eigen
274 {
275  namespace internal
276  {
277  template<>
278  struct reinterpreter<Packet8i>
279  {
280  EIGEN_STRONG_INLINE Packet8f to_float(const Packet8i& x)
281  {
282  return _mm256_castsi256_ps(x);
283  }
284 
285  EIGEN_STRONG_INLINE Packet4d to_double(const Packet8i& x)
286  {
287  return _mm256_castsi256_pd(x);
288  }
289 
290  EIGEN_STRONG_INLINE Packet8i to_int(const Packet8i& x)
291  {
292  return x;
293  }
294  };
295 
296  template<>
297  struct reinterpreter<Packet8f>
298  {
299  EIGEN_STRONG_INLINE Packet8f to_float(const Packet8f& x)
300  {
301  return x;
302  }
303 
304  EIGEN_STRONG_INLINE Packet4d to_double(const Packet8f& x)
305  {
306  return _mm256_castps_pd(x);
307  }
308 
309  EIGEN_STRONG_INLINE Packet8i to_int(const Packet8f& x)
310  {
311  return _mm256_castps_si256(x);
312  }
313  };
314 
315  template<>
316  struct reinterpreter<Packet4d>
317  {
318  EIGEN_STRONG_INLINE Packet8f to_float(const Packet4d& x)
319  {
320  return _mm256_castpd_ps(x);
321  }
322 
323  EIGEN_STRONG_INLINE Packet4d to_double(const Packet4d& x)
324  {
325  return x;
326  }
327 
328  EIGEN_STRONG_INLINE Packet8i to_int(const Packet4d& x)
329  {
330  return _mm256_castpd_si256(x);
331  }
332  };
333 
334  EIGEN_STRONG_INLINE void split_two(const Packet8i& x, Packet4i& a, Packet4i& b)
335  {
336  a = _mm256_extractf128_si256(x, 0);
337  b = _mm256_extractf128_si256(x, 1);
338  }
339 
340  EIGEN_STRONG_INLINE Packet8i combine_two(const Packet4i& a, const Packet4i& b)
341  {
342  return _mm256_insertf128_si256(_mm256_castsi128_si256(a), b, 1);
343  }
344 
345  EIGEN_STRONG_INLINE void split_two(const Packet8f& x, Packet4f& a, Packet4f& b)
346  {
347  a = _mm256_extractf128_ps(x, 0);
348  b = _mm256_extractf128_ps(x, 1);
349  }
350 
351  EIGEN_STRONG_INLINE Packet8f combine_two(const Packet4f& a, const Packet4f& b)
352  {
353  return _mm256_insertf128_ps(_mm256_castps128_ps256(a), b, 1);
354  }
355 
356 
357  EIGEN_STRONG_INLINE Packet4i combine_low32(const Packet8i& a)
358  {
359 #ifdef EIGEN_VECTORIZE_AVX2
360  return _mm256_castsi256_si128(_mm256_permutevar8x32_epi32(a, _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7)));
361 #else
362  auto sc = _mm256_permutevar_ps(_mm256_castsi256_ps(a), _mm256_setr_epi32(0, 2, 1, 3, 1, 3, 0, 2));
363  return _mm_castps_si128(_mm_blend_ps(_mm256_extractf128_ps(sc, 0), _mm256_extractf128_ps(sc, 1), 0b1100));
364 #endif
365  }
366 
367  template<>
368  EIGEN_STRONG_INLINE Packet8i pseti64<Packet8i>(uint64_t a)
369  {
370  return _mm256_set1_epi64x(a);
371  }
372 
373  template<>
374  EIGEN_STRONG_INLINE Packet8i pcmpeq<Packet8i>(const Packet8i& a, const Packet8i& b)
375  {
376 #ifdef EIGEN_VECTORIZE_AVX2
377  return _mm256_cmpeq_epi32(a, b);
378 #else
379  Packet4i a1, a2, b1, b2;
380  split_two(a, a1, a2);
381  split_two(b, b1, b2);
382  return combine_two(_mm_cmpeq_epi32(a1, b1), _mm_cmpeq_epi32(a2, b2));
383 #endif
384  }
385 
386  template<>
387  EIGEN_STRONG_INLINE Packet8i psll<Packet8i>(const Packet8i& a, int b)
388  {
389 #ifdef EIGEN_VECTORIZE_AVX2
390  return _mm256_slli_epi32(a, b);
391 #else
392  Packet4i a1, a2;
393  split_two(a, a1, a2);
394  return combine_two(_mm_slli_epi32(a1, b), _mm_slli_epi32(a2, b));
395 #endif
396  }
397 
398  template<>
399  EIGEN_STRONG_INLINE Packet8i psrl<Packet8i>(const Packet8i& a, int b)
400  {
401 #ifdef EIGEN_VECTORIZE_AVX2
402  return _mm256_srli_epi32(a, b);
403 #else
404  Packet4i a1, a2;
405  split_two(a, a1, a2);
406  return combine_two(_mm_srli_epi32(a1, b), _mm_srli_epi32(a2, b));
407 #endif
408  }
409 
410  template<>
411  EIGEN_STRONG_INLINE Packet8i psll64<Packet8i>(const Packet8i& a, int b)
412  {
413 #ifdef EIGEN_VECTORIZE_AVX2
414  return _mm256_slli_epi64(a, b);
415 #else
416  Packet4i a1, a2;
417  split_two(a, a1, a2);
418  return combine_two(_mm_slli_epi64(a1, b), _mm_slli_epi64(a2, b));
419 #endif
420  }
421 
422  template<>
423  EIGEN_STRONG_INLINE Packet8i psrl64<Packet8i>(const Packet8i& a, int b)
424  {
425 #ifdef EIGEN_VECTORIZE_AVX2
426  return _mm256_srli_epi64(a, b);
427 #else
428  Packet4i a1, a2;
429  split_two(a, a1, a2);
430  return combine_two(_mm_srli_epi64(a1, b), _mm_srli_epi64(a2, b));
431 #endif
432  }
433 
434  template<> EIGEN_STRONG_INLINE Packet8i padd<Packet8i>(const Packet8i& a, const Packet8i& b)
435  {
436 #ifdef EIGEN_VECTORIZE_AVX2
437  return _mm256_add_epi32(a, b);
438 #else
439  Packet4i a1, a2, b1, b2;
440  split_two(a, a1, a2);
441  split_two(b, b1, b2);
442  return combine_two(_mm_add_epi32(a1, b1), _mm_add_epi32(a2, b2));
443 #endif
444  }
445 
446  template<> EIGEN_STRONG_INLINE Packet8i psub<Packet8i>(const Packet8i& a, const Packet8i& b)
447  {
448 #ifdef EIGEN_VECTORIZE_AVX2
449  return _mm256_sub_epi32(a, b);
450 #else
451  Packet4i a1, a2, b1, b2;
452  split_two(a, a1, a2);
453  split_two(b, b1, b2);
454  return combine_two(_mm_sub_epi32(a1, b1), _mm_sub_epi32(a2, b2));
455 #endif
456  }
457 
458  template<> EIGEN_STRONG_INLINE Packet8i pand<Packet8i>(const Packet8i& a, const Packet8i& b)
459  {
460 #ifdef EIGEN_VECTORIZE_AVX2
461  return _mm256_and_si256(a, b);
462 #else
463  return reinterpret_to_int(_mm256_and_ps(reinterpret_to_float(a), reinterpret_to_float(b)));
464 #endif
465  }
466 
467  template<> EIGEN_STRONG_INLINE Packet8i pandnot<Packet8i>(const Packet8i& a, const Packet8i& b)
468  {
469 #ifdef EIGEN_VECTORIZE_AVX2
470  return _mm256_andnot_si256(a, b);
471 #else
472  return reinterpret_to_int(_mm256_andnot_ps(reinterpret_to_float(a), reinterpret_to_float(b)));
473 #endif
474  }
475 
476  template<> EIGEN_STRONG_INLINE Packet8i por<Packet8i>(const Packet8i& a, const Packet8i& b)
477  {
478 #ifdef EIGEN_VECTORIZE_AVX2
479  return _mm256_or_si256(a, b);
480 #else
481  return reinterpret_to_int(_mm256_or_ps(reinterpret_to_float(a), reinterpret_to_float(b)));
482 #endif
483  }
484 
485  template<> EIGEN_STRONG_INLINE Packet8i pxor<Packet8i>(const Packet8i& a, const Packet8i& b)
486  {
487 #ifdef EIGEN_VECTORIZE_AVX2
488  return _mm256_xor_si256(a, b);
489 #else
490  return reinterpret_to_int(_mm256_xor_ps(reinterpret_to_float(a), reinterpret_to_float(b)));
491 #endif
492  }
493 
494  template<>
495  EIGEN_STRONG_INLINE Packet8i pcmplt<Packet8i>(const Packet8i& a, const Packet8i& b)
496  {
497  return _mm256_cmpgt_epi32(b, a);
498  }
499 
500  template<>
501  EIGEN_STRONG_INLINE Packet8f pcmplt<Packet8f>(const Packet8f& a, const Packet8f& b)
502  {
503  return _mm256_cmp_ps(a, b, _CMP_LT_OQ);
504  }
505 
506  template<>
507  EIGEN_STRONG_INLINE Packet8f pcmple<Packet8f>(const Packet8f& a, const Packet8f& b)
508  {
509  return _mm256_cmp_ps(a, b, _CMP_LE_OQ);
510  }
511 
512  template<>
513  EIGEN_STRONG_INLINE Packet4d pcmplt<Packet4d>(const Packet4d& a, const Packet4d& b)
514  {
515  return _mm256_cmp_pd(a, b, _CMP_LT_OQ);
516  }
517 
518  template<>
519  EIGEN_STRONG_INLINE Packet4d pcmple<Packet4d>(const Packet4d& a, const Packet4d& b)
520  {
521  return _mm256_cmp_pd(a, b, _CMP_LE_OQ);
522  }
523 
524  template<>
525  EIGEN_STRONG_INLINE Packet8f pblendv(const Packet8f& ifPacket, const Packet8f& thenPacket, const Packet8f& elsePacket)
526  {
527  return _mm256_blendv_ps(elsePacket, thenPacket, ifPacket);
528  }
529 
530  template<>
531  EIGEN_STRONG_INLINE Packet8i pblendv(const Packet8i& ifPacket, const Packet8i& thenPacket, const Packet8i& elsePacket)
532  {
533  return _mm256_castps_si256(_mm256_blendv_ps(
534  _mm256_castsi256_ps(elsePacket),
535  _mm256_castsi256_ps(thenPacket),
536  _mm256_castsi256_ps(ifPacket)
537  ));
538  }
539 
540  template<>
541  EIGEN_STRONG_INLINE Packet4d pblendv(const Packet4d& ifPacket, const Packet4d& thenPacket, const Packet4d& elsePacket)
542  {
543  return _mm256_blendv_pd(elsePacket, thenPacket, ifPacket);
544  }
545 
546  template<>
547  EIGEN_STRONG_INLINE Packet8i pgather<Packet8i>(const int* addr, const Packet8i& index)
548  {
549 #ifdef EIGEN_VECTORIZE_AVX2
550  return _mm256_i32gather_epi32(addr, index, 4);
551 #else
552  uint32_t u[8];
553  _mm256_storeu_si256((Packet8i*)u, index);
554  return _mm256_setr_epi32(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]],
555  addr[u[4]], addr[u[5]], addr[u[6]], addr[u[7]]);
556 #endif
557  }
558 
559  template<>
560  EIGEN_STRONG_INLINE Packet8f pgather<Packet8i>(const float *addr, const Packet8i& index)
561  {
562 #ifdef EIGEN_VECTORIZE_AVX2
563  return _mm256_i32gather_ps(addr, index, 4);
564 #else
565  uint32_t u[8];
566  _mm256_storeu_si256((Packet8i*)u, index);
567  return _mm256_setr_ps(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]],
568  addr[u[4]], addr[u[5]], addr[u[6]], addr[u[7]]);
569 #endif
570  }
571 
572  template<>
573  EIGEN_STRONG_INLINE Packet4d pgather<Packet8i>(const double *addr, const Packet8i& index, bool upperhalf)
574  {
575 #ifdef EIGEN_VECTORIZE_AVX2
576  return _mm256_i32gather_pd(addr, _mm256_castsi256_si128(index), 8);
577 #else
578  uint32_t u[8];
579  _mm256_storeu_si256((Packet8i*)u, index);
580  if (upperhalf)
581  {
582  return _mm256_setr_pd(addr[u[4]], addr[u[5]], addr[u[6]], addr[u[7]]);
583  }
584  else
585  {
586  return _mm256_setr_pd(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]]);
587  }
588 #endif
589  }
590 
591  template<>
592  EIGEN_STRONG_INLINE int pmovemask<Packet8f>(const Packet8f& a)
593  {
594  return _mm256_movemask_ps(a);
595  }
596 
597  template<>
598  EIGEN_STRONG_INLINE int pmovemask<Packet4d>(const Packet4d& a)
599  {
600  return _mm256_movemask_pd(a);
601  }
602 
603  template<>
604  EIGEN_STRONG_INLINE int pmovemask<Packet8i>(const Packet8i& a)
605  {
606  return pmovemask(_mm256_castsi256_ps(a));
607  }
608 
609  template<>
610  EIGEN_STRONG_INLINE Packet8f ptruncate<Packet8f>(const Packet8f& a)
611  {
612  return _mm256_round_ps(a, _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC);
613  }
614 
615  template<>
616  EIGEN_STRONG_INLINE Packet4d ptruncate<Packet4d>(const Packet4d& a)
617  {
618  return _mm256_round_pd(a, _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC);
619  }
620  }
621 }
622 #endif
623 
624 #ifdef EIGEN_VECTORIZE_SSE2
625 #include <xmmintrin.h>
626 
627 namespace Eigen
628 {
629  namespace internal
630  {
631  template<>
632  struct reinterpreter<Packet4i>
633  {
634  EIGEN_STRONG_INLINE Packet4f to_float(const Packet4i& x)
635  {
636  return _mm_castsi128_ps(x);
637  }
638 
639  EIGEN_STRONG_INLINE Packet2d to_double(const Packet4i& x)
640  {
641  return _mm_castsi128_pd(x);
642  }
643 
644  EIGEN_STRONG_INLINE Packet4i to_int(const Packet4i& x)
645  {
646  return x;
647  }
648  };
649 
650  template<>
651  struct reinterpreter<Packet4f>
652  {
653  EIGEN_STRONG_INLINE Packet4f to_float(const Packet4f& x)
654  {
655  return x;
656  }
657 
658  EIGEN_STRONG_INLINE Packet2d to_double(const Packet4f& x)
659  {
660  return _mm_castps_pd(x);
661  }
662 
663  EIGEN_STRONG_INLINE Packet4i to_int(const Packet4f& x)
664  {
665  return _mm_castps_si128(x);
666  }
667  };
668 
669  template<>
670  struct reinterpreter<Packet2d>
671  {
672  EIGEN_STRONG_INLINE Packet4f to_float(const Packet2d& x)
673  {
674  return _mm_castpd_ps(x);
675  }
676 
677  EIGEN_STRONG_INLINE Packet2d to_double(const Packet2d& x)
678  {
679  return x;
680  }
681 
682  EIGEN_STRONG_INLINE Packet4i to_int(const Packet2d& x)
683  {
684  return _mm_castpd_si128(x);
685  }
686  };
687 
688  EIGEN_STRONG_INLINE void split_two(const Packet4i& x, uint64_t& a, uint64_t& b)
689  {
690 #ifdef EIGEN_VECTORIZE_SSE4_1
691  a = _mm_extract_epi64(x, 0);
692  b = _mm_extract_epi64(x, 1);
693 #else
694  uint64_t u[2];
695  _mm_storeu_si128((__m128i*)u, x);
696  a = u[0];
697  b = u[1];
698 #endif
699  }
700 
701  EIGEN_STRONG_INLINE Packet4i combine_low32(const Packet4i& a, const Packet4i& b)
702  {
703  auto sa = _mm_shuffle_epi32(a, _MM_SHUFFLE(3, 1, 2, 0));
704  auto sb = _mm_shuffle_epi32(b, _MM_SHUFFLE(2, 0, 3, 1));
705  sa = _mm_and_si128(sa, _mm_setr_epi32(-1, -1, 0, 0));
706  sb = _mm_and_si128(sb, _mm_setr_epi32(0, 0, -1, -1));
707  return _mm_or_si128(sa, sb);
708  }
709 
710  template<>
711  EIGEN_STRONG_INLINE Packet4i pseti64<Packet4i>(uint64_t a)
712  {
713  return _mm_set1_epi64x(a);
714  }
715 
716  template<>
717  EIGEN_STRONG_INLINE Packet4i pcmpeq<Packet4i>(const Packet4i& a, const Packet4i& b)
718  {
719  return _mm_cmpeq_epi32(a, b);
720  }
721 
722  template<>
723  EIGEN_STRONG_INLINE Packet4i psll<Packet4i>(const Packet4i& a, int b)
724  {
725  return _mm_slli_epi32(a, b);
726  }
727 
728  template<>
729  EIGEN_STRONG_INLINE Packet4i psrl<Packet4i>(const Packet4i& a, int b)
730  {
731  return _mm_srli_epi32(a, b);
732  }
733 
734 
735  template<>
736  EIGEN_STRONG_INLINE Packet4i psll64<Packet4i>(const Packet4i& a, int b)
737  {
738  return _mm_slli_epi64(a, b);
739  }
740 
741  template<>
742  EIGEN_STRONG_INLINE Packet4i psrl64<Packet4i>(const Packet4i& a, int b)
743  {
744  return _mm_srli_epi64(a, b);
745  }
746 
747  template<>
748  EIGEN_STRONG_INLINE Packet4i pcmplt<Packet4i>(const Packet4i& a, const Packet4i& b)
749  {
750  return _mm_cmplt_epi32(a, b);
751  }
752 
753  template<>
754  EIGEN_STRONG_INLINE Packet4f pcmplt<Packet4f>(const Packet4f& a, const Packet4f& b)
755  {
756  return _mm_cmplt_ps(a, b);
757  }
758 
759  template<>
760  EIGEN_STRONG_INLINE Packet4f pcmple<Packet4f>(const Packet4f& a, const Packet4f& b)
761  {
762  return _mm_cmple_ps(a, b);
763  }
764 
765  template<>
766  EIGEN_STRONG_INLINE Packet2d pcmplt<Packet2d>(const Packet2d& a, const Packet2d& b)
767  {
768  return _mm_cmplt_pd(a, b);
769  }
770 
771  template<>
772  EIGEN_STRONG_INLINE Packet2d pcmple<Packet2d>(const Packet2d& a, const Packet2d& b)
773  {
774  return _mm_cmple_pd(a, b);
775  }
776 
777  template<>
778  EIGEN_STRONG_INLINE Packet4f pblendv(const Packet4f& ifPacket, const Packet4f& thenPacket, const Packet4f& elsePacket)
779  {
780 #ifdef EIGEN_VECTORIZE_SSE4_1
781  return _mm_blendv_ps(elsePacket, thenPacket, ifPacket);
782 #else
783  return _mm_or_ps(_mm_and_ps(ifPacket, thenPacket), _mm_andnot_ps(ifPacket, elsePacket));
784 #endif
785  }
786 
787  template<>
788  EIGEN_STRONG_INLINE Packet4i pblendv(const Packet4i& ifPacket, const Packet4i& thenPacket, const Packet4i& elsePacket)
789  {
790 #ifdef EIGEN_VECTORIZE_SSE4_1
791  return _mm_castps_si128(_mm_blendv_ps(_mm_castsi128_ps(elsePacket), _mm_castsi128_ps(thenPacket), _mm_castsi128_ps(ifPacket)));
792 #else
793  return _mm_or_si128(_mm_and_si128(ifPacket, thenPacket), _mm_andnot_si128(ifPacket, elsePacket));
794 #endif
795  }
796 
797  template<>
798  EIGEN_STRONG_INLINE Packet2d pblendv(const Packet2d& ifPacket, const Packet2d& thenPacket, const Packet2d& elsePacket)
799  {
800 #ifdef EIGEN_VECTORIZE_SSE4_1
801  return _mm_blendv_pd(elsePacket, thenPacket, ifPacket);
802 #else
803  return _mm_or_pd(_mm_and_pd(ifPacket, thenPacket), _mm_andnot_pd(ifPacket, elsePacket));
804 #endif
805  }
806 
807  template<>
808  EIGEN_STRONG_INLINE Packet4i pgather<Packet4i>(const int* addr, const Packet4i& index)
809  {
810 #ifdef EIGEN_VECTORIZE_AVX2
811  return _mm_i32gather_epi32(addr, index, 4);
812 #else
813  uint32_t u[4];
814  _mm_storeu_si128((Packet4i*)u, index);
815  return _mm_setr_epi32(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]]);
816 #endif
817  }
818 
819  template<>
820  EIGEN_STRONG_INLINE Packet4f pgather<Packet4i>(const float* addr, const Packet4i& index)
821  {
822 #ifdef EIGEN_VECTORIZE_AVX2
823  return _mm_i32gather_ps(addr, index, 4);
824 #else
825  uint32_t u[4];
826  _mm_storeu_si128((Packet4i*)u, index);
827  return _mm_setr_ps(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]]);
828 #endif
829  }
830 
831  template<>
832  EIGEN_STRONG_INLINE Packet2d pgather<Packet4i>(const double* addr, const Packet4i& index, bool upperhalf)
833  {
834 #ifdef EIGEN_VECTORIZE_AVX2
835  return _mm_i32gather_pd(addr, index, 8);
836 #else
837  uint32_t u[4];
838  _mm_storeu_si128((Packet4i*)u, index);
839  if (upperhalf)
840  {
841  return _mm_setr_pd(addr[u[2]], addr[u[3]]);
842  }
843  else
844  {
845  return _mm_setr_pd(addr[u[0]], addr[u[1]]);
846  }
847 #endif
848  }
849 
850  template<>
851  EIGEN_STRONG_INLINE int pmovemask<Packet4f>(const Packet4f& a)
852  {
853  return _mm_movemask_ps(a);
854  }
855 
856  template<>
857  EIGEN_STRONG_INLINE int pmovemask<Packet2d>(const Packet2d& a)
858  {
859  return _mm_movemask_pd(a);
860  }
861 
862  template<>
863  EIGEN_STRONG_INLINE int pmovemask<Packet4i>(const Packet4i& a)
864  {
865  return pmovemask(_mm_castsi128_ps(a));
866  }
867 
868  template<>
869  EIGEN_STRONG_INLINE Packet4f ptruncate<Packet4f>(const Packet4f& a)
870  {
871 #ifdef EIGEN_VECTORIZE_SSE4_1
872  return _mm_round_ps(a, _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC);
873 #else
874  auto round = _MM_GET_ROUNDING_MODE();
875  _MM_SET_ROUNDING_MODE(_MM_ROUND_TOWARD_ZERO);
876  auto ret = _mm_cvtepi32_ps(_mm_cvtps_epi32(a));
877  _MM_SET_ROUNDING_MODE(round);
878  return ret;
879 #endif
880  }
881 
882  template<>
883  EIGEN_STRONG_INLINE Packet2d ptruncate<Packet2d>(const Packet2d& a)
884  {
885 #ifdef EIGEN_VECTORIZE_SSE4_1
886  return _mm_round_pd(a, _MM_FROUND_TO_ZERO | _MM_FROUND_NO_EXC);
887 #else
888  auto round = _MM_GET_ROUNDING_MODE();
889  _MM_SET_ROUNDING_MODE(_MM_ROUND_TOWARD_ZERO);
890  auto ret = _mm_cvtepi32_pd(_mm_cvtpd_epi32(a));
891  _MM_SET_ROUNDING_MODE(round);
892  return ret;
893 #endif
894  }
895  }
896 }
897 #endif
898 
899 #endif