eigen.hpp 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248
  1. ///////////////////////////////////////////////////////////////////////////////
  2. // Copyright 2018 John Maddock. Distributed under the Boost
  3. // Software License, Version 1.0. (See accompanying file
  4. // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  5. #ifndef BOOST_MP_EIGEN_HPP
  6. #define BOOST_MP_EIGEN_HPP
  7. #include <boost/multiprecision/number.hpp>
  8. #include <Eigen/Core>
  9. //
  10. // Generic Eigen support code:
  11. //
  12. namespace Eigen {
  13. template <class B1, class B2>
  14. struct NumTraitsImp;
  15. template <class B1>
  16. struct NumTraitsImp<B1, B1>
  17. {
  18. using self_type = B1;
  19. using Real = typename boost::multiprecision::scalar_result_from_possible_complex<self_type>::type;
  20. using NonInteger = self_type; // Not correct but we can't do much better??
  21. using Literal = double;
  22. using Nested = self_type;
  23. enum
  24. {
  25. IsComplex = boost::multiprecision::number_category<self_type>::value == boost::multiprecision::number_kind_complex,
  26. IsInteger = boost::multiprecision::number_category<self_type>::value == boost::multiprecision::number_kind_integer,
  27. ReadCost = 1,
  28. AddCost = 4,
  29. MulCost = 8,
  30. IsSigned = std::numeric_limits<self_type>::is_specialized ? std::numeric_limits<self_type>::is_signed : true,
  31. RequireInitialization = 1,
  32. };
  33. static Real epsilon()
  34. {
  35. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  36. return std::numeric_limits<Real>::epsilon();
  37. }
  38. static Real dummy_precision()
  39. {
  40. return 1000 * epsilon();
  41. }
  42. static Real highest()
  43. {
  44. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  45. return (std::numeric_limits<Real>::max)();
  46. }
  47. static Real lowest()
  48. {
  49. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  50. return (std::numeric_limits<Real>::min)();
  51. }
  52. static int digits10_imp(const std::integral_constant<bool, true>&)
  53. {
  54. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  55. return std::numeric_limits<Real>::digits10;
  56. }
  57. template <bool B>
  58. static int digits10_imp(const std::integral_constant<bool, B>&)
  59. {
  60. return Real::thread_default_precision();
  61. }
  62. static int digits10()
  63. {
  64. return digits10_imp(std::integral_constant < bool, std::numeric_limits<Real>::digits10 && (std::numeric_limits<Real>::digits10 != INT_MAX) ? true : false > ());
  65. }
  66. static int digits()
  67. {
  68. // return the number of digits in the component type in case Real is complex
  69. // and we have no numeric_limits specialization.
  70. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  71. return std::numeric_limits<Real>::digits;
  72. }
  73. static int min_exponent()
  74. {
  75. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  76. return std::numeric_limits<Real>::min_exponent;
  77. }
  78. static int max_exponent()
  79. {
  80. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  81. return std::numeric_limits<Real>::max_exponent;
  82. }
  83. static Real infinity()
  84. {
  85. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  86. return std::numeric_limits<Real>::infinity();
  87. }
  88. static Real quiet_NaN()
  89. {
  90. static_assert(std::numeric_limits<Real>::is_specialized, "Eigen's NumTraits instantiated on a type with no numeric_limits support. Are you using a variable precision type?");
  91. return std::numeric_limits<Real>::quiet_NaN();
  92. }
  93. };
  94. template <class B1, class B2>
  95. struct NumTraitsImp : public NumTraitsImp<B2, B2>
  96. {
  97. //
  98. // This version is instantiated when B1 and B2 are different types, this happens for rational/complex/interval
  99. // types, in which case many methods defer to those of the "component type" B2.
  100. //
  101. using self_type = B1;
  102. using Real = typename boost::multiprecision::scalar_result_from_possible_complex<self_type>::type;
  103. using NonInteger = self_type; // Not correct but we can't do much better??
  104. using Literal = double;
  105. using Nested = self_type;
  106. enum
  107. {
  108. IsComplex = boost::multiprecision::number_category<self_type>::value == boost::multiprecision::number_kind_complex,
  109. IsInteger = boost::multiprecision::number_category<self_type>::value == boost::multiprecision::number_kind_integer,
  110. ReadCost = 1,
  111. AddCost = 4,
  112. MulCost = 8,
  113. IsSigned = std::numeric_limits<self_type>::is_specialized ? std::numeric_limits<self_type>::is_signed : true,
  114. RequireInitialization = 1,
  115. };
  116. static B2 epsilon()
  117. {
  118. return NumTraitsImp<B2, B2>::epsilon();
  119. }
  120. static B2 dummy_precision()
  121. {
  122. return 1000 * epsilon();
  123. }
  124. static B2 highest()
  125. {
  126. return NumTraitsImp<B2, B2>::highest();
  127. }
  128. static B2 lowest()
  129. {
  130. return NumTraitsImp<B2, B2>::lowest();
  131. }
  132. static int digits10()
  133. {
  134. return NumTraitsImp<B2, B2>::digits10();
  135. }
  136. static int digits()
  137. {
  138. return NumTraitsImp<B2, B2>::digits();
  139. }
  140. static int min_exponent()
  141. {
  142. return NumTraitsImp<B2, B2>::min_exponent();
  143. }
  144. static int max_exponent()
  145. {
  146. return NumTraitsImp<B2, B2>::max_exponent();
  147. }
  148. static B2 infinity()
  149. {
  150. return NumTraitsImp<B2, B2>::infinity();
  151. }
  152. static B2 quiet_NaN()
  153. {
  154. return NumTraitsImp<B2, B2>::quiet_NaN();
  155. }
  156. };
  157. template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates>
  158. struct NumTraits<boost::multiprecision::number<Backend, ExpressionTemplates> > : public NumTraitsImp<boost::multiprecision::number<Backend, ExpressionTemplates>, typename boost::multiprecision::number<Backend, ExpressionTemplates>::value_type>
  159. {};
  160. template <class tag, class Arg1, class Arg2, class Arg3, class Arg4>
  161. struct NumTraits<boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4> > : public NumTraits<typename boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type>
  162. {};
  163. #define BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(A) \
  164. template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, typename BinaryOp> \
  165. struct ScalarBinaryOpTraits<boost::multiprecision::number<Backend, ExpressionTemplates>, A, BinaryOp> \
  166. { \
  167. /*static_assert(boost::multiprecision::is_compatible_arithmetic_type<A, boost::multiprecision::number<Backend, ExpressionTemplates> >::value, "Interoperability with this arithmetic type is not supported.");*/ \
  168. using ReturnType = boost::multiprecision::number<Backend, ExpressionTemplates>; \
  169. }; \
  170. template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, typename BinaryOp> \
  171. struct ScalarBinaryOpTraits<A, boost::multiprecision::number<Backend, ExpressionTemplates>, BinaryOp> \
  172. { \
  173. /*static_assert(boost::multiprecision::is_compatible_arithmetic_type<A, boost::multiprecision::number<Backend, ExpressionTemplates> >::value, "Interoperability with this arithmetic type is not supported.");*/ \
  174. using ReturnType = boost::multiprecision::number<Backend, ExpressionTemplates>; \
  175. };
  176. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(float)
  177. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(double)
  178. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(long double)
  179. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(char)
  180. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned char)
  181. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(signed char)
  182. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(short)
  183. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned short)
  184. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(int)
  185. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned int)
  186. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(long)
  187. BOOST_MP_EIGEN_SCALAR_TRAITS_DECL(unsigned long)
  188. #if 0
  189. template<class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, class Backend2, boost::multiprecision::expression_template_option ExpressionTemplates2, typename BinaryOp>
  190. struct ScalarBinaryOpTraits<boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::number<Backend2, ExpressionTemplates2>, BinaryOp>
  191. {
  192. static_assert(
  193. boost::multiprecision::is_compatible_arithmetic_type<boost::multiprecision::number<Backend2, ExpressionTemplates2>, boost::multiprecision::number<Backend, ExpressionTemplates> >::value
  194. || boost::multiprecision::is_compatible_arithmetic_type<boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::number<Backend2, ExpressionTemplates2> >::value, "Interoperability with this arithmetic type is not supported.");
  195. using ReturnType = typename std::conditional<std::is_convertible<boost::multiprecision::number<Backend2, ExpressionTemplates2>, boost::multiprecision::number<Backend, ExpressionTemplates> >::value,
  196. boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::number<Backend2, ExpressionTemplates2> >::type;
  197. };
  198. template<unsigned D, typename BinaryOp>
  199. struct ScalarBinaryOpTraits<boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<D>, boost::multiprecision::et_on>, boost::multiprecision::mpfr_float, BinaryOp>
  200. {
  201. using ReturnType = boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<D>, boost::multiprecision::et_on>;
  202. };
  203. template<typename BinaryOp>
  204. struct ScalarBinaryOpTraits<boost::multiprecision::mpfr_float, boost::multiprecision::mpc_complex, BinaryOp>
  205. {
  206. using ReturnType = boost::multiprecision::number<boost::multiprecision::backends::mpc_complex_backend<0>, boost::multiprecision::et_on>;
  207. };
  208. template<class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, typename BinaryOp>
  209. struct ScalarBinaryOpTraits<boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::number<Backend, ExpressionTemplates>, BinaryOp>
  210. {
  211. using ReturnType = boost::multiprecision::number<Backend, ExpressionTemplates>;
  212. };
  213. #endif
  214. template <class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, class tag, class Arg1, class Arg2, class Arg3, class Arg4, typename BinaryOp>
  215. struct ScalarBinaryOpTraits<boost::multiprecision::number<Backend, ExpressionTemplates>, boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, BinaryOp>
  216. {
  217. static_assert(std::is_convertible<typename boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type, boost::multiprecision::number<Backend, ExpressionTemplates> >::value, "Interoperability with this arithmetic type is not supported.");
  218. using ReturnType = boost::multiprecision::number<Backend, ExpressionTemplates>;
  219. };
  220. template <class tag, class Arg1, class Arg2, class Arg3, class Arg4, class Backend, boost::multiprecision::expression_template_option ExpressionTemplates, typename BinaryOp>
  221. struct ScalarBinaryOpTraits<boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>, boost::multiprecision::number<Backend, ExpressionTemplates>, BinaryOp>
  222. {
  223. static_assert(std::is_convertible<typename boost::multiprecision::detail::expression<tag, Arg1, Arg2, Arg3, Arg4>::result_type, boost::multiprecision::number<Backend, ExpressionTemplates> >::value, "Interoperability with this arithmetic type is not supported.");
  224. using ReturnType = boost::multiprecision::number<Backend, ExpressionTemplates>;
  225. };
  226. namespace numext {
  227. using boost::multiprecision::conj;
  228. }
  229. } // namespace Eigen
  230. #endif