00001
00002
00003
00004
00005
00006 #ifndef mmx_numerics_rounding_mode_hpp
00007 #define mmx_numerics_rounding_mode_hpp
00008
00009 #include <fenv.h>
00010
00011 namespace mmx {
00012
00013 namespace numerics
00014 {
00015 inline int get_cw() { return fegetround(); };
00016 inline int get_rnd() { return fegetround(); };
00017 inline void set_cw( int cw ) { fesetround(cw);};
00018 inline static int rnd_up() { return FE_UPWARD;};
00019 inline static int rnd_dw() { return FE_DOWNWARD;};
00020 inline static int rnd_nr() { return FE_TONEAREST;};
00021 inline static int rnd_z() { return FE_TOWARDZERO;};
00022 template< typename T > struct LongVersion { typedef T result_t; };
00023 template<> struct LongVersion<float> { typedef long double result_t; };
00024 template<> struct LongVersion<double> { typedef long double result_t; };
00025 template< class T > struct fpu_rounding {
00026 static const int id = 0;
00027 typedef int rnd_t;
00028 inline static rnd_t getrnd() { return 0; };
00029 inline static void setrnd( rnd_t r ) { };
00030 inline static int rnd_up() { return 0; };
00031 inline static int rnd_dw() { return 0; };
00032 inline static int rnd_nr() { return 0; };
00033 inline static int rnd_z() { return 0; };
00034 fpu_rounding( int ){};
00035 };
00036 template<> struct fpu_rounding<long double>
00037 {
00038 typedef int rnd_t;
00039 inline static rnd_t getrnd() { return fegetround(); };
00040 inline static void setrnd( rnd_t r ) { fesetround(r); };
00041 inline static int rnd_up() { return FE_UPWARD;};
00042 inline static int rnd_dw() { return FE_DOWNWARD;};
00043 inline static int rnd_nr() { return FE_TONEAREST;};
00044 inline static int rnd_z() { return FE_TOWARDZERO;};
00045 static const int id = 1;
00046 rnd_t m_prev;
00047 bool m_chg;
00048 fpu_rounding( rnd_t rnd )
00049 {
00050 m_chg = false;
00051 m_prev = getrnd();
00052 if ( m_prev != rnd )
00053 {
00054 m_chg = true;
00055 setrnd(rnd);
00056 };
00057 };
00058 ~fpu_rounding() {
00059 if ( m_chg ) setrnd(m_prev);
00060 };
00061 };
00062
00063
00064 template< typename T >
00065 struct rounding : fpu_rounding< typename LongVersion< T >::result_t >
00066 { rounding( typename fpu_rounding<typename LongVersion<T>::result_t >::rnd_t rnd ) : fpu_rounding< typename LongVersion< T >::result_t >( rnd ) {}; };
00067
00068 template< typename T > inline T div_dw(const T& a, const T& b) { rounding<T> r(rounding<T>::rnd_dw()); return a/b; };
00069 template< typename T > inline T div_up(const T& a, const T& b) { rounding<T> r(rounding<T>::rnd_up()); return a/b; };
00070 template< typename T > inline T add_dw(const T& a, const T& b) { rounding<T> r(rounding<T>::rnd_dw()); return a+b; };
00071 template< typename T > inline T add_up(const T& a, const T& b) { rounding<T> r(rounding<T>::rnd_up()); return a+b; };
00072 template< typename T > inline T mul_dw(const T& a, const T& b) { rounding<T> r(rounding<T>::rnd_dw()); return a*b; };
00073 template< typename T > inline T mul_up(const T& a, const T& b) { rounding<T> r(rounding<T>::rnd_up()); return a*b; };
00074 template< typename T > inline T sub_dw(const T& a, const T& b) { rounding<T> r(rounding<T>::rnd_dw()); return a-b; };
00075 template< typename T > inline T sub_up(const T& a, const T& b) { rounding<T> r(rounding<T>::rnd_up()); return a-b; };
00076
00077 namespace rdw
00078 {
00079 template< typename T > inline
00080 T upadd( const T& a, const T& b ) { T tmp = -a; tmp -= b; return -tmp; };
00081 template< typename T > inline
00082 T dwadd( const T& a, const T& b ) { return a+b; };
00083 template< typename T > inline
00084 T upsub( const T& a, const T& b ) { T tmp = b; tmp -= a; return -tmp; };
00085 template< typename T > inline
00086 T dwsub( const T& a, const T& b ) { return a-b; };
00087 template< typename T > inline
00088 T upmul( const T& a, const T& b ) { T tmp = -a; tmp *= b; return -tmp; };
00089 template< typename T > inline
00090 T dwmul( const T& a, const T& b ) { return a*b; };
00091 template< typename T > inline
00092 T updiv( const T& a, const T& b ) { T tmp = -a; tmp /= b; return -tmp; };
00093 template< typename T > inline
00094 T dwdiv( const T& a, const T& b ) { return a/b; };
00095 template< typename T > inline
00096 T updet( const T& a, const T& b, const T& c, const T& d )
00097 { return upsub(upmul(a,d),dwmul(b,c)); };
00098 template< typename T > inline
00099 T dwdet( const T& a, const T& b, const T& c, const T& d )
00100 { return dwsub(dwmul(a,d),upmul(b,c)); };
00101 };
00102
00103 namespace rup
00104 {
00105 template< typename T > inline
00106 T upadd( const T& a, const T& b ) { return a+b; };
00107 template< typename T > inline
00108 T dwadd( const T& a, const T& b ) { T tmp = -a; tmp -= b; return -tmp; };
00109 template< typename T > inline
00110 T upsub( const T& a, const T& b ) { return a-b; };
00111 template< typename T > inline
00112 T dwsub( const T& a, const T& b ) { T tmp = b; tmp -= a; return -tmp; };
00113 template< typename T > inline
00114 T upmul( const T& a, const T& b ) { return a*b; };
00115 template< typename T > inline
00116 T dwmul( const T& a, const T& b ) { T tmp = -a; tmp *= b; return -tmp; };
00117 template< typename T > inline
00118 T updiv( const T& a, const T& b ) { return a/b; };
00119 template< typename T > inline
00120 T dwdiv( const T& a, const T& b ) { T tmp = -a; tmp /= b; return -tmp; };
00121 template< typename T > inline
00122 T updet( const T& a, const T& b, const T& c, const T& d )
00123 { return upsub(upmul(a,d),dwmul(b,c)); };
00124 template< typename T > inline
00125 T dwdet( const T& a, const T& b, const T& c, const T& d )
00126 { return dwsub(dwmul(a,d),upmul(b,c)); };
00127 }
00128 };
00129
00130 }
00131
00132 #endif //
