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

Go to the documentation of this file.

/******************************************************************************
* MODULE     : rounded.hpp
* DESCRIPTION: Rounding management
* COPYRIGHT  : (C) 2008  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__ROUNDED_HPP
#define __MMX__ROUNDED_HPP
#include <basix/operators.hpp>
#include <basix/double.hpp>
#include <numerix/floating.hpp>
namespace mmx {
#define TMPL template<typename C>

/******************************************************************************
* Default upward rounding
******************************************************************************/

struct rounded_none_up {
  template<typename Op, typename C>
  static inline C
  rounded_as () {
    typedef Return_op(Op,C) Rop;
    C r  = Rop::op ();
    C err= rounding_error (r);
    return r + (err + err);
  }

  template<typename Op, typename C, typename X1>
  static inline C
  rounded_as (const X1& x1) {
    typedef Return_op(Op,C) Rop;
    C r  = Rop::op (x1);
    C err= rounding_error (r);
    return r + (err + err);
  }

  template<typename Op, typename C, typename X1, typename X2>
  static inline C
  rounded_as (const X1& x1, const X2& x2) {
    typedef Return_op(Op,C) Rop;
    C r  = Rop::op (x1, x2);
    C err= rounding_error (r);
    return r + (err + err);
  }
};

/******************************************************************************
* Default downward rounding
******************************************************************************/

struct rounded_none_down {
  template<typename Op, typename C>
  static inline C
  rounded_as () {
    typedef Return_op(Op,C) Rop;
    C r  = Rop::op ();
    C err= rounding_error (r);
    return r - (err + err);
  }

  template<typename Op, typename C, typename X1>
  static inline C
  rounded_as (const X1& x1) {
    typedef Return_op(Op,C) Rop;
    C r  = Rop::op (x1);
    C err= rounding_error (r);
    return r - (err + err);
  }

  template<typename Op, typename C, typename X1, typename X2>
  static inline C
  rounded_as (const X1& x1, const X2& x2) {
    typedef Return_op(Op,C) Rop;
    C r  = Rop::op (x1, x2);
    C err= rounding_error (r);
    return r - (err + err);
  }
};

/******************************************************************************
* Global modification of rounding mode
******************************************************************************/

template<mmx_rounding Mode>
struct rounded_global {
  template<typename Op, typename C>
  static inline C
  rounded_as () {
    typedef Return_op(Op,C) Rop;
    typedef rounding_helper<C> H;
    H::set_rounding (H::template translate<Mode>::val);
    return Rop::op ();
  }

  template<typename Op, typename C, typename X1>
  static inline C
  rounded_as (const X1& x1) {
    typedef Return_op(Op,C) Rop;
    typedef rounding_helper<C> H;
    H::set_rounding (H::template translate<Mode>::val);
    return Rop::op (x1);
  }

  template<typename Op, typename C, typename X1, typename X2>
  static inline C
  rounded_as (const X1& x1, const X2& x2) {
    typedef Return_op(Op,C) Rop;
    typedef rounding_helper<C> H;
    H::set_rounding (H::template translate<Mode>::val);
    return Rop::op (x1, x2);
  }
};

/******************************************************************************
* Local modification of rounding mode
******************************************************************************/

template<mmx_rounding Mode>
struct rounded_local {
  template<typename Op, typename C>
  static inline C
  rounded_as () {
    typedef Return_op(Op,C) Rop;
    typedef rounding_helper<C> H;
    typename H::rounding_mode old= H::get_rounding ();
    H::set_rounding (H::template translate<Mode>::val);
    C r= Rop::op ();
    H::set_rounding (old);
    return r;
  }

  template<typename Op, typename C, typename X1>
  static inline C
  rounded_as (const X1& x1) {
    typedef Return_op(Op,C) Rop;
    typedef rounding_helper<C> H;
    typename H::rounding_mode old= H::get_rounding ();
    H::set_rounding (H::template translate<Mode>::val);
    C r= Rop::op (x1);
    H::set_rounding (old);
    return r;
  }

  template<typename Op, typename C, typename X1, typename X2>
  static inline C
  rounded_as (const X1& x1, const X2& x2) {
    typedef Return_op(Op,C) Rop;
    typedef rounding_helper<C> H;
    typename H::rounding old= H::get_rounding ();
    H::set_rounding (H::template translate<Mode>::val);
    C r= Rop::op (x1, x2);
    H::set_rounding (old);
    return r;
  }
};

/******************************************************************************
* Rounded operations on floating types
******************************************************************************/

template<typename W, typename V> inline floating<W>
view (const floating<V>& x) {
  return *((floating<W>*) ((void*) &x));
}

template<mmx_rounding Mode>
struct rounded_floating {
  template<typename Op, typename C>
  struct helper {};

  template<typename Op, typename V>
  struct helper<Op,floating<V> > {
    typedef rounding_helper<floating<V> > H;
    typedef rnd_floating<V,H::template translate<Mode>::val> RV;
    typedef Return_op (Op, floating<RV> ) Rop;

    static inline floating<V>
    op () {
      return view<V> (Rop::op ()); }

    template<typename X1>
    static inline floating<V>
    op (const X1& x1) {
      return view<V> (Rop::op (x1)); }

    template<typename X1, typename X2>
    static inline floating<V>
    op (const X1& x1, const X2& x2) {
      return view<V> (Rop::op (x1, x2)); }

    template<typename W1>
    static inline floating<V>
    op (const floating<W1>& x1) {
      return view<V> (Rop::op (view<RV> (x1))); }

    template<typename W1, typename W2>
    static inline floating<V>
    op (const floating<W1>& x1, const floating<W2>& x2) {
      return view<V> (Rop::op (view<RV> (x1), view<RV> (x2))); }
  };

  template<typename Op, typename C>
  static inline C rounded_as () {
    return helper<Op,C>::op ();
  };

  template<typename Op, typename C, typename X1>
  static inline C rounded_as (const X1& x1) {
    return helper<Op,C>::op (x1);
  };

  template<typename Op, typename C, typename X1, typename X2>
  static inline C rounded_as (const X1& x1, const X2& x2) {
    return helper<Op,C>::op (x1, x2);
  };
};

/******************************************************************************
* Standard rounded operations
******************************************************************************/

template<typename V>
struct rounded_operations: public V {

  template<typename Op, typename C>
  static inline C
  rounded () {
    return V::template rounded_as<Op,C> (); }

  template<typename Op, typename C>
  static inline C
  rounded (const C& x1) {
    return V::template rounded_as<Op,C> (x1); }

  template<typename Op, typename C>
  static inline C
  rounded (const C& x1, const C& x2) {
    return V::template rounded_as<Op,C> (x1, x2); }

  TMPL static inline C log2 () {
    return rounded<log2_as_op,C> (); }
  TMPL static inline C pi () {
    return rounded<pi_as_op,C> (); }
  TMPL static inline C euler () {
    return rounded<euler_as_op,C> (); }
  TMPL static inline C catalan () {
    return rounded<catalan_as_op,C> (); }
  TMPL static inline C smallest () {
    return rounded<smallest_as_op,C> (); }
  TMPL static inline C largest () {
    return rounded<largest_as_op,C> (); }
  TMPL static inline C accuracy () {
    return rounded<accuracy_as_op,C> (); }

  TMPL static inline C add (const C& x1, const C& x2) {
    return rounded<add_op> (x1, x2); }
  TMPL static inline C sub (const C& x1, const C& x2) {
    return rounded<sub_op> (x1, x2); }
  TMPL static inline C mul (const C& x1, const C& x2) {
    return rounded<mul_op> (x1, x2); }
  TMPL static inline C div (const C& x1, const C& x2) {
    return rounded<div_op> (x1, x2); }
  TMPL static inline C square (const C& x) {
    return rounded<square_op> (x); }

  TMPL static inline C sqrt (const C& x1) {
    return rounded<sqrt_op> (x1); }
  TMPL static inline C cbrt (const C& x1) {
    return rounded<cbrt_op> (x1); }
  TMPL static inline C hypot (const C& x1, const C& x2) {
    return rounded<hypot_op> (x1, x2); }
  TMPL static inline C exp (const C& x1) {
    return rounded<exp_op> (x1); }
  TMPL static inline C exp2 (const C& x1) {
    return rounded<exp2_op> (x1); }
  TMPL static inline C log (const C& x1) {
    return rounded<log_op> (x1); }
  TMPL static inline C log2 (const C& x1) {
    return rounded<log2_op> (x1); }
  TMPL static inline C log10 (const C& x1) {
    return rounded<log10_op> (x1); }
  TMPL static inline C pow (const C& x1, const C& x2) {
    return rounded<pow_op> (x1, x2); }
  TMPL static inline C cos (const C& x1) {
    return rounded<cos_op> (x1); }
  TMPL static inline C sin (const C& x1) {
    return rounded<sin_op> (x1); }
  TMPL static inline C tan (const C& x1) {
    return rounded<tan_op> (x1); }
  TMPL static inline C acos (const C& x1) {
    return rounded<acos_op> (x1); }
  TMPL static inline C asin (const C& x1) {
    return rounded<asin_op> (x1); }
  TMPL static inline C atan (const C& x1) {
    return rounded<atan_op> (x1); }
  TMPL static inline C atan2 (const C& x1, const C& x2) {
    return rounded<atan2_op> (x1, x2); }
  TMPL static inline C cosh (const C& x1) {
    return rounded<cosh_op> (x1); }
  TMPL static inline C sinh (const C& x1) {
    return rounded<sinh_op> (x1); }
  TMPL static inline C tanh (const C& x1) {
    return rounded<tanh_op> (x1); }
  TMPL static inline C acosh (const C& x1) {
    return rounded<acosh_op> (x1); }
  TMPL static inline C asinh (const C& x1) {
    return rounded<asinh_op> (x1); }
  TMPL static inline C atanh (const C& x1) {
    return rounded<atanh_op> (x1); }

  TMPL static inline C gamma (const C& x1) {
    return rounded<gamma> (x1); }
  TMPL static inline C zeta (const C& x1) {
    return rounded<zeta> (x1); }
  TMPL static inline C erf (const C& x1) {
    return rounded<erf> (x1); }

  TMPL static inline C floor (const C& x1) {
    return rounded<floor_op> (x1); }
  TMPL static inline C trunc (const C& x1) {
    return rounded<trunc_op> (x1); }
  TMPL static inline C ceil (const C& x1) {
    return rounded<ceil_op> (x1); }
  TMPL static inline C round (const C& x1) {
    return rounded<round_op> (x1); }

  template<typename C, typename X1>
  static inline C
  as (const X1& x1) {
    return V::template rounded_as<as_op,C> (x1); }

  template<typename C, typename X1>
  static inline C
  abs_as (const X1& x1) {
    return V::template rounded_as<abs_as_op,C> (x1); }

  template<typename C, typename X1, typename X2>
  static inline C
  add_as (const X1& x1, const X2& x2) {
    return V::template rounded_as<add_op,C> (x1, x2); }

  template<typename C, typename X1, typename X2>
  static inline C
  sub_as (const X1& x1, const X2& x2) {
    return V::template rounded_as<sub_op,C> (x1, x2); }

};

/******************************************************************************
* Standard rounded arithmetic for an opposite rounding mode
******************************************************************************/

template<typename V>
struct rounded_opposite: public V {
  TMPL static inline C add (const C& x1, const C& x2) {
    C tmp= -x1; tmp -= x2; return -tmp; }
  TMPL static inline C sub (const C& x1, const C& x2) {
    C tmp=  x2; tmp -= x1; return -tmp; }
  TMPL static inline C mul (const C& x1, const C& x2) {
    C tmp= -x1; tmp *= x2; return -tmp; }
  TMPL static inline C div (const C& x1, const C& x2) {
    C tmp= -x1; tmp /= x2; return -tmp; }
  TMPL static inline C square (const C& x1) {
    C tmp= -x1; tmp *= x1; return -tmp; }
};

/******************************************************************************
* Type dependent rounding modes
******************************************************************************/

STMPL
struct rounding_helper<double> {
  typedef int rounding_mode;
  template<mmx_rounding Mode, typename K= nat> struct translate {
    static const rounding_mode val= FE_TONEAREST; };
  template<typename K> struct translate<MMX_ROUND_NEAR,K> {
    static const rounding_mode val= FE_TONEAREST; };
  template<typename K> struct translate<MMX_ROUND_ZERO,K> {
    static const rounding_mode val= FE_TOWARDZERO; };
  template<typename K> struct translate<MMX_ROUND_UP,K> {
    static const rounding_mode val= FE_UPWARD; };
  template<typename K> struct translate<MMX_ROUND_DOWN,K> {
    static const rounding_mode val= FE_DOWNWARD; };
  static inline rounding_mode get_rounding () {
    return fegetround (); }
  static inline void set_rounding (rounding_mode r) {
    fesetround (r); }
  static const bool exact_type = false;
  typedef rounded_operations<rounded_global<MMX_ROUND_UP> > UV;
  typedef rounded_operations<rounded_global<MMX_ROUND_DOWN> > DV;
};

template<typename V>
struct rounding_helper<floating<V> > {
  typedef mpfr_rnd_t rounding_mode;
  template<mmx_rounding Mode, typename K= nat> struct translate {
    static const rounding_mode val= GMP_RNDN; };
  template<typename K> struct translate<MMX_ROUND_NEAR,K> {
    static const rounding_mode val= GMP_RNDN; };
  template<typename K> struct translate<MMX_ROUND_ZERO,K> {
    static const rounding_mode val= GMP_RNDZ; };
  template<typename K> struct translate<MMX_ROUND_UP,K> {
    static const rounding_mode val= GMP_RNDU; };
  template<typename K> struct translate<MMX_ROUND_DOWN,K> {
    static const rounding_mode val= GMP_RNDD; };
  static inline rounding_mode get_rounding () {
    return mmx_rounding_mode; }
  static inline void set_rounding (rounding_mode r) {
    mmx_rounding_mode= r; }
  static const bool exact_type = false;
  typedef rounded_operations<rounded_floating<MMX_ROUND_UP> > UV;
  typedef rounded_operations<rounded_floating<MMX_ROUND_DOWN> > DV;
};

#undef TMPL
}
#endif // __MMX__ROUNDED_HPP