zmath.h 6.5 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251
  1. #pragma once
  2. // Complex number math operations that act as no-ops for other dtypes.
  3. #include <c10/util/complex.h>
  4. #include <c10/util/math_compat.h>
  5. #include <c10/util/MathConstants.h>
  6. #include<ATen/NumericUtils.h>
  7. namespace at { namespace native {
  8. inline namespace CPU_CAPABILITY {
  9. template <typename SCALAR_TYPE, typename VALUE_TYPE=SCALAR_TYPE>
  10. inline VALUE_TYPE zabs (SCALAR_TYPE z) {
  11. return z;
  12. }
  13. template<>
  14. inline c10::complex<float> zabs <c10::complex<float>> (c10::complex<float> z) {
  15. return c10::complex<float>(std::abs(z));
  16. }
  17. template<>
  18. inline float zabs <c10::complex<float>, float> (c10::complex<float> z) {
  19. return std::abs(z);
  20. }
  21. template<>
  22. inline c10::complex<double> zabs <c10::complex<double>> (c10::complex<double> z) {
  23. return c10::complex<double>(std::abs(z));
  24. }
  25. template<>
  26. inline double zabs <c10::complex<double>, double> (c10::complex<double> z) {
  27. return std::abs(z);
  28. }
  29. // This overload corresponds to non-complex dtypes.
  30. // The function is consistent with its NumPy equivalent
  31. // for non-complex dtypes where `pi` is returned for
  32. // negative real numbers and `0` is returned for 0 or positive
  33. // real numbers.
  34. // Note: `nan` is propagated.
  35. template <typename SCALAR_TYPE, typename VALUE_TYPE=SCALAR_TYPE>
  36. inline VALUE_TYPE angle_impl (SCALAR_TYPE z) {
  37. if (at::_isnan(z)) {
  38. return z;
  39. }
  40. return z < 0 ? c10::pi<double> : 0;
  41. }
  42. template<>
  43. inline c10::complex<float> angle_impl <c10::complex<float>> (c10::complex<float> z) {
  44. return c10::complex<float>(std::arg(z), 0.0);
  45. }
  46. template<>
  47. inline float angle_impl <c10::complex<float>, float> (c10::complex<float> z) {
  48. return std::arg(z);
  49. }
  50. template<>
  51. inline c10::complex<double> angle_impl <c10::complex<double>> (c10::complex<double> z) {
  52. return c10::complex<double>(std::arg(z), 0.0);
  53. }
  54. template<>
  55. inline double angle_impl <c10::complex<double>, double> (c10::complex<double> z) {
  56. return std::arg(z);
  57. }
  58. template <typename SCALAR_TYPE, typename VALUE_TYPE=SCALAR_TYPE>
  59. constexpr VALUE_TYPE real_impl (SCALAR_TYPE z) {
  60. return z; //No-Op
  61. }
  62. template<>
  63. constexpr c10::complex<float> real_impl <c10::complex<float>> (c10::complex<float> z) {
  64. return c10::complex<float>(z.real(), 0.0);
  65. }
  66. template<>
  67. constexpr float real_impl <c10::complex<float>, float> (c10::complex<float> z) {
  68. return z.real();
  69. }
  70. template<>
  71. constexpr c10::complex<double> real_impl <c10::complex<double>> (c10::complex<double> z) {
  72. return c10::complex<double>(z.real(), 0.0);
  73. }
  74. template<>
  75. constexpr double real_impl <c10::complex<double>, double> (c10::complex<double> z) {
  76. return z.real();
  77. }
  78. template <typename SCALAR_TYPE, typename VALUE_TYPE=SCALAR_TYPE>
  79. constexpr VALUE_TYPE imag_impl (SCALAR_TYPE /*z*/) {
  80. return 0;
  81. }
  82. template<>
  83. constexpr c10::complex<float> imag_impl <c10::complex<float>> (c10::complex<float> z) {
  84. return c10::complex<float>(z.imag(), 0.0);
  85. }
  86. template<>
  87. constexpr float imag_impl <c10::complex<float>, float> (c10::complex<float> z) {
  88. return z.imag();
  89. }
  90. template<>
  91. constexpr c10::complex<double> imag_impl <c10::complex<double>> (c10::complex<double> z) {
  92. return c10::complex<double>(z.imag(), 0.0);
  93. }
  94. template<>
  95. constexpr double imag_impl <c10::complex<double>, double> (c10::complex<double> z) {
  96. return z.imag();
  97. }
  98. template <typename TYPE>
  99. inline TYPE conj_impl (TYPE z) {
  100. return z; //No-Op
  101. }
  102. template<>
  103. inline c10::complex<at::Half> conj_impl <c10::complex<at::Half>> (c10::complex<at::Half> z) {
  104. return c10::complex<at::Half>{z.real(), -z.imag()};
  105. }
  106. template<>
  107. inline c10::complex<float> conj_impl <c10::complex<float>> (c10::complex<float> z) {
  108. return c10::complex<float>(z.real(), -z.imag());
  109. }
  110. template<>
  111. inline c10::complex<double> conj_impl <c10::complex<double>> (c10::complex<double> z) {
  112. return c10::complex<double>(z.real(), -z.imag());
  113. }
  114. template <typename TYPE>
  115. inline TYPE ceil_impl (TYPE z) {
  116. return std::ceil(z);
  117. }
  118. template <>
  119. inline c10::complex<float> ceil_impl (c10::complex<float> z) {
  120. return c10::complex<float>(std::ceil(z.real()), std::ceil(z.imag()));
  121. }
  122. template <>
  123. inline c10::complex<double> ceil_impl (c10::complex<double> z) {
  124. return c10::complex<double>(std::ceil(z.real()), std::ceil(z.imag()));
  125. }
  126. template<typename T>
  127. inline c10::complex<T> sgn_impl (c10::complex<T> z) {
  128. if (z == c10::complex<T>(0, 0)) {
  129. return c10::complex<T>(0, 0);
  130. } else {
  131. return z / zabs(z);
  132. }
  133. }
  134. template <typename TYPE>
  135. inline TYPE floor_impl (TYPE z) {
  136. return std::floor(z);
  137. }
  138. template <>
  139. inline c10::complex<float> floor_impl (c10::complex<float> z) {
  140. return c10::complex<float>(std::floor(z.real()), std::floor(z.imag()));
  141. }
  142. template <>
  143. inline c10::complex<double> floor_impl (c10::complex<double> z) {
  144. return c10::complex<double>(std::floor(z.real()), std::floor(z.imag()));
  145. }
  146. template <typename TYPE>
  147. inline TYPE round_impl (TYPE z) {
  148. return std::nearbyint(z);
  149. }
  150. template <>
  151. inline c10::complex<float> round_impl (c10::complex<float> z) {
  152. return c10::complex<float>(std::nearbyint(z.real()), std::nearbyint(z.imag()));
  153. }
  154. template <>
  155. inline c10::complex<double> round_impl (c10::complex<double> z) {
  156. return c10::complex<double>(std::nearbyint(z.real()), std::nearbyint(z.imag()));
  157. }
  158. template <typename TYPE>
  159. inline TYPE trunc_impl (TYPE z) {
  160. return std::trunc(z);
  161. }
  162. template <>
  163. inline c10::complex<float> trunc_impl (c10::complex<float> z) {
  164. return c10::complex<float>(std::trunc(z.real()), std::trunc(z.imag()));
  165. }
  166. template <>
  167. inline c10::complex<double> trunc_impl (c10::complex<double> z) {
  168. return c10::complex<double>(std::trunc(z.real()), std::trunc(z.imag()));
  169. }
  170. template <typename TYPE, std::enable_if_t<!c10::is_complex<TYPE>::value, int> = 0>
  171. inline TYPE max_impl (TYPE a, TYPE b) {
  172. if (_isnan<TYPE>(a) || _isnan<TYPE>(b)) {
  173. return std::numeric_limits<TYPE>::quiet_NaN();
  174. } else {
  175. return std::max(a, b);
  176. }
  177. }
  178. template <typename TYPE, std::enable_if_t<c10::is_complex<TYPE>::value, int> = 0>
  179. inline TYPE max_impl (TYPE a, TYPE b) {
  180. if (_isnan<TYPE>(a)) {
  181. return a;
  182. } else if (_isnan<TYPE>(b)) {
  183. return b;
  184. } else {
  185. return std::abs(a) > std::abs(b) ? a : b;
  186. }
  187. }
  188. template <typename TYPE, std::enable_if_t<!c10::is_complex<TYPE>::value, int> = 0>
  189. inline TYPE min_impl (TYPE a, TYPE b) {
  190. if (_isnan<TYPE>(a) || _isnan<TYPE>(b)) {
  191. return std::numeric_limits<TYPE>::quiet_NaN();
  192. } else {
  193. return std::min(a, b);
  194. }
  195. }
  196. template <typename TYPE, std::enable_if_t<c10::is_complex<TYPE>::value, int> = 0>
  197. inline TYPE min_impl (TYPE a, TYPE b) {
  198. if (_isnan<TYPE>(a)) {
  199. return a;
  200. } else if (_isnan<TYPE>(b)) {
  201. return b;
  202. } else {
  203. return std::abs(a) < std::abs(b) ? a : b;
  204. }
  205. }
  206. } // end namespace
  207. }} //end at::native