1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2009 Gael Guennebaud <[email protected]>
5 // Copyright (C) 2006-2008 Benoit Jacob <[email protected]>
6 //
7 // This Source Code Form is subject to the terms of the Mozilla
8 // Public License v. 2.0. If a copy of the MPL was not distributed
9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 #include "packetmath_test_shared.h"
12 #include "random_without_cast_overflow.h"
13
14 template <typename T>
REF_ADD(const T & a,const T & b)15 inline T REF_ADD(const T& a, const T& b) {
16 return a + b;
17 }
18 template <typename T>
REF_SUB(const T & a,const T & b)19 inline T REF_SUB(const T& a, const T& b) {
20 return a - b;
21 }
22 template <typename T>
REF_MUL(const T & a,const T & b)23 inline T REF_MUL(const T& a, const T& b) {
24 return a * b;
25 }
26 template <typename T>
REF_DIV(const T & a,const T & b)27 inline T REF_DIV(const T& a, const T& b) {
28 return a / b;
29 }
30 template <typename T>
REF_ABS_DIFF(const T & a,const T & b)31 inline T REF_ABS_DIFF(const T& a, const T& b) {
32 return a > b ? a - b : b - a;
33 }
34
35 // Specializations for bool.
36 template <>
REF_ADD(const bool & a,const bool & b)37 inline bool REF_ADD(const bool& a, const bool& b) {
38 return a || b;
39 }
40 template <>
REF_SUB(const bool & a,const bool & b)41 inline bool REF_SUB(const bool& a, const bool& b) {
42 return a ^ b;
43 }
44 template <>
REF_MUL(const bool & a,const bool & b)45 inline bool REF_MUL(const bool& a, const bool& b) {
46 return a && b;
47 }
48
49 template <typename T>
REF_FREXP(const T & x,T & exp)50 inline T REF_FREXP(const T& x, T& exp) {
51 int iexp;
52 EIGEN_USING_STD(frexp)
53 const T out = static_cast<T>(frexp(x, &iexp));
54 exp = static_cast<T>(iexp);
55 return out;
56 }
57
58 template <typename T>
REF_LDEXP(const T & x,const T & exp)59 inline T REF_LDEXP(const T& x, const T& exp) {
60 EIGEN_USING_STD(ldexp)
61 return static_cast<T>(ldexp(x, static_cast<int>(exp)));
62 }
63
64 // Uses pcast to cast from one array to another.
65 template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
66 struct pcast_array;
67
68 template <typename SrcPacket, typename TgtPacket, int TgtCoeffRatio>
69 struct pcast_array<SrcPacket, TgtPacket, 1, TgtCoeffRatio> {
70 typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar;
71 typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar;
castpcast_array72 static void cast(const SrcScalar* src, size_t size, TgtScalar* dst) {
73 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
74 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
75 size_t i;
76 for (i = 0; i < size && i + SrcPacketSize <= size; i += TgtPacketSize) {
77 internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(internal::ploadu<SrcPacket>(src + i)));
78 }
79 // Leftovers that cannot be loaded into a packet.
80 for (; i < size; ++i) {
81 dst[i] = static_cast<TgtScalar>(src[i]);
82 }
83 }
84 };
85
86 template <typename SrcPacket, typename TgtPacket>
87 struct pcast_array<SrcPacket, TgtPacket, 2, 1> {
castpcast_array88 static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size,
89 typename internal::unpacket_traits<TgtPacket>::type* dst) {
90 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
91 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
92 for (size_t i = 0; i < size; i += TgtPacketSize) {
93 SrcPacket a = internal::ploadu<SrcPacket>(src + i);
94 SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize);
95 internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b));
96 }
97 }
98 };
99
100 template <typename SrcPacket, typename TgtPacket>
101 struct pcast_array<SrcPacket, TgtPacket, 4, 1> {
castpcast_array102 static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size,
103 typename internal::unpacket_traits<TgtPacket>::type* dst) {
104 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
105 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
106 for (size_t i = 0; i < size; i += TgtPacketSize) {
107 SrcPacket a = internal::ploadu<SrcPacket>(src + i);
108 SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize);
109 SrcPacket c = internal::ploadu<SrcPacket>(src + i + 2 * SrcPacketSize);
110 SrcPacket d = internal::ploadu<SrcPacket>(src + i + 3 * SrcPacketSize);
111 internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b, c, d));
112 }
113 }
114 };
115
116 template <typename SrcPacket, typename TgtPacket>
117 struct pcast_array<SrcPacket, TgtPacket, 8, 1> {
castpcast_array118 static void cast(const typename internal::unpacket_traits<SrcPacket>::type* src, size_t size,
119 typename internal::unpacket_traits<TgtPacket>::type* dst) {
120 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
121 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
122 for (size_t i = 0; i < size; i += TgtPacketSize) {
123 SrcPacket a = internal::ploadu<SrcPacket>(src + i);
124 SrcPacket b = internal::ploadu<SrcPacket>(src + i + SrcPacketSize);
125 SrcPacket c = internal::ploadu<SrcPacket>(src + i + 2 * SrcPacketSize);
126 SrcPacket d = internal::ploadu<SrcPacket>(src + i + 3 * SrcPacketSize);
127 SrcPacket e = internal::ploadu<SrcPacket>(src + i + 4 * SrcPacketSize);
128 SrcPacket f = internal::ploadu<SrcPacket>(src + i + 5 * SrcPacketSize);
129 SrcPacket g = internal::ploadu<SrcPacket>(src + i + 6 * SrcPacketSize);
130 SrcPacket h = internal::ploadu<SrcPacket>(src + i + 7 * SrcPacketSize);
131 internal::pstoreu(dst + i, internal::pcast<SrcPacket, TgtPacket>(a, b, c, d, e, f, g, h));
132 }
133 }
134 };
135
136 template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio, bool CanCast = false>
137 struct test_cast_helper;
138
139 template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
140 struct test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, false> {
runtest_cast_helper141 static void run() {}
142 };
143
144 template <typename SrcPacket, typename TgtPacket, int SrcCoeffRatio, int TgtCoeffRatio>
145 struct test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, true> {
runtest_cast_helper146 static void run() {
147 typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar;
148 typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar;
149 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
150 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
151 static const int BlockSize = SrcPacketSize * SrcCoeffRatio;
152 eigen_assert(BlockSize == TgtPacketSize * TgtCoeffRatio && "Packet sizes and cast ratios are mismatched.");
153
154 static const int DataSize = 10 * BlockSize;
155 EIGEN_ALIGN_MAX SrcScalar data1[DataSize];
156 EIGEN_ALIGN_MAX TgtScalar data2[DataSize];
157 EIGEN_ALIGN_MAX TgtScalar ref[DataSize];
158
159 // Construct a packet of scalars that will not overflow when casting
160 for (int i = 0; i < DataSize; ++i) {
161 data1[i] = internal::random_without_cast_overflow<SrcScalar, TgtScalar>::value();
162 }
163
164 for (int i = 0; i < DataSize; ++i) {
165 ref[i] = static_cast<const TgtScalar>(data1[i]);
166 }
167
168 pcast_array<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio>::cast(data1, DataSize, data2);
169
170 VERIFY(test::areApprox(ref, data2, DataSize) && "internal::pcast<>");
171 }
172 };
173
174 template <typename SrcPacket, typename TgtPacket>
175 struct test_cast {
runtest_cast176 static void run() {
177 typedef typename internal::unpacket_traits<SrcPacket>::type SrcScalar;
178 typedef typename internal::unpacket_traits<TgtPacket>::type TgtScalar;
179 typedef typename internal::type_casting_traits<SrcScalar, TgtScalar> TypeCastingTraits;
180 static const int SrcCoeffRatio = TypeCastingTraits::SrcCoeffRatio;
181 static const int TgtCoeffRatio = TypeCastingTraits::TgtCoeffRatio;
182 static const int SrcPacketSize = internal::unpacket_traits<SrcPacket>::size;
183 static const int TgtPacketSize = internal::unpacket_traits<TgtPacket>::size;
184 static const bool HasCast =
185 internal::unpacket_traits<SrcPacket>::vectorizable && internal::unpacket_traits<TgtPacket>::vectorizable &&
186 TypeCastingTraits::VectorizedCast && (SrcPacketSize * SrcCoeffRatio == TgtPacketSize * TgtCoeffRatio);
187 test_cast_helper<SrcPacket, TgtPacket, SrcCoeffRatio, TgtCoeffRatio, HasCast>::run();
188 }
189 };
190
191 template <typename SrcPacket, typename TgtScalar,
192 typename TgtPacket = typename internal::packet_traits<TgtScalar>::type,
193 bool Vectorized = internal::packet_traits<TgtScalar>::Vectorizable,
194 bool HasHalf = !internal::is_same<typename internal::unpacket_traits<TgtPacket>::half, TgtPacket>::value>
195 struct test_cast_runner;
196
197 template <typename SrcPacket, typename TgtScalar, typename TgtPacket>
198 struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, true, false> {
runtest_cast_runner199 static void run() { test_cast<SrcPacket, TgtPacket>::run(); }
200 };
201
202 template <typename SrcPacket, typename TgtScalar, typename TgtPacket>
203 struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, true, true> {
runtest_cast_runner204 static void run() {
205 test_cast<SrcPacket, TgtPacket>::run();
206 test_cast_runner<SrcPacket, TgtScalar, typename internal::unpacket_traits<TgtPacket>::half>::run();
207 }
208 };
209
210 template <typename SrcPacket, typename TgtScalar, typename TgtPacket>
211 struct test_cast_runner<SrcPacket, TgtScalar, TgtPacket, false, false> {
runtest_cast_runner212 static void run() {}
213 };
214
215 template <typename Scalar, typename Packet, typename EnableIf = void>
216 struct packetmath_pcast_ops_runner {
runpacketmath_pcast_ops_runner217 static void run() {
218 test_cast_runner<Packet, float>::run();
219 test_cast_runner<Packet, double>::run();
220 test_cast_runner<Packet, int8_t>::run();
221 test_cast_runner<Packet, uint8_t>::run();
222 test_cast_runner<Packet, int16_t>::run();
223 test_cast_runner<Packet, uint16_t>::run();
224 test_cast_runner<Packet, int32_t>::run();
225 test_cast_runner<Packet, uint32_t>::run();
226 test_cast_runner<Packet, int64_t>::run();
227 test_cast_runner<Packet, uint64_t>::run();
228 test_cast_runner<Packet, bool>::run();
229 test_cast_runner<Packet, std::complex<float> >::run();
230 test_cast_runner<Packet, std::complex<double> >::run();
231 test_cast_runner<Packet, half>::run();
232 test_cast_runner<Packet, bfloat16>::run();
233 }
234 };
235
236 // Only some types support cast from std::complex<>.
237 template <typename Scalar, typename Packet>
238 struct packetmath_pcast_ops_runner<Scalar, Packet, typename internal::enable_if<NumTraits<Scalar>::IsComplex>::type> {
runpacketmath_pcast_ops_runner239 static void run() {
240 test_cast_runner<Packet, std::complex<float> >::run();
241 test_cast_runner<Packet, std::complex<double> >::run();
242 test_cast_runner<Packet, half>::run();
243 test_cast_runner<Packet, bfloat16>::run();
244 }
245 };
246
247 template <typename Scalar, typename Packet>
packetmath_boolean_mask_ops()248 void packetmath_boolean_mask_ops() {
249 const int PacketSize = internal::unpacket_traits<Packet>::size;
250 const int size = 2 * PacketSize;
251 EIGEN_ALIGN_MAX Scalar data1[size];
252 EIGEN_ALIGN_MAX Scalar data2[size];
253 EIGEN_ALIGN_MAX Scalar ref[size];
254
255 for (int i = 0; i < size; ++i) {
256 data1[i] = internal::random<Scalar>();
257 }
258 CHECK_CWISE1(internal::ptrue, internal::ptrue);
259 CHECK_CWISE2_IF(true, internal::pandnot, internal::pandnot);
260 for (int i = 0; i < PacketSize; ++i) {
261 data1[i] = Scalar(i);
262 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
263 }
264
265 CHECK_CWISE2_IF(true, internal::pcmp_eq, internal::pcmp_eq);
266
267 //Test (-0) == (0) for signed operations
268 for (int i = 0; i < PacketSize; ++i) {
269 data1[i] = Scalar(-0.0);
270 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
271 }
272 CHECK_CWISE2_IF(true, internal::pcmp_eq, internal::pcmp_eq);
273
274 //Test NaN
275 for (int i = 0; i < PacketSize; ++i) {
276 data1[i] = NumTraits<Scalar>::quiet_NaN();
277 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
278 }
279 CHECK_CWISE2_IF(true, internal::pcmp_eq, internal::pcmp_eq);
280 }
281
282 template <typename Scalar, typename Packet>
packetmath_boolean_mask_ops_real()283 void packetmath_boolean_mask_ops_real() {
284 const int PacketSize = internal::unpacket_traits<Packet>::size;
285 const int size = 2 * PacketSize;
286 EIGEN_ALIGN_MAX Scalar data1[size];
287 EIGEN_ALIGN_MAX Scalar data2[size];
288 EIGEN_ALIGN_MAX Scalar ref[size];
289
290 for (int i = 0; i < PacketSize; ++i) {
291 data1[i] = internal::random<Scalar>();
292 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
293 }
294
295 CHECK_CWISE2_IF(true, internal::pcmp_lt_or_nan, internal::pcmp_lt_or_nan);
296
297 //Test (-0) <=/< (0) for signed operations
298 for (int i = 0; i < PacketSize; ++i) {
299 data1[i] = Scalar(-0.0);
300 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
301 }
302 CHECK_CWISE2_IF(true, internal::pcmp_lt_or_nan, internal::pcmp_lt_or_nan);
303
304 //Test NaN
305 for (int i = 0; i < PacketSize; ++i) {
306 data1[i] = NumTraits<Scalar>::quiet_NaN();
307 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
308 }
309 CHECK_CWISE2_IF(true, internal::pcmp_lt_or_nan, internal::pcmp_lt_or_nan);
310 }
311
312 template <typename Scalar, typename Packet>
packetmath_boolean_mask_ops_notcomplex()313 void packetmath_boolean_mask_ops_notcomplex() {
314 const int PacketSize = internal::unpacket_traits<Packet>::size;
315 const int size = 2 * PacketSize;
316 EIGEN_ALIGN_MAX Scalar data1[size];
317 EIGEN_ALIGN_MAX Scalar data2[size];
318 EIGEN_ALIGN_MAX Scalar ref[size];
319
320 for (int i = 0; i < PacketSize; ++i) {
321 data1[i] = internal::random<Scalar>();
322 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
323 }
324
325 CHECK_CWISE2_IF(true, internal::pcmp_le, internal::pcmp_le);
326 CHECK_CWISE2_IF(true, internal::pcmp_lt, internal::pcmp_lt);
327
328 //Test (-0) <=/< (0) for signed operations
329 for (int i = 0; i < PacketSize; ++i) {
330 data1[i] = Scalar(-0.0);
331 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
332 }
333 CHECK_CWISE2_IF(true, internal::pcmp_le, internal::pcmp_le);
334 CHECK_CWISE2_IF(true, internal::pcmp_lt, internal::pcmp_lt);
335
336 //Test NaN
337 for (int i = 0; i < PacketSize; ++i) {
338 data1[i] = NumTraits<Scalar>::quiet_NaN();
339 data1[i + PacketSize] = internal::random<bool>() ? data1[i] : Scalar(0);
340 }
341 CHECK_CWISE2_IF(true, internal::pcmp_le, internal::pcmp_le);
342 CHECK_CWISE2_IF(true, internal::pcmp_lt, internal::pcmp_lt);
343 }
344
345 // Packet16b representing bool does not support ptrue, pandnot or pcmp_eq, since the scalar path
346 // (for some compilers) compute the bitwise and with 0x1 of the results to keep the value in [0,1].
347 template<>
packetmath_boolean_mask_ops()348 void packetmath_boolean_mask_ops<bool, internal::packet_traits<bool>::type>() {}
349 template<>
packetmath_boolean_mask_ops_notcomplex()350 void packetmath_boolean_mask_ops_notcomplex<bool, internal::packet_traits<bool>::type>() {}
351
352 template <typename Scalar, typename Packet>
packetmath_minus_zero_add()353 void packetmath_minus_zero_add() {
354 const int PacketSize = internal::unpacket_traits<Packet>::size;
355 const int size = 2 * PacketSize;
356 EIGEN_ALIGN_MAX Scalar data1[size];
357 EIGEN_ALIGN_MAX Scalar data2[size];
358 EIGEN_ALIGN_MAX Scalar ref[size];
359
360 for (int i = 0; i < PacketSize; ++i) {
361 data1[i] = Scalar(-0.0);
362 data1[i + PacketSize] = Scalar(-0.0);
363 }
364 CHECK_CWISE2_IF(internal::packet_traits<Scalar>::HasAdd, REF_ADD, internal::padd);
365 }
366
367 // Ensure optimization barrier compiles and doesn't modify contents.
368 // Only applies to raw types, so will not work for std::complex, Eigen::half
369 // or Eigen::bfloat16. For those you would need to refer to an underlying
370 // storage element.
371 template<typename Packet, typename EnableIf = void>
372 struct eigen_optimization_barrier_test {
runeigen_optimization_barrier_test373 static void run() {}
374 };
375
376 template<typename Packet>
377 struct eigen_optimization_barrier_test<Packet, typename internal::enable_if<
378 !NumTraits<Packet>::IsComplex &&
379 !internal::is_same<Packet, Eigen::half>::value &&
380 !internal::is_same<Packet, Eigen::bfloat16>::value
381 >::type> {
runeigen_optimization_barrier_test382 static void run() {
383 typedef typename internal::unpacket_traits<Packet>::type Scalar;
384 Scalar s = internal::random<Scalar>();
385 Packet barrier = internal::pset1<Packet>(s);
386 EIGEN_OPTIMIZATION_BARRIER(barrier);
387 eigen_assert(s == internal::pfirst(barrier) && "EIGEN_OPTIMIZATION_BARRIER");
388 }
389 };
390
391 template <typename Scalar, typename Packet>
packetmath()392 void packetmath() {
393 typedef internal::packet_traits<Scalar> PacketTraits;
394 const int PacketSize = internal::unpacket_traits<Packet>::size;
395 typedef typename NumTraits<Scalar>::Real RealScalar;
396
397 if (g_first_pass)
398 std::cerr << "=== Testing packet of type '" << typeid(Packet).name() << "' and scalar type '"
399 << typeid(Scalar).name() << "' and size '" << PacketSize << "' ===\n";
400
401 const int max_size = PacketSize > 4 ? PacketSize : 4;
402 const int size = PacketSize * max_size;
403 EIGEN_ALIGN_MAX Scalar data1[size];
404 EIGEN_ALIGN_MAX Scalar data2[size];
405 EIGEN_ALIGN_MAX Scalar data3[size];
406 EIGEN_ALIGN_MAX Scalar ref[size];
407 RealScalar refvalue = RealScalar(0);
408
409 eigen_optimization_barrier_test<Packet>::run();
410 eigen_optimization_barrier_test<Scalar>::run();
411
412 for (int i = 0; i < size; ++i) {
413 data1[i] = internal::random<Scalar>() / RealScalar(PacketSize);
414 data2[i] = internal::random<Scalar>() / RealScalar(PacketSize);
415 refvalue = (std::max)(refvalue, numext::abs(data1[i]));
416 }
417
418 internal::pstore(data2, internal::pload<Packet>(data1));
419 VERIFY(test::areApprox(data1, data2, PacketSize) && "aligned load/store");
420
421 for (int offset = 0; offset < PacketSize; ++offset) {
422 internal::pstore(data2, internal::ploadu<Packet>(data1 + offset));
423 VERIFY(test::areApprox(data1 + offset, data2, PacketSize) && "internal::ploadu");
424 }
425
426 for (int offset = 0; offset < PacketSize; ++offset) {
427 internal::pstoreu(data2 + offset, internal::pload<Packet>(data1));
428 VERIFY(test::areApprox(data1, data2 + offset, PacketSize) && "internal::pstoreu");
429 }
430
431 if (internal::unpacket_traits<Packet>::masked_load_available) {
432 test::packet_helper<internal::unpacket_traits<Packet>::masked_load_available, Packet> h;
433 unsigned long long max_umask = (0x1ull << PacketSize);
434
435 for (int offset = 0; offset < PacketSize; ++offset) {
436 for (unsigned long long umask = 0; umask < max_umask; ++umask) {
437 h.store(data2, h.load(data1 + offset, umask));
438 for (int k = 0; k < PacketSize; ++k) data3[k] = ((umask & (0x1ull << k)) >> k) ? data1[k + offset] : Scalar(0);
439 VERIFY(test::areApprox(data3, data2, PacketSize) && "internal::ploadu masked");
440 }
441 }
442 }
443
444 if (internal::unpacket_traits<Packet>::masked_store_available) {
445 test::packet_helper<internal::unpacket_traits<Packet>::masked_store_available, Packet> h;
446 unsigned long long max_umask = (0x1ull << PacketSize);
447
448 for (int offset = 0; offset < PacketSize; ++offset) {
449 for (unsigned long long umask = 0; umask < max_umask; ++umask) {
450 internal::pstore(data2, internal::pset1<Packet>(Scalar(0)));
451 h.store(data2, h.loadu(data1 + offset), umask);
452 for (int k = 0; k < PacketSize; ++k) data3[k] = ((umask & (0x1ull << k)) >> k) ? data1[k + offset] : Scalar(0);
453 VERIFY(test::areApprox(data3, data2, PacketSize) && "internal::pstoreu masked");
454 }
455 }
456 }
457
458 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasAdd);
459 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasSub);
460 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasMul);
461
462 CHECK_CWISE2_IF(PacketTraits::HasAdd, REF_ADD, internal::padd);
463 CHECK_CWISE2_IF(PacketTraits::HasSub, REF_SUB, internal::psub);
464 CHECK_CWISE2_IF(PacketTraits::HasMul, REF_MUL, internal::pmul);
465 CHECK_CWISE2_IF(PacketTraits::HasDiv, REF_DIV, internal::pdiv);
466
467 if (PacketTraits::HasNegate) CHECK_CWISE1(internal::negate, internal::pnegate);
468 CHECK_CWISE1(numext::conj, internal::pconj);
469
470 for (int offset = 0; offset < 3; ++offset) {
471 for (int i = 0; i < PacketSize; ++i) ref[i] = data1[offset];
472 internal::pstore(data2, internal::pset1<Packet>(data1[offset]));
473 VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::pset1");
474 }
475
476 {
477 for (int i = 0; i < PacketSize * 4; ++i) ref[i] = data1[i / PacketSize];
478 Packet A0, A1, A2, A3;
479 internal::pbroadcast4<Packet>(data1, A0, A1, A2, A3);
480 internal::pstore(data2 + 0 * PacketSize, A0);
481 internal::pstore(data2 + 1 * PacketSize, A1);
482 internal::pstore(data2 + 2 * PacketSize, A2);
483 internal::pstore(data2 + 3 * PacketSize, A3);
484 VERIFY(test::areApprox(ref, data2, 4 * PacketSize) && "internal::pbroadcast4");
485 }
486
487 {
488 for (int i = 0; i < PacketSize * 2; ++i) ref[i] = data1[i / PacketSize];
489 Packet A0, A1;
490 internal::pbroadcast2<Packet>(data1, A0, A1);
491 internal::pstore(data2 + 0 * PacketSize, A0);
492 internal::pstore(data2 + 1 * PacketSize, A1);
493 VERIFY(test::areApprox(ref, data2, 2 * PacketSize) && "internal::pbroadcast2");
494 }
495
496 VERIFY(internal::isApprox(data1[0], internal::pfirst(internal::pload<Packet>(data1))) && "internal::pfirst");
497
498 if (PacketSize > 1) {
499 // apply different offsets to check that ploaddup is robust to unaligned inputs
500 for (int offset = 0; offset < 4; ++offset) {
501 for (int i = 0; i < PacketSize / 2; ++i) ref[2 * i + 0] = ref[2 * i + 1] = data1[offset + i];
502 internal::pstore(data2, internal::ploaddup<Packet>(data1 + offset));
503 VERIFY(test::areApprox(ref, data2, PacketSize) && "ploaddup");
504 }
505 }
506
507 if (PacketSize > 2) {
508 // apply different offsets to check that ploadquad is robust to unaligned inputs
509 for (int offset = 0; offset < 4; ++offset) {
510 for (int i = 0; i < PacketSize / 4; ++i)
511 ref[4 * i + 0] = ref[4 * i + 1] = ref[4 * i + 2] = ref[4 * i + 3] = data1[offset + i];
512 internal::pstore(data2, internal::ploadquad<Packet>(data1 + offset));
513 VERIFY(test::areApprox(ref, data2, PacketSize) && "ploadquad");
514 }
515 }
516
517 ref[0] = Scalar(0);
518 for (int i = 0; i < PacketSize; ++i) ref[0] += data1[i];
519 VERIFY(test::isApproxAbs(ref[0], internal::predux(internal::pload<Packet>(data1)), refvalue) && "internal::predux");
520
521 if (!internal::is_same<Packet, typename internal::unpacket_traits<Packet>::half>::value) {
522 int HalfPacketSize = PacketSize > 4 ? PacketSize / 2 : PacketSize;
523 for (int i = 0; i < HalfPacketSize; ++i) ref[i] = Scalar(0);
524 for (int i = 0; i < PacketSize; ++i) ref[i % HalfPacketSize] += data1[i];
525 internal::pstore(data2, internal::predux_half_dowto4(internal::pload<Packet>(data1)));
526 VERIFY(test::areApprox(ref, data2, HalfPacketSize) && "internal::predux_half_dowto4");
527 }
528
529 ref[0] = Scalar(1);
530 for (int i = 0; i < PacketSize; ++i) ref[0] = REF_MUL(ref[0], data1[i]);
531 VERIFY(internal::isApprox(ref[0], internal::predux_mul(internal::pload<Packet>(data1))) && "internal::predux_mul");
532
533 for (int i = 0; i < PacketSize; ++i) ref[i] = data1[PacketSize - i - 1];
534 internal::pstore(data2, internal::preverse(internal::pload<Packet>(data1)));
535 VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::preverse");
536
537 internal::PacketBlock<Packet> kernel;
538 for (int i = 0; i < PacketSize; ++i) {
539 kernel.packet[i] = internal::pload<Packet>(data1 + i * PacketSize);
540 }
541 ptranspose(kernel);
542 for (int i = 0; i < PacketSize; ++i) {
543 internal::pstore(data2, kernel.packet[i]);
544 for (int j = 0; j < PacketSize; ++j) {
545 VERIFY(test::isApproxAbs(data2[j], data1[i + j * PacketSize], refvalue) && "ptranspose");
546 }
547 }
548
549 // GeneralBlockPanelKernel also checks PacketBlock<Packet,(PacketSize%4)==0?4:PacketSize>;
550 if (PacketSize > 4 && PacketSize % 4 == 0) {
551 internal::PacketBlock<Packet, PacketSize%4==0?4:PacketSize> kernel2;
552 for (int i = 0; i < 4; ++i) {
553 kernel2.packet[i] = internal::pload<Packet>(data1 + i * PacketSize);
554 }
555 ptranspose(kernel2);
556 int data_counter = 0;
557 for (int i = 0; i < PacketSize; ++i) {
558 for (int j = 0; j < 4; ++j) {
559 data2[data_counter++] = data1[j*PacketSize + i];
560 }
561 }
562 for (int i = 0; i < 4; ++i) {
563 internal::pstore(data3, kernel2.packet[i]);
564 for (int j = 0; j < PacketSize; ++j) {
565 VERIFY(test::isApproxAbs(data3[j], data2[i*PacketSize + j], refvalue) && "ptranspose");
566 }
567 }
568 }
569
570 if (PacketTraits::HasBlend) {
571 Packet thenPacket = internal::pload<Packet>(data1);
572 Packet elsePacket = internal::pload<Packet>(data2);
573 EIGEN_ALIGN_MAX internal::Selector<PacketSize> selector;
574 for (int i = 0; i < PacketSize; ++i) {
575 selector.select[i] = i;
576 }
577
578 Packet blend = internal::pblend(selector, thenPacket, elsePacket);
579 EIGEN_ALIGN_MAX Scalar result[size];
580 internal::pstore(result, blend);
581 for (int i = 0; i < PacketSize; ++i) {
582 VERIFY(test::isApproxAbs(result[i], (selector.select[i] ? data1[i] : data2[i]), refvalue));
583 }
584 }
585
586 {
587 for (int i = 0; i < PacketSize; ++i) {
588 // "if" mask
589 unsigned char v = internal::random<bool>() ? 0xff : 0;
590 char* bytes = (char*)(data1 + i);
591 for (int k = 0; k < int(sizeof(Scalar)); ++k) {
592 bytes[k] = v;
593 }
594 // "then" packet
595 data1[i + PacketSize] = internal::random<Scalar>();
596 // "else" packet
597 data1[i + 2 * PacketSize] = internal::random<Scalar>();
598 }
599 CHECK_CWISE3_IF(true, internal::pselect, internal::pselect);
600 }
601
602 for (int i = 0; i < size; ++i) {
603 data1[i] = internal::random<Scalar>();
604 }
605 CHECK_CWISE1(internal::pzero, internal::pzero);
606 CHECK_CWISE2_IF(true, internal::por, internal::por);
607 CHECK_CWISE2_IF(true, internal::pxor, internal::pxor);
608 CHECK_CWISE2_IF(true, internal::pand, internal::pand);
609
610 packetmath_boolean_mask_ops<Scalar, Packet>();
611 packetmath_pcast_ops_runner<Scalar, Packet>::run();
612 packetmath_minus_zero_add<Scalar, Packet>();
613
614 for (int i = 0; i < size; ++i) {
615 data1[i] = numext::abs(internal::random<Scalar>());
616 }
617 CHECK_CWISE1_IF(PacketTraits::HasSqrt, numext::sqrt, internal::psqrt);
618 CHECK_CWISE1_IF(PacketTraits::HasRsqrt, numext::rsqrt, internal::prsqrt);
619 }
620
621 // Notice that this definition works for complex types as well.
622 // c++11 has std::log2 for real, but not for complex types.
623 template <typename Scalar>
log2(Scalar x)624 Scalar log2(Scalar x) {
625 return Scalar(EIGEN_LOG2E) * std::log(x);
626 }
627
628 template <typename Scalar, typename Packet>
packetmath_real()629 void packetmath_real() {
630 typedef internal::packet_traits<Scalar> PacketTraits;
631 const int PacketSize = internal::unpacket_traits<Packet>::size;
632
633 const int size = PacketSize * 4;
634 EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4];
635 EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4];
636 EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4];
637
638 for (int i = 0; i < size; ++i) {
639 data1[i] = Scalar(internal::random<double>(0, 1) * std::pow(10., internal::random<double>(-6, 6)));
640 data2[i] = Scalar(internal::random<double>(0, 1) * std::pow(10., internal::random<double>(-6, 6)));
641 }
642
643 if (internal::random<float>(0, 1) < 0.1f) data1[internal::random<int>(0, PacketSize)] = Scalar(0);
644
645 CHECK_CWISE1_IF(PacketTraits::HasLog, std::log, internal::plog);
646 CHECK_CWISE1_IF(PacketTraits::HasLog, log2, internal::plog2);
647 CHECK_CWISE1_IF(PacketTraits::HasRsqrt, numext::rsqrt, internal::prsqrt);
648
649 for (int i = 0; i < size; ++i) {
650 data1[i] = Scalar(internal::random<double>(-1, 1) * std::pow(10., internal::random<double>(-3, 3)));
651 data2[i] = Scalar(internal::random<double>(-1, 1) * std::pow(10., internal::random<double>(-3, 3)));
652 }
653 CHECK_CWISE1_IF(PacketTraits::HasSin, std::sin, internal::psin);
654 CHECK_CWISE1_IF(PacketTraits::HasCos, std::cos, internal::pcos);
655 CHECK_CWISE1_IF(PacketTraits::HasTan, std::tan, internal::ptan);
656
657 CHECK_CWISE1_EXACT_IF(PacketTraits::HasRound, numext::round, internal::pround);
658 CHECK_CWISE1_EXACT_IF(PacketTraits::HasCeil, numext::ceil, internal::pceil);
659 CHECK_CWISE1_EXACT_IF(PacketTraits::HasFloor, numext::floor, internal::pfloor);
660 CHECK_CWISE1_EXACT_IF(PacketTraits::HasRint, numext::rint, internal::print);
661
662 packetmath_boolean_mask_ops_real<Scalar,Packet>();
663
664 // Rounding edge cases.
665 if (PacketTraits::HasRound || PacketTraits::HasCeil || PacketTraits::HasFloor || PacketTraits::HasRint) {
666 typedef typename internal::make_integer<Scalar>::type IntType;
667 // Start with values that cannot fit inside an integer, work down to less than one.
668 Scalar val = numext::mini(
669 Scalar(2) * static_cast<Scalar>(NumTraits<IntType>::highest()),
670 NumTraits<Scalar>::highest());
671 std::vector<Scalar> values;
672 while (val > Scalar(0.25)) {
673 // Cover both even and odd, positive and negative cases.
674 values.push_back(val);
675 values.push_back(val + Scalar(0.3));
676 values.push_back(val + Scalar(0.5));
677 values.push_back(val + Scalar(0.8));
678 values.push_back(val + Scalar(1));
679 values.push_back(val + Scalar(1.3));
680 values.push_back(val + Scalar(1.5));
681 values.push_back(val + Scalar(1.8));
682 values.push_back(-val);
683 values.push_back(-val - Scalar(0.3));
684 values.push_back(-val - Scalar(0.5));
685 values.push_back(-val - Scalar(0.8));
686 values.push_back(-val - Scalar(1));
687 values.push_back(-val - Scalar(1.3));
688 values.push_back(-val - Scalar(1.5));
689 values.push_back(-val - Scalar(1.8));
690 values.push_back(Scalar(-1.5) + val); // Bug 1785.
691 val = val / Scalar(2);
692 }
693 values.push_back(NumTraits<Scalar>::infinity());
694 values.push_back(-NumTraits<Scalar>::infinity());
695 values.push_back(NumTraits<Scalar>::quiet_NaN());
696
697 for (size_t k=0; k<values.size(); ++k) {
698 data1[0] = values[k];
699 CHECK_CWISE1_EXACT_IF(PacketTraits::HasRound, numext::round, internal::pround);
700 CHECK_CWISE1_EXACT_IF(PacketTraits::HasCeil, numext::ceil, internal::pceil);
701 CHECK_CWISE1_EXACT_IF(PacketTraits::HasFloor, numext::floor, internal::pfloor);
702 CHECK_CWISE1_EXACT_IF(PacketTraits::HasRint, numext::rint, internal::print);
703 }
704 }
705
706 for (int i = 0; i < size; ++i) {
707 data1[i] = Scalar(internal::random<double>(-1, 1));
708 data2[i] = Scalar(internal::random<double>(-1, 1));
709 }
710 CHECK_CWISE1_IF(PacketTraits::HasASin, std::asin, internal::pasin);
711 CHECK_CWISE1_IF(PacketTraits::HasACos, std::acos, internal::pacos);
712
713 for (int i = 0; i < size; ++i) {
714 data1[i] = Scalar(internal::random<double>(-87, 88));
715 data2[i] = Scalar(internal::random<double>(-87, 88));
716 }
717 CHECK_CWISE1_IF(PacketTraits::HasExp, std::exp, internal::pexp);
718
719 CHECK_CWISE1_BYREF1_IF(PacketTraits::HasExp, REF_FREXP, internal::pfrexp);
720 if (PacketTraits::HasExp) {
721 // Check denormals:
722 for (int j=0; j<3; ++j) {
723 data1[0] = Scalar(std::ldexp(1, NumTraits<Scalar>::min_exponent()-j));
724 CHECK_CWISE1_BYREF1_IF(PacketTraits::HasExp, REF_FREXP, internal::pfrexp);
725 data1[0] = -data1[0];
726 CHECK_CWISE1_BYREF1_IF(PacketTraits::HasExp, REF_FREXP, internal::pfrexp);
727 }
728
729 // zero
730 data1[0] = Scalar(0);
731 CHECK_CWISE1_BYREF1_IF(PacketTraits::HasExp, REF_FREXP, internal::pfrexp);
732
733 // inf and NaN only compare output fraction, not exponent.
734 test::packet_helper<PacketTraits::HasExp,Packet> h;
735 Packet pout;
736 Scalar sout;
737 Scalar special[] = { NumTraits<Scalar>::infinity(),
738 -NumTraits<Scalar>::infinity(),
739 NumTraits<Scalar>::quiet_NaN()};
740 for (int i=0; i<3; ++i) {
741 data1[0] = special[i];
742 ref[0] = Scalar(REF_FREXP(data1[0], ref[PacketSize]));
743 h.store(data2, internal::pfrexp(h.load(data1), h.forward_reference(pout, sout)));
744 VERIFY(test::areApprox(ref, data2, 1) && "internal::pfrexp");
745 }
746 }
747
748 for (int i = 0; i < PacketSize; ++i) {
749 data1[i] = Scalar(internal::random<double>(-1, 1));
750 data2[i] = Scalar(internal::random<double>(-1, 1));
751 }
752 for (int i = 0; i < PacketSize; ++i) {
753 data1[i+PacketSize] = Scalar(internal::random<int>(-4, 4));
754 data2[i+PacketSize] = Scalar(internal::random<double>(-4, 4));
755 }
756 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp);
757 if (PacketTraits::HasExp) {
758 data1[0] = Scalar(-1);
759 // underflow to zero
760 data1[PacketSize] = Scalar(NumTraits<Scalar>::min_exponent()-55);
761 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp);
762 // overflow to inf
763 data1[PacketSize] = Scalar(NumTraits<Scalar>::max_exponent()+10);
764 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp);
765 // NaN stays NaN
766 data1[0] = NumTraits<Scalar>::quiet_NaN();
767 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp);
768 VERIFY((numext::isnan)(data2[0]));
769 // inf stays inf
770 data1[0] = NumTraits<Scalar>::infinity();
771 data1[PacketSize] = Scalar(NumTraits<Scalar>::min_exponent()-10);
772 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp);
773 // zero stays zero
774 data1[0] = Scalar(0);
775 data1[PacketSize] = Scalar(NumTraits<Scalar>::max_exponent()+10);
776 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp);
777 // Small number big exponent.
778 data1[0] = Scalar(std::ldexp(Scalar(1.0), NumTraits<Scalar>::min_exponent()-1));
779 data1[PacketSize] = Scalar(-NumTraits<Scalar>::min_exponent()
780 +NumTraits<Scalar>::max_exponent());
781 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp);
782 // Big number small exponent.
783 data1[0] = Scalar(std::ldexp(Scalar(1.0), NumTraits<Scalar>::max_exponent()-1));
784 data1[PacketSize] = Scalar(+NumTraits<Scalar>::min_exponent()
785 -NumTraits<Scalar>::max_exponent());
786 CHECK_CWISE2_IF(PacketTraits::HasExp, REF_LDEXP, internal::pldexp);
787 }
788
789 for (int i = 0; i < size; ++i) {
790 data1[i] = Scalar(internal::random<double>(-1, 1) * std::pow(10., internal::random<double>(-6, 6)));
791 data2[i] = Scalar(internal::random<double>(-1, 1) * std::pow(10., internal::random<double>(-6, 6)));
792 }
793 data1[0] = Scalar(1e-20);
794 CHECK_CWISE1_IF(PacketTraits::HasTanh, std::tanh, internal::ptanh);
795 if (PacketTraits::HasExp && PacketSize >= 2) {
796 const Scalar small = NumTraits<Scalar>::epsilon();
797 data1[0] = NumTraits<Scalar>::quiet_NaN();
798 data1[1] = small;
799 test::packet_helper<PacketTraits::HasExp, Packet> h;
800 h.store(data2, internal::pexp(h.load(data1)));
801 VERIFY((numext::isnan)(data2[0]));
802 // TODO(rmlarsen): Re-enable for bfloat16.
803 if (!internal::is_same<Scalar, bfloat16>::value) {
804 VERIFY_IS_APPROX(std::exp(small), data2[1]);
805 }
806
807 data1[0] = -small;
808 data1[1] = Scalar(0);
809 h.store(data2, internal::pexp(h.load(data1)));
810 // TODO(rmlarsen): Re-enable for bfloat16.
811 if (!internal::is_same<Scalar, bfloat16>::value) {
812 VERIFY_IS_APPROX(std::exp(-small), data2[0]);
813 }
814 VERIFY_IS_EQUAL(std::exp(Scalar(0)), data2[1]);
815
816 data1[0] = (std::numeric_limits<Scalar>::min)();
817 data1[1] = -(std::numeric_limits<Scalar>::min)();
818 h.store(data2, internal::pexp(h.load(data1)));
819 VERIFY_IS_APPROX(std::exp((std::numeric_limits<Scalar>::min)()), data2[0]);
820 VERIFY_IS_APPROX(std::exp(-(std::numeric_limits<Scalar>::min)()), data2[1]);
821
822 data1[0] = std::numeric_limits<Scalar>::denorm_min();
823 data1[1] = -std::numeric_limits<Scalar>::denorm_min();
824 h.store(data2, internal::pexp(h.load(data1)));
825 VERIFY_IS_APPROX(std::exp(std::numeric_limits<Scalar>::denorm_min()), data2[0]);
826 VERIFY_IS_APPROX(std::exp(-std::numeric_limits<Scalar>::denorm_min()), data2[1]);
827 }
828
829 if (PacketTraits::HasTanh) {
830 // NOTE this test migh fail with GCC prior to 6.3, see MathFunctionsImpl.h for details.
831 data1[0] = NumTraits<Scalar>::quiet_NaN();
832 test::packet_helper<internal::packet_traits<Scalar>::HasTanh, Packet> h;
833 h.store(data2, internal::ptanh(h.load(data1)));
834 VERIFY((numext::isnan)(data2[0]));
835 }
836
837 if (PacketTraits::HasExp) {
838 internal::scalar_logistic_op<Scalar> logistic;
839 for (int i = 0; i < size; ++i) {
840 data1[i] = Scalar(internal::random<double>(-20, 20));
841 }
842
843 test::packet_helper<PacketTraits::HasExp, Packet> h;
844 h.store(data2, logistic.packetOp(h.load(data1)));
845 for (int i = 0; i < PacketSize; ++i) {
846 VERIFY_IS_APPROX(data2[i], logistic(data1[i]));
847 }
848 }
849
850 #if EIGEN_HAS_C99_MATH && (EIGEN_COMP_CXXVER >= 11)
851 data1[0] = NumTraits<Scalar>::infinity();
852 data1[1] = Scalar(-1);
853 CHECK_CWISE1_IF(PacketTraits::HasLog1p, std::log1p, internal::plog1p);
854 data1[0] = NumTraits<Scalar>::infinity();
855 data1[1] = -NumTraits<Scalar>::infinity();
856 CHECK_CWISE1_IF(PacketTraits::HasExpm1, std::expm1, internal::pexpm1);
857 #endif
858
859 if (PacketSize >= 2) {
860 data1[0] = NumTraits<Scalar>::quiet_NaN();
861 data1[1] = NumTraits<Scalar>::epsilon();
862 if (PacketTraits::HasLog) {
863 test::packet_helper<PacketTraits::HasLog, Packet> h;
864 h.store(data2, internal::plog(h.load(data1)));
865 VERIFY((numext::isnan)(data2[0]));
866 // TODO(cantonios): Re-enable for bfloat16.
867 if (!internal::is_same<Scalar, bfloat16>::value) {
868 VERIFY_IS_APPROX(std::log(data1[1]), data2[1]);
869 }
870
871 data1[0] = -NumTraits<Scalar>::epsilon();
872 data1[1] = Scalar(0);
873 h.store(data2, internal::plog(h.load(data1)));
874 VERIFY((numext::isnan)(data2[0]));
875 VERIFY_IS_EQUAL(std::log(Scalar(0)), data2[1]);
876
877 data1[0] = (std::numeric_limits<Scalar>::min)();
878 data1[1] = -(std::numeric_limits<Scalar>::min)();
879 h.store(data2, internal::plog(h.load(data1)));
880 // TODO(cantonios): Re-enable for bfloat16.
881 if (!internal::is_same<Scalar, bfloat16>::value) {
882 VERIFY_IS_APPROX(std::log((std::numeric_limits<Scalar>::min)()), data2[0]);
883 }
884 VERIFY((numext::isnan)(data2[1]));
885
886 // Note: 32-bit arm always flushes denorms to zero.
887 #if !EIGEN_ARCH_ARM
888 if (std::numeric_limits<Scalar>::has_denorm == std::denorm_present) {
889 data1[0] = std::numeric_limits<Scalar>::denorm_min();
890 data1[1] = -std::numeric_limits<Scalar>::denorm_min();
891 h.store(data2, internal::plog(h.load(data1)));
892 // TODO(rmlarsen): Reenable.
893 // VERIFY_IS_EQUAL(std::log(std::numeric_limits<Scalar>::denorm_min()), data2[0]);
894 VERIFY((numext::isnan)(data2[1]));
895 }
896 #endif
897
898 data1[0] = Scalar(-1.0f);
899 h.store(data2, internal::plog(h.load(data1)));
900 VERIFY((numext::isnan)(data2[0]));
901
902 data1[0] = NumTraits<Scalar>::infinity();
903 h.store(data2, internal::plog(h.load(data1)));
904 VERIFY((numext::isinf)(data2[0]));
905 }
906 if (PacketTraits::HasLog1p) {
907 test::packet_helper<PacketTraits::HasLog1p, Packet> h;
908 data1[0] = Scalar(-2);
909 data1[1] = -NumTraits<Scalar>::infinity();
910 h.store(data2, internal::plog1p(h.load(data1)));
911 VERIFY((numext::isnan)(data2[0]));
912 VERIFY((numext::isnan)(data2[1]));
913 }
914 if (PacketTraits::HasSqrt) {
915 test::packet_helper<PacketTraits::HasSqrt, Packet> h;
916 data1[0] = Scalar(-1.0f);
917 if (std::numeric_limits<Scalar>::has_denorm == std::denorm_present) {
918 data1[1] = -std::numeric_limits<Scalar>::denorm_min();
919 } else {
920 data1[1] = -NumTraits<Scalar>::epsilon();
921 }
922 h.store(data2, internal::psqrt(h.load(data1)));
923 VERIFY((numext::isnan)(data2[0]));
924 VERIFY((numext::isnan)(data2[1]));
925 }
926 // TODO(rmlarsen): Re-enable for half and bfloat16.
927 if (PacketTraits::HasCos
928 && !internal::is_same<Scalar, half>::value
929 && !internal::is_same<Scalar, bfloat16>::value) {
930 test::packet_helper<PacketTraits::HasCos, Packet> h;
931 for (Scalar k = Scalar(1); k < Scalar(10000) / NumTraits<Scalar>::epsilon(); k *= Scalar(2)) {
932 for (int k1 = 0; k1 <= 1; ++k1) {
933 data1[0] = Scalar((2 * double(k) + k1) * double(EIGEN_PI) / 2 * internal::random<double>(0.8, 1.2));
934 data1[1] = Scalar((2 * double(k) + 2 + k1) * double(EIGEN_PI) / 2 * internal::random<double>(0.8, 1.2));
935 h.store(data2, internal::pcos(h.load(data1)));
936 h.store(data2 + PacketSize, internal::psin(h.load(data1)));
937 VERIFY(data2[0] <= Scalar(1.) && data2[0] >= Scalar(-1.));
938 VERIFY(data2[1] <= Scalar(1.) && data2[1] >= Scalar(-1.));
939 VERIFY(data2[PacketSize + 0] <= Scalar(1.) && data2[PacketSize + 0] >= Scalar(-1.));
940 VERIFY(data2[PacketSize + 1] <= Scalar(1.) && data2[PacketSize + 1] >= Scalar(-1.));
941
942 VERIFY_IS_APPROX(data2[0], std::cos(data1[0]));
943 VERIFY_IS_APPROX(data2[1], std::cos(data1[1]));
944 VERIFY_IS_APPROX(data2[PacketSize + 0], std::sin(data1[0]));
945 VERIFY_IS_APPROX(data2[PacketSize + 1], std::sin(data1[1]));
946
947 VERIFY_IS_APPROX(numext::abs2(data2[0]) + numext::abs2(data2[PacketSize + 0]), Scalar(1));
948 VERIFY_IS_APPROX(numext::abs2(data2[1]) + numext::abs2(data2[PacketSize + 1]), Scalar(1));
949 }
950 }
951
952 data1[0] = NumTraits<Scalar>::infinity();
953 data1[1] = -NumTraits<Scalar>::infinity();
954 h.store(data2, internal::psin(h.load(data1)));
955 VERIFY((numext::isnan)(data2[0]));
956 VERIFY((numext::isnan)(data2[1]));
957
958 h.store(data2, internal::pcos(h.load(data1)));
959 VERIFY((numext::isnan)(data2[0]));
960 VERIFY((numext::isnan)(data2[1]));
961
962 data1[0] = NumTraits<Scalar>::quiet_NaN();
963 h.store(data2, internal::psin(h.load(data1)));
964 VERIFY((numext::isnan)(data2[0]));
965 h.store(data2, internal::pcos(h.load(data1)));
966 VERIFY((numext::isnan)(data2[0]));
967
968 data1[0] = -Scalar(0.);
969 h.store(data2, internal::psin(h.load(data1)));
970 VERIFY(internal::biteq(data2[0], data1[0]));
971 h.store(data2, internal::pcos(h.load(data1)));
972 VERIFY_IS_EQUAL(data2[0], Scalar(1));
973 }
974 }
975 }
976
977 #define CAST_CHECK_CWISE1_IF(COND, REFOP, POP, SCALAR, REFTYPE) if(COND) { \
978 test::packet_helper<COND,Packet> h; \
979 for (int i=0; i<PacketSize; ++i) \
980 ref[i] = SCALAR(REFOP(static_cast<REFTYPE>(data1[i]))); \
981 h.store(data2, POP(h.load(data1))); \
982 VERIFY(test::areApprox(ref, data2, PacketSize) && #POP); \
983 }
984
985 template <typename Scalar>
propagate_nan_max(const Scalar & a,const Scalar & b)986 Scalar propagate_nan_max(const Scalar& a, const Scalar& b) {
987 if ((numext::isnan)(a)) return a;
988 if ((numext::isnan)(b)) return b;
989 return (numext::maxi)(a,b);
990 }
991
992 template <typename Scalar>
propagate_nan_min(const Scalar & a,const Scalar & b)993 Scalar propagate_nan_min(const Scalar& a, const Scalar& b) {
994 if ((numext::isnan)(a)) return a;
995 if ((numext::isnan)(b)) return b;
996 return (numext::mini)(a,b);
997 }
998
999 template <typename Scalar>
propagate_number_max(const Scalar & a,const Scalar & b)1000 Scalar propagate_number_max(const Scalar& a, const Scalar& b) {
1001 if ((numext::isnan)(a)) return b;
1002 if ((numext::isnan)(b)) return a;
1003 return (numext::maxi)(a,b);
1004 }
1005
1006 template <typename Scalar>
propagate_number_min(const Scalar & a,const Scalar & b)1007 Scalar propagate_number_min(const Scalar& a, const Scalar& b) {
1008 if ((numext::isnan)(a)) return b;
1009 if ((numext::isnan)(b)) return a;
1010 return (numext::mini)(a,b);
1011 }
1012
1013 template <typename Scalar, typename Packet>
packetmath_notcomplex()1014 void packetmath_notcomplex() {
1015 typedef internal::packet_traits<Scalar> PacketTraits;
1016 const int PacketSize = internal::unpacket_traits<Packet>::size;
1017
1018 EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4];
1019 EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4];
1020 EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4];
1021
1022 Array<Scalar, Dynamic, 1>::Map(data1, PacketSize * 4).setRandom();
1023
1024 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasMin);
1025 VERIFY((!PacketTraits::Vectorizable) || PacketTraits::HasMax);
1026
1027 CHECK_CWISE2_IF(PacketTraits::HasMin, (std::min), internal::pmin);
1028 CHECK_CWISE2_IF(PacketTraits::HasMax, (std::max), internal::pmax);
1029
1030 CHECK_CWISE2_IF(PacketTraits::HasMin, propagate_number_min, internal::pmin<PropagateNumbers>);
1031 CHECK_CWISE2_IF(PacketTraits::HasMax, propagate_number_max, internal::pmax<PropagateNumbers>);
1032 CHECK_CWISE1(numext::abs, internal::pabs);
1033 CHECK_CWISE2_IF(PacketTraits::HasAbsDiff, REF_ABS_DIFF, internal::pabsdiff);
1034
1035 ref[0] = data1[0];
1036 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmin(ref[0], data1[i]);
1037 VERIFY(internal::isApprox(ref[0], internal::predux_min(internal::pload<Packet>(data1))) && "internal::predux_min");
1038 ref[0] = data1[0];
1039 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmax(ref[0], data1[i]);
1040 VERIFY(internal::isApprox(ref[0], internal::predux_max(internal::pload<Packet>(data1))) && "internal::predux_max");
1041
1042 for (int i = 0; i < PacketSize; ++i) ref[i] = data1[0] + Scalar(i);
1043 internal::pstore(data2, internal::plset<Packet>(data1[0]));
1044 VERIFY(test::areApprox(ref, data2, PacketSize) && "internal::plset");
1045
1046 {
1047 unsigned char* data1_bits = reinterpret_cast<unsigned char*>(data1);
1048 // predux_all - not needed yet
1049 // for (unsigned int i=0; i<PacketSize*sizeof(Scalar); ++i) data1_bits[i] = 0xff;
1050 // VERIFY(internal::predux_all(internal::pload<Packet>(data1)) && "internal::predux_all(1111)");
1051 // for(int k=0; k<PacketSize; ++k)
1052 // {
1053 // for (unsigned int i=0; i<sizeof(Scalar); ++i) data1_bits[k*sizeof(Scalar)+i] = 0x0;
1054 // VERIFY( (!internal::predux_all(internal::pload<Packet>(data1))) && "internal::predux_all(0101)");
1055 // for (unsigned int i=0; i<sizeof(Scalar); ++i) data1_bits[k*sizeof(Scalar)+i] = 0xff;
1056 // }
1057
1058 // predux_any
1059 for (unsigned int i = 0; i < PacketSize * sizeof(Scalar); ++i) data1_bits[i] = 0x0;
1060 VERIFY((!internal::predux_any(internal::pload<Packet>(data1))) && "internal::predux_any(0000)");
1061 for (int k = 0; k < PacketSize; ++k) {
1062 for (unsigned int i = 0; i < sizeof(Scalar); ++i) data1_bits[k * sizeof(Scalar) + i] = 0xff;
1063 VERIFY(internal::predux_any(internal::pload<Packet>(data1)) && "internal::predux_any(0101)");
1064 for (unsigned int i = 0; i < sizeof(Scalar); ++i) data1_bits[k * sizeof(Scalar) + i] = 0x00;
1065 }
1066 }
1067
1068
1069 // Test NaN propagation.
1070 if (!NumTraits<Scalar>::IsInteger) {
1071 // Test reductions with no NaNs.
1072 ref[0] = data1[0];
1073 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmin<PropagateNumbers>(ref[0], data1[i]);
1074 VERIFY(internal::isApprox(ref[0], internal::predux_min<PropagateNumbers>(internal::pload<Packet>(data1))) && "internal::predux_min<PropagateNumbers>");
1075 ref[0] = data1[0];
1076 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmin<PropagateNaN>(ref[0], data1[i]);
1077 VERIFY(internal::isApprox(ref[0], internal::predux_min<PropagateNaN>(internal::pload<Packet>(data1))) && "internal::predux_min<PropagateNaN>");
1078 ref[0] = data1[0];
1079 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmax<PropagateNumbers>(ref[0], data1[i]);
1080 VERIFY(internal::isApprox(ref[0], internal::predux_max<PropagateNumbers>(internal::pload<Packet>(data1))) && "internal::predux_max<PropagateNumbers>");
1081 ref[0] = data1[0];
1082 for (int i = 0; i < PacketSize; ++i) ref[0] = internal::pmax<PropagateNaN>(ref[0], data1[i]);
1083 VERIFY(internal::isApprox(ref[0], internal::predux_max<PropagateNaN>(internal::pload<Packet>(data1))) && "internal::predux_max<PropagateNumbers>");
1084 // A single NaN.
1085 const size_t index = std::numeric_limits<size_t>::quiet_NaN() % PacketSize;
1086 data1[index] = NumTraits<Scalar>::quiet_NaN();
1087 VERIFY(PacketSize==1 || !(numext::isnan)(internal::predux_min<PropagateNumbers>(internal::pload<Packet>(data1))));
1088 VERIFY((numext::isnan)(internal::predux_min<PropagateNaN>(internal::pload<Packet>(data1))));
1089 VERIFY(PacketSize==1 || !(numext::isnan)(internal::predux_max<PropagateNumbers>(internal::pload<Packet>(data1))));
1090 VERIFY((numext::isnan)(internal::predux_max<PropagateNaN>(internal::pload<Packet>(data1))));
1091 // All NaNs.
1092 for (int i = 0; i < 4 * PacketSize; ++i) data1[i] = NumTraits<Scalar>::quiet_NaN();
1093 VERIFY((numext::isnan)(internal::predux_min<PropagateNumbers>(internal::pload<Packet>(data1))));
1094 VERIFY((numext::isnan)(internal::predux_min<PropagateNaN>(internal::pload<Packet>(data1))));
1095 VERIFY((numext::isnan)(internal::predux_max<PropagateNumbers>(internal::pload<Packet>(data1))));
1096 VERIFY((numext::isnan)(internal::predux_max<PropagateNaN>(internal::pload<Packet>(data1))));
1097
1098 // Test NaN propagation for coefficient-wise min and max.
1099 for (int i = 0; i < PacketSize; ++i) {
1100 data1[i] = internal::random<bool>() ? NumTraits<Scalar>::quiet_NaN() : Scalar(0);
1101 data1[i + PacketSize] = internal::random<bool>() ? NumTraits<Scalar>::quiet_NaN() : Scalar(0);
1102 }
1103 // Note: NaN propagation is implementation defined for pmin/pmax, so we do not test it here.
1104 CHECK_CWISE2_IF(PacketTraits::HasMin, propagate_number_min, (internal::pmin<PropagateNumbers>));
1105 CHECK_CWISE2_IF(PacketTraits::HasMax, propagate_number_max, internal::pmax<PropagateNumbers>);
1106 CHECK_CWISE2_IF(PacketTraits::HasMin, propagate_nan_min, (internal::pmin<PropagateNaN>));
1107 CHECK_CWISE2_IF(PacketTraits::HasMax, propagate_nan_max, internal::pmax<PropagateNaN>);
1108 }
1109
1110 packetmath_boolean_mask_ops_notcomplex<Scalar, Packet>();
1111 }
1112
1113 template <typename Scalar, typename Packet, bool ConjLhs, bool ConjRhs>
test_conj_helper(Scalar * data1,Scalar * data2,Scalar * ref,Scalar * pval)1114 void test_conj_helper(Scalar* data1, Scalar* data2, Scalar* ref, Scalar* pval) {
1115 const int PacketSize = internal::unpacket_traits<Packet>::size;
1116
1117 internal::conj_if<ConjLhs> cj0;
1118 internal::conj_if<ConjRhs> cj1;
1119 internal::conj_helper<Scalar, Scalar, ConjLhs, ConjRhs> cj;
1120 internal::conj_helper<Packet, Packet, ConjLhs, ConjRhs> pcj;
1121
1122 for (int i = 0; i < PacketSize; ++i) {
1123 ref[i] = cj0(data1[i]) * cj1(data2[i]);
1124 VERIFY(internal::isApprox(ref[i], cj.pmul(data1[i], data2[i])) && "conj_helper pmul");
1125 }
1126 internal::pstore(pval, pcj.pmul(internal::pload<Packet>(data1), internal::pload<Packet>(data2)));
1127 VERIFY(test::areApprox(ref, pval, PacketSize) && "conj_helper pmul");
1128
1129 for (int i = 0; i < PacketSize; ++i) {
1130 Scalar tmp = ref[i];
1131 ref[i] += cj0(data1[i]) * cj1(data2[i]);
1132 VERIFY(internal::isApprox(ref[i], cj.pmadd(data1[i], data2[i], tmp)) && "conj_helper pmadd");
1133 }
1134 internal::pstore(
1135 pval, pcj.pmadd(internal::pload<Packet>(data1), internal::pload<Packet>(data2), internal::pload<Packet>(pval)));
1136 VERIFY(test::areApprox(ref, pval, PacketSize) && "conj_helper pmadd");
1137 }
1138
1139 template <typename Scalar, typename Packet>
packetmath_complex()1140 void packetmath_complex() {
1141 typedef internal::packet_traits<Scalar> PacketTraits;
1142 typedef typename Scalar::value_type RealScalar;
1143 const int PacketSize = internal::unpacket_traits<Packet>::size;
1144
1145 const int size = PacketSize * 4;
1146 EIGEN_ALIGN_MAX Scalar data1[PacketSize * 4];
1147 EIGEN_ALIGN_MAX Scalar data2[PacketSize * 4];
1148 EIGEN_ALIGN_MAX Scalar ref[PacketSize * 4];
1149 EIGEN_ALIGN_MAX Scalar pval[PacketSize * 4];
1150
1151 for (int i = 0; i < size; ++i) {
1152 data1[i] = internal::random<Scalar>() * Scalar(1e2);
1153 data2[i] = internal::random<Scalar>() * Scalar(1e2);
1154 }
1155
1156 test_conj_helper<Scalar, Packet, false, false>(data1, data2, ref, pval);
1157 test_conj_helper<Scalar, Packet, false, true>(data1, data2, ref, pval);
1158 test_conj_helper<Scalar, Packet, true, false>(data1, data2, ref, pval);
1159 test_conj_helper<Scalar, Packet, true, true>(data1, data2, ref, pval);
1160
1161 // Test pcplxflip.
1162 {
1163 for (int i = 0; i < PacketSize; ++i) ref[i] = Scalar(std::imag(data1[i]), std::real(data1[i]));
1164 internal::pstore(pval, internal::pcplxflip(internal::pload<Packet>(data1)));
1165 VERIFY(test::areApprox(ref, pval, PacketSize) && "pcplxflip");
1166 }
1167
1168 if (PacketTraits::HasSqrt) {
1169 for (int i = 0; i < size; ++i) {
1170 data1[i] = Scalar(internal::random<RealScalar>(), internal::random<RealScalar>());
1171 }
1172 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, size);
1173
1174 // Test misc. corner cases.
1175 const RealScalar zero = RealScalar(0);
1176 const RealScalar one = RealScalar(1);
1177 const RealScalar inf = std::numeric_limits<RealScalar>::infinity();
1178 const RealScalar nan = std::numeric_limits<RealScalar>::quiet_NaN();
1179 data1[0] = Scalar(zero, zero);
1180 data1[1] = Scalar(-zero, zero);
1181 data1[2] = Scalar(one, zero);
1182 data1[3] = Scalar(zero, one);
1183 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4);
1184 data1[0] = Scalar(-one, zero);
1185 data1[1] = Scalar(zero, -one);
1186 data1[2] = Scalar(one, one);
1187 data1[3] = Scalar(-one, -one);
1188 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4);
1189 data1[0] = Scalar(inf, zero);
1190 data1[1] = Scalar(zero, inf);
1191 data1[2] = Scalar(-inf, zero);
1192 data1[3] = Scalar(zero, -inf);
1193 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4);
1194 data1[0] = Scalar(inf, inf);
1195 data1[1] = Scalar(-inf, inf);
1196 data1[2] = Scalar(inf, -inf);
1197 data1[3] = Scalar(-inf, -inf);
1198 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4);
1199 data1[0] = Scalar(nan, zero);
1200 data1[1] = Scalar(zero, nan);
1201 data1[2] = Scalar(nan, one);
1202 data1[3] = Scalar(one, nan);
1203 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4);
1204 data1[0] = Scalar(nan, nan);
1205 data1[1] = Scalar(inf, nan);
1206 data1[2] = Scalar(nan, inf);
1207 data1[3] = Scalar(-inf, nan);
1208 CHECK_CWISE1_N(numext::sqrt, internal::psqrt, 4);
1209 }
1210 }
1211
1212 template <typename Scalar, typename Packet>
packetmath_scatter_gather()1213 void packetmath_scatter_gather() {
1214 typedef typename NumTraits<Scalar>::Real RealScalar;
1215 const int PacketSize = internal::unpacket_traits<Packet>::size;
1216 EIGEN_ALIGN_MAX Scalar data1[PacketSize];
1217 RealScalar refvalue = RealScalar(0);
1218 for (int i = 0; i < PacketSize; ++i) {
1219 data1[i] = internal::random<Scalar>() / RealScalar(PacketSize);
1220 }
1221
1222 int stride = internal::random<int>(1, 20);
1223
1224 // Buffer of zeros.
1225 EIGEN_ALIGN_MAX Scalar buffer[PacketSize * 20] = {};
1226
1227 Packet packet = internal::pload<Packet>(data1);
1228 internal::pscatter<Scalar, Packet>(buffer, packet, stride);
1229
1230 for (int i = 0; i < PacketSize * 20; ++i) {
1231 if ((i % stride) == 0 && i < stride * PacketSize) {
1232 VERIFY(test::isApproxAbs(buffer[i], data1[i / stride], refvalue) && "pscatter");
1233 } else {
1234 VERIFY(test::isApproxAbs(buffer[i], Scalar(0), refvalue) && "pscatter");
1235 }
1236 }
1237
1238 for (int i = 0; i < PacketSize * 7; ++i) {
1239 buffer[i] = internal::random<Scalar>() / RealScalar(PacketSize);
1240 }
1241 packet = internal::pgather<Scalar, Packet>(buffer, 7);
1242 internal::pstore(data1, packet);
1243 for (int i = 0; i < PacketSize; ++i) {
1244 VERIFY(test::isApproxAbs(data1[i], buffer[i * 7], refvalue) && "pgather");
1245 }
1246 }
1247
1248 namespace Eigen {
1249 namespace test {
1250
1251 template <typename Scalar, typename PacketType>
1252 struct runall<Scalar, PacketType, false, false> { // i.e. float or double
runEigen::test::runall1253 static void run() {
1254 packetmath<Scalar, PacketType>();
1255 packetmath_scatter_gather<Scalar, PacketType>();
1256 packetmath_notcomplex<Scalar, PacketType>();
1257 packetmath_real<Scalar, PacketType>();
1258 }
1259 };
1260
1261 template <typename Scalar, typename PacketType>
1262 struct runall<Scalar, PacketType, false, true> { // i.e. int
runEigen::test::runall1263 static void run() {
1264 packetmath<Scalar, PacketType>();
1265 packetmath_scatter_gather<Scalar, PacketType>();
1266 packetmath_notcomplex<Scalar, PacketType>();
1267 }
1268 };
1269
1270 template <typename Scalar, typename PacketType>
1271 struct runall<Scalar, PacketType, true, false> { // i.e. complex
runEigen::test::runall1272 static void run() {
1273 packetmath<Scalar, PacketType>();
1274 packetmath_scatter_gather<Scalar, PacketType>();
1275 packetmath_complex<Scalar, PacketType>();
1276 }
1277 };
1278
1279 } // namespace test
1280 } // namespace Eigen
1281
EIGEN_DECLARE_TEST(packetmath)1282 EIGEN_DECLARE_TEST(packetmath) {
1283 g_first_pass = true;
1284 for (int i = 0; i < g_repeat; i++) {
1285 CALL_SUBTEST_1(test::runner<float>::run());
1286 CALL_SUBTEST_2(test::runner<double>::run());
1287 CALL_SUBTEST_3(test::runner<int8_t>::run());
1288 CALL_SUBTEST_4(test::runner<uint8_t>::run());
1289 CALL_SUBTEST_5(test::runner<int16_t>::run());
1290 CALL_SUBTEST_6(test::runner<uint16_t>::run());
1291 CALL_SUBTEST_7(test::runner<int32_t>::run());
1292 CALL_SUBTEST_8(test::runner<uint32_t>::run());
1293 CALL_SUBTEST_9(test::runner<int64_t>::run());
1294 CALL_SUBTEST_10(test::runner<uint64_t>::run());
1295 CALL_SUBTEST_11(test::runner<std::complex<float> >::run());
1296 CALL_SUBTEST_12(test::runner<std::complex<double> >::run());
1297 CALL_SUBTEST_13(test::runner<half>::run());
1298 CALL_SUBTEST_14((packetmath<bool, internal::packet_traits<bool>::type>()));
1299 CALL_SUBTEST_15(test::runner<bfloat16>::run());
1300 g_first_pass = false;
1301 }
1302 }
1303