numerix: include/numerix/interval.hpp Source File

00001 
00002 /******************************************************************************
00003 * MODULE     : interval.hpp
00004 * DESCRIPTION: Intervals
00005 * COPYRIGHT  : (C) 2006  Joris van der Hoeven
00006 *******************************************************************************
00007 * This software falls under the GNU general public license and comes WITHOUT
00008 * ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
00009 * If you don't have this file, write to the Free Software Foundation, Inc.,
00010 * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00011 ******************************************************************************/
00012 
00013 #ifndef __MMX_INTERVAL_HPP
00014 #define __MMX_INTERVAL_HPP
00015 #include <basix/function.hpp>
00016 #include <numerix/rounded.hpp>
00017 namespace mmx {
00018 #define TMPL_DEF template<typename C,typename V=std_interval>
00019 #define TMPL template<typename C,typename V>
00020 #define Interval interval<C,V>
00021 #define Up V::template rnd_helper<C>::UV
00022 #define Down V::template rnd_helper<C>::DV
00023 TMPL class interval;
00024 TMPL inline C lower (const Interval& z);
00025 TMPL inline C upper (const Interval& z);
00026 
00027 /******************************************************************************
00028 * Interval variants
00029 ******************************************************************************/
00030 
00031 struct std_interval {
00032   template<typename C>
00033   struct rnd_helper {
00034     typedef Round_up(C) UV;
00035     typedef Round_down(C) DV;
00036   };
00037 };
00038 
00039 /******************************************************************************
00040 * The interval type
00041 ******************************************************************************/
00042 
00043 TMPL_DEF
00044 class interval {
00045 MMX_ALLOCATORS
00046 private:
00047   C l;
00048   C r;
00049 
00050 public:      
00051   typedef C boundary_type;
00052 
00053   inline interval (): l (0), r (0) {}
00054   template<typename T> inline interval (const T& x):
00055     l (Down::template as<C> (x)),
00056     r (Up::template as<C> (x)) {}
00057   template<typename T> inline interval (const interval<T,V>& x):
00058     l (Down::template as<C> (lower (x))),
00059     r (Up::template as<C> (upper (x))) {}
00060   template<typename LT,typename UT> inline
00061   interval (const LT& l2, const UT& u2):
00062     l (Down::template as<C> (l2)),
00063     r (Up::template as<C> (u2)) {}
00064   inline interval (const C& x): l (x), r (x) {}
00065   inline interval (const C& l2, const C& r2): l (l2), r (r2) {}
00066   friend C lower LESSGTR (const Interval& x);
00067   friend C upper LESSGTR (const Interval& x);
00068 
00069   static inline Interval nan ();
00070   static inline Interval maximal ();
00071   static inline Interval minimal ();
00072   static inline Interval infinity ();
00073   static inline Interval fuzzy ();
00074   static inline Interval log2 ();
00075   static inline Interval pi ();
00076   static inline Interval euler ();
00077 
00078   inline Interval& operator += (const Interval& x);
00079   inline Interval& operator -= (const Interval& x);
00080   inline Interval& operator <<= (const xint& shift);
00081   inline Interval& operator >>= (const xint& shift);
00082 };
00083 DEFINE_UNARY_FORMAT_1 (interval)
00084 
00085 UNARY_RETURN_TYPE(TMPL,center_op,Interval,C);
00086 UNARY_RETURN_TYPE(TMPL,radius_op,Interval,C);
00087 TMPL struct rounding_helper<Interval >: public rounding_helper<C> {};
00088 
00089 /******************************************************************************
00090 * Standard constructors and accessors
00091 ******************************************************************************/
00092 
00093 template<typename C, typename V, typename C2, typename R2>
00094 struct make_ball_helper<Interval,C2,R2> {
00095   static inline Interval val (const C2& c, const R2& r) {
00096     typedef Round_up(C2) Rnd_Up;
00097     typedef Round_down(C2) Rnd_Down;
00098     return Interval (Rnd_Down::template sub_as<C2> (c, r),
00099                      Rnd_Up::template add_as<C2> (c, r));
00100   }
00101 };
00102 
00103 template<typename C, typename V, typename C2>
00104 struct make_interval_helper<Interval,C2> {
00105   static inline Interval val (const C2& l, const C2& r) {
00106     return Interval (l, r); }
00107 };
00108 
00109 TMPL inline C lower (const Interval& x) { return x.l; }
00110 TMPL inline C upper (const Interval& x) { return x.r; }
00111 TMPL inline C center (const Interval& x) {
00112   return decexp2 (lower (x) + upper (x)); }
00113 TMPL inline C radius (const Interval& x) {
00114   if (is_infinite (x)) return 0;
00115   return Up::add (decexp2 (Up::sub (upper (x), lower (x))),
00116                   multiplicative_error (center (x))); }
00117 TMPL inline format<C> CF (const Interval& x) {
00118   return get_format (lower (x)); }
00119 
00120 TMPL inline Interval abs (const Interval& x);
00121 TMPL inline C abs_down (const Interval& x) { return lower (abs (x)); }
00122 TMPL inline C abs_up (const Interval& x) { return upper (abs (x)); }
00123 TMPL inline C bnd_down (const Interval& x) { return lower (x); }
00124 TMPL inline C bnd_up (const Interval& x) { return upper (x); }
00125 
00126 /******************************************************************************
00127 * Basic routines
00128 ******************************************************************************/
00129 
00130 TMPL inline Interval copy (const Interval& x) {
00131   return Interval (copy (lower (x)), copy (upper (x))); }
00132 TMPL inline Interval duplicate (const Interval& x) {
00133   return Interval (duplicate (lower (x)), duplicate (upper (x))); }
00134 
00135 TMPL inline bool is_exact_zero (const Interval& x) {
00136   return is_exact_zero (lower (x)) && is_exact_zero (upper (x)); }
00137 
00138 template<typename Op, typename C, typename V> nat
00139 unary_hash (const Interval& x) {
00140   nat h= Op::op (lower (x));
00141   return (h<<3) ^ (h<<15) ^ (h>>29) ^ Op::op (upper (x));
00142 }
00143 
00144 template<typename Op, typename C, typename V> bool
00145 binary_test (const Interval& x1, const Interval& x2) {
00146   return Op::op (lower (x1), lower (x2)) && Op::op (upper (x1), upper (x2));
00147 }
00148 
00149 template<typename S1, typename D> interval<D>
00150 map (const function_1<D,Argument(S1) >& fun, const interval<S1>& x) {
00151   return interval<D> (fun (lower (x)), fun (upper (x)));
00152 }
00153 
00154 template<typename S1, typename D> interval<D>
00155 map (D (*fun) (const S1&), const interval<S1>& x) {
00156   return interval<D> (fun (lower (x)), fun (upper (x)));
00157 }
00158 
00159 EXACT_IDENTITY_OP_SUGAR(TMPL,Interval)
00160 HARD_IDENTITY_OP_SUGAR(TMPL,Interval)
00161 TMPL nat hash (const Interval& x) { return 12721; }
00162 
00163 TMPL syntactic
00164 flatten (const Interval& x) {
00165   if (is_nan (x)) return Nan (syntactic);
00166   if (lower (x) == 0 && upper (x) == 0) return 0;
00167   return flatten_range (lower (x), upper (x));
00168   // return flatten ("interval", lower (x), upper (x));
00169 }
00170 
00171 TMPL syntactic
00172 flatten (const Interval& x, xnat dd) {
00173   if (is_nan (x)) return Nan (syntactic);
00174   if (lower (x) == 0 && upper (x) == 0) return 0;
00175   return flatten_range (lower (x), upper (x), dd);
00176 }
00177 
00178 TMPL
00179 struct binary_helper<Interval >: public void_binary_helper<Interval > {
00180   static inline string short_type_name () {
00181     return "I" * Short_type_name (C); }
00182   static inline generic full_type_name () {
00183     return gen ("Interval", Full_type_name (C)); }
00184   static inline generic disassemble (const Interval& x) {
00185     return gen_vec (as<generic> (lower (x)), as<generic> (upper (x))); }
00186   static inline Interval assemble (const generic& v) {
00187     return Interval (as<C> (vector_access (v, 0)),
00188                      as<C> (vector_access (v, 1))); }
00189   static inline void write (const port& out, const Interval& x) {
00190     binary_write<C> (out, lower (x));
00191     binary_write<C> (out, upper (x)); }
00192   static inline Interval read (const port& in) {
00193     C l= binary_read<C> (in);
00194     C r= binary_read<C> (in);
00195     return Interval (l, r); }
00196 };
00197 
00198 /******************************************************************************
00199 * Comparisons
00200 ******************************************************************************/
00201 
00202 TMPL inline bool
00203 intersect (Interval& r, const Interval &a, const Interval &b) {
00204   if ((upper (a) < lower (b)) || (lower (a) > upper (b))) return false;
00205   r = Interval (max (lower (a), lower (b)), min (upper (a), upper (b)));
00206   return true;
00207 }
00208 
00209 TMPL inline bool
00210 is_zero (const Interval& x) {
00211   return is_finite (x) && sign (lower (x)) * sign (upper (x)) <= 0;
00212 }
00213 
00214 TMPL inline bool
00215 is_non_zero (const Interval& x) {
00216   return !is_finite (x) || sign (lower (x)) * sign (upper (x)) > 0;
00217 }
00218 
00219 TMPL inline bool
00220 operator == (const Interval& x1, const Interval& x2) {
00221   if (is_nan (x1) || is_nan (x2)) return is_nan (x1) && is_nan (x2);
00222   return lower (x1) <= upper (x2) && lower (x2) <= upper (x1);
00223 }
00224 
00225 TMPL inline bool
00226 operator != (const Interval& x1, const Interval& x2) {
00227   if (is_nan (x1) || is_nan (x2)) return !(is_nan (x1) && is_nan (x2));
00228   return lower (x1) > upper (x2) || lower (x2) > upper (x1);
00229 }
00230 
00231 TMPL inline bool
00232 operator <= (const Interval& x1, const Interval& x2) {
00233   if (is_nan (x1) || is_nan (x2)) return false;
00234   return lower (x1) <= upper (x2);
00235 }
00236 
00237 TMPL inline bool
00238 operator < (const Interval& x1, const Interval& x2) {
00239   if (is_nan (x1) || is_nan (x2)) return false;
00240   return lower (x2) > upper (x1);
00241 }
00242 
00243 TMPL inline bool
00244 operator >= (const Interval& x1, const Interval& x2) {
00245   if (is_nan (x1) || is_nan (x2)) return false;
00246   return lower (x2) <= upper (x1);
00247 }
00248 
00249 TMPL inline bool
00250 operator > (const Interval& x1, const Interval& x2) {
00251   if (is_nan (x1) || is_nan (x2)) return false;
00252   return lower (x1) > upper (x2);
00253 }
00254 
00255 EQUAL_INT_SUGAR(TMPL,Interval);
00256 COMPARE_INT_SUGAR(TMPL,Interval);
00257 EQUAL_SCALAR_SUGAR(TMPL,Interval,C);
00258 COMPARE_SCALAR_SUGAR(TMPL,Interval,C);
00259 EQUAL_SCALAR_SUGAR_BIS(TMPL,Interval,C);
00260 COMPARE_SCALAR_SUGAR_BIS(TMPL,Interval,C);
00261 
00262 /******************************************************************************
00263 * Other order related functions
00264 ******************************************************************************/
00265 
00266 TMPL inline Interval
00267 abs (const Interval& x) {
00268   if (is_nan (x)) return x;
00269   C r= max (abs (lower (x)), abs (upper (x)));
00270   if (sign (lower (x)) * sign (upper (x)) <= 0) return Interval (C(0), r);
00271   return Interval (min (abs (lower (x)), abs (upper (x))), r);
00272 }
00273 
00274 TMPL Interval
00275 min (const Interval& x1, const Interval& x2) {
00276   return Interval (min (lower (x1), lower (x2)), min (upper (x1), upper (x2)));
00277 }
00278 
00279 TMPL Interval
00280 max (const Interval& x1, const Interval& x2) {
00281   return Interval (max (lower (x1), lower (x2)), max (upper (x1), upper (x2)));
00282 }
00283 
00284 TMPL Interval
00285 floor (const Interval& x) {
00286   return Interval (Down::floor (lower (x)), Up::floor (upper (x)));
00287 }
00288 
00289 TMPL Interval
00290 trunc (const Interval& x) {
00291   return Interval (Down::trunc (lower (x)), Up::trunc (upper (x)));
00292 }
00293 
00294 TMPL Interval
00295 ceil (const Interval& x) {
00296   return Interval (Down::ceil (lower (x)), Up::ceil (upper (x)));
00297 }
00298 
00299 TMPL Interval
00300 round (const Interval& x) {
00301   return Interval (Down::round (lower (x)), Up::round (upper (x)));
00302 }
00303 
00304 /******************************************************************************
00305 * Constants
00306 ******************************************************************************/
00307 
00308 TMPL inline void set_default (Interval& x) {
00309   x= Interval (Default (C), Default (C)); }
00310 TMPL inline void set_nan (Interval& x) {
00311   x= Interval (Nan (C), Nan (C)); }
00312 TMPL inline void set_maximal (Interval& x) {
00313   x= Interval (Maximal (C), Maximal (C)); }
00314 TMPL inline void set_minimal (Interval& x) {
00315   x= Interval (Minimal (C), Minimal (C)); }
00316 TMPL inline void set_infinity (Interval& x) {
00317   x= Interval (Infinity (C), Infinity (C)); }
00318 TMPL inline void set_fuzz (Interval& x) {
00319   x= Interval (Minimal (C), Maximal (C)); }
00320 TMPL inline void set_smallest (Interval& x) {
00321   x= Interval (Smallest (C), Smallest (C)); }
00322 TMPL inline void set_largest (Interval& x) {
00323   x= Interval (Largest (C), Largest (C)); }
00324 TMPL inline void set_accuracy (Interval& x) {
00325   x= Interval (Accuracy (C), Accuracy (C)); }
00326 TMPL inline void set_log2 (Interval& x) {
00327   x= Interval (Down::template log2<C>(), Up::template log2<C>()); }
00328 TMPL inline void set_pi (Interval& x) {
00329   x= Interval (Down::template pi<C>(), Up::template pi<C>()); }
00330 TMPL inline void set_euler (Interval& x) {
00331   x= Interval (Down::template euler<C>(), Up::template euler<C>()); }
00332 TMPL inline void set_catalan (Interval& x) {
00333   x= Interval (Down::template catalan<C>(), Up::template catalan<C>()); }
00334 TMPL inline Interval times_infinity (const Interval& x) {
00335   return x * infinity_cst<Interval > (); }
00336 
00337 /******************************************************************************
00338 * Basic arithmetic
00339 ******************************************************************************/
00340 
00341 TMPL inline Interval
00342 operator - (const Interval& x) {
00343   return Interval (-upper (x), -lower (x));
00344 }
00345 
00346 TMPL inline Interval
00347 operator + (const Interval& x1, const Interval& x2) {
00348   return Interval (Down::add (lower (x1), lower (x2)),
00349                    Up::add (upper (x1), upper (x2)));
00350 }
00351 
00352 TMPL inline Interval
00353 operator + (const Interval& x1, const C& x2) {
00354   return Interval (Down::add (lower (x1), x2), Up::add (upper (x1), x2));
00355 }
00356 
00357 TMPL inline Interval
00358 operator + (const C& x1, const Interval& x2) {
00359   return Interval (Down::add (x1, lower (x2)), Up::add (x1, upper (x2)));
00360 }
00361 
00362 TMPL inline Interval
00363 operator - (const Interval& x1, const Interval& x2) {
00364   return Interval (Down::sub (lower (x1), upper (x2)),
00365                    Up::sub (upper (x1), lower (x2)));
00366 }
00367 
00368 TMPL inline Interval
00369 operator - (const Interval& x1, const C& x2) {
00370   return Interval (Down::add (lower (x1), x2), Up::add (upper (x1), x2));
00371 }
00372 
00373 TMPL inline Interval
00374 operator - (const C& x1, const Interval& x2) {
00375   return Interval (Down::add (x1, upper (x2)), Up::add (x1, lower (x2)));
00376 }
00377 
00378 TMPL Interval
00379 operator * (const Interval& x1, const Interval& x2) {
00380   if (lower (x1) >= 0)
00381     return Interval (lower (x2) >= 0?
00382                      Down::mul (lower (x1), lower (x2)):
00383                      Down::mul (upper (x1), lower (x2)),
00384                      upper (x2) >= 0?
00385                      Up::mul (upper (x1), upper (x2)):
00386                      Up::mul (lower (x1), upper (x2)));
00387   else if (upper (x1) <= 0)
00388     return Interval (upper (x2) >= 0?
00389                      Down::mul (lower (x1), upper (x2)):
00390                      Down::mul (upper (x1), upper (x2)),
00391                      lower (x2) >= 0?
00392                      Up::mul (upper (x1), lower (x2)):
00393                      Up::mul (lower (x1), lower (x2)));
00394   else if (lower (x2) >= 0)
00395     return Interval (Down::mul (lower (x1), upper (x2)),
00396                      Up::mul (upper (x1), upper (x2)));
00397   else if (upper (x2) <= 0)
00398     return Interval (Down::mul (upper (x1), lower (x2)),
00399                      Up::mul (lower (x1), lower (x2)));
00400   else
00401     return Interval (min (Down::mul (lower (x1), upper (x2)),
00402                           Down::mul (upper (x1), lower (x2))),
00403                      max (Up::mul (lower (x1), lower (x2)),
00404                           Up::mul (upper (x1), upper (x2))));
00405 }
00406 
00407 TMPL inline Interval operator * (const Interval& x1, const C& x2) {
00408   return x1 * Interval (x2); }
00409 TMPL inline Interval operator * (const C& x1, const Interval& x2) {
00410   return Interval (x1) * x2; }
00411 
00412 TMPL inline Interval
00413 square (const Interval& x) {
00414   if (lower (x) >= 0)
00415     return Interval (Down::mul (lower (x), lower (x)),
00416                      Up::mul (upper (x), upper (x)));
00417   else if (upper (x) <= 0)
00418     return Interval (Down::mul (upper (x), upper (x)),
00419                      Up::mul (lower (x), lower (x)));
00420   else return Interval (Down::mul (lower (x), upper (x)),
00421                         max (Up::mul (lower (x), lower (x)),
00422                              Up::mul (upper (x), upper (x))));
00423 }
00424 
00425 TMPL Interval
00426 operator / (const Interval& x1, const Interval& x2) {
00427   if (is_zero (x2)) return Nan (Interval);
00428   else if (lower (x2) > 0)
00429     return Interval (lower (x1) > 0?
00430                      Down::div (lower (x1), upper (x2)):
00431                      Down::div (lower (x1), lower (x2)),
00432                      upper (x1) > 0?
00433                      Up::div (upper (x1), lower (x2)):
00434                      Up::div (upper (x1), upper (x2)));
00435   else if (upper (x2) < 0)
00436     return Interval (upper (x1) > 0?
00437                      Down::div (upper (x1), upper (x2)):
00438                      Down::div (upper (x1), lower (x2)),
00439                      lower (x1) > 0?
00440                      Up::div (lower (x1), lower (x2)):
00441                      Up::div (lower (x1), upper (x2)));
00442   else return Nan (Interval);
00443 }
00444 
00445 TMPL inline Interval operator / (const Interval& x1, const C& x2) {
00446   return x1 / Interval (x2); }
00447 TMPL inline Interval operator / (const C& x1, const Interval& x2) {
00448   return Interval (x1) / x2; }
00449 
00450 ARITH_INT_SUGAR(TMPL,Interval);
00451 
00452 /******************************************************************************
00453 * In place arithmetic operations
00454 ******************************************************************************/
00455 
00456 TMPL inline Interval&
00457 Interval::operator += (const Interval& x) {
00458   l= Down::add (l, lower (x));
00459   r= Up::add (r, upper (x));
00460   return *this;
00461 }
00462 
00463 TMPL inline Interval&
00464 Interval::operator -= (const Interval& x) {
00465   l= Down::sub (l, upper (x));
00466   r= Up::sub (r, lower (x));
00467   return *this;
00468 }
00469 
00470 /******************************************************************************
00471 * Elementary functions
00472 ******************************************************************************/
00473 
00474 TMPL Interval
00475 sqrt (const Interval& x) {
00476   return Interval (Down::sqrt (lower (x)), Up::sqrt (upper (x)));
00477 }
00478 
00479 TMPL Interval
00480 hypot (const Interval& x, const Interval& y) {
00481   Interval z= square (x) + square (y);
00482   if (z > 0) return sqrt (z);
00483   else return Interval (0, Up::sqrt (upper (z)));
00484 }
00485 
00486 TMPL Interval
00487 exp (const Interval& x) {
00488   return Interval (Down::exp (lower (x)), Up::exp (upper (x)));
00489 }
00490 
00491 TMPL Interval
00492 log (const Interval& x) {
00493   return Interval (Down::log (lower (x)), Up::log (upper (x)));
00494 }
00495 
00496 TMPL inline Interval
00497 pow (const Interval& x, const Interval& y) {
00498   return exp (y * log (x));
00499 }
00500 
00501 TMPL Interval
00502 cos (const Interval& x) {
00503   C pi_lo= promote (314159, lower (x)) / promote (100000, lower (x));
00504   if (Up::sub (upper (x), lower (x)) >= pi_lo)
00505     return Interval (promote (-1, lower (x)), promote (1, lower (x)));
00506   else {
00507     int s1= -sign (sin (lower (x)));
00508     int s2= -sign (sin (upper (x)));
00509     if (s1 >= 0 && s2 >= 0)
00510       return Interval (Down::cos (lower (x)), Up::cos (upper (x)));
00511     else if (s1 <= 0 && s2 <= 0)
00512       return Interval (Down::cos (upper (x)), Down::cos (lower (x)));
00513     else if (s1 >= 0)
00514       return Interval (min (Down::cos (lower (x)), Down::cos (upper (x))),
00515                        promote (1, lower (x)));
00516     else
00517       return Interval (promote (-1, lower (x)),
00518                        max (Down::cos (lower (x)), Down::cos (upper (x))));
00519   }
00520 }
00521 
00522 TMPL Interval
00523 sin (const Interval& x) {
00524   C pi_lo= promote (314159, lower (x)) / promote (100000, lower (x));
00525   if (Up::sub (upper (x), lower (x)) >= pi_lo)
00526     return Interval (promote (-1, lower (x)), promote (1, lower (x)));
00527   else {
00528     int s1= sign (cos (lower (x)));
00529     int s2= sign (cos (upper (x)));
00530     if (s1 >= 0 && s2 >= 0)
00531       return Interval (Down::sin (lower (x)), Up::sin (upper (x)));
00532     else if (s1 <= 0 && s2 <= 0)
00533       return Interval (Down::sin (upper (x)), Down::sin (lower (x)));
00534     else if (s1 >= 0)
00535       return Interval (min (Down::sin (lower (x)), Down::sin (upper (x))),
00536                        promote (1, lower (x)));
00537     else
00538       return Interval (promote (-1, lower (x)),
00539                        max (Down::sin (lower (x)), Down::sin (upper (x))));
00540   }
00541 }
00542 
00543 TMPL Interval
00544 tan (const Interval& x) {
00545   return sin (x) / cos (x);
00546 }
00547 
00548 TMPL inline Interval
00549 cosh (const Interval& x) {
00550   if (is_nan (x)) return Nan (Interval);
00551   else if (lower (x) >= 0)
00552     return Interval (Down::cosh (lower (x)), Up::cosh (upper (x)));
00553   else if (upper (x) <= 0)
00554     return Interval (Down::cosh (upper (x)), Up::cosh (lower (x)));
00555   else return Interval (1, max (Up::cosh (lower (x)), Up::cosh (upper (x))));
00556 }
00557 
00558 TMPL inline Interval
00559 sinh (const Interval& x) {
00560   return Interval (Down::sinh (lower (x)), Up::sinh (upper (x)));
00561 }
00562 
00563 TMPL Interval
00564 tanh (const Interval& x) {
00565   return sinh (x) / cosh (x);
00566 }
00567 
00568 TMPL Interval
00569 acos (const Interval& x) {
00570   return Interval (Down::acos (upper (x)), Up::acos (lower (x)));
00571 }
00572 
00573 TMPL Interval
00574 asin (const Interval& x) {
00575   return Interval (Down::asin (lower (x)), Up::asin (upper (x)));
00576 }
00577 
00578 TMPL Interval
00579 atan (const Interval& x) {
00580   return Interval (Down::atan (lower (x)), Up::atan (upper (x)));
00581 }
00582 
00583 TMPL Interval
00584 atan2 (const Interval& y, const Interval& x) {
00585   if (sign (lower (x)) > 0) return atan (y/x);
00586   if (sign (lower (y)) > 0) return atan (-x/y) + (Pi (Interval) >> 1);
00587   if (sign (upper (y)) < 0) return atan (-x/y) - (Pi (Interval) >> 1);
00588   if (sign (upper (x)) < 0) return atan (y/x) + Pi (Interval);
00589   C pi= Up::template pi<C> ();
00590   return Interval (-pi, pi);
00591 }
00592 
00593 INV_TRIGO_SUGAR(TMPL,Interval)
00594 INV_HYPER_SUGAR(TMPL,Interval)
00595 ARG_HYPER_SUGAR(TMPL,Interval)
00596 
00597 /******************************************************************************
00598 * Floating point related functions
00599 ******************************************************************************/
00600 
00601 TMPL inline bool is_finite (const Interval& x) {
00602   return is_finite (lower (x)) && is_finite (upper (x)); }
00603 TMPL inline bool is_infinite (const Interval& x) {
00604   return is_infinite (lower (x)) && upper (x) == lower (x); }
00605 TMPL inline bool is_fuzz (const Interval& x) {
00606   return is_infinite (lower (x)) && is_infinite (upper (x)) &&
00607     lower (x) != upper (x); }
00608 TMPL inline bool is_nan (const Interval& x) {
00609   return is_nan (lower (x)) || is_nan (upper (x)); }
00610 TMPL inline bool is_reliable (const Interval& x) {
00611   (void) x; return true; }
00612 
00613 TMPL inline Interval change_precision (const Interval& x, xnat prec) {
00614   return Interval (change_precision (lower (x), prec),
00615                    change_precision (upper (x), prec)); }
00616 TMPL inline xnat precision (const Interval& x) {
00617   return min (precision (lower (x)), precision (upper (x))); }
00618 
00619 TMPL inline xint exponent (const Interval& x) {
00620   return max (exponent (lower (x)), exponent (upper (x))); }
00621 TMPL inline double magnitude (const Interval& x) {
00622   return max (magnitude (lower (x)), magnitude (upper (x))); }
00623 TMPL inline Interval operator << (const Interval& x, const xint& shift) {
00624   return Interval (incexp2 (lower (x), shift), incexp2 (upper (x), shift)); }
00625 TMPL inline Interval operator >> (const Interval& x, const xint& shift) {
00626   return Interval (decexp2 (lower (x), shift), decexp2 (upper (x), shift)); }
00627 
00628 TMPL inline Interval&
00629 Interval::operator <<= (const xint& shift) {
00630   incexp2_assign (l, shift);
00631   incexp2_assign (r, shift);
00632   return *this;
00633 }
00634 
00635 TMPL inline Interval&
00636 Interval::operator >>= (const xint& shift) {
00637   decexp2_assign (l, shift);
00638   decexp2_assign (r, shift);
00639   return *this;
00640 }
00641 
00642 TMPL inline Interval sharpen (const Interval& x) {
00643   return Interval (center (x)); }
00644 TMPL inline Interval blur (const Interval& x, const C& r) {
00645   return x + Interval (-r, r); }
00646 TMPL inline Interval blur (const Interval& x, const Interval& r) {
00647   return blur (x, abs_up (r)); }
00648 
00649 #undef TMPL_DEF
00650 #undef TMPL
00651 #undef Interval
00652 #undef Up
00653 #undef Down
00654 } // namespace mmx
00655 #endif // __MMX_INTERVAL_HPP