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
27 namespace boost {
28 namespace numeric {
29
30 template<class T, class Policies> inline
fmod(const interval<T,Policies> & x,const interval<T,Policies> & y)31 interval<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
43 template<class T, class Policies> inline
fmod(const interval<T,Policies> & x,const T & y)44 interval<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
54 template<class T, class Policies> inline
fmod(const T & x,const interval<T,Policies> & y)55 interval<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
66 namespace interval_lib {
67
68 template<class T, class Policies> inline
division_part1(const interval<T,Policies> & x,const interval<T,Policies> & y,bool & b)69 interval<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
91 template<class T, class Policies> inline
division_part2(const interval<T,Policies> & x,const interval<T,Policies> & y,bool b=true)92 interval<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
99 template<class T, class Policies> inline
multiplicative_inverse(const interval<T,Policies> & x)100 interval<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
123 namespace detail {
124
125 template<class T, class Rounding> inline
pow_aux(const T & x_,int pwr,Rounding & rnd)126 T 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
142 template<class T, class Policies> inline
pow(const interval<T,Policies> & x,int pwr)143 interval<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
181 template<class T, class Policies> inline
sqrt(const interval<T,Policies> & x)182 interval<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
192 template<class T, class Policies> inline
square(const interval<T,Policies> & x)193 interval<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
213