numeric/interval/transc.hpp

0001 /* Boost interval/transc.hpp template implementation file
0002  *
0003  * Copyright 2000 Jens Maurer
0004  * Copyright 2002 Hervé Brönnimann, Guillaume Melquiond, Sylvain Pion
0005  *
0006  * Distributed under the Boost Software License, Version 1.0.
0007  * (See accompanying file LICENSE_1_0.txt or
0008  * copy at http://www.boost.org/LICENSE_1_0.txt)
0009  */
0010
0011 #ifndef BOOST_NUMERIC_INTERVAL_TRANSC_HPP
0012 #define BOOST_NUMERIC_INTERVAL_TRANSC_HPP
0013
0014 #include <boost/config.hpp>
0015 #include <boost/numeric/interval/detail/interval_prototype.hpp>
0016 #include <boost/numeric/interval/detail/bugs.hpp>
0017 #include <boost/numeric/interval/detail/test_input.hpp>
0018 #include <boost/numeric/interval/rounding.hpp>
0019 #include <boost/numeric/interval/constants.hpp>
0020 #include <boost/numeric/interval/arith.hpp>
0021 #include <boost/numeric/interval/arith2.hpp>
0022 #include <algorithm>
0023
0024 namespace boost {
0025 namespace numeric {
0026
0027 template<class T, class Policies> inline
0028 interval<T, Policies> exp(const interval<T, Policies>& x)
0029 {
0030   typedef interval<T, Policies> I;
0031   if (interval_lib::detail::test_input(x))
0032     return I::empty();
0033   typename Policies::rounding rnd;
0034   return I(rnd.exp_down(x.lower()), rnd.exp_up(x.upper()), true);
0035 }
0036
0037 template<class T, class Policies> inline
0038 interval<T, Policies> log(const interval<T, Policies>& x)
0039 {
0040   typedef interval<T, Policies> I;
0041   if (interval_lib::detail::test_input(x) ||
0042       !interval_lib::user::is_pos(x.upper()))
0043     return I::empty();
0044   typename Policies::rounding rnd;
0045   typedef typename Policies::checking checking;
0046   T l = !interval_lib::user::is_pos(x.lower())
0047              ? checking::neg_inf() : rnd.log_down(x.lower());
0048   return I(l, rnd.log_up(x.upper()), true);
0049 }
0050
0051 template<class T, class Policies> inline
0052 interval<T, Policies> cos(const interval<T, Policies>& x)
0053 {
0054   if (interval_lib::detail::test_input(x))
0055     return interval<T, Policies>::empty();
0056   typename Policies::rounding rnd;
0057   typedef interval<T, Policies> I;
0058   typedef typename interval_lib::unprotect<I>::type R;
0059
0060   // get lower bound within [0, pi]
0061   const R pi2 = interval_lib::pi_twice<R>();
0062   R tmp = fmod((const R&)x, pi2);
0063   if (width(tmp) >= pi2.lower())
0064     return I(static_cast<T>(-1), static_cast<T>(1), true); // we are covering a full period
0065   if (tmp.lower() >= interval_lib::constants::pi_upper<T>())
0066     return -cos(tmp - interval_lib::pi<R>());
0067   T l = tmp.lower();
0068   T u = tmp.upper();
0069
0070   BOOST_USING_STD_MIN();
0071   // separate into monotone subintervals
0072   if (u <= interval_lib::constants::pi_lower<T>())
0073     return I(rnd.cos_down(u), rnd.cos_up(l), true);
0074   else if (u <= pi2.lower())
0075     return I(static_cast<T>(-1), rnd.cos_up(min BOOST_PREVENT_MACRO_SUBSTITUTION(rnd.sub_down(pi2.lower(), u), l)), true);
0076   else
0077     return I(static_cast<T>(-1), static_cast<T>(1), true);
0078 }
0079
0080 template<class T, class Policies> inline
0081 interval<T, Policies> sin(const interval<T, Policies>& x)
0082 {
0083   typedef interval<T, Policies> I;
0084   if (interval_lib::detail::test_input(x))
0085     return I::empty();
0086   typename Policies::rounding rnd;
0087   typedef typename interval_lib::unprotect<I>::type R;
0088   I r = cos((const R&)x - interval_lib::pi_half<R>());
0089   (void)&rnd;
0090   return r;
0091 }
0092
0093 template<class T, class Policies> inline
0094 interval<T, Policies> tan(const interval<T, Policies>& x)
0095 {
0096   typedef interval<T, Policies> I;
0097   if (interval_lib::detail::test_input(x))
0098     return I::empty();
0099   typename Policies::rounding rnd;
0100   typedef typename interval_lib::unprotect<I>::type R;
0101
0102   // get lower bound within [-pi/2, pi/2]
0103   const R pi = interval_lib::pi<R>();
0104   R tmp = fmod((const R&)x, pi);
0105   const T pi_half_d = interval_lib::constants::pi_half_lower<T>();
0106   if (tmp.lower() >= pi_half_d)
0107     tmp -= pi;
0108   if (tmp.lower() <= -pi_half_d || tmp.upper() >= pi_half_d)
0109     return I::whole();
0110   return I(rnd.tan_down(tmp.lower()), rnd.tan_up(tmp.upper()), true);
0111 }
0112
0113 template<class T, class Policies> inline
0114 interval<T, Policies> asin(const interval<T, Policies>& x)
0115 {
0116   typedef interval<T, Policies> I;
0117   if (interval_lib::detail::test_input(x)
0118      || x.upper() < static_cast<T>(-1) || x.lower() > static_cast<T>(1))
0119     return I::empty();
0120   typename Policies::rounding rnd;
0121   T l = (x.lower() <= static_cast<T>(-1))
0122              ? -interval_lib::constants::pi_half_upper<T>()
0123              : rnd.asin_down(x.lower());
0124   T u = (x.upper() >= static_cast<T>(1) )
0125              ?  interval_lib::constants::pi_half_upper<T>()
0126              : rnd.asin_up  (x.upper());
0127   return I(l, u, true);
0128 }
0129
0130 template<class T, class Policies> inline
0131 interval<T, Policies> acos(const interval<T, Policies>& x)
0132 {
0133   typedef interval<T, Policies> I;
0134   if (interval_lib::detail::test_input(x)
0135      || x.upper() < static_cast<T>(-1) || x.lower() > static_cast<T>(1))
0136     return I::empty();
0137   typename Policies::rounding rnd;
0138   T l = (x.upper() >= static_cast<T>(1) )
0139           ? static_cast<T>(0)
0140           : rnd.acos_down(x.upper());
0141   T u = (x.lower() <= static_cast<T>(-1))
0142           ? interval_lib::constants::pi_upper<T>()
0143           : rnd.acos_up  (x.lower());
0144   return I(l, u, true);
0145 }
0146
0147 template<class T, class Policies> inline
0148 interval<T, Policies> atan(const interval<T, Policies>& x)
0149 {
0150   typedef interval<T, Policies> I;
0151   if (interval_lib::detail::test_input(x))
0152     return I::empty();
0153   typename Policies::rounding rnd;
0154   return I(rnd.atan_down(x.lower()), rnd.atan_up(x.upper()), true);
0155 }
0156
0157 template<class T, class Policies> inline
0158 interval<T, Policies> sinh(const interval<T, Policies>& x)
0159 {
0160   typedef interval<T, Policies> I;
0161   if (interval_lib::detail::test_input(x))
0162     return I::empty();
0163   typename Policies::rounding rnd;
0164   return I(rnd.sinh_down(x.lower()), rnd.sinh_up(x.upper()), true);
0165 }
0166
0167 template<class T, class Policies> inline
0168 interval<T, Policies> cosh(const interval<T, Policies>& x)
0169 {
0170   typedef interval<T, Policies> I;
0171   if (interval_lib::detail::test_input(x))
0172     return I::empty();
0173   typename Policies::rounding rnd;
0174   if (interval_lib::user::is_neg(x.upper()))
0175     return I(rnd.cosh_down(x.upper()), rnd.cosh_up(x.lower()), true);
0176   else if (!interval_lib::user::is_neg(x.lower()))
0177     return I(rnd.cosh_down(x.lower()), rnd.cosh_up(x.upper()), true);
0178   else
0179     return I(static_cast<T>(1), rnd.cosh_up(-x.lower() > x.upper() ? x.lower() : x.upper()), true);
0180 }
0181
0182 template<class T, class Policies> inline
0183 interval<T, Policies> tanh(const interval<T, Policies>& x)
0184 {
0185   typedef interval<T, Policies> I;
0186   if (interval_lib::detail::test_input(x))
0187     return I::empty();
0188   typename Policies::rounding rnd;
0189   return I(rnd.tanh_down(x.lower()), rnd.tanh_up(x.upper()), true);
0190 }
0191
0192 template<class T, class Policies> inline
0193 interval<T, Policies> asinh(const interval<T, Policies>& x)
0194 {
0195   typedef interval<T, Policies> I;
0196   if (interval_lib::detail::test_input(x))
0197     return I::empty();
0198   typename Policies::rounding rnd;
0199   return I(rnd.asinh_down(x.lower()), rnd.asinh_up(x.upper()), true);
0200 }
0201
0202 template<class T, class Policies> inline
0203 interval<T, Policies> acosh(const interval<T, Policies>& x)
0204 {
0205   typedef interval<T, Policies> I;
0206   if (interval_lib::detail::test_input(x) || x.upper() < static_cast<T>(1))
0207     return I::empty();
0208   typename Policies::rounding rnd;
0209   T l = x.lower() <= static_cast<T>(1) ? static_cast<T>(0) : rnd.acosh_down(x.lower());
0210   return I(l, rnd.acosh_up(x.upper()), true);
0211 }
0212
0213 template<class T, class Policies> inline
0214 interval<T, Policies> atanh(const interval<T, Policies>& x)
0215 {
0216   typedef interval<T, Policies> I;
0217   if (interval_lib::detail::test_input(x)
0218       || x.upper() < static_cast<T>(-1) || x.lower() > static_cast<T>(1))
0219     return I::empty();
0220   typename Policies::rounding rnd;
0221   typedef typename Policies::checking checking;
0222   T l = (x.lower() <= static_cast<T>(-1))
0223              ? checking::neg_inf() : rnd.atanh_down(x.lower());
0224   T u = (x.upper() >= static_cast<T>(1) )
0225              ? checking::pos_inf() : rnd.atanh_up  (x.upper());
0226   return I(l, u, true);
0227 }
0228
0229 } // namespace numeric
0230 } // namespace boost
0231
0232 #endif // BOOST_NUMERIC_INTERVAL_TRANSC_HPP