special_packetmath.cpp 6.2 KB

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