arith.hpp 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305
  1. /* Boost interval/arith.hpp template implementation file
  2. *
  3. * Copyright 2000 Jens Maurer
  4. * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
  5. *
  6. * Distributed under the Boost Software License, Version 1.0.
  7. * (See accompanying file LICENSE_1_0.txt or
  8. * copy at http://www.boost.org/LICENSE_1_0.txt)
  9. */
  10. #ifndef BOOST_NUMERIC_INTERVAL_ARITH_HPP
  11. #define BOOST_NUMERIC_INTERVAL_ARITH_HPP
  12. #include <boost/config.hpp>
  13. #include <boost/numeric/interval/interval.hpp>
  14. #include <boost/numeric/interval/detail/bugs.hpp>
  15. #include <boost/numeric/interval/detail/test_input.hpp>
  16. #include <boost/numeric/interval/detail/division.hpp>
  17. #include <algorithm>
  18. namespace boost {
  19. namespace numeric {
  20. /*
  21. * Basic arithmetic operators
  22. */
  23. template<class T, class Policies> inline
  24. const interval<T, Policies>& operator+(const interval<T, Policies>& x)
  25. {
  26. return x;
  27. }
  28. template<class T, class Policies> inline
  29. interval<T, Policies> operator-(const interval<T, Policies>& x)
  30. {
  31. if (interval_lib::detail::test_input(x))
  32. return interval<T, Policies>::empty();
  33. return interval<T, Policies>(-x.upper(), -x.lower(), true);
  34. }
  35. template<class T, class Policies> inline
  36. interval<T, Policies>& interval<T, Policies>::operator+=(const interval<T, Policies>& r)
  37. {
  38. if (interval_lib::detail::test_input(*this, r))
  39. set_empty();
  40. else {
  41. typename Policies::rounding rnd;
  42. set(rnd.add_down(low, r.low), rnd.add_up(up, r.up));
  43. }
  44. return *this;
  45. }
  46. template<class T, class Policies> inline
  47. interval<T, Policies>& interval<T, Policies>::operator+=(const T& r)
  48. {
  49. if (interval_lib::detail::test_input(*this, r))
  50. set_empty();
  51. else {
  52. typename Policies::rounding rnd;
  53. set(rnd.add_down(low, r), rnd.add_up(up, r));
  54. }
  55. return *this;
  56. }
  57. template<class T, class Policies> inline
  58. interval<T, Policies>& interval<T, Policies>::operator-=(const interval<T, Policies>& r)
  59. {
  60. if (interval_lib::detail::test_input(*this, r))
  61. set_empty();
  62. else {
  63. typename Policies::rounding rnd;
  64. set(rnd.sub_down(low, r.up), rnd.sub_up(up, r.low));
  65. }
  66. return *this;
  67. }
  68. template<class T, class Policies> inline
  69. interval<T, Policies>& interval<T, Policies>::operator-=(const T& r)
  70. {
  71. if (interval_lib::detail::test_input(*this, r))
  72. set_empty();
  73. else {
  74. typename Policies::rounding rnd;
  75. set(rnd.sub_down(low, r), rnd.sub_up(up, r));
  76. }
  77. return *this;
  78. }
  79. template<class T, class Policies> inline
  80. interval<T, Policies>& interval<T, Policies>::operator*=(const interval<T, Policies>& r)
  81. {
  82. return *this = *this * r;
  83. }
  84. template<class T, class Policies> inline
  85. interval<T, Policies>& interval<T, Policies>::operator*=(const T& r)
  86. {
  87. return *this = r * *this;
  88. }
  89. template<class T, class Policies> inline
  90. interval<T, Policies>& interval<T, Policies>::operator/=(const interval<T, Policies>& r)
  91. {
  92. return *this = *this / r;
  93. }
  94. template<class T, class Policies> inline
  95. interval<T, Policies>& interval<T, Policies>::operator/=(const T& r)
  96. {
  97. return *this = *this / r;
  98. }
  99. template<class T, class Policies> inline
  100. interval<T, Policies> operator+(const interval<T, Policies>& x,
  101. const interval<T, Policies>& y)
  102. {
  103. if (interval_lib::detail::test_input(x, y))
  104. return interval<T, Policies>::empty();
  105. typename Policies::rounding rnd;
  106. return interval<T,Policies>(rnd.add_down(x.lower(), y.lower()),
  107. rnd.add_up (x.upper(), y.upper()), true);
  108. }
  109. template<class T, class Policies> inline
  110. interval<T, Policies> operator+(const T& x, const interval<T, Policies>& y)
  111. {
  112. if (interval_lib::detail::test_input(x, y))
  113. return interval<T, Policies>::empty();
  114. typename Policies::rounding rnd;
  115. return interval<T,Policies>(rnd.add_down(x, y.lower()),
  116. rnd.add_up (x, y.upper()), true);
  117. }
  118. template<class T, class Policies> inline
  119. interval<T, Policies> operator+(const interval<T, Policies>& x, const T& y)
  120. { return y + x; }
  121. template<class T, class Policies> inline
  122. interval<T, Policies> operator-(const interval<T, Policies>& x,
  123. const interval<T, Policies>& y)
  124. {
  125. if (interval_lib::detail::test_input(x, y))
  126. return interval<T, Policies>::empty();
  127. typename Policies::rounding rnd;
  128. return interval<T,Policies>(rnd.sub_down(x.lower(), y.upper()),
  129. rnd.sub_up (x.upper(), y.lower()), true);
  130. }
  131. template<class T, class Policies> inline
  132. interval<T, Policies> operator-(const T& x, const interval<T, Policies>& y)
  133. {
  134. if (interval_lib::detail::test_input(x, y))
  135. return interval<T, Policies>::empty();
  136. typename Policies::rounding rnd;
  137. return interval<T,Policies>(rnd.sub_down(x, y.upper()),
  138. rnd.sub_up (x, y.lower()), true);
  139. }
  140. template<class T, class Policies> inline
  141. interval<T, Policies> operator-(const interval<T, Policies>& x, const T& y)
  142. {
  143. if (interval_lib::detail::test_input(x, y))
  144. return interval<T, Policies>::empty();
  145. typename Policies::rounding rnd;
  146. return interval<T,Policies>(rnd.sub_down(x.lower(), y),
  147. rnd.sub_up (x.upper(), y), true);
  148. }
  149. template<class T, class Policies> inline
  150. interval<T, Policies> operator*(const interval<T, Policies>& x,
  151. const interval<T, Policies>& y)
  152. {
  153. BOOST_USING_STD_MIN();
  154. BOOST_USING_STD_MAX();
  155. typedef interval<T, Policies> I;
  156. if (interval_lib::detail::test_input(x, y))
  157. return I::empty();
  158. typename Policies::rounding rnd;
  159. const T& xl = x.lower();
  160. const T& xu = x.upper();
  161. const T& yl = y.lower();
  162. const T& yu = y.upper();
  163. if (interval_lib::user::is_neg(xl))
  164. if (interval_lib::user::is_pos(xu))
  165. if (interval_lib::user::is_neg(yl))
  166. if (interval_lib::user::is_pos(yu)) // M * M
  167. return I(min BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_down(xl, yu), rnd.mul_down(xu, yl)),
  168. max BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_up (xl, yl), rnd.mul_up (xu, yu)), true);
  169. else // M * N
  170. return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yl), true);
  171. else
  172. if (interval_lib::user::is_pos(yu)) // M * P
  173. return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yu), true);
  174. else // M * Z
  175. return I(static_cast<T>(0), static_cast<T>(0), true);
  176. else
  177. if (interval_lib::user::is_neg(yl))
  178. if (interval_lib::user::is_pos(yu)) // N * M
  179. return I(rnd.mul_down(xl, yu), rnd.mul_up(xl, yl), true);
  180. else // N * N
  181. return I(rnd.mul_down(xu, yu), rnd.mul_up(xl, yl), true);
  182. else
  183. if (interval_lib::user::is_pos(yu)) // N * P
  184. return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yl), true);
  185. else // N * Z
  186. return I(static_cast<T>(0), static_cast<T>(0), true);
  187. else
  188. if (interval_lib::user::is_pos(xu))
  189. if (interval_lib::user::is_neg(yl))
  190. if (interval_lib::user::is_pos(yu)) // P * M
  191. return I(rnd.mul_down(xu, yl), rnd.mul_up(xu, yu), true);
  192. else // P * N
  193. return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yu), true);
  194. else
  195. if (interval_lib::user::is_pos(yu)) // P * P
  196. return I(rnd.mul_down(xl, yl), rnd.mul_up(xu, yu), true);
  197. else // P * Z
  198. return I(static_cast<T>(0), static_cast<T>(0), true);
  199. else // Z * ?
  200. return I(static_cast<T>(0), static_cast<T>(0), true);
  201. }
  202. template<class T, class Policies> inline
  203. interval<T, Policies> operator*(const T& x, const interval<T, Policies>& y)
  204. {
  205. typedef interval<T, Policies> I;
  206. if (interval_lib::detail::test_input(x, y))
  207. return I::empty();
  208. typename Policies::rounding rnd;
  209. const T& yl = y.lower();
  210. const T& yu = y.upper();
  211. // x is supposed not to be infinite
  212. if (interval_lib::user::is_neg(x))
  213. return I(rnd.mul_down(x, yu), rnd.mul_up(x, yl), true);
  214. else if (interval_lib::user::is_zero(x))
  215. return I(static_cast<T>(0), static_cast<T>(0), true);
  216. else
  217. return I(rnd.mul_down(x, yl), rnd.mul_up(x, yu), true);
  218. }
  219. template<class T, class Policies> inline
  220. interval<T, Policies> operator*(const interval<T, Policies>& x, const T& y)
  221. { return y * x; }
  222. template<class T, class Policies> inline
  223. interval<T, Policies> operator/(const interval<T, Policies>& x,
  224. const interval<T, Policies>& y)
  225. {
  226. if (interval_lib::detail::test_input(x, y))
  227. return interval<T, Policies>::empty();
  228. if (zero_in(y))
  229. if (!interval_lib::user::is_zero(y.lower()))
  230. if (!interval_lib::user::is_zero(y.upper()))
  231. return interval_lib::detail::div_zero(x);
  232. else
  233. return interval_lib::detail::div_negative(x, y.lower());
  234. else
  235. if (!interval_lib::user::is_zero(y.upper()))
  236. return interval_lib::detail::div_positive(x, y.upper());
  237. else
  238. return interval<T, Policies>::empty();
  239. else
  240. return interval_lib::detail::div_non_zero(x, y);
  241. }
  242. template<class T, class Policies> inline
  243. interval<T, Policies> operator/(const T& x, const interval<T, Policies>& y)
  244. {
  245. if (interval_lib::detail::test_input(x, y))
  246. return interval<T, Policies>::empty();
  247. if (zero_in(y))
  248. if (!interval_lib::user::is_zero(y.lower()))
  249. if (!interval_lib::user::is_zero(y.upper()))
  250. return interval_lib::detail::div_zero<T, Policies>(x);
  251. else
  252. return interval_lib::detail::div_negative<T, Policies>(x, y.lower());
  253. else
  254. if (!interval_lib::user::is_zero(y.upper()))
  255. return interval_lib::detail::div_positive<T, Policies>(x, y.upper());
  256. else
  257. return interval<T, Policies>::empty();
  258. else
  259. return interval_lib::detail::div_non_zero(x, y);
  260. }
  261. template<class T, class Policies> inline
  262. interval<T, Policies> operator/(const interval<T, Policies>& x, const T& y)
  263. {
  264. if (interval_lib::detail::test_input(x, y) || interval_lib::user::is_zero(y))
  265. return interval<T, Policies>::empty();
  266. typename Policies::rounding rnd;
  267. const T& xl = x.lower();
  268. const T& xu = x.upper();
  269. if (interval_lib::user::is_neg(y))
  270. return interval<T, Policies>(rnd.div_down(xu, y), rnd.div_up(xl, y), true);
  271. else
  272. return interval<T, Policies>(rnd.div_down(xl, y), rnd.div_up(xu, y), true);
  273. }
  274. } // namespace numeric
  275. } // namespace boost
  276. #endif // BOOST_NUMERIC_INTERVAL_ARITH_HPP