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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : ball_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_ROUNDED_HPP
#define __MMX_BALL_ROUNDED_HPP
#include <numerix/rounded.hpp>
namespace mmx {
#define TMPL template<typename C,typename R,typename V>
#define Ball ball<C,R,V>

TMPL class ball;
TMPL inline C center (const Ball& z);
TMPL inline R radius (const Ball& z);
TMPL inline C& center (Ball& z);
TMPL inline R& radius (Ball& z);

/******************************************************************************
* The ball_rounded variant
******************************************************************************/

struct ball_rounded {};

template<typename C>
struct ball_variant_helper {
  typedef ball_rounded BV;
};

/******************************************************************************
* Rounding and ball limit related functions
******************************************************************************/

struct ball_rounding {};

template<typename W>
struct implementation<ball_rounding,W,ball_rounded> {

TMPL static inline void
add_rounding_error (Ball& d) {
  typedef Round_up(R) Up;
  R err= Up::template as<R> (rounding_error (center (d)));
  radius (d)= Up::add (radius (d), err);
}

TMPL static inline void
add_additive_error (Ball& d) {
  typedef Round_up(R) Up;
  R err= Up::template as<R> (additive_error (center (d)));
  radius (d)= Up::add (radius (d), err);
}

TMPL static inline void
add_multiplicative_error (Ball& d) {
  typedef Round_up(R) Up;
  R err= Up::template as<R> (multiplicative_error (center (d)));
  radius (d)= Up::add (radius (d), err);
}

TMPL static inline void
add_elementary_error (Ball& d) {
  typedef Round_up(R) Up;
  R err= Up::template as<R> (elementary_error (center (d)));
  radius (d)= Up::add (radius (d), err);
}

TMPL static inline C
lower (const Ball& z) {
  typedef Round_down(C) Down;
  return Down::template sub_as<C> (center (z), radius (z));
}

TMPL static inline C
upper (const Ball& z) {
  typedef Round_up(C) Up;
  return Up::template add_as<C> (center (z), radius (z));
}

TMPL static inline void
make_interval (Ball& d, const C& l, const C& r) {
  typedef Round_up(C) Up;
  C c= decexp2 (l + r, 1);
  center (d)= c;
  radius (d)= Up::template as<R> (max (Up::sub (r, c), Up::sub (c, l)));
}

TMPL static inline R
abs_down (const Ball& z) {
  typedef Round_down(R) Down;
  R a= Down::template abs_as<R> (center (z));
  if (a <= radius (z)) return R(0);
  else return Down::sub (a, radius (z));
}

TMPL static inline R
abs_up (const Ball& z) {
  typedef Round_up(R) Up;
  return Up::add (Up::template abs_as<R> (center(z)), radius(z));
}

TMPL static inline R
bnd_down (const Ball& z) {
  typedef Round_down(R) Down;
  return Down::sub (Down::template as<R> (center(z)), radius(z));
}

TMPL static inline R
bnd_up (const Ball& z) {
  typedef Round_up(R) Up;
  return Up::add (Up::template as<R> (center(z)), radius(z));
}

}; // implementation<ball_rounding,W,ball_rounded>

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

struct ball_additive {};

template<typename W>
struct implementation<ball_additive,W,ball_rounded> {
  typedef implementation<ball_rounding,W> Rnd;

TMPL static inline void
neg (Ball& d, const Ball& z) {
  center (d)= -center (z);
  radius (d)= radius (z);
}

TMPL static inline void
add (Ball& d, const Ball& z1, const Ball& z2) {
  typedef Round_up(R) Up;
  center (d)= center (z1) + center (z2);
  radius (d)= Up::add (radius (z1), radius (z2));
  Rnd::add_additive_error (d);
}

TMPL static inline void
add (Ball& d, const C& z1, const Ball& z2) {
  center (d)= z1 + center (z2);
  radius (d)= radius (z2);
  Rnd::add_additive_error (d);
}

TMPL static inline void
add (Ball& d, const Ball& z1, const C& z2) {
  center (d)= center (z1) + z2;
  radius (d)= radius (z1);
  Rnd::add_additive_error (d);
}

TMPL static inline void
sub (Ball& d, const Ball& z1, const Ball& z2) {
  typedef Round_up(R) Up;
  center (d)= center (z1) - center (z2);
  radius (d)= Up::add (radius (z1), radius (z2));
  Rnd::add_additive_error (d);
}

TMPL static inline void
sub (Ball& d, const C& z1, const Ball& z2) {
  center (d)= z1 - center (z2);
  radius (d)= radius (z2);
  Rnd::add_additive_error (d);
}

TMPL static inline void
sub (Ball& d, const Ball& z1, const C& z2) {
  center (d)= center (z1) - z2;
  radius (d)= radius (z1);
  Rnd::add_additive_error (d);
}

}; // implementation<ball_additive,W,ball_rounded>

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

struct ball_multiplicative {};

template<typename W>
struct implementation<ball_multiplicative,W,ball_rounded> {
  typedef implementation<ball_rounding,W> Rnd;

TMPL static inline void
mul (Ball& d, const Ball& z1, const Ball& z2) {
  typedef Round_up(R) Up;
  R a1= Rnd::abs_up (z1);
  R a2= Up::template abs_as<R> (center (z2));
  center (d)= center (z1) * center (z2);
  radius (d)= Up::add (Up::mul (a1, radius (z2)), Up::mul (radius (z1), a2));
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
mul (Ball& d, const C& z1, const Ball& z2) {
  typedef Round_up(R) Up;
  R a1= Up::template abs_as<R> (z1);
  center (d)= z1 * center (z2);
  radius (d)= Up::mul (a1, radius (z2));
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
mul (Ball& d, const Ball& z1, const C& z2) {
  typedef Round_up(R) Up;
  R a2= Up::template abs_as<R> (z2);
  center (d)= center (z1) * z2;
  radius (d)= Up::mul (radius (z1), a2);
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
square (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  R a= Up::template abs_as<R> (center (z));
  center (d)= square_op::op (center (z));
  radius (d)= Up::mul (Up::add (a + a, radius (z)), radius (z));
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
invert (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  typedef Round_down(R) Down;
  center (d)= invert_op::op (center (z));
  radius (d)= Up::div (radius (z), Down::square (Rnd::abs_down (z)));
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
invert (Ball& d, const C& z) {
  center (d)= invert_op::op (z);
  radius (d)= multiplicative_error (center (d));
}

TMPL static inline void
div (Ball& d, const Ball& z1, const Ball& z2) {
  d= z1 * invert_op::op (z2);
}

TMPL static inline void
div (Ball& d, const Ball& z1, const C& z2) {
  d= z1 * invert_op::op (z2);
}

TMPL static inline void
div (Ball& d, const C& z1, const Ball& z2) {
  d= z1 * invert_op::op (z2);
}

}; // implementation<ball_multiplicative,W,ball_rounded>

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

struct ball_root {};

template<typename W>
struct implementation<ball_root,W,ball_rounded> {
  typedef implementation<ball_rounding,W> Rnd;

TMPL static void
sqrt (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  typedef Round_down(R) Down;
  if (is_negative_or_zero (z)) d= Nan (Ball);
  else {
    R aux= Down::sqrt (Rnd::abs_down (z));
    center (d)= sqrt_op::op (center (z));
    radius (d)= Up::div (radius (z), aux + aux);
    Rnd::add_multiplicative_error (d);
  }
}

TMPL static void
hypot (Ball& d, const Ball& x, const Ball& y) {
  typedef Round_up(R) Up;
  center (d)= hypot_op::op (center (x), center (y));
  radius (d)= Up::hypot (radius (x), radius (y));
  Rnd::add_multiplicative_error (d);
}

}; // implementation<ball_root,W,ball_rounded>

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

struct ball_elementary {};

template<typename W>
struct implementation<ball_elementary,W,ball_rounded> {
  typedef implementation<ball_rounding,W> 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 (Rnd::bnd_up (z)), radius (z));
  Rnd::add_elementary_error (d);
}

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

TMPL static inline void
pow (Ball& d, const Ball& z, const Ball& y) {
  d= exp_op::op (y * log_op::op (z));
}

TMPL static void
cos (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  typedef Round_down(R) Down;
  if (radius (z) >= Down::template pi<R> ()) {
    center (d)= promote (0, center (z));
    radius (d)= promote (1, radius (z));
  }
  else {
    center (d)= cos_op::op (center (z));
    radius (d)= radius (z);
    Rnd::add_elementary_error (d);
    if (radius (z) < abs (center (d))) {
      R u1= Up::sin (Rnd::bnd_down (z));
      R u2= Up::sin (Rnd::bnd_up (z));
      radius (d)= Up::mul (max (abs (u1), abs (u2)), radius (z));
      Rnd::add_elementary_error (d);
    }
  }
}

TMPL static void
sin (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  typedef Round_down(R) Down;
  if (radius (z) >= Down::template pi<R> ()) {
    center (d)= promote (0, center (z));
    radius (d)= promote (1, radius (z));
  }
  else {
    center (d)= sin_op::op (center (z));
    radius (d)= radius (z);
    Rnd::add_elementary_error (d);
    if (radius (z) < abs (center (d))) {
      R u1= Up::cos (Rnd::bnd_down (z));
      R u2= Up::cos (Rnd::bnd_up (z));
      radius (d)= Up::mul (max (abs (u1), abs (u2)), radius (z));
      Rnd::add_elementary_error (d);
    }
  }
}

TMPL static inline 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 (Rnd::bnd_up (abs (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 (Rnd::bnd_up (abs (z))), radius (z));
  Rnd::add_elementary_error (d);
}

TMPL static inline 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) {
  typedef Round_up(R) Up;
  typedef Round_down(R) Down;
  Ball u= promote (1, z) - square_op::op (z);
  R v= Rnd::bnd_down (u);
  if (v <= 0) d= Nan (Ball);
  else {
    center (d)= acos_op::op (center (z));
    radius (d)= Up::div (radius (z), Down::sqrt (v));
    Rnd::add_elementary_error (d);
  }
}

TMPL static void
asin (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  typedef Round_down(R) Down;
  Ball u= promote (1, z) - square_op::op (z);
  R v= Rnd::bnd_down (u);
  if (v <= 0) d= Nan (Ball);
  else {
    center (d)= asin_op::op (center (z));
    radius (d)= Up::div (radius (z), Down::sqrt (v));
    Rnd::add_elementary_error (d);
  }
}

TMPL static void
atan (Ball& d, const Ball& z) {
  typedef Round_up(R) Up;
  typedef Round_down(R) Down;
  center (d)= atan_op::op (center (z));
  R u= Rnd::bnd_down (z);
  R v= Down::add (promote (1, radius (z)), Down::square (u));
  if (u <= 0) radius (d)= radius (z);
  else radius (d)= Up::div (radius (z), v);
  Rnd::add_elementary_error (d);
}

TMPL static void
atan2 (Ball& d, const Ball& y, const Ball& x) {
  typedef Round_up(R) Up;
  Ball z= square_op::op (x) + square_op::op (y);
  if (z <= 0) d= Nan (Ball);
  else {
    R b = Up::div (promote (1, radius (y)), Rnd::bnd_down (z));
    R dx= Up::mul (Rnd::bnd_up (abs (y)), b);
    R dy= Up::mul (Rnd::bnd_up (abs (x)), b);
    center (d)= atan2_op::op (center (y), center (x));
    radius (d)= Up::add (Up::mul (dx, radius (x)), Up::mul (dy, radius (y)));
    Rnd::add_elementary_error (d);
  }
}

}; // implementation<ball_elementary,W,ball_rounded>

/******************************************************************************
* Operations on ordered center types
******************************************************************************/

struct ball_ordered {};

template<typename W>
struct implementation<ball_ordered,W,ball_rounded> {
  typedef implementation<ball_rounding,W> Rnd;

TMPL static void
min (Ball& d, const Ball& z1, const Ball& z2) {
  if (z1 < z2) d= z1;
  else if (z2 < z1) d= z2;
  else {
    center (d)= min_op::op (center (z1), center (z2));
    radius (d)= max_op::op (radius (z1), radius (z2));
  }
}

TMPL static void
max (Ball& d, const Ball& z1, const Ball& z2) {
  if (z1 > z2) d= z1;
  else if (z2 > z1) d= z2;
  else {
    center (d)= max_op::op (center (z1), center (z2));
    radius (d)= max_op::op (radius (z1), radius (z2));
  }
}

TMPL static void
floor (Ball& d, const Ball& z) {
  typedef Round_up(C) Up;
  typedef Round_down(C) Down;
  C l= Down::floor (lower (z));
  C r= Up  ::floor (upper (z));
  Rnd::make_interval (d, l, r);
}

TMPL static void
trunc (Ball& d, const Ball& z) {
  typedef Round_up(C) Up;
  typedef Round_down(C) Down;
  C l= Down::trunc (lower (z));
  C r= Up  ::trunc (upper (z));
  Rnd::make_interval (d, l, r);
}

TMPL static void
ceil (Ball& d, const Ball& z) {
  typedef Round_up(C) Up;
  typedef Round_down(C) Down;
  C l= Down::ceil (lower (z));
  C r= Up  ::ceil (upper (z));
  Rnd::make_interval (d, l, r);
}

TMPL static void
round (Ball& d, const Ball& z) {
  typedef Round_up(C) Up;
  typedef Round_down(C) Down;
  C l= Down::round (lower (z));
  C r= Up  ::round (upper (z));
  Rnd::make_interval (d, l, r);
}

}; // implementation<ball_ordered,W,ball_rounded>

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

struct ball_abs {};

template<typename W>
struct implementation<ball_abs,W,ball_rounded> {
  typedef implementation<ball_rounding,W> Rnd;

template<typename C ,typename R ,typename V,
         typename C2,typename V2> static inline void
abs (ball<C2,R,V2>& d, const Ball& z) {
  center (d)= abs_op::op (center (z));
  radius (d)= radius (z);
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
abs (Ball& d, const Ball& z) {
  center (d)= abs_op::op (center (z));
  radius (d)= radius (z);
}

}; // implementation<ball_abs,W,ball_rounded>

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

struct ball_shift {};

template<typename W>
struct implementation<ball_shift,W,ball_rounded> {
  typedef implementation<ball_rounding,W> Rnd;

TMPL static inline void
shiftl (Ball& d, const xint& shift) {
  center (d)= incexp2 (center (d), shift);
  radius (d)= incexp2 (radius (d), shift);
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
shiftl (Ball& d, const Ball& z, const xint& shift) {
  center (d)= incexp2 (center (z), shift);
  radius (d)= incexp2 (radius (z), shift);
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
shiftr (Ball& d, const xint& shift) {
  center (d)= decexp2 (center (d), shift);
  radius (d)= decexp2 (radius (d), shift);
  Rnd::add_multiplicative_error (d);
}

TMPL static inline void
shiftr (Ball& d, const Ball& z, const xint& shift) {
  center (d)= decexp2 (center (z), shift);
  radius (d)= decexp2 (radius (z), shift);
  Rnd::add_multiplicative_error (d);
}

}; // implementation<ball_shift,W,ball_rounded>

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