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 <limits>
12 #include "packetmath_test_shared.h"
13 #include "../Eigen/SpecialFunctions"
14
packetmath_real()15 template<typename Scalar,typename Packet> void packetmath_real()
16 {
17 using std::abs;
18 typedef internal::packet_traits<Scalar> PacketTraits;
19 const int PacketSize = internal::unpacket_traits<Packet>::size;
20
21 const int size = PacketSize*4;
22 EIGEN_ALIGN_MAX Scalar data1[PacketSize*4];
23 EIGEN_ALIGN_MAX Scalar data2[PacketSize*4];
24 EIGEN_ALIGN_MAX Scalar ref[PacketSize*4];
25
26 #if EIGEN_HAS_C99_MATH
27 {
28 data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
29 test::packet_helper<internal::packet_traits<Scalar>::HasLGamma,Packet> h;
30 h.store(data2, internal::plgamma(h.load(data1)));
31 VERIFY((numext::isnan)(data2[0]));
32 }
33 if (internal::packet_traits<Scalar>::HasErf) {
34 data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
35 test::packet_helper<internal::packet_traits<Scalar>::HasErf,Packet> h;
36 h.store(data2, internal::perf(h.load(data1)));
37 VERIFY((numext::isnan)(data2[0]));
38 }
39 {
40 data1[0] = std::numeric_limits<Scalar>::quiet_NaN();
41 test::packet_helper<internal::packet_traits<Scalar>::HasErfc,Packet> h;
42 h.store(data2, internal::perfc(h.load(data1)));
43 VERIFY((numext::isnan)(data2[0]));
44 }
45 {
46 for (int i=0; i<size; ++i) {
47 data1[i] = internal::random<Scalar>(Scalar(0),Scalar(1));
48 }
49 CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasNdtri, numext::ndtri, internal::pndtri);
50 }
51 #endif // EIGEN_HAS_C99_MATH
52
53 // For bessel_i*e and bessel_j*, the valid range is negative reals.
54 {
55 const int max_exponent = numext::mini(std::numeric_limits<Scalar>::max_exponent10-1, 6);
56 for (int i=0; i<size; ++i)
57 {
58 data1[i] = internal::random<Scalar>(Scalar(-1),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-max_exponent),Scalar(max_exponent))));
59 data2[i] = internal::random<Scalar>(Scalar(-1),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-max_exponent),Scalar(max_exponent))));
60 }
61
62 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0e, internal::pbessel_i0e);
63 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1e, internal::pbessel_i1e);
64 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j0, internal::pbessel_j0);
65 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_j1, internal::pbessel_j1);
66 }
67
68 // Use a smaller data range for the bessel_i* as these can become very large.
69 // Following #1693, we also restrict this range further to avoid inf's due to
70 // differences in pexp and exp.
71 for (int i=0; i<size; ++i) {
72 data1[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
73 Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1),Scalar(2))));
74 data2[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
75 Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1),Scalar(2))));
76 }
77 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i0, internal::pbessel_i0);
78 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_i1, internal::pbessel_i1);
79
80
81 // y_i, and k_i are valid for x > 0.
82 {
83 const int max_exponent = numext::mini(std::numeric_limits<Scalar>::max_exponent10-1, 5);
84 for (int i=0; i<size; ++i)
85 {
86 data1[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-2),Scalar(max_exponent))));
87 data2[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) * Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-2),Scalar(max_exponent))));
88 }
89 }
90
91 // TODO(srvasude): Re-enable this test once properly investigated why the
92 // scalar and vector paths differ.
93 // CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y0, internal::pbessel_y0);
94 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_y1, internal::pbessel_y1);
95 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0e, internal::pbessel_k0e);
96 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1e, internal::pbessel_k1e);
97
98 // Following #1693, we restrict the range for exp to avoid zeroing out too
99 // fast.
100 for (int i=0; i<size; ++i) {
101 data1[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
102 Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1),Scalar(2))));
103 data2[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
104 Scalar(std::pow(Scalar(9), internal::random<Scalar>(Scalar(-1),Scalar(2))));
105 }
106 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k0, internal::pbessel_k0);
107 CHECK_CWISE1_IF(PacketTraits::HasBessel, numext::bessel_k1, internal::pbessel_k1);
108
109
110 for (int i=0; i<size; ++i) {
111 data1[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
112 Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-1),Scalar(2))));
113 data2[i] = internal::random<Scalar>(Scalar(0.01),Scalar(1)) *
114 Scalar(std::pow(Scalar(10), internal::random<Scalar>(Scalar(-1),Scalar(2))));
115 }
116
117 #if EIGEN_HAS_C99_MATH && (EIGEN_COMP_CXXVER >= 11)
118 CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasLGamma, std::lgamma, internal::plgamma);
119 CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErf, std::erf, internal::perf);
120 CHECK_CWISE1_IF(internal::packet_traits<Scalar>::HasErfc, std::erfc, internal::perfc);
121 #endif
122
123 }
124
125 namespace Eigen {
126 namespace test {
127
128 template<typename Scalar,typename PacketType, bool IsComplex, bool IsInteger>
129 struct runall {
runEigen::test::runall130 static void run() {
131 packetmath_real<Scalar,PacketType>();
132 }
133 };
134
135 }
136 }
137
EIGEN_DECLARE_TEST(special_packetmath)138 EIGEN_DECLARE_TEST(special_packetmath)
139 {
140 g_first_pass = true;
141 for(int i = 0; i < g_repeat; i++) {
142
143 CALL_SUBTEST_1( test::runner<float>::run() );
144 CALL_SUBTEST_2( test::runner<double>::run() );
145 CALL_SUBTEST_3( test::runner<Eigen::half>::run() );
146 CALL_SUBTEST_4( test::runner<Eigen::bfloat16>::run() );
147 g_first_pass = false;
148 }
149 }
150