EigenRand  0.4.0-alpha
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 #define EIGENRAND_PRINT_PACKET(p) do { using _MTy = typename std::remove_const<typename std::remove_reference<decltype(p)>::type>::type; typename std::conditional<Eigen::internal::IsFloatPacket<_MTy>::value, float, typename std::conditional<Eigen::internal::IsDoublePacket<_MTy>::value, double, int>::type>::type f[4]; Eigen::internal::pstore(f, p); std::cout << #p " " << f[0] << " " << f[1] << " " << f[2] << " " << f[3] << std::endl; } while(0)
18 
19 namespace Eigen
20 {
21  namespace internal
22  {
23  template<typename Ty>
24  struct IsIntPacket : std::false_type {};
25 
26  template<typename Ty>
27  struct IsFloatPacket : std::false_type {};
28 
29  template<typename Ty>
30  struct IsDoublePacket : std::false_type {};
31 
32  template<typename Ty>
33  struct HalfPacket;
34 
35  template<typename Packet>
36  struct reinterpreter{};
37 
38  template<typename Packet>
39  inline auto reinterpret_to_float(const Packet& x)
40  -> decltype(reinterpreter<Packet>{}.to_float(x))
41  {
42  return reinterpreter<Packet>{}.to_float(x);
43  }
44 
45  template<typename Packet>
46  inline auto reinterpret_to_double(const Packet& x)
47  -> decltype(reinterpreter<Packet>{}.to_double(x))
48  {
49  return reinterpreter<Packet>{}.to_double(x);
50  }
51 
52  template<typename Packet>
53  inline auto reinterpret_to_int(const Packet& x)
54  -> decltype(reinterpreter<Packet>{}.to_int(x))
55  {
56  return reinterpreter<Packet>{}.to_int(x);
57  }
58 
59  template<typename Packet>
60  EIGEN_STRONG_INLINE void split_two(const Packet& p, typename HalfPacket<Packet>::type& a, typename HalfPacket<Packet>::type& b);
61 
62  template<typename Packet>
63  EIGEN_STRONG_INLINE Packet pseti64(uint64_t a);
64 
65  template<typename Packet>
66  EIGEN_STRONG_INLINE Packet padd64(const Packet& a, const Packet& b);
67 
68  template<typename Packet>
69  EIGEN_STRONG_INLINE Packet psub64(const Packet& a, const Packet& b);
70 
71  template <typename SrcPacket, typename TgtPacket>
72  EIGEN_STRONG_INLINE TgtPacket pcast64(const SrcPacket& a);
73 
74  template<typename Packet>
75  EIGEN_STRONG_INLINE Packet pcmpeq(const Packet& a, const Packet& b);
76 
77  template<typename Packet>
78  struct BitShifter {};
79 
80  template<int b, typename Packet>
81  EIGEN_STRONG_INLINE Packet psll(const Packet& a);
82 
83  template<int _b, typename Packet>
84  EIGEN_STRONG_INLINE Packet psrl(const Packet& a, int b = _b);
85 
86  template<int b, typename Packet>
87  EIGEN_STRONG_INLINE Packet psll64(const Packet& a);
88 
89  template<int b, typename Packet>
90  EIGEN_STRONG_INLINE Packet psrl64(const Packet& a);
91 
92  /*template<typename Packet>
93  EIGEN_STRONG_INLINE Packet psll(const Packet& a, int b);
94 
95  template<typename Packet>
96  EIGEN_STRONG_INLINE Packet psrl(const Packet& a, int b);
97 
98  template<typename Packet>
99  EIGEN_STRONG_INLINE Packet psll64(const Packet& a, int b);
100 
101  template<typename Packet>
102  EIGEN_STRONG_INLINE Packet psrl64(const Packet& a, int b);*/
103 
104  template<typename Packet>
105  EIGEN_STRONG_INLINE int pmovemask(const Packet& a);
106 
107  template<typename Packet>
108  EIGEN_STRONG_INLINE typename std::enable_if<
109  IsFloatPacket<Packet>::value, Packet
110  >::type pext_sign(const Packet& a)
111  {
112  using IntPacket = decltype(reinterpret_to_int(a));
113  return reinterpret_to_float(
114  pand(reinterpret_to_int(a), pset1<IntPacket>(0x80000000))
115  );
116  }
117 
118  template<typename Packet>
119  EIGEN_STRONG_INLINE typename std::enable_if<
120  IsDoublePacket<Packet>::value, Packet
121  >::type pext_sign(const Packet& a)
122  {
123  using IntPacket = decltype(reinterpret_to_int(a));
124  return reinterpret_to_double(
125  pand(reinterpret_to_int(a), pseti64<IntPacket>(0x8000000000000000))
126  );
127  }
128 
129  /*template<>
130  EIGEN_STRONG_INLINE uint64_t psll64<uint64_t>(const uint64_t& a, int b)
131  {
132  return a << b;
133  }
134 
135  template<>
136  EIGEN_STRONG_INLINE uint64_t psrl64<uint64_t>(const uint64_t& a, int b)
137  {
138  return a >> b;
139  }*/
140 
141  // 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))
142  template<typename Packet>
143  EIGEN_STRONG_INLINE Packet plgamma_approx(const Packet& x)
144  {
145  auto x_3 = padd(x, pset1<Packet>(3));
146  auto ret = pmul(padd(x_3, pset1<Packet>(-0.5)), plog(x_3));
147  ret = psub(ret, x_3);
148  ret = padd(ret, pset1<Packet>(0.9189385332046727));
149  ret = padd(ret, pdiv(pset1<Packet>(1 / 12.), x_3));
150  ret = psub(ret, plog(pmul(
151  pmul(psub(x_3, pset1<Packet>(1)), psub(x_3, pset1<Packet>(2))), x)));
152  return ret;
153  }
154 
155  template<typename Packet>
156  EIGEN_STRONG_INLINE Packet pcmplt(const Packet& a, const Packet& b);
157 
158  template<typename Packet>
159  EIGEN_STRONG_INLINE Packet pcmple(const Packet& a, const Packet& b);
160 
161  template<typename Packet>
162  EIGEN_STRONG_INLINE Packet pbitnot(const Packet& a);
163 
164  template<typename PacketIf, typename Packet>
165  EIGEN_STRONG_INLINE Packet pblendv(const PacketIf& ifPacket, const Packet& thenPacket, const Packet& elsePacket);
166 
167  template<typename Packet>
168  EIGEN_STRONG_INLINE Packet pgather(const int* addr, const Packet& index);
169 
170  template<typename Packet>
171  EIGEN_STRONG_INLINE auto pgather(const float* addr, const Packet& index) -> decltype(reinterpret_to_float(std::declval<Packet>()));
172 
173  template<typename Packet>
174  EIGEN_STRONG_INLINE auto pgather(const double* addr, const Packet& index, bool upperhalf = false) -> decltype(reinterpret_to_double(std::declval<Packet>()));
175 
176  template<typename Packet>
177  EIGEN_STRONG_INLINE Packet ptruncate(const Packet& a);
178 
179  template<typename Packet>
180  EIGEN_STRONG_INLINE Packet pcmpeq64(const Packet& a, const Packet& b);
181 
182  template<typename Packet>
183  EIGEN_STRONG_INLINE Packet pcmplt64(const Packet& a, const Packet& b);
184 
185  template<typename Packet>
186  EIGEN_STRONG_INLINE Packet pmuluadd64(const Packet& a, uint64_t b, uint64_t c);
187 
188  template<typename IntPacket>
189  EIGEN_STRONG_INLINE auto bit_to_ur_float(const IntPacket& x) -> decltype(reinterpret_to_float(x))
190  {
191  using FloatPacket = decltype(reinterpret_to_float(x));
192 
193  const IntPacket lower = pset1<IntPacket>(0x7FFFFF),
194  upper = pset1<IntPacket>(127 << 23);
195  const FloatPacket one = pset1<FloatPacket>(1);
196 
197  return psub(reinterpret_to_float(por(pand(x, lower), upper)), one);
198  }
199 
200  template<typename IntPacket>
201  EIGEN_STRONG_INLINE auto bit_to_ur_double(const IntPacket& x) -> decltype(reinterpret_to_double(x))
202  {
203  using DoublePacket = decltype(reinterpret_to_double(x));
204 
205  const IntPacket lower = pseti64<IntPacket>(0xFFFFFFFFFFFFFull),
206  upper = pseti64<IntPacket>(1023ull << 52);
207  const DoublePacket one = pset1<DoublePacket>(1);
208 
209  return psub(reinterpret_to_double(por(pand(x, lower), upper)), one);
210  }
211 
212  template<typename _Scalar>
213  struct BitScalar;
214 
215  template<>
216  struct BitScalar<float>
217  {
218  float to_ur(uint32_t x)
219  {
220  union
221  {
222  uint32_t u;
223  float f;
224  };
225  u = (x & 0x7FFFFF) | (127 << 23);
226  return f - 1.f;
227  }
228 
229  float to_nzur(uint32_t x)
230  {
231  return to_ur(x) + std::numeric_limits<float>::epsilon() / 8;
232  }
233  };
234 
235  template<>
236  struct BitScalar<double>
237  {
238  double to_ur(uint64_t x)
239  {
240  union
241  {
242  uint64_t u;
243  double f;
244  };
245  u = (x & 0xFFFFFFFFFFFFFull) | (1023ull << 52);
246  return f - 1.;
247  }
248 
249  double to_nzur(uint64_t x)
250  {
251  return to_ur(x) + std::numeric_limits<double>::epsilon() / 8;
252  }
253  };
254 
255 
256  struct float2
257  {
258  float f[2];
259  };
260 
261  EIGEN_STRONG_INLINE float2 bit_to_ur_float(uint64_t x)
262  {
263  BitScalar<float> bs;
264  float2 ret;
265  ret.f[0] = bs.to_ur(x & 0xFFFFFFFF);
266  ret.f[1] = bs.to_ur(x >> 32);
267  return ret;
268  }
269 
270  template<typename Packet>
271  EIGEN_STRONG_INLINE typename std::enable_if<
272  IsFloatPacket<Packet>::value
273  >::type psincos(Packet x, Packet& s, Packet& c)
274  {
275  Packet xmm1, xmm2, xmm3 = pset1<Packet>(0), sign_bit_sin, y;
276  using IntPacket = decltype(reinterpret_to_int(x));
277  IntPacket emm0, emm2, emm4;
278 
279  sign_bit_sin = x;
280  /* take the absolute value */
281  x = pabs(x);
282  /* extract the sign bit (upper one) */
283  sign_bit_sin = pext_sign(sign_bit_sin);
284 
285  /* scale by 4/Pi */
286  y = pmul(x, pset1<Packet>(1.27323954473516));
287 
288  /* store the integer part of y in emm2 */
289  emm2 = pcast<Packet, IntPacket>(y);
290 
291  /* j=(j+1) & (~1) (see the cephes sources) */
292  emm2 = padd(emm2, pset1<IntPacket>(1));
293  emm2 = pand(emm2, pset1<IntPacket>(~1));
294  y = pcast<IntPacket, Packet>(emm2);
295 
296  emm4 = emm2;
297 
298  /* get the swap sign flag for the sine */
299  emm0 = pand(emm2, pset1<IntPacket>(4));
300  emm0 = psll<29>(emm0);
301  Packet swap_sign_bit_sin = reinterpret_to_float(emm0);
302 
303  /* get the polynom selection mask for the sine*/
304  emm2 = pand(emm2, pset1<IntPacket>(2));
305 
306  emm2 = pcmpeq(emm2, pset1<IntPacket>(0));
307  Packet poly_mask = reinterpret_to_float(emm2);
308 
309  /* The magic pass: "Extended precision modular arithmetic"
310  x = ((x - y * DP1) - y * DP2) - y * DP3; */
311  xmm1 = pset1<Packet>(-0.78515625);
312  xmm2 = pset1<Packet>(-2.4187564849853515625e-4);
313  xmm3 = pset1<Packet>(-3.77489497744594108e-8);
314  xmm1 = pmul(y, xmm1);
315  xmm2 = pmul(y, xmm2);
316  xmm3 = pmul(y, xmm3);
317  x = padd(x, xmm1);
318  x = padd(x, xmm2);
319  x = padd(x, xmm3);
320 
321  emm4 = psub(emm4, pset1<IntPacket>(2));
322  #if defined(EIGEN_VECTORIZE_NEON) || defined(EIGENRAND_EIGEN_34_MODE)
323  emm4 = pandnot(pset1<IntPacket>(4), emm4);
324  #else
325  emm4 = pandnot(emm4, pset1<IntPacket>(4));
326  #endif
327  emm4 = psll<29>(emm4);
328  Packet sign_bit_cos = reinterpret_to_float(emm4);
329  sign_bit_sin = pxor(sign_bit_sin, swap_sign_bit_sin);
330 
331 
332  /* Evaluate the first polynom (0 <= x <= Pi/4) */
333  Packet z = pmul(x, x);
334  y = pset1<Packet>(2.443315711809948E-005);
335 
336  y = pmul(y, z);
337  y = padd(y, pset1<Packet>(-1.388731625493765E-003));
338  y = pmul(y, z);
339  y = padd(y, pset1<Packet>(4.166664568298827E-002));
340  y = pmul(y, z);
341  y = pmul(y, z);
342  Packet tmp = pmul(z, pset1<Packet>(0.5));
343  y = psub(y, tmp);
344  y = padd(y, pset1<Packet>(1));
345 
346  /* Evaluate the second polynom (Pi/4 <= x <= 0) */
347 
348  Packet y2 = pset1<Packet>(-1.9515295891E-4);
349  y2 = pmul(y2, z);
350  y2 = padd(y2, pset1<Packet>(8.3321608736E-3));
351  y2 = pmul(y2, z);
352  y2 = padd(y2, pset1<Packet>(-1.6666654611E-1));
353  y2 = pmul(y2, z);
354  y2 = pmul(y2, x);
355  y2 = padd(y2, x);
356 
357  /* select the correct result from the two polynoms */
358  xmm3 = poly_mask;
359  Packet ysin2 = pand(xmm3, y2);
360  #if defined(EIGEN_VECTORIZE_NEON) || defined(EIGENRAND_EIGEN_34_MODE)
361  Packet ysin1 = pandnot(y, xmm3);
362  #else
363  Packet ysin1 = pandnot(xmm3, y);
364  #endif
365  y2 = psub(y2, ysin2);
366  y = psub(y, ysin1);
367 
368  xmm1 = padd(ysin1, ysin2);
369  xmm2 = padd(y, y2);
370 
371  /* update the sign */
372  s = pxor(xmm1, sign_bit_sin);
373  c = pxor(xmm2, sign_bit_cos);
374  }
375 
376  template<typename Packet>
377  EIGEN_STRONG_INLINE typename std::enable_if<
378  IsDoublePacket<Packet>::value
379  >::type psincos(Packet x, Packet& s, Packet& c)
380  {
381  Packet xmm1, xmm2, xmm3 = pset1<Packet>(0), sign_bit_sin, y;
382  using IntPacket = decltype(reinterpret_to_int(x));
383  IntPacket emm0, emm2, emm4;
384 
385  sign_bit_sin = x;
386  /* take the absolute value */
387  x = pabs(x);
388  /* extract the sign bit (upper one) */
389  sign_bit_sin = pext_sign(sign_bit_sin);
390 
391  /* scale by 4/Pi */
392  y = pmul(x, pset1<Packet>(1.27323954473516));
393 
394  /* store the integer part of y in emm2 */
395  emm2 = pcast64<Packet, IntPacket>(y);
396 
397  /* j=(j+1) & (~1) (see the cephes sources) */
398  emm2 = padd64(emm2, pseti64<IntPacket>(1));
399  emm2 = pand(emm2, pseti64<IntPacket>(~1ll));
400  y = pcast64<IntPacket, Packet>(emm2);
401 
402  emm4 = emm2;
403 
404  /* get the swap sign flag for the sine */
405  emm0 = pand(emm2, pseti64<IntPacket>(4));
406  emm0 = psll64<61>(emm0);
407  Packet swap_sign_bit_sin = reinterpret_to_double(emm0);
408 
409  /* get the polynom selection mask for the sine*/
410  emm2 = pand(emm2, pseti64<IntPacket>(2));
411 
412  emm2 = pcmpeq64(emm2, pseti64<IntPacket>(0));
413  Packet poly_mask = reinterpret_to_double(emm2);
414 
415  /* The magic pass: "Extended precision modular arithmetic"
416  x = ((x - y * DP1) - y * DP2) - y * DP3; */
417  xmm1 = pset1<Packet>(-0.78515625);
418  xmm2 = pset1<Packet>(-2.4187564849853515625e-4);
419  xmm3 = pset1<Packet>(-3.77489497744594108e-8);
420  xmm1 = pmul(y, xmm1);
421  xmm2 = pmul(y, xmm2);
422  xmm3 = pmul(y, xmm3);
423  x = padd(x, xmm1);
424  x = padd(x, xmm2);
425  x = padd(x, xmm3);
426 
427  emm4 = psub64(emm4, pseti64<IntPacket>(2));
428  #if defined(EIGEN_VECTORIZE_NEON) || defined(EIGENRAND_EIGEN_34_MODE)
429  emm4 = pandnot(pseti64<IntPacket>(4), emm4);
430  #else
431  emm4 = pandnot(emm4, pseti64<IntPacket>(4));
432  #endif
433  emm4 = psll64<61>(emm4);
434  Packet sign_bit_cos = reinterpret_to_double(emm4);
435  sign_bit_sin = pxor(sign_bit_sin, swap_sign_bit_sin);
436 
437 
438  /* Evaluate the first polynom (0 <= x <= Pi/4) */
439  Packet z = pmul(x, x);
440  y = pset1<Packet>(2.443315711809948E-005);
441 
442  y = pmul(y, z);
443  y = padd(y, pset1<Packet>(-1.388731625493765E-003));
444  y = pmul(y, z);
445  y = padd(y, pset1<Packet>(4.166664568298827E-002));
446  y = pmul(y, z);
447  y = pmul(y, z);
448  Packet tmp = pmul(z, pset1<Packet>(0.5));
449  y = psub(y, tmp);
450  y = padd(y, pset1<Packet>(1));
451 
452  /* Evaluate the second polynom (Pi/4 <= x <= 0) */
453 
454  Packet y2 = pset1<Packet>(-1.9515295891E-4);
455  y2 = pmul(y2, z);
456  y2 = padd(y2, pset1<Packet>(8.3321608736E-3));
457  y2 = pmul(y2, z);
458  y2 = padd(y2, pset1<Packet>(-1.6666654611E-1));
459  y2 = pmul(y2, z);
460  y2 = pmul(y2, x);
461  y2 = padd(y2, x);
462 
463  /* select the correct result from the two polynoms */
464  xmm3 = poly_mask;
465  Packet ysin2 = pand(xmm3, y2);
466  #if defined(EIGEN_VECTORIZE_NEON) || defined(EIGENRAND_EIGEN_34_MODE)
467  Packet ysin1 = pandnot(y, xmm3);
468  #else
469  Packet ysin1 = pandnot(xmm3, y);
470  #endif
471  y2 = psub(y2, ysin2);
472  y = psub(y, ysin1);
473 
474  xmm1 = padd(ysin1, ysin2);
475  xmm2 = padd(y, y2);
476 
477  /* update the sign */
478  s = pxor(xmm1, sign_bit_sin);
479  c = pxor(xmm2, sign_bit_cos);
480  }
481 
482  template<typename Packet>
483  EIGEN_STRONG_INLINE typename std::enable_if<
484  IsDoublePacket<Packet>::value, Packet
485  >::type _psin(Packet x)
486  {
487  Packet xmm1, xmm2, xmm3 = pset1<Packet>(0), sign_bit_sin, y;
488  using IntPacket = decltype(reinterpret_to_int(x));
489  IntPacket emm0, emm2;
490 
491  sign_bit_sin = x;
492  /* take the absolute value */
493  x = pabs(x);
494  /* extract the sign bit (upper one) */
495  sign_bit_sin = pext_sign(sign_bit_sin);
496 
497  /* scale by 4/Pi */
498  y = pmul(x, pset1<Packet>(1.27323954473516));
499 
500  /* store the integer part of y in emm2 */
501  emm2 = pcast64<Packet, IntPacket>(y);
502 
503  /* j=(j+1) & (~1) (see the cephes sources) */
504  emm2 = padd64(emm2, pseti64<IntPacket>(1));
505  emm2 = pand(emm2, pseti64<IntPacket>(~1ll));
506  y = pcast64<IntPacket, Packet>(emm2);
507 
508  /* get the swap sign flag for the sine */
509  emm0 = pand(emm2, pseti64<IntPacket>(4));
510  emm0 = psll64<61>(emm0);
511  Packet swap_sign_bit_sin = reinterpret_to_double(emm0);
512 
513  /* get the polynom selection mask for the sine*/
514  emm2 = pand(emm2, pseti64<IntPacket>(2));
515 
516  emm2 = pcmpeq64(emm2, pseti64<IntPacket>(0));
517  Packet poly_mask = reinterpret_to_double(emm2);
518 
519  /* The magic pass: "Extended precision modular arithmetic"
520  x = ((x - y * DP1) - y * DP2) - y * DP3; */
521  xmm1 = pset1<Packet>(-0.78515625);
522  xmm2 = pset1<Packet>(-2.4187564849853515625e-4);
523  xmm3 = pset1<Packet>(-3.77489497744594108e-8);
524  xmm1 = pmul(y, xmm1);
525  xmm2 = pmul(y, xmm2);
526  xmm3 = pmul(y, xmm3);
527  x = padd(x, xmm1);
528  x = padd(x, xmm2);
529  x = padd(x, xmm3);
530 
531  sign_bit_sin = pxor(sign_bit_sin, swap_sign_bit_sin);
532 
533 
534  /* Evaluate the first polynom (0 <= x <= Pi/4) */
535  Packet z = pmul(x, x);
536  y = pset1<Packet>(2.443315711809948E-005);
537 
538  y = pmul(y, z);
539  y = padd(y, pset1<Packet>(-1.388731625493765E-003));
540  y = pmul(y, z);
541  y = padd(y, pset1<Packet>(4.166664568298827E-002));
542  y = pmul(y, z);
543  y = pmul(y, z);
544  Packet tmp = pmul(z, pset1<Packet>(0.5));
545  y = psub(y, tmp);
546  y = padd(y, pset1<Packet>(1));
547 
548  /* Evaluate the second polynom (Pi/4 <= x <= 0) */
549 
550  Packet y2 = pset1<Packet>(-1.9515295891E-4);
551  y2 = pmul(y2, z);
552  y2 = padd(y2, pset1<Packet>(8.3321608736E-3));
553  y2 = pmul(y2, z);
554  y2 = padd(y2, pset1<Packet>(-1.6666654611E-1));
555  y2 = pmul(y2, z);
556  y2 = pmul(y2, x);
557  y2 = padd(y2, x);
558 
559  /* select the correct result from the two polynoms */
560  xmm3 = poly_mask;
561  Packet ysin2 = pand(xmm3, y2);
562  #if defined(EIGEN_VECTORIZE_NEON) || defined(EIGENRAND_EIGEN_34_MODE)
563  Packet ysin1 = pandnot(y, xmm3);
564  #else
565  Packet ysin1 = pandnot(xmm3, y);
566  #endif
567 
568  xmm1 = padd(ysin1, ysin2);
569 
570  /* update the sign */
571  return pxor(xmm1, sign_bit_sin);
572  }
573  }
574 }
575 
576 #ifdef EIGEN_VECTORIZE_AVX
577 #include "arch/AVX/MorePacketMath.h"
578 #endif
579 
580 #ifdef EIGEN_VECTORIZE_SSE2
581 #include "arch/SSE/MorePacketMath.h"
582 #endif
583 
584 #ifdef EIGEN_VECTORIZE_NEON
585 #include "arch/NEON/MorePacketMath.h"
586 #endif
587 
588 namespace Eigen
589 {
590  namespace internal
591  {
592  template<int b, typename Packet>
593  EIGEN_STRONG_INLINE Packet psll(const Packet& a)
594  {
595  return BitShifter<Packet>{}.template sll<b>(a);
596  }
597 
598  template<int _b, typename Packet>
599  EIGEN_STRONG_INLINE Packet psrl(const Packet& a, int b)
600  {
601  return BitShifter<Packet>{}.template srl<_b>(a, b);
602  }
603 
604  template<int b, typename Packet>
605  EIGEN_STRONG_INLINE Packet psll64(const Packet& a)
606  {
607  return BitShifter<Packet>{}.template sll64<b>(a);
608  }
609 
610  template<int b, typename Packet>
611  EIGEN_STRONG_INLINE Packet psrl64(const Packet& a)
612  {
613  return BitShifter<Packet>{}.template srl64<b>(a);
614  }
615  }
616 }
617 
618 #endif