source: NonGTP/Boost/boost/numeric/interval/arith2.hpp @ 857

Revision 857, 7.0 KB checked in by igarcia, 19 years ago (diff)
Line 
1/* Boost interval/arith2.hpp template implementation file
2 *
3 * This header provides some auxiliary arithmetic
4 * functions: fmod, sqrt, square, pov, inverse and
5 * a multi-interval division.
6 *
7 * Copyright 2002-2003 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
8 *
9 * Distributed under the Boost Software License, Version 1.0.
10 * (See accompanying file LICENSE_1_0.txt or
11 * copy at http://www.boost.org/LICENSE_1_0.txt)
12 */
13
14#ifndef BOOST_NUMERIC_INTERVAL_ARITH2_HPP
15#define BOOST_NUMERIC_INTERVAL_ARITH2_HPP
16
17#include <boost/config.hpp>
18#include <boost/numeric/interval/detail/interval_prototype.hpp>
19#include <boost/numeric/interval/detail/test_input.hpp>
20#include <boost/numeric/interval/detail/bugs.hpp>
21#include <boost/numeric/interval/detail/division.hpp>
22#include <boost/numeric/interval/arith.hpp>
23#include <boost/numeric/interval/policies.hpp>
24#include <algorithm>
25#include <cmath>
26
27namespace boost {
28namespace numeric {
29
30template<class T, class Policies> inline
31interval<T, Policies> fmod(const interval<T, Policies>& x,
32                           const interval<T, Policies>& y)
33{
34  if (interval_lib::detail::test_input(x, y))
35    return interval<T, Policies>::empty();
36  typename Policies::rounding rnd;
37  typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
38  T const &yb = interval_lib::user::is_neg(x.lower()) ? y.lower() : y.upper();
39  T n = rnd.int_down(rnd.div_down(x.lower(), yb));
40  return (const I&)x - n * (const I&)y;
41}
42
43template<class T, class Policies> inline
44interval<T, Policies> fmod(const interval<T, Policies>& x, const T& y)
45{
46  if (interval_lib::detail::test_input(x, y))
47    return interval<T, Policies>::empty();
48  typename Policies::rounding rnd;
49  typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
50  T n = rnd.int_down(rnd.div_down(x.lower(), y));
51  return (const I&)x - n * I(y);
52}
53
54template<class T, class Policies> inline
55interval<T, Policies> fmod(const T& x, const interval<T, Policies>& y)
56{
57  if (interval_lib::detail::test_input(x, y))
58    return interval<T, Policies>::empty();
59  typename Policies::rounding rnd;
60  typedef typename interval_lib::unprotect<interval<T, Policies> >::type I;
61  T const &yb = interval_lib::user::is_neg(x) ? y.lower() : y.upper();
62  T n = rnd.int_down(rnd.div_down(x, yb));
63  return x - n * (const I&)y;
64}
65
66namespace interval_lib {
67
68template<class T, class Policies> inline
69interval<T, Policies> division_part1(const interval<T, Policies>& x,
70                                     const interval<T, Policies>& y, bool& b)
71{
72  typedef interval<T, Policies> I;
73  b = false;
74  if (detail::test_input(x, y))
75    return I::empty();
76  if (in_zero(y))
77    if (!user::is_zero(y.lower()))
78      if (!user::is_zero(y.upper()))
79        return detail::div_zero_part1(x, y, b);
80      else
81        return detail::div_negative(x, y.lower());
82    else
83      if (!user::is_zero(y.upper()))
84        return detail::div_positive(x, y.upper());
85      else
86        return I::empty();
87  else
88    return detail::div_non_zero(x, y);
89}
90
91template<class T, class Policies> inline
92interval<T, Policies> division_part2(const interval<T, Policies>& x,
93                                     const interval<T, Policies>& y, bool b = true)
94{
95  if (!b) return interval<T, Policies>::empty();
96  return detail::div_zero_part2(x, y);
97}
98
99template<class T, class Policies> inline
100interval<T, Policies> multiplicative_inverse(const interval<T, Policies>& x)
101{
102  typedef interval<T, Policies> I;
103  if (detail::test_input(x))
104    return I::empty();
105  T one = static_cast<T>(1);
106  typename Policies::rounding rnd;
107  if (in_zero(x)) {
108    typedef typename Policies::checking checking;
109    if (!user::is_zero(x.lower()))
110      if (!user::is_zero(x.upper()))
111        return I::whole();
112      else
113        return I(checking::neg_inf(), rnd.div_up(one, x.lower()), true);
114    else
115      if (!user::is_zero(x.upper()))
116        return I(rnd.div_down(one, x.upper()), checking::pos_inf(), true);
117      else
118        return I::empty();
119  } else
120    return I(rnd.div_down(one, x.upper()), rnd.div_up(one, x.lower()), true);
121}
122
123namespace detail {
124
125template<class T, class Rounding> inline
126T pow_aux(const T& x_, int pwr, Rounding& rnd) // x and pwr are positive
127{
128  T x = x_;
129  T y = (pwr & 1) ? x_ : static_cast<T>(1);
130  pwr >>= 1;
131  while (pwr > 0) {
132    x = rnd.mul_up(x, x);
133    if (pwr & 1) y = rnd.mul_up(x, y);
134    pwr >>= 1;
135  }
136  return y;
137}
138
139} // namespace detail
140} // namespace interval_lib
141
142template<class T, class Policies> inline
143interval<T, Policies> pow(const interval<T, Policies>& x, int pwr)
144{
145  BOOST_USING_STD_MAX();
146  using interval_lib::detail::pow_aux;
147  typedef interval<T, Policies> I;
148
149  if (interval_lib::detail::test_input(x))
150    return I::empty();
151
152  if (pwr == 0)
153    if (interval_lib::user::is_zero(x.lower())
154        && interval_lib::user::is_zero(x.upper()))
155      return I::empty();
156    else
157      return I(static_cast<T>(1));
158  else if (pwr < 0)
159    return interval_lib::multiplicative_inverse(pow(x, -pwr));
160
161  typename Policies::rounding rnd;
162 
163  if (interval_lib::user::is_neg(x.upper())) {        // [-2,-1]
164    T yl = pow_aux(-x.upper(), pwr, rnd);
165    T yu = pow_aux(-x.lower(), pwr, rnd);
166    if (pwr & 1)     // [-2,-1]^1
167      return I(-yu, -yl, true);
168    else             // [-2,-1]^2
169      return I(yl, yu, true);
170  } else if (interval_lib::user::is_neg(x.lower())) { // [-1,1]
171    if (pwr & 1) {   // [-1,1]^1
172      return I(-pow_aux(-x.lower(), pwr, rnd), pow_aux(x.upper(), pwr, rnd), true);
173    } else {         // [-1,1]^2
174      return I(static_cast<T>(0), pow_aux(max BOOST_PREVENT_MACRO_SUBSTITUTION(-x.lower(), x.upper()), pwr, rnd), true);
175    }
176  } else {                                // [1,2]
177    return I(pow_aux(x.lower(), pwr, rnd), pow_aux(x.upper(), pwr, rnd), true);
178  }
179}
180
181template<class T, class Policies> inline
182interval<T, Policies> sqrt(const interval<T, Policies>& x)
183{
184  typedef interval<T, Policies> I;
185  if (interval_lib::detail::test_input(x) || interval_lib::user::is_neg(x.upper()))
186    return I::empty();
187  typename Policies::rounding rnd;
188  T l = !interval_lib::user::is_pos(x.lower()) ? static_cast<T>(0) : rnd.sqrt_down(x.lower());
189  return I(l, rnd.sqrt_up(x.upper()), true);
190}
191
192template<class T, class Policies> inline
193interval<T, Policies> square(const interval<T, Policies>& x)
194{
195  typedef interval<T, Policies> I;
196  if (interval_lib::detail::test_input(x))
197    return I::empty();
198  typename Policies::rounding rnd;
199  const T& xl = x.lower();
200  const T& xu = x.upper();
201  if (interval_lib::user::is_neg(xu))
202    return I(rnd.mul_down(xu, xu), rnd.mul_up(xl, xl), true);
203  else if (interval_lib::user::is_pos(x.lower()))
204    return I(rnd.mul_down(xl, xl), rnd.mul_up(xu, xu), true);
205  else
206    return I(static_cast<T>(0), (-xl > xu ? rnd.mul_up(xl, xl) : rnd.mul_up(xu, xu)), true);
207}
208
209} // namespace numeric
210} // namespace boost
211
212#endif // BOOST_NUMERIC_INTERVAL_ARITH2_HPP
Note: See TracBrowser for help on using the repository browser.