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 void psincos(Packet x, Packet &s, Packet &c)
67  {
68  Packet xmm1, xmm2, xmm3 = pset1<Packet>(0), sign_bit_sin, y;
69  using IntPacket = decltype(reinterpret_to_int(x));
70  IntPacket emm0, emm2, emm4;
71 
72  sign_bit_sin = x;
73  /* take the absolute value */
74  x = pabs(x);
75  /* extract the sign bit (upper one) */
76  sign_bit_sin = reinterpret_to_float(
77  pand(reinterpret_to_int(sign_bit_sin), pset1<IntPacket>(0x80000000))
78  );
79 
80  /* scale by 4/Pi */
81  y = pmul(x, pset1<Packet>(1.27323954473516));
82 
83  /* store the integer part of y in emm2 */
84  emm2 = pcast<Packet, IntPacket>(y);
85 
86  /* j=(j+1) & (~1) (see the cephes sources) */
87  emm2 = padd(emm2, pset1<IntPacket>(1));
88  emm2 = pand(emm2, pset1<IntPacket>(~1));
89  y = pcast<IntPacket, Packet>(emm2);
90 
91  emm4 = emm2;
92 
93  /* get the swap sign flag for the sine */
94  emm0 = pand(emm2, pset1<IntPacket>(4));
95  emm0 = psll(emm0, 29);
96  Packet swap_sign_bit_sin = reinterpret_to_float(emm0);
97 
98  /* get the polynom selection mask for the sine*/
99  emm2 = pand(emm2, pset1<IntPacket>(2));
100 
101  emm2 = pcmpeq(emm2, pset1<IntPacket>(0));
102  Packet poly_mask = reinterpret_to_float(emm2);
103 
104  /* The magic pass: "Extended precision modular arithmetic"
105  x = ((x - y * DP1) - y * DP2) - y * DP3; */
106  xmm1 = pset1<Packet>(-0.78515625);
107  xmm2 = pset1<Packet>(-2.4187564849853515625e-4);
108  xmm3 = pset1<Packet>(-3.77489497744594108e-8);
109  xmm1 = pmul(y, xmm1);
110  xmm2 = pmul(y, xmm2);
111  xmm3 = pmul(y, xmm3);
112  x = padd(x, xmm1);
113  x = padd(x, xmm2);
114  x = padd(x, xmm3);
115 
116  emm4 = psub(emm4, pset1<IntPacket>(2));
117  emm4 = pandnot(emm4, pset1<IntPacket>(4));
118  emm4 = psll(emm4, 29);
119  Packet sign_bit_cos = reinterpret_to_float(emm4);
120  sign_bit_sin = pxor(sign_bit_sin, swap_sign_bit_sin);
121 
122 
123  /* Evaluate the first polynom (0 <= x <= Pi/4) */
124  Packet z = pmul(x, x);
125  y = pset1<Packet>(2.443315711809948E-005);
126 
127  y = pmul(y, z);
128  y = padd(y, pset1<Packet>(-1.388731625493765E-003));
129  y = pmul(y, z);
130  y = padd(y, pset1<Packet>(4.166664568298827E-002));
131  y = pmul(y, z);
132  y = pmul(y, z);
133  Packet tmp = pmul(z, pset1<Packet>(0.5));
134  y = psub(y, tmp);
135  y = padd(y, pset1<Packet>(1));
136 
137  /* Evaluate the second polynom (Pi/4 <= x <= 0) */
138 
139  Packet y2 = pset1<Packet>(-1.9515295891E-4);
140  y2 = pmul(y2, z);
141  y2 = padd(y2, pset1<Packet>(8.3321608736E-3));
142  y2 = pmul(y2, z);
143  y2 = padd(y2, pset1<Packet>(-1.6666654611E-1));
144  y2 = pmul(y2, z);
145  y2 = pmul(y2, x);
146  y2 = padd(y2, x);
147 
148  /* select the correct result from the two polynoms */
149  xmm3 = poly_mask;
150  Packet ysin2 = pand(xmm3, y2);
151  Packet ysin1 = pandnot(xmm3, y);
152  y2 = psub(y2, ysin2);
153  y = psub(y, ysin1);
154 
155  xmm1 = padd(ysin1, ysin2);
156  xmm2 = padd(y, y2);
157 
158  /* update the sign */
159  s = pxor(xmm1, sign_bit_sin);
160  c = pxor(xmm2, sign_bit_cos);
161  }
162 
163  template<typename Packet>
164  EIGEN_STRONG_INLINE Packet pcmplt(const Packet& a, const Packet& b);
165 
166  template<typename Packet>
167  EIGEN_STRONG_INLINE Packet pcmple(const Packet& a, const Packet& b);
168 
169  template<typename Packet>
170  EIGEN_STRONG_INLINE Packet pblendv(const Packet& ifPacket, const Packet& thenPacket, const Packet& elsePacket);
171 
172  template<typename Packet>
173  EIGEN_STRONG_INLINE Packet pgather(const int* addr, const Packet& index);
174 
175  template<typename Packet>
176  EIGEN_STRONG_INLINE auto pgather(const float* addr, const Packet& index) -> decltype(reinterpret_to_float(std::declval<Packet>()));
177 
178  template<typename Packet>
179  EIGEN_STRONG_INLINE auto pgather(const double* addr, const Packet& index, bool upperhalf = false) -> decltype(reinterpret_to_double(std::declval<Packet>()));
180  }
181 }
182 
183 #ifdef EIGEN_VECTORIZE_AVX
184 #include <immintrin.h>
185 
186 namespace Eigen
187 {
188  namespace internal
189  {
190  template<>
191  struct reinterpreter<Packet8i>
192  {
193  EIGEN_STRONG_INLINE Packet8f to_float(const Packet8i& x)
194  {
195  return _mm256_castsi256_ps(x);
196  }
197 
198  EIGEN_STRONG_INLINE Packet4d to_double(const Packet8i& x)
199  {
200  return _mm256_castsi256_pd(x);
201  }
202 
203  EIGEN_STRONG_INLINE Packet8i to_int(const Packet8i& x)
204  {
205  return x;
206  }
207  };
208 
209  template<>
210  struct reinterpreter<Packet8f>
211  {
212  EIGEN_STRONG_INLINE Packet8f to_float(const Packet8f& x)
213  {
214  return x;
215  }
216 
217  EIGEN_STRONG_INLINE Packet4d to_double(const Packet8f& x)
218  {
219  return _mm256_castps_pd(x);
220  }
221 
222  EIGEN_STRONG_INLINE Packet8i to_int(const Packet8f& x)
223  {
224  return _mm256_castps_si256(x);
225  }
226  };
227 
228  template<>
229  struct reinterpreter<Packet4d>
230  {
231  EIGEN_STRONG_INLINE Packet8f to_float(const Packet4d& x)
232  {
233  return _mm256_castpd_ps(x);
234  }
235 
236  EIGEN_STRONG_INLINE Packet4d to_double(const Packet4d& x)
237  {
238  return x;
239  }
240 
241  EIGEN_STRONG_INLINE Packet8i to_int(const Packet4d& x)
242  {
243  return _mm256_castpd_si256(x);
244  }
245  };
246 
247  EIGEN_STRONG_INLINE void split_two(const Packet8i& x, Packet4i& a, Packet4i& b)
248  {
249  a = _mm256_extractf128_si256(x, 0);
250  b = _mm256_extractf128_si256(x, 1);
251  }
252 
253  EIGEN_STRONG_INLINE Packet8i combine_two(const Packet4i& a, const Packet4i& b)
254  {
255  return _mm256_insertf128_si256(_mm256_castsi128_si256(a), b, 1);
256  }
257 
258  EIGEN_STRONG_INLINE void split_two(const Packet8f& x, Packet4f& a, Packet4f& b)
259  {
260  a = _mm256_extractf128_ps(x, 0);
261  b = _mm256_extractf128_ps(x, 1);
262  }
263 
264  EIGEN_STRONG_INLINE Packet8f combine_two(const Packet4f& a, const Packet4f& b)
265  {
266  return _mm256_insertf128_ps(_mm256_castps128_ps256(a), b, 1);
267  }
268 
269 
270  EIGEN_STRONG_INLINE Packet4i combine_low32(const Packet8i& a)
271  {
272 #ifdef EIGEN_VECTORIZE_AVX2
273  return _mm256_castsi256_si128(_mm256_permutevar8x32_epi32(a, _mm256_setr_epi32(0, 2, 4, 6, 1, 3, 5, 7)));
274 #else
275  auto sc = _mm256_permutevar_ps(_mm256_castsi256_ps(a), _mm256_setr_epi32(0, 2, 1, 3, 1, 3, 0, 2));
276  return _mm_castps_si128(_mm_blend_ps(_mm256_extractf128_ps(sc, 0), _mm256_extractf128_ps(sc, 1), 0b1100));
277 #endif
278  }
279 
280  template<>
281  EIGEN_STRONG_INLINE Packet8i pseti64<Packet8i>(uint64_t a)
282  {
283  return _mm256_set1_epi64x(a);
284  }
285 
286  template<>
287  EIGEN_STRONG_INLINE Packet8i pcmpeq<Packet8i>(const Packet8i& a, const Packet8i& b)
288  {
289 #ifdef EIGEN_VECTORIZE_AVX2
290  return _mm256_cmpeq_epi32(a, b);
291 #else
292  Packet4i a1, a2, b1, b2;
293  split_two(a, a1, a2);
294  split_two(b, b1, b2);
295  return combine_two(_mm_cmpeq_epi32(a1, b1), _mm_cmpeq_epi32(a2, b2));
296 #endif
297  }
298 
299  template<>
300  EIGEN_STRONG_INLINE Packet8i psll<Packet8i>(const Packet8i& a, int b)
301  {
302 #ifdef EIGEN_VECTORIZE_AVX2
303  return _mm256_slli_epi32(a, b);
304 #else
305  Packet4i a1, a2;
306  split_two(a, a1, a2);
307  return combine_two(_mm_slli_epi32(a1, b), _mm_slli_epi32(a2, b));
308 #endif
309  }
310 
311  template<>
312  EIGEN_STRONG_INLINE Packet8i psrl<Packet8i>(const Packet8i& a, int b)
313  {
314 #ifdef EIGEN_VECTORIZE_AVX2
315  return _mm256_srli_epi32(a, b);
316 #else
317  Packet4i a1, a2;
318  split_two(a, a1, a2);
319  return combine_two(_mm_srli_epi32(a1, b), _mm_srli_epi32(a2, b));
320 #endif
321  }
322 
323  template<>
324  EIGEN_STRONG_INLINE Packet8i psll64<Packet8i>(const Packet8i& a, int b)
325  {
326 #ifdef EIGEN_VECTORIZE_AVX2
327  return _mm256_slli_epi64(a, b);
328 #else
329  Packet4i a1, a2;
330  split_two(a, a1, a2);
331  return combine_two(_mm_slli_epi64(a1, b), _mm_slli_epi64(a2, b));
332 #endif
333  }
334 
335  template<>
336  EIGEN_STRONG_INLINE Packet8i psrl64<Packet8i>(const Packet8i& a, int b)
337  {
338 #ifdef EIGEN_VECTORIZE_AVX2
339  return _mm256_srli_epi64(a, b);
340 #else
341  Packet4i a1, a2;
342  split_two(a, a1, a2);
343  return combine_two(_mm_srli_epi64(a1, b), _mm_srli_epi64(a2, b));
344 #endif
345  }
346 
347  template<> EIGEN_STRONG_INLINE Packet8i padd<Packet8i>(const Packet8i& a, const Packet8i& b)
348  {
349 #ifdef EIGEN_VECTORIZE_AVX2
350  return _mm256_add_epi32(a, b);
351 #else
352  Packet4i a1, a2, b1, b2;
353  split_two(a, a1, a2);
354  split_two(b, b1, b2);
355  return combine_two(_mm_add_epi32(a1, b1), _mm_add_epi32(a2, b2));
356 #endif
357  }
358 
359  template<> EIGEN_STRONG_INLINE Packet8i psub<Packet8i>(const Packet8i& a, const Packet8i& b)
360  {
361 #ifdef EIGEN_VECTORIZE_AVX2
362  return _mm256_sub_epi32(a, b);
363 #else
364  Packet4i a1, a2, b1, b2;
365  split_two(a, a1, a2);
366  split_two(b, b1, b2);
367  return combine_two(_mm_sub_epi32(a1, b1), _mm_sub_epi32(a2, b2));
368 #endif
369  }
370 
371  template<> EIGEN_STRONG_INLINE Packet8i pand<Packet8i>(const Packet8i& a, const Packet8i& b)
372  {
373 #ifdef EIGEN_VECTORIZE_AVX2
374  return _mm256_and_si256(a, b);
375 #else
376  return reinterpret_to_int(_mm256_and_ps(reinterpret_to_float(a), reinterpret_to_float(b)));
377 #endif
378  }
379 
380  template<> EIGEN_STRONG_INLINE Packet8i pandnot<Packet8i>(const Packet8i& a, const Packet8i& b)
381  {
382 #ifdef EIGEN_VECTORIZE_AVX2
383  return _mm256_andnot_si256(a, b);
384 #else
385  return reinterpret_to_int(_mm256_andnot_ps(reinterpret_to_float(a), reinterpret_to_float(b)));
386 #endif
387  }
388 
389  template<> EIGEN_STRONG_INLINE Packet8i por<Packet8i>(const Packet8i& a, const Packet8i& b)
390  {
391 #ifdef EIGEN_VECTORIZE_AVX2
392  return _mm256_or_si256(a, b);
393 #else
394  return reinterpret_to_int(_mm256_or_ps(reinterpret_to_float(a), reinterpret_to_float(b)));
395 #endif
396  }
397 
398  template<> EIGEN_STRONG_INLINE Packet8i pxor<Packet8i>(const Packet8i& a, const Packet8i& b)
399  {
400 #ifdef EIGEN_VECTORIZE_AVX2
401  return _mm256_xor_si256(a, b);
402 #else
403  return reinterpret_to_int(_mm256_xor_ps(reinterpret_to_float(a), reinterpret_to_float(b)));
404 #endif
405  }
406 
407  template<>
408  EIGEN_STRONG_INLINE Packet8i pcmplt<Packet8i>(const Packet8i& a, const Packet8i& b)
409  {
410  return _mm256_cmpgt_epi32(b, a);
411  }
412 
413  template<>
414  EIGEN_STRONG_INLINE Packet8f pcmplt<Packet8f>(const Packet8f& a, const Packet8f& b)
415  {
416  return _mm256_cmp_ps(a, b, _CMP_LT_OQ);
417  }
418 
419  template<>
420  EIGEN_STRONG_INLINE Packet8f pcmple<Packet8f>(const Packet8f& a, const Packet8f& b)
421  {
422  return _mm256_cmp_ps(a, b, _CMP_LE_OQ);
423  }
424 
425  template<>
426  EIGEN_STRONG_INLINE Packet4d pcmplt<Packet4d>(const Packet4d& a, const Packet4d& b)
427  {
428  return _mm256_cmp_pd(a, b, _CMP_LT_OQ);
429  }
430 
431  template<>
432  EIGEN_STRONG_INLINE Packet4d pcmple<Packet4d>(const Packet4d& a, const Packet4d& b)
433  {
434  return _mm256_cmp_pd(a, b, _CMP_LE_OQ);
435  }
436 
437  template<>
438  EIGEN_STRONG_INLINE Packet8f pblendv(const Packet8f& ifPacket, const Packet8f& thenPacket, const Packet8f& elsePacket)
439  {
440  return _mm256_blendv_ps(elsePacket, thenPacket, ifPacket);
441  }
442 
443  template<>
444  EIGEN_STRONG_INLINE Packet8i pblendv(const Packet8i& ifPacket, const Packet8i& thenPacket, const Packet8i& elsePacket)
445  {
446  return _mm256_castps_si256(_mm256_blendv_ps(
447  _mm256_castsi256_ps(elsePacket),
448  _mm256_castsi256_ps(thenPacket),
449  _mm256_castsi256_ps(ifPacket)
450  ));
451  }
452 
453  template<>
454  EIGEN_STRONG_INLINE Packet4d pblendv(const Packet4d& ifPacket, const Packet4d& thenPacket, const Packet4d& elsePacket)
455  {
456  return _mm256_blendv_pd(elsePacket, thenPacket, ifPacket);
457  }
458 
459  template<>
460  EIGEN_STRONG_INLINE Packet8i pgather<Packet8i>(const int* addr, const Packet8i& index)
461  {
462 #ifdef EIGEN_VECTORIZE_AVX2
463  return _mm256_i32gather_epi32(addr, index, 4);
464 #else
465  uint32_t u[8];
466  _mm256_storeu_si256((Packet8i*)u, index);
467  return _mm256_setr_epi32(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]],
468  addr[u[4]], addr[u[5]], addr[u[6]], addr[u[7]]);
469 #endif
470  }
471 
472  template<>
473  EIGEN_STRONG_INLINE Packet8f pgather<Packet8i>(const float *addr, const Packet8i& index)
474  {
475 #ifdef EIGEN_VECTORIZE_AVX2
476  return _mm256_i32gather_ps(addr, index, 4);
477 #else
478  uint32_t u[8];
479  _mm256_storeu_si256((Packet8i*)u, index);
480  return _mm256_setr_ps(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]],
481  addr[u[4]], addr[u[5]], addr[u[6]], addr[u[7]]);
482 #endif
483  }
484 
485  template<>
486  EIGEN_STRONG_INLINE Packet4d pgather<Packet8i>(const double *addr, const Packet8i& index, bool upperhalf)
487  {
488 #ifdef EIGEN_VECTORIZE_AVX2
489  return _mm256_i32gather_pd(addr, _mm256_castsi256_si128(index), 8);
490 #else
491  uint32_t u[8];
492  _mm256_storeu_si256((Packet8i*)u, index);
493  if (upperhalf)
494  {
495  return _mm256_setr_pd(addr[u[4]], addr[u[5]], addr[u[6]], addr[u[7]]);
496  }
497  else
498  {
499  return _mm256_setr_pd(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]]);
500  }
501 #endif
502  }
503  }
504 }
505 #endif
506 
507 #ifdef EIGEN_VECTORIZE_SSE2
508 #include <xmmintrin.h>
509 
510 namespace Eigen
511 {
512  namespace internal
513  {
514  template<>
515  struct reinterpreter<Packet4i>
516  {
517  EIGEN_STRONG_INLINE Packet4f to_float(const Packet4i& x)
518  {
519  return _mm_castsi128_ps(x);
520  }
521 
522  EIGEN_STRONG_INLINE Packet2d to_double(const Packet4i& x)
523  {
524  return _mm_castsi128_pd(x);
525  }
526 
527  EIGEN_STRONG_INLINE Packet4i to_int(const Packet4i& x)
528  {
529  return x;
530  }
531  };
532 
533  template<>
534  struct reinterpreter<Packet4f>
535  {
536  EIGEN_STRONG_INLINE Packet4f to_float(const Packet4f& x)
537  {
538  return x;
539  }
540 
541  EIGEN_STRONG_INLINE Packet2d to_double(const Packet4f& x)
542  {
543  return _mm_castps_pd(x);
544  }
545 
546  EIGEN_STRONG_INLINE Packet4i to_int(const Packet4f& x)
547  {
548  return _mm_castps_si128(x);
549  }
550  };
551 
552  template<>
553  struct reinterpreter<Packet2d>
554  {
555  EIGEN_STRONG_INLINE Packet4f to_float(const Packet2d& x)
556  {
557  return _mm_castpd_ps(x);
558  }
559 
560  EIGEN_STRONG_INLINE Packet2d to_double(const Packet2d& x)
561  {
562  return x;
563  }
564 
565  EIGEN_STRONG_INLINE Packet4i to_int(const Packet2d& x)
566  {
567  return _mm_castpd_si128(x);
568  }
569  };
570 
571  EIGEN_STRONG_INLINE void split_two(const Packet4i& x, uint64_t& a, uint64_t& b)
572  {
573 #ifdef EIGEN_VECTORIZE_SSE4_1
574  a = _mm_extract_epi64(x, 0);
575  b = _mm_extract_epi64(x, 1);
576 #else
577  uint64_t u[2];
578  _mm_storeu_si128((__m128i*)u, x);
579  a = u[0];
580  b = u[1];
581 #endif
582  }
583 
584  EIGEN_STRONG_INLINE Packet4i combine_low32(const Packet4i& a, const Packet4i& b)
585  {
586  auto sa = _mm_shuffle_epi32(a, _MM_SHUFFLE(3, 1, 2, 0));
587  auto sb = _mm_shuffle_epi32(b, _MM_SHUFFLE(2, 0, 3, 1));
588  sa = _mm_and_si128(sa, _mm_setr_epi32(-1, -1, 0, 0));
589  sb = _mm_and_si128(sb, _mm_setr_epi32(0, 0, -1, -1));
590  return _mm_or_si128(sa, sb);
591  }
592 
593  template<>
594  EIGEN_STRONG_INLINE Packet4i pseti64<Packet4i>(uint64_t a)
595  {
596  return _mm_set1_epi64x(a);
597  }
598 
599  template<>
600  Packet4i pcmpeq<Packet4i>(const Packet4i& a, const Packet4i& b)
601  {
602  return _mm_cmpeq_epi32(a, b);
603  }
604 
605  template<>
606  Packet4i psll<Packet4i>(const Packet4i& a, int b)
607  {
608  return _mm_slli_epi32(a, b);
609  }
610 
611  template<>
612  Packet4i psrl<Packet4i>(const Packet4i& a, int b)
613  {
614  return _mm_srli_epi32(a, b);
615  }
616 
617 
618  template<>
619  Packet4i psll64<Packet4i>(const Packet4i& a, int b)
620  {
621  return _mm_slli_epi64(a, b);
622  }
623 
624  template<>
625  Packet4i psrl64<Packet4i>(const Packet4i& a, int b)
626  {
627  return _mm_srli_epi64(a, b);
628  }
629 
630  template<>
631  EIGEN_STRONG_INLINE Packet4i pcmplt<Packet4i>(const Packet4i& a, const Packet4i& b)
632  {
633  return _mm_cmplt_epi32(a, b);
634  }
635 
636  template<>
637  EIGEN_STRONG_INLINE Packet4f pcmplt<Packet4f>(const Packet4f& a, const Packet4f& b)
638  {
639  return _mm_cmplt_ps(a, b);
640  }
641 
642  template<>
643  EIGEN_STRONG_INLINE Packet4f pcmple<Packet4f>(const Packet4f& a, const Packet4f& b)
644  {
645  return _mm_cmple_ps(a, b);
646  }
647 
648  template<>
649  EIGEN_STRONG_INLINE Packet2d pcmplt<Packet2d>(const Packet2d& a, const Packet2d& b)
650  {
651  return _mm_cmplt_pd(a, b);
652  }
653 
654  template<>
655  EIGEN_STRONG_INLINE Packet2d pcmple<Packet2d>(const Packet2d& a, const Packet2d& b)
656  {
657  return _mm_cmple_pd(a, b);
658  }
659 
660  template<>
661  EIGEN_STRONG_INLINE Packet4f pblendv(const Packet4f& ifPacket, const Packet4f& thenPacket, const Packet4f& elsePacket)
662  {
663 #ifdef XXX_EIGEN_VECTORIZE_SSE4_1
664  return _mm_blendv_ps(elsePacket, thenPacket, ifPacket);
665 #else
666  return _mm_or_ps(_mm_and_ps(ifPacket, thenPacket), _mm_andnot_ps(ifPacket, elsePacket));
667 #endif
668  }
669 
670  template<>
671  EIGEN_STRONG_INLINE Packet4i pblendv(const Packet4i& ifPacket, const Packet4i& thenPacket, const Packet4i& elsePacket)
672  {
673 #ifdef XXX_EIGEN_VECTORIZE_SSE4_1
674  return _mm_castps_si128(_mm_blendv_ps(_mm_castsi128_ps(elsePacket), _mm_castsi128_ps(thenPacket), _mm_castsi128_ps(ifPacket)));
675 #else
676  return _mm_or_si128(_mm_and_si128(ifPacket, thenPacket), _mm_andnot_si128(ifPacket, elsePacket));
677 #endif
678  }
679 
680  template<>
681  EIGEN_STRONG_INLINE Packet2d pblendv(const Packet2d& ifPacket, const Packet2d& thenPacket, const Packet2d& elsePacket)
682  {
683 #ifdef XXX_EIGEN_VECTORIZE_SSE4_1
684  return _mm_blendv_pd(elsePacket, thenPacket, ifPacket);
685 #else
686  return _mm_or_pd(_mm_and_pd(ifPacket, thenPacket), _mm_andnot_pd(ifPacket, elsePacket));
687 #endif
688  }
689 
690  template<>
691  EIGEN_STRONG_INLINE Packet4i pgather<Packet4i>(const int* addr, const Packet4i& index)
692  {
693 #ifdef EIGEN_VECTORIZE_AVX2
694  return _mm_i32gather_epi32(addr, index, 4);
695 #else
696  uint32_t u[4];
697  _mm_storeu_si128((Packet4i*)u, index);
698  return _mm_setr_epi32(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]]);
699 #endif
700  }
701 
702  template<>
703  EIGEN_STRONG_INLINE Packet4f pgather<Packet4i>(const float* addr, const Packet4i& index)
704  {
705 #ifdef EIGEN_VECTORIZE_AVX2
706  return _mm_i32gather_ps(addr, index, 4);
707 #else
708  uint32_t u[4];
709  _mm_storeu_si128((Packet4i*)u, index);
710  return _mm_setr_ps(addr[u[0]], addr[u[1]], addr[u[2]], addr[u[3]]);
711 #endif
712  }
713 
714  template<>
715  EIGEN_STRONG_INLINE Packet2d pgather<Packet4i>(const double* addr, const Packet4i& index, bool upperhalf)
716  {
717 #ifdef EIGEN_VECTORIZE_AVX2
718  return _mm_i32gather_pd(addr, index, 8);
719 #else
720  uint32_t u[4];
721  _mm_storeu_si128((Packet4i*)u, index);
722  if (upperhalf)
723  {
724  return _mm_setr_pd(addr[u[2]], addr[u[3]]);
725  }
726  else
727  {
728  return _mm_setr_pd(addr[u[0]], addr[u[1]]);
729  }
730 #endif
731  }
732  }
733 }
734 #endif
735 
736 #endif