source: XIOS/dev/dev_olga/src/extern/boost/include/boost/numeric/interval/arith.hpp @ 1022

Last change on this file since 1022 was 1022, checked in by mhnguyen, 7 years ago
File size: 10.1 KB
Line 
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
11#ifndef BOOST_NUMERIC_INTERVAL_ARITH_HPP
12#define BOOST_NUMERIC_INTERVAL_ARITH_HPP
13
14#include <boost/config.hpp>
15#include <boost/numeric/interval/interval.hpp>
16#include <boost/numeric/interval/detail/bugs.hpp>
17#include <boost/numeric/interval/detail/test_input.hpp>
18#include <boost/numeric/interval/detail/division.hpp>
19#include <algorithm>
20
21namespace boost {
22namespace numeric {
23
24/*
25 * Basic arithmetic operators
26 */
27
28template<class T, class Policies> inline
29const interval<T, Policies>& operator+(const interval<T, Policies>& x)
30{
31  return x;
32}
33
34template<class T, class Policies> inline
35interval<T, Policies> operator-(const interval<T, Policies>& x)
36{
37  if (interval_lib::detail::test_input(x))
38    return interval<T, Policies>::empty();
39  return interval<T, Policies>(-x.upper(), -x.lower(), true);
40}
41
42template<class T, class Policies> inline
43interval<T, Policies>& interval<T, Policies>::operator+=(const interval<T, Policies>& r)
44{
45  if (interval_lib::detail::test_input(*this, r))
46    set_empty();
47  else {
48    typename Policies::rounding rnd;
49    set(rnd.add_down(low, r.low), rnd.add_up(up, r.up));
50  }
51  return *this;
52}
53
54template<class T, class Policies> inline
55interval<T, Policies>& interval<T, Policies>::operator+=(const T& r)
56{
57  if (interval_lib::detail::test_input(*this, r))
58    set_empty();
59  else {
60    typename Policies::rounding rnd;
61    set(rnd.add_down(low, r), rnd.add_up(up, r));
62  }
63  return *this;
64}
65
66template<class T, class Policies> inline
67interval<T, Policies>& interval<T, Policies>::operator-=(const interval<T, Policies>& r)
68{
69  if (interval_lib::detail::test_input(*this, r))
70    set_empty();
71  else {
72    typename Policies::rounding rnd;
73    set(rnd.sub_down(low, r.up), rnd.sub_up(up, r.low));
74  }
75  return *this;
76}
77
78template<class T, class Policies> inline
79interval<T, Policies>& interval<T, Policies>::operator-=(const T& r)
80{
81  if (interval_lib::detail::test_input(*this, r))
82    set_empty();
83  else {
84    typename Policies::rounding rnd;
85    set(rnd.sub_down(low, r), rnd.sub_up(up, r));
86  }
87  return *this;
88}
89
90template<class T, class Policies> inline
91interval<T, Policies>& interval<T, Policies>::operator*=(const interval<T, Policies>& r)
92{
93  return *this = *this * r;
94}
95
96template<class T, class Policies> inline
97interval<T, Policies>& interval<T, Policies>::operator*=(const T& r)
98{
99  return *this = r * *this;
100}
101
102template<class T, class Policies> inline
103interval<T, Policies>& interval<T, Policies>::operator/=(const interval<T, Policies>& r)
104{
105  return *this = *this / r;
106}
107
108template<class T, class Policies> inline
109interval<T, Policies>& interval<T, Policies>::operator/=(const T& r)
110{
111  return *this = *this / r;
112}
113
114template<class T, class Policies> inline
115interval<T, Policies> operator+(const interval<T, Policies>& x,
116                                const interval<T, Policies>& y)
117{
118  if (interval_lib::detail::test_input(x, y))
119    return interval<T, Policies>::empty();
120  typename Policies::rounding rnd;
121  return interval<T,Policies>(rnd.add_down(x.lower(), y.lower()),
122                              rnd.add_up  (x.upper(), y.upper()), true);
123}
124
125template<class T, class Policies> inline
126interval<T, Policies> operator+(const T& x, const interval<T, Policies>& y)
127{
128  if (interval_lib::detail::test_input(x, y))
129    return interval<T, Policies>::empty();
130  typename Policies::rounding rnd;
131  return interval<T,Policies>(rnd.add_down(x, y.lower()),
132                              rnd.add_up  (x, y.upper()), true);
133}
134
135template<class T, class Policies> inline
136interval<T, Policies> operator+(const interval<T, Policies>& x, const T& y)
137{ return y + x; }
138
139template<class T, class Policies> inline
140interval<T, Policies> operator-(const interval<T, Policies>& x,
141                                const interval<T, Policies>& 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.upper()),
147                              rnd.sub_up  (x.upper(), y.lower()), true);
148}
149
150template<class T, class Policies> inline
151interval<T, Policies> operator-(const T& x, const interval<T, Policies>& y)
152{
153  if (interval_lib::detail::test_input(x, y))
154    return interval<T, Policies>::empty();
155  typename Policies::rounding rnd;
156  return interval<T,Policies>(rnd.sub_down(x, y.upper()),
157                              rnd.sub_up  (x, y.lower()), true);
158}
159
160template<class T, class Policies> inline
161interval<T, Policies> operator-(const interval<T, Policies>& x, const T& y)
162{
163  if (interval_lib::detail::test_input(x, y))
164    return interval<T, Policies>::empty();
165  typename Policies::rounding rnd;
166  return interval<T,Policies>(rnd.sub_down(x.lower(), y),
167                              rnd.sub_up  (x.upper(), y), true);
168}
169
170template<class T, class Policies> inline
171interval<T, Policies> operator*(const interval<T, Policies>& x,
172                                const interval<T, Policies>& y)
173{
174  BOOST_USING_STD_MIN();
175  BOOST_USING_STD_MAX();
176  typedef interval<T, Policies> I;
177  if (interval_lib::detail::test_input(x, y))
178    return I::empty();
179  typename Policies::rounding rnd;
180  const T& xl = x.lower();
181  const T& xu = x.upper();
182  const T& yl = y.lower();
183  const T& yu = y.upper();
184
185  if (interval_lib::user::is_neg(xl))
186    if (interval_lib::user::is_pos(xu))
187      if (interval_lib::user::is_neg(yl))
188        if (interval_lib::user::is_pos(yu)) // M * M
189          return I(min BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_down(xl, yu), rnd.mul_down(xu, yl)),
190                   max BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.mul_up  (xl, yl), rnd.mul_up  (xu, yu)), true);
191        else                    // M * N
192          return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yl), true);
193      else
194        if (interval_lib::user::is_pos(yu)) // M * P
195          return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yu), true);
196        else                    // M * Z
197          return I(static_cast<T>(0), static_cast<T>(0), true);
198    else
199      if (interval_lib::user::is_neg(yl))
200        if (interval_lib::user::is_pos(yu)) // N * M
201          return I(rnd.mul_down(xl, yu), rnd.mul_up(xl, yl), true);
202        else                    // N * N
203          return I(rnd.mul_down(xu, yu), rnd.mul_up(xl, yl), true);
204      else
205        if (interval_lib::user::is_pos(yu)) // N * P
206          return I(rnd.mul_down(xl, yu), rnd.mul_up(xu, yl), true);
207        else                    // N * Z
208          return I(static_cast<T>(0), static_cast<T>(0), true);
209  else
210    if (interval_lib::user::is_pos(xu))
211      if (interval_lib::user::is_neg(yl))
212        if (interval_lib::user::is_pos(yu)) // P * M
213          return I(rnd.mul_down(xu, yl), rnd.mul_up(xu, yu), true);
214        else                    // P * N
215          return I(rnd.mul_down(xu, yl), rnd.mul_up(xl, yu), true);
216      else
217        if (interval_lib::user::is_pos(yu)) // P * P
218          return I(rnd.mul_down(xl, yl), rnd.mul_up(xu, yu), true);
219        else                    // P * Z
220          return I(static_cast<T>(0), static_cast<T>(0), true);
221    else                        // Z * ?
222      return I(static_cast<T>(0), static_cast<T>(0), true);
223}
224
225template<class T, class Policies> inline
226interval<T, Policies> operator*(const T& x, const interval<T, Policies>& y)
227{ 
228  typedef interval<T, Policies> I;
229  if (interval_lib::detail::test_input(x, y))
230    return I::empty();
231  typename Policies::rounding rnd;
232  const T& yl = y.lower();
233  const T& yu = y.upper();
234  // x is supposed not to be infinite
235  if (interval_lib::user::is_neg(x))
236    return I(rnd.mul_down(x, yu), rnd.mul_up(x, yl), true);
237  else if (interval_lib::user::is_zero(x))
238    return I(static_cast<T>(0), static_cast<T>(0), true);
239  else
240    return I(rnd.mul_down(x, yl), rnd.mul_up(x, yu), true);
241}
242
243template<class T, class Policies> inline
244interval<T, Policies> operator*(const interval<T, Policies>& x, const T& y)
245{ return y * x; }
246
247template<class T, class Policies> inline
248interval<T, Policies> operator/(const interval<T, Policies>& x,
249                                const interval<T, Policies>& y)
250{
251  if (interval_lib::detail::test_input(x, y))
252    return interval<T, Policies>::empty();
253  if (zero_in(y))
254    if (!interval_lib::user::is_zero(y.lower()))
255      if (!interval_lib::user::is_zero(y.upper()))
256        return interval_lib::detail::div_zero(x);
257      else
258        return interval_lib::detail::div_negative(x, y.lower());
259    else
260      if (!interval_lib::user::is_zero(y.upper()))
261        return interval_lib::detail::div_positive(x, y.upper());
262      else
263        return interval<T, Policies>::empty();
264  else
265    return interval_lib::detail::div_non_zero(x, y);
266}
267
268template<class T, class Policies> inline
269interval<T, Policies> operator/(const T& x, const interval<T, Policies>& y)
270{
271  if (interval_lib::detail::test_input(x, y))
272    return interval<T, Policies>::empty();
273  if (zero_in(y))
274    if (!interval_lib::user::is_zero(y.lower()))
275      if (!interval_lib::user::is_zero(y.upper()))
276        return interval_lib::detail::div_zero<T, Policies>(x);
277      else
278        return interval_lib::detail::div_negative<T, Policies>(x, y.lower());
279    else
280      if (!interval_lib::user::is_zero(y.upper()))
281        return interval_lib::detail::div_positive<T, Policies>(x, y.upper());
282      else
283        return interval<T, Policies>::empty();
284  else
285    return interval_lib::detail::div_non_zero(x, y);
286}
287
288template<class T, class Policies> inline
289interval<T, Policies> operator/(const interval<T, Policies>& x, const T& y)
290{
291  if (interval_lib::detail::test_input(x, y) || interval_lib::user::is_zero(y))
292    return interval<T, Policies>::empty();
293  typename Policies::rounding rnd;
294  const T& xl = x.lower();
295  const T& xu = x.upper();
296  if (interval_lib::user::is_neg(y))
297    return interval<T, Policies>(rnd.div_down(xu, y), rnd.div_up(xl, y), true);
298  else
299    return interval<T, Policies>(rnd.div_down(xl, y), rnd.div_up(xu, y), true);
300}
301
302} // namespace numeric
303} // namespace boost
304
305#endif // BOOST_NUMERIC_INTERVAL_ARITH_HPP
Note: See TracBrowser for help on using the repository browser.