eop_aux.hpp 21 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286
  1. // Copyright 2008-2016 Conrad Sanderson (http://conradsanderson.id.au)
  2. // Copyright 2008-2016 National ICT Australia (NICTA)
  3. //
  4. // Licensed under the Apache License, Version 2.0 (the "License");
  5. // you may not use this file except in compliance with the License.
  6. // You may obtain a copy of the License at
  7. // http://www.apache.org/licenses/LICENSE-2.0
  8. //
  9. // Unless required by applicable law or agreed to in writing, software
  10. // distributed under the License is distributed on an "AS IS" BASIS,
  11. // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  12. // See the License for the specific language governing permissions and
  13. // limitations under the License.
  14. // ------------------------------------------------------------------------
  15. //! \addtogroup eop_aux
  16. //! @{
  17. //! use of the SFINAE approach to work around compiler limitations
  18. //! http://en.wikipedia.org/wiki/SFINAE
  19. class eop_aux
  20. {
  21. public:
  22. template<typename eT> arma_inline static typename arma_integral_only<eT>::result acos (const eT x) { return eT( std::acos(double(x)) ); }
  23. template<typename eT> arma_inline static typename arma_integral_only<eT>::result asin (const eT x) { return eT( std::asin(double(x)) ); }
  24. template<typename eT> arma_inline static typename arma_integral_only<eT>::result atan (const eT x) { return eT( std::atan(double(x)) ); }
  25. template<typename eT> arma_inline static typename arma_real_only<eT>::result acos (const eT x) { return std::acos(x); }
  26. template<typename eT> arma_inline static typename arma_real_only<eT>::result asin (const eT x) { return std::asin(x); }
  27. template<typename eT> arma_inline static typename arma_real_only<eT>::result atan (const eT x) { return std::atan(x); }
  28. template<typename eT> arma_inline static typename arma_cx_only<eT>::result acos (const eT x) { return arma_acos(x); }
  29. template<typename eT> arma_inline static typename arma_cx_only<eT>::result asin (const eT x) { return arma_asin(x); }
  30. template<typename eT> arma_inline static typename arma_cx_only<eT>::result atan (const eT x) { return arma_atan(x); }
  31. template<typename eT> arma_inline static typename arma_integral_only<eT>::result acosh (const eT x) { return eT( arma_acosh(double(x)) ); }
  32. template<typename eT> arma_inline static typename arma_integral_only<eT>::result asinh (const eT x) { return eT( arma_asinh(double(x)) ); }
  33. template<typename eT> arma_inline static typename arma_integral_only<eT>::result atanh (const eT x) { return eT( arma_atanh(double(x)) ); }
  34. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result acosh (const eT x) { return arma_acosh(x); }
  35. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result asinh (const eT x) { return arma_asinh(x); }
  36. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result atanh (const eT x) { return arma_atanh(x); }
  37. template<typename eT> arma_inline static typename arma_not_cx<eT>::result conj(const eT x) { return x; }
  38. template<typename T> arma_inline static std::complex<T> conj(const std::complex<T>& x) { return std::conj(x); }
  39. template<typename eT> arma_inline static typename arma_integral_only<eT>::result sqrt (const eT x) { return eT( std::sqrt (double(x)) ); }
  40. template<typename eT> arma_inline static typename arma_integral_only<eT>::result log10 (const eT x) { return eT( std::log10(double(x)) ); }
  41. template<typename eT> arma_inline static typename arma_integral_only<eT>::result log (const eT x) { return eT( std::log (double(x)) ); }
  42. template<typename eT> arma_inline static typename arma_integral_only<eT>::result exp (const eT x) { return eT( std::exp (double(x)) ); }
  43. template<typename eT> arma_inline static typename arma_integral_only<eT>::result cos (const eT x) { return eT( std::cos (double(x)) ); }
  44. template<typename eT> arma_inline static typename arma_integral_only<eT>::result sin (const eT x) { return eT( std::sin (double(x)) ); }
  45. template<typename eT> arma_inline static typename arma_integral_only<eT>::result tan (const eT x) { return eT( std::tan (double(x)) ); }
  46. template<typename eT> arma_inline static typename arma_integral_only<eT>::result cosh (const eT x) { return eT( std::cosh (double(x)) ); }
  47. template<typename eT> arma_inline static typename arma_integral_only<eT>::result sinh (const eT x) { return eT( std::sinh (double(x)) ); }
  48. template<typename eT> arma_inline static typename arma_integral_only<eT>::result tanh (const eT x) { return eT( std::tanh (double(x)) ); }
  49. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result sqrt (const eT x) { return std::sqrt (x); }
  50. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result log10 (const eT x) { return std::log10(x); }
  51. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result log (const eT x) { return std::log (x); }
  52. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result exp (const eT x) { return std::exp (x); }
  53. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result cos (const eT x) { return std::cos (x); }
  54. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result sin (const eT x) { return std::sin (x); }
  55. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result tan (const eT x) { return std::tan (x); }
  56. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result cosh (const eT x) { return std::cosh (x); }
  57. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result sinh (const eT x) { return std::sinh (x); }
  58. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result tanh (const eT x) { return std::tanh (x); }
  59. template<typename eT> arma_inline static typename arma_unsigned_integral_only<eT>::result neg (const eT x) { return x; }
  60. template<typename eT> arma_inline static typename arma_signed_only<eT>::result neg (const eT x) { return -x; }
  61. template<typename eT> arma_inline static typename arma_integral_only<eT>::result floor (const eT x) { return x; }
  62. template<typename eT> arma_inline static typename arma_real_only<eT>::result floor (const eT x) { return std::floor(x); }
  63. template<typename eT> arma_inline static typename arma_cx_only<eT>::result floor (const eT& x) { return eT( std::floor(x.real()), std::floor(x.imag()) ); }
  64. template<typename eT> arma_inline static typename arma_integral_only<eT>::result ceil (const eT x) { return x; }
  65. template<typename eT> arma_inline static typename arma_real_only<eT>::result ceil (const eT x) { return std::ceil(x); }
  66. template<typename eT> arma_inline static typename arma_cx_only<eT>::result ceil (const eT& x) { return eT( std::ceil(x.real()), std::ceil(x.imag()) ); }
  67. #if defined(ARMA_USE_CXX11)
  68. template<typename eT> arma_inline static typename arma_integral_only<eT>::result round (const eT x) { return x; }
  69. template<typename eT> arma_inline static typename arma_real_only<eT>::result round (const eT x) { return std::round(x); }
  70. template<typename eT> arma_inline static typename arma_cx_only<eT>::result round (const eT& x) { return eT( std::round(x.real()), std::round(x.imag()) ); }
  71. #else
  72. template<typename eT> arma_inline static typename arma_integral_only<eT>::result round (const eT x) { return x; }
  73. template<typename eT> arma_inline static typename arma_real_only<eT>::result round (const eT x) { return (x >= eT(0)) ? std::floor(x+0.5) : std::ceil(x-0.5); }
  74. template<typename eT> arma_inline static typename arma_cx_only<eT>::result round (const eT& x) { return eT( eop_aux::round(x.real()), eop_aux::round(x.imag()) ); }
  75. #endif
  76. #if defined(ARMA_USE_CXX11)
  77. template<typename eT> arma_inline static typename arma_integral_only<eT>::result trunc (const eT x) { return x; }
  78. template<typename eT> arma_inline static typename arma_real_only<eT>::result trunc (const eT x) { return std::trunc(x); }
  79. template<typename eT> arma_inline static typename arma_cx_only<eT>::result trunc (const eT& x) { return eT( std::trunc(x.real()), std::trunc(x.imag()) ); }
  80. #else
  81. template<typename eT> arma_inline static typename arma_integral_only<eT>::result trunc (const eT x) { return x; }
  82. template<typename eT> arma_inline static typename arma_real_only<eT>::result trunc (const eT x) { return (x >= eT(0)) ? std::floor(x) : std::ceil(x); }
  83. template<typename eT> arma_inline static typename arma_cx_only<eT>::result trunc (const eT& x) { return eT( eop_aux::trunc(x.real()), eop_aux::trunc(x.imag()) ); }
  84. #endif
  85. #if defined(ARMA_USE_CXX11)
  86. template<typename eT> arma_inline static typename arma_integral_only<eT>::result log2 (const eT x) { return eT( std::log(double(x))/ double(0.69314718055994530942) ); }
  87. template<typename eT> arma_inline static typename arma_real_only<eT>::result log2 (const eT x) { return std::log2(x); }
  88. template<typename eT> arma_inline static typename arma_cx_only<eT>::result log2 (const eT& x) { typedef typename get_pod_type<eT>::result T; return std::log(x) / T(0.69314718055994530942); }
  89. #else
  90. template<typename eT> arma_inline static typename arma_integral_only<eT>::result log2 (const eT x) { return eT( std::log(double(x))/ double(0.69314718055994530942) ); }
  91. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result log2 (const eT x) { typedef typename get_pod_type<eT>::result T; return std::log(x) / T(0.69314718055994530942); }
  92. #endif
  93. #if defined(ARMA_USE_CXX11)
  94. template<typename eT> arma_inline static typename arma_integral_only<eT>::result log1p (const eT x) { return eT( std::log1p(double(x)) ); }
  95. template<typename eT> arma_inline static typename arma_real_only<eT>::result log1p (const eT x) { return std::log1p(x); }
  96. template<typename eT> arma_inline static typename arma_cx_only<eT>::result log1p (const eT& x) { arma_ignore(x); return eT(0); }
  97. #elif defined(ARMA_HAVE_TR1)
  98. template<typename eT> arma_inline static typename arma_integral_only<eT>::result log1p (const eT x) { return eT( std::tr1::log1p(double(x)) ); }
  99. template<typename eT> arma_inline static typename arma_real_only<eT>::result log1p (const eT x) { return std::tr1::log1p(x); }
  100. template<typename eT> arma_inline static typename arma_cx_only<eT>::result log1p (const eT& x) { arma_ignore(x); return eT(0); }
  101. #else
  102. template<typename eT> arma_inline static eT log1p (const eT x) { arma_ignore(x); arma_stop_logic_error("log1p(): C++11 compiler required"); return eT(0); }
  103. #endif
  104. #if defined(ARMA_USE_CXX11)
  105. template<typename eT> arma_inline static typename arma_integral_only<eT>::result exp2 (const eT x) { return eT( std::pow(double(2), double(x)) ); }
  106. template<typename eT> arma_inline static typename arma_real_only<eT>::result exp2 (const eT x) { return std::exp2(x); }
  107. template<typename eT> arma_inline static typename arma_cx_only<eT>::result exp2 (const eT& x) { typedef typename get_pod_type<eT>::result T; return std::pow( T(2), x); }
  108. #else
  109. template<typename eT> arma_inline static typename arma_integral_only<eT>::result exp2 (const eT x) { return eT( std::pow(double(2), double(x)) ); }
  110. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result exp2 (const eT x) { typedef typename get_pod_type<eT>::result T; return std::pow( T(2), x); }
  111. #endif
  112. template<typename eT> arma_inline static typename arma_integral_only<eT>::result exp10 (const eT x) { return eT( std::pow(double(10), double(x)) ); }
  113. template<typename eT> arma_inline static typename arma_real_or_cx_only<eT>::result exp10 (const eT x) { typedef typename get_pod_type<eT>::result T; return std::pow( T(10), x); }
  114. #if defined(ARMA_USE_CXX11)
  115. template<typename eT> arma_inline static typename arma_integral_only<eT>::result expm1 (const eT x) { return eT( std::expm1(double(x)) ); }
  116. template<typename eT> arma_inline static typename arma_real_only<eT>::result expm1 (const eT x) { return std::expm1(x); }
  117. template<typename eT> arma_inline static typename arma_cx_only<eT>::result expm1 (const eT& x) { arma_ignore(x); return eT(0); }
  118. #elif defined(ARMA_HAVE_TR1)
  119. template<typename eT> arma_inline static typename arma_integral_only<eT>::result expm1 (const eT x) { return eT( std::tr1::expm1(double(x)) ); }
  120. template<typename eT> arma_inline static typename arma_real_only<eT>::result expm1 (const eT x) { return std::tr1::expm1(x); }
  121. template<typename eT> arma_inline static typename arma_cx_only<eT>::result expm1 (const eT& x) { arma_ignore(x); return eT(0); }
  122. #else
  123. template<typename eT> arma_inline static eT expm1 (const eT x) { arma_ignore(x); arma_stop_logic_error("expm1(): C++11 compiler required"); return eT(0); }
  124. #endif
  125. template<typename eT> arma_inline static typename arma_unsigned_integral_only<eT>::result arma_abs (const eT x) { return x; }
  126. template<typename eT> arma_inline static typename arma_signed_integral_only<eT>::result arma_abs (const eT x) { return std::abs(x); }
  127. template<typename eT> arma_inline static typename arma_real_only<eT>::result arma_abs (const eT x) { return std::abs(x); }
  128. template<typename T> arma_inline static typename arma_real_only< T>::result arma_abs (const std::complex<T>& x) { return std::abs(x); }
  129. #if defined(ARMA_USE_CXX11)
  130. template<typename eT> arma_inline static typename arma_integral_only<eT>::result erf (const eT x) { return eT( std::erf(double(x)) ); }
  131. template<typename eT> arma_inline static typename arma_real_only<eT>::result erf (const eT x) { return std::erf(x); }
  132. template<typename eT> arma_inline static typename arma_cx_only<eT>::result erf (const eT& x) { arma_ignore(x); return eT(0); }
  133. #elif defined(ARMA_HAVE_TR1)
  134. template<typename eT> arma_inline static typename arma_integral_only<eT>::result erf (const eT x) { return eT( std::tr1::erf(double(x)) ); }
  135. template<typename eT> arma_inline static typename arma_real_only<eT>::result erf (const eT x) { return std::tr1::erf(x); }
  136. template<typename eT> arma_inline static typename arma_cx_only<eT>::result erf (const eT& x) { arma_ignore(x); return eT(0); }
  137. #else
  138. template<typename eT> arma_inline static eT erf (const eT x) { arma_ignore(x); arma_stop_logic_error("erf(): C++11 compiler required"); return eT(0); }
  139. #endif
  140. #if defined(ARMA_USE_CXX11)
  141. template<typename eT> arma_inline static typename arma_integral_only<eT>::result erfc (const eT x) { return eT( std::erfc(double(x)) ); }
  142. template<typename eT> arma_inline static typename arma_real_only<eT>::result erfc (const eT x) { return std::erfc(x); }
  143. template<typename eT> arma_inline static typename arma_cx_only<eT>::result erfc (const eT& x) { arma_ignore(x); return eT(0); }
  144. #elif defined(ARMA_HAVE_TR1)
  145. template<typename eT> arma_inline static typename arma_integral_only<eT>::result erfc (const eT x) { return eT( std::tr1::erfc(double(x)) ); }
  146. template<typename eT> arma_inline static typename arma_real_only<eT>::result erfc (const eT x) { return std::tr1::erfc(x); }
  147. template<typename eT> arma_inline static typename arma_cx_only<eT>::result erfc (const eT& x) { arma_ignore(x); return eT(0); }
  148. #else
  149. template<typename eT> arma_inline static eT erfc (const eT x) { arma_ignore(x); arma_stop_logic_error("erfc(): C++11 compiler required"); return eT(0); }
  150. #endif
  151. #if defined(ARMA_USE_CXX11)
  152. template<typename eT> arma_inline static typename arma_integral_only<eT>::result lgamma (const eT x) { return eT( std::lgamma(double(x)) ); }
  153. template<typename eT> arma_inline static typename arma_real_only<eT>::result lgamma (const eT x) { return std::lgamma(x); }
  154. template<typename eT> arma_inline static typename arma_cx_only<eT>::result lgamma (const eT& x) { arma_ignore(x); return eT(0); }
  155. #elif defined(ARMA_HAVE_TR1)
  156. template<typename eT> arma_inline static typename arma_integral_only<eT>::result lgamma (const eT x) { return eT( std::tr1::lgamma(double(x)) ); }
  157. template<typename eT> arma_inline static typename arma_real_only<eT>::result lgamma (const eT x) { return std::tr1::lgamma(x); }
  158. template<typename eT> arma_inline static typename arma_cx_only<eT>::result lgamma (const eT& x) { arma_ignore(x); return eT(0); }
  159. #else
  160. template<typename eT> arma_inline static eT lgamma (const eT x) { arma_ignore(x); arma_stop_logic_error("lgamma(): C++11 compiler required"); return eT(0); }
  161. #endif
  162. #ifndef LUOYC20220706
  163. static double pow_new (double base, double exponent)
  164. {
  165. double cc;
  166. if(isnan(base))
  167. {
  168. cc = 0*log((double)0);
  169. }else
  170. {
  171. cc = std::pow(base, exponent);
  172. }
  173. return cc;
  174. }
  175. template<typename T1, typename T2>
  176. static T1 pow_new(T1 base, T2 exponent)
  177. {
  178. T1 cc;
  179. if(isnan(base))
  180. {
  181. cc = 0*log((double)0);
  182. }else
  183. {
  184. cc = std::pow(base, exponent);
  185. }
  186. return cc;
  187. }
  188. template<typename T1, typename T2> arma_inline static typename arma_integral_only<T1>::result pow (const T1 base, const T2 exponent) { return T1( pow_new( double(base), double(exponent) ) ); }
  189. template<typename T1, typename T2> arma_inline static typename arma_real_or_cx_only<T1>::result pow (const T1 base, const T2 exponent) { return T1( pow_new( base, exponent ) ); }
  190. #else
  191. template<typename T1, typename T2> arma_inline static typename arma_integral_only<T1>::result pow (const T1 base, const T2 exponent) { return T1( std::pow( double(base), double(exponent) ) ); }
  192. template<typename T1, typename T2> arma_inline static typename arma_real_or_cx_only<T1>::result pow (const T1 base, const T2 exponent) { return T1( std::pow( base, exponent ) ); }
  193. #endif
  194. template<typename eT>
  195. arma_inline
  196. static
  197. typename arma_integral_only<eT>::result
  198. direct_eps(const eT)
  199. {
  200. return eT(0);
  201. }
  202. template<typename eT>
  203. inline
  204. static
  205. typename arma_real_only<eT>::result
  206. direct_eps(const eT x)
  207. {
  208. //arma_extra_debug_sigprint();
  209. // acording to IEEE Standard for Floating-Point Arithmetic (IEEE 754)
  210. // the mantissa length for double is 53 bits = std::numeric_limits<double>::digits
  211. // the mantissa length for float is 24 bits = std::numeric_limits<float >::digits
  212. //return std::pow( std::numeric_limits<eT>::radix, (std::floor(std::log10(std::abs(x))/std::log10(std::numeric_limits<eT>::radix))-(std::numeric_limits<eT>::digits-1)) );
  213. const eT radix_eT = eT(std::numeric_limits<eT>::radix);
  214. const eT digits_m1_eT = eT(std::numeric_limits<eT>::digits - 1);
  215. // return std::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
  216. return eop_aux::pow( radix_eT, eT(std::floor(std::log10(std::abs(x))/std::log10(radix_eT)) - digits_m1_eT) );
  217. }
  218. template<typename T>
  219. inline
  220. static
  221. typename arma_real_only<T>::result
  222. direct_eps(const std::complex<T>& x)
  223. {
  224. //arma_extra_debug_sigprint();
  225. //return std::pow( std::numeric_limits<T>::radix, (std::floor(std::log10(std::abs(x))/std::log10(std::numeric_limits<T>::radix))-(std::numeric_limits<T>::digits-1)) );
  226. const T radix_T = T(std::numeric_limits<T>::radix);
  227. const T digits_m1_T = T(std::numeric_limits<T>::digits - 1);
  228. return std::pow( radix_T, T(std::floor(std::log10(std::abs(x))/std::log10(radix_T)) - digits_m1_T) );
  229. }
  230. };
  231. //! @}