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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : ball_rough.hpp
* DESCRIPTION: More efficient but rough ball arithmetic
* COPYRIGHT  : (C) 2010  Joris van der Hoeven
*******************************************************************************
* This file implements a more efficient implementation of ball arithmetic,
* which works for all rounding modes and never changes the rounding mode.
* Nevertheless, it does assume IEEE 754 compliant floating point types.
* The implementation does not attempt to provide correctness for
* infinite ball argnuments. Furthermore, some operations which should
* return a fuzzy ball with an infinite radius are allowed to only return
* a ball with a largest possible finite radius Largest(C).
*******************************************************************************
* 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_ROUGH_HPP
#define __MMX_BALL_ROUGH_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_rough variant
******************************************************************************/

struct ball_rough {};

template<typename F, typename V>
struct implementation<F,V,ball_rough>:
  public implementation<F,V,ball_rounded> {};

STMPL
struct ball_variant_helper<double> {
  typedef ball_rough BV;
};

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

TMPL inline void
add_rough_additive_error (Ball& d) {
  radius (d) += additive_error (center (d));
  radius (d) *= (1 + incexp2 (Accuracy (R), 2));
}

template<typename V> inline void
add_rough_additive_error (ball<double,double,V>& d) {
  const double alpha= 1.000000000000008882;
  radius (d)= alpha * (radius (d) + ldexp (abs (center (d)), -50));
}

template<typename W>
struct implementation<ball_additive,W,ball_rough> {

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) {
  center (d)= center (z1) + center (z2);
  radius (d)= radius (z1) + radius (z2);
  add_rough_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);
  add_rough_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);
  add_rough_additive_error (d);
}

TMPL static inline void
sub (Ball& d, const Ball& z1, const Ball& z2) {
  center (d)= center (z1) - center (z2);
  radius (d)= radius (z1) + radius (z2);
  add_rough_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);
  add_rough_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);
  add_rough_additive_error (d);
}

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

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

TMPL inline void
add_rough_multiplicative_error (Ball& d) {
  radius (d) += multiplicative_error (center (d));
  radius (d) *= (1 + incexp2 (Accuracy (R), 4));
  radius (d) += incexp2 (Smallest (R), 4);
}

template<typename V> inline void
add_rough_multiplicative_error (ball<double,double,V>& d) {
  const double alpha= 1.000000000000001777;
  const double beta = 3.952525166729972354e-323;
  radius (d)= alpha * (radius (d) + ldexp (abs (center (d)), -50)) + beta;
}

template<typename C> inline C
rough_next_above (const C& z) {
  return (1 + Accuracy (C)) * z;
}

inline double
rough_next_above (const double& z) {
  const double alpha= 1.000000000000000223;
  return alpha * z;
}

template<typename W>
struct implementation<ball_multiplicative,W,ball_rough> {

TMPL static inline void
mul (Ball& d, const Ball& z1, const Ball& z2) {
  R a1= as<R> (abs (center (z1))) + radius (z1);
  R a2= as<R> (abs (center (z2)));
  center (d)= center (z1) * center (z2);
  radius (d)= a1 * radius (z2) + radius (z1) * a2;
  add_rough_multiplicative_error (d);
}

TMPL static inline void
mul (Ball& d, const C& z1, const Ball& z2) {
  R a1= as<R> (abs (z1));
  center (d)= z1 * center (z2);
  radius (d)= a1 * radius (z2);
  add_rough_multiplicative_error (d);
}

TMPL static inline void
mul (Ball& d, const Ball& z1, const C& z2) {
  R a2= as<R> (abs (z2));
  center (d)= center (z1) * z2;
  radius (d)= radius (z1) * a2;
  add_rough_multiplicative_error (d);
}

TMPL static inline void
square (Ball& d, const Ball& z) {
  R a= as<R> (abs (center (z)));
  center (d)= square_op::op (center (z));
  radius (d)= (a + a + radius (z)) * radius (z);
  add_rough_multiplicative_error (d);
}

TMPL static inline void
invert (Ball& d, const Ball& z) {
  R a= max (as<R> (abs (center (z))) - rough_next_above (radius (z)), R(0));
  center (d)= invert_op::op (center (z));
  radius (d)= radius (z) / square_op::op (a);
  add_rough_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_rough>

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

template<typename W>
struct implementation<ball_root,W,ball_rough> {

TMPL static void
sqrt (Ball& d, const Ball& z) {
  R l= as<R> (center (z)) - rough_next_above (radius (z));
  if (l <= 0) d= Nan (Ball);
  else {
    center (d)= sqrt_op::op (center (z));
    radius (d)= decexp2 (radius (z)) / sqrt_op::op (l);
    add_rough_multiplicative_error (d);
  }
}

TMPL static void
hypot (Ball& d, const Ball& x, const Ball& y) {
  center (d)= hypot_op::op (center (x), center (y));
  radius (d)= hypot_op::op (radius (x), radius (y));
  add_rough_multiplicative_error (d);
}

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

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

template<typename W>
struct implementation<ball_shift,W,ball_rough> {

TMPL static inline void
shiftl (Ball& d, const xint& shift) {
  center (d)= incexp2 (center (d), shift);
  radius (d)= incexp2 (radius (d), shift) + Smallest (R);
}

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) + Smallest (R);
}

TMPL static inline void
shiftr (Ball& d, const xint& shift) {
  center (d)= decexp2 (center (d), shift);
  radius (d)= decexp2 (radius (d), shift) + Smallest (R);
}

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) + Smallest (R);
}

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

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