special_packetmath.cpp (6372B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com> 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 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 { 130 static void run() { 131 packetmath_real<Scalar,PacketType>(); 132 } 133 }; 134 135 } 136 } 137 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 }