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