/Users/mourrain/Devel/mmx/numerix/include/numerix/ball_infinities.hpp

Go to the documentation of this file.

/******************************************************************************
* MODULE     : ball_infinities.hpp
* DESCRIPTION: Refined error estimates when computing with infinities
* COPYRIGHT  : (C) 2006  Joris van der Hoeven
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#ifndef __MMX_BALL_INFINITIES_HPP
#define __MMX_BALL_INFINITIES_HPP
#include <numerix/ball_rounded.hpp>
namespace mmx {
#define TMPL template<typename C,typename R,typename V>
#define Ball ball<C,R,V>

/******************************************************************************
* The ball_infinities<V> variant
******************************************************************************/

template<typename V>
struct ball_infinities: public V {};

template<typename F, typename V, typename W>
struct implementation<F,V,ball_infinities<W> >:
  public implementation<F,V,W> {};

template<typename V>
struct ball_variant_helper<floating<V> > {
  typedef ball_infinities<ball_rounded> BV;
};

/******************************************************************************
* Additive operations
******************************************************************************/

template<typename IV, typename BV>
struct implementation<ball_additive,IV,ball_infinities<BV> > {
  typedef implementation<ball_additive,BV> Fallback;
  typedef implementation<ball_rounding,BV> Rnd;
  
TMPL static void
neg (Ball& d, const Ball& z) {
  if (is_finite (z) || is_infinite (z)) Fallback::neg (d, z);
  else d= z;
}

TMPL static void
add (Ball& d, const Ball& z1, const Ball& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::add (d, z1, z2);
  else if (is_nan (z1) || is_nan (z2)) d= Nan (Ball);
  else if (is_fuzz (z1) || is_fuzz (z2)) d= Fuzz (Ball);
  else if (is_infinite (z1) && is_infinite (z2)) {
    if (Projective (C) || center (z1) != center (z2)) d= Fuzz (Ball);
    else d= z1;
  }
  else d= (is_infinite (z1)? z1: z2);
}

TMPL static void
add (Ball& d, const C& z1, const Ball& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::add (d, z1, z2);
  else add (d, Ball (z1), z2);
}

TMPL static void
add (Ball& d, const Ball& z1, const C& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::add (d, z1, z2);
  else add (d, z1, Ball (z2));
}

TMPL static void
sub (Ball& d, const Ball& z1, const Ball& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::sub (d, z1, z2);
  else if (is_nan (z1) || is_nan (z2)) d= Nan (Ball);
  else if (is_fuzz (z1) || is_fuzz (z2)) d= Fuzz (Ball);
  else if (is_infinite (z1) && is_infinite (z2)) {
    if (Projective (C) || center (z1) == center (z2)) d= Fuzz (Ball);
    else d= z1;
  }
  else d= (is_infinite (z1)? z1: -z2);
}

TMPL static void
sub (Ball& d, const C& z1, const Ball& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::sub (d, z1, z2);
  else sub (d, Ball (z1), z2);
}

TMPL static void
sub (Ball& d, const Ball& z1, const C& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::sub (d, z1, z2);
  else sub (d, z1, Ball (z2));
}

}; // implementation<ball_additive,IV,ball_infinities<BV> >

/******************************************************************************
* Multiplicative operations
******************************************************************************/

template<typename IV, typename BV>
struct implementation<ball_multiplicative,IV,ball_infinities<BV> > {
  typedef implementation<ball_multiplicative,BV> Fallback;
  typedef implementation<ball_rounding,BV> Rnd;

TMPL static void
mul (Ball& d, const Ball& z1, const Ball& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::mul (d, z1, z2);
  else if (is_nan (z1) || is_nan (z2)) d= Nan (Ball);
  else if (is_fuzz (z1) || is_fuzz (z2)) d= Fuzz (Ball);
  else if (is_infinite (z1)) {
    if (is_zero (z2)) d= Fuzz (Ball);
    else if (Projective (C)) d= Infinity (Ball);
    else {
      center (d)= center (z1) * center (z2);
      radius (d)= 0;
    }
  }
  else {
    if (is_zero (z1)) d= Fuzz (Ball);
    else if (Projective (C)) d= Infinity (Ball);
    else {
      center (d)= center (z1) * center (z2);
      radius (d)= 0;
    }
  }
}

TMPL static void
mul (Ball& d, const C& z1, const Ball& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::mul (d, z1, z2);
  else mul (d, Ball (z1), z2);
}

TMPL static void
mul (Ball& d, const Ball& z1, const C& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::mul (d, z1, z2);
  else mul (d, z1, Ball (z2));
}

TMPL static void
square (Ball& d, const Ball& z) {
  if (is_finite (z)) Fallback::square (d, z);
  else if (is_infinite (z)) d= Infinity (Ball);
  else d= z;
}

TMPL static void
invert (Ball& d, const Ball& z) {
  if (is_zero (z)) d= Nan (Ball);
  else if (is_finite (z)) Fallback::invert (d, z);
  else if (is_nan (z)) d= Nan (Ball);
  else d= 0;
}

TMPL static void
invert (Ball& d, const C& z) {
  if (is_zero (z)) d= Nan (Ball);
  else if (is_finite (z)) Fallback::invert (d, z);
  else if (is_nan (z)) d= Nan (Ball);
  else d= 0;
}

TMPL static void
div (Ball& d, const Ball& z1, const Ball& z2) {
  if (is_zero (z2)) d= Nan (Ball);
  else if (is_finite (z1) && is_finite (z2)) Fallback::div (d, z1, z2);
  else if (is_nan (z1) || is_nan (z2)) d= Nan (Ball);
  else if (is_fuzz (z1)) d= Fuzz (Ball);
  else if (is_infinite (z1)) {
    if (is_infinite (z2)) d= Fuzz (Ball);
    else if (Projective (C)) d= Infinity (Ball);
    else {
      center (d)= center (z1) / center (z2);
      radius (d)= 0;
    }
  }
  else d= 0;
}

TMPL static void
div (Ball& d, const C& z1, const Ball& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::div (d, z1, z2);
  else div (d, Ball (z1), z2);
}

TMPL static void
div (Ball& d, const Ball& z1, const C& z2) {
  if (is_finite (z1) && is_finite (z2)) Fallback::div (d, z1, z2);
  else div (d, z1, Ball (z2));
}

}; // implementation<ball_multiplicative,IV,ball_infinities<BV> >

/******************************************************************************
* Square root related functions
******************************************************************************/

template<typename IV, typename BV>
struct implementation<ball_root,IV,ball_infinities<BV> > {
  typedef implementation<ball_root,BV> Fallback;
  typedef implementation<ball_rounding,BV> Rnd;

TMPL static void
sqrt (Ball& d, const Ball& z) {
  if (is_negative_or_zero (z)) d= Nan (Ball);
  else if (is_finite (z)) Fallback::sqrt (d, z);
  else if (is_infinite (z)) d= Maximal (Ball);
  else d= Nan (Ball);
}

TMPL static void
hypot (Ball& d, const Ball& x, const Ball& y) {
  if (is_finite (x) && is_finite (y)) Fallback::hypot (d, x, y);
  else if (is_nan (x) || is_nan (y)) d= Nan (Ball);
  else if (is_fuzz (x) || is_fuzz (y)) d= Fuzz (Ball);
  else d= Maximal (Ball);
}

}; // implementation<ball_root,IV,ball_infinities<BV> >

/******************************************************************************
* Absolute value
******************************************************************************/

#define ATMPL template<typename C,typename R,typename V,typename C2,typename V2>
#define ABall ball<C2,R,V2>

template<typename IV, typename BV>
struct implementation<ball_abs,IV,ball_infinities<BV> > {
  typedef implementation<ball_abs,BV> Fallback;
  typedef implementation<ball_rounding,BV> Rnd;

ATMPL static inline void
abs (ABall& d, const Ball& z) {
  if (is_finite (z)) Fallback::abs (d, z);
  else if (is_infinite (z)) d= Infinity (ABall);
  else if (is_fuzz (z)) d= Fuzz (ABall);
  else d= Nan (ABall);
}

TMPL static inline void
abs (Ball& d, const Ball& z) {
  Fallback::abs (d, z);
}

}; // implementation<ball_abs,IV,ball_infinities<BV> >

#undef ATMPL
#undef ABall

/******************************************************************************
* Shift operations
******************************************************************************/

template<typename IV, typename BV>
struct implementation<ball_shift,IV,ball_infinities<BV> > {
  typedef implementation<ball_shift,BV> Fallback;
  typedef implementation<ball_rounding,BV> Rnd;

TMPL static inline void
shiftl (Ball& d, const xint& shift) {
  if (is_finite (d)) Fallback::shiftl (d, shift);
}

TMPL static inline void
shiftl (Ball& d, const Ball& z, const xint& shift) {
  if (is_finite (z)) Fallback::shiftl (d, z, shift);
  else d= z;
}

TMPL static inline void
shiftr (Ball& d, const xint& shift) {
  if (is_finite (d)) Fallback::shiftr (d, shift);
}

TMPL static inline void
shiftr (Ball& d, const Ball& z, const xint& shift) {
  if (is_finite (z)) Fallback::shiftr (d, z, shift);
  else d= z;
}

}; // implementation<ball_shift,IV,ball_infinities<BV> >

#undef TMPL
#undef Ball
} // namespace mmx
#endif // __MMX_BALL_INFINITIES_HPP