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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : ball_complex_rounded.hpp
* DESCRIPTION: Balls over a center type with exactly rounded arithmetic
* 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_COMPLEX_ROUNDED_HPP
#define __MMX_BALL_COMPLEX_ROUNDED_HPP
#include <numerix/ball.hpp>
#include <numerix/complex.hpp>
namespace mmx {
#define TMPL template<typename C,typename R,typename V>
#define CTMPL template<typename C,typename R,typename V,typename CC,typename CV>
#define Ball ball<C,R,V>
#define CBall ball<CC,R,CV>

/******************************************************************************
* The ball_complex<V> variant
******************************************************************************/

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

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

template<typename C>
struct ball_variant_helper<complex<C> > {
  typedef Ball_variant(C) BV;
};

/******************************************************************************
* Complex ball constructors
******************************************************************************/

struct ball_complex_construct {};

template<typename IV, typename BV>
struct implementation<ball_complex_construct,IV,ball_complex<BV> > {
  typedef implementation<ball_rounding,BV> Rnd;

CTMPL static void
gauss (CBall& d, const Ball& re, const Ball& im) {
  typedef Round_up(R) Up;
  gaussian_as (center (d), center (re), center (im));
  radius (d)= Up::hypot (radius (re), radius (im));
  Rnd::add_elementary_error (d);
}

CTMPL static void
polar (CBall& d, const Ball& r, const Ball& t) {
  gauss (d, r * cos (t), r * sin (t));
}

}; // implementation<ball_complex_construct,IV,ball_complex<BV> >

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

template<typename IV, typename BV>
struct implementation<ball_root,IV,ball_complex<BV> >:
  public implementation<ball_root,BV>
{
  typedef implementation<ball_rounding,BV> Rnd;

TMPL static void
sqrt (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  typedef Round_down(R) Down;
  if (is_zero (z)) d= Fuzz (Ball);
  else if (is_infinite (z)) d= Infinity (Ball);
  else {
    center (d)= sqrt_op::op (center (z));
    radius (d)= Up::div (decexp2 (radius (z)), Down::sqrt (abs_down (z)));
    Rnd::add_multiplicative_error (d);
  }
}

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

/******************************************************************************
* Elementary functions
******************************************************************************/

template<typename IV, typename BV>
struct implementation<ball_elementary,IV,ball_complex<BV> >:
  public implementation<ball_elementary,BV>
{
  typedef implementation<ball_rounding,BV> Rnd;

TMPL static void
exp (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  center (d)= exp_op::op (center (z));
  radius (d)= Up::mul (Up::exp (bnd_up (Re (z))), radius (z));
  Rnd::add_elementary_error (d);
}

TMPL static void
log (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  if (is_zero (z)) d= Nan (Ball);
  else {
    center (d)= log_op::op (center (z));
    radius (d)= Up::div (radius (z), abs_down (z));
    Rnd::add_elementary_error (d);
  }
}

TMPL static void
cos (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  center (d)= cos_op::op (center (z));
  radius (d)= Up::mul (Up::sinh (bnd_up (Im (z))), radius (z));
  Rnd::add_elementary_error (d);
}

TMPL static void
sin (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  center (d)= sin_op::op (center (z));
  radius (d)= Up::mul (Up::cosh (bnd_up (Im (z))), radius (z));
  Rnd::add_elementary_error (d);
}

TMPL static void
tan (Ball& d, const Ball& z) {
  d= sin_op::op (z) / cos_op::op (z);
}

TMPL static inline void
cosh (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  center (d)= cosh_op::op (center (z));
  radius (d)= Up::mul (Up::sinh (bnd_up (abs (Re (z)))), radius (z));
  Rnd::add_elementary_error (d);
}

TMPL static inline void
sinh (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  center (d)= sinh_op::op (center (z));
  radius (d)= Up::mul (Up::cosh (bnd_up (abs (Re (z)))), radius (z));
  Rnd::add_elementary_error (d);
}

TMPL static void
tanh (Ball& d, const Ball& z) {
  d= sinh_op::op (z) / cosh_op::op (z);
}

TMPL static void
acos (Ball& d, const Ball& z) {
  d= over_i (acosh_op::op (z));
}

TMPL static void
asin (Ball& d, const Ball& z) {
  d= over_i (asinh_op::op (times_i (z)));
}

TMPL static void
atan (Ball& d, const Ball& z) {
  d= over_i (atanh_op::op (times_i (z)));
}

}; // implementation<ball_elementary,IV,ball_complex<BV> >

#undef TMPL
#undef CTMPL
#undef Ball
#undef CBall
} // namespace mmx
#endif // __MMX_BALL_COMPLEX_ROUNDED_HPP