00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX__ROUNDED_HPP
00014 #define __MMX__ROUNDED_HPP
00015 #include <basix/operators.hpp>
00016 #include <basix/double.hpp>
00017 #include <numerix/floating.hpp>
00018 namespace mmx {
00019 #define TMPL template<typename C>
00020
00021
00022
00023
00024
00025 struct rounded_none_up {
00026 template<typename Op, typename C>
00027 static inline C
00028 rounded_as () {
00029 typedef Return_op(Op,C) Rop;
00030 C r = Rop::op ();
00031 C err= rounding_error (r);
00032 return r + (err + err);
00033 }
00034
00035 template<typename Op, typename C, typename X1>
00036 static inline C
00037 rounded_as (const X1& x1) {
00038 typedef Return_op(Op,C) Rop;
00039 C r = Rop::op (x1);
00040 C err= rounding_error (r);
00041 return r + (err + err);
00042 }
00043
00044 template<typename Op, typename C, typename X1, typename X2>
00045 static inline C
00046 rounded_as (const X1& x1, const X2& x2) {
00047 typedef Return_op(Op,C) Rop;
00048 C r = Rop::op (x1, x2);
00049 C err= rounding_error (r);
00050 return r + (err + err);
00051 }
00052 };
00053
00054
00055
00056
00057
00058 struct rounded_none_down {
00059 template<typename Op, typename C>
00060 static inline C
00061 rounded_as () {
00062 typedef Return_op(Op,C) Rop;
00063 C r = Rop::op ();
00064 C err= rounding_error (r);
00065 return r - (err + err);
00066 }
00067
00068 template<typename Op, typename C, typename X1>
00069 static inline C
00070 rounded_as (const X1& x1) {
00071 typedef Return_op(Op,C) Rop;
00072 C r = Rop::op (x1);
00073 C err= rounding_error (r);
00074 return r - (err + err);
00075 }
00076
00077 template<typename Op, typename C, typename X1, typename X2>
00078 static inline C
00079 rounded_as (const X1& x1, const X2& x2) {
00080 typedef Return_op(Op,C) Rop;
00081 C r = Rop::op (x1, x2);
00082 C err= rounding_error (r);
00083 return r - (err + err);
00084 }
00085 };
00086
00087
00088
00089
00090
00091 template<mmx_rounding Mode>
00092 struct rounded_global {
00093 template<typename Op, typename C>
00094 static inline C
00095 rounded_as () {
00096 typedef Return_op(Op,C) Rop;
00097 typedef rounding_helper<C> H;
00098 H::set_rounding (H::template translate<Mode>::val);
00099 return Rop::op ();
00100 }
00101
00102 template<typename Op, typename C, typename X1>
00103 static inline C
00104 rounded_as (const X1& x1) {
00105 typedef Return_op(Op,C) Rop;
00106 typedef rounding_helper<C> H;
00107 H::set_rounding (H::template translate<Mode>::val);
00108 return Rop::op (x1);
00109 }
00110
00111 template<typename Op, typename C, typename X1, typename X2>
00112 static inline C
00113 rounded_as (const X1& x1, const X2& x2) {
00114 typedef Return_op(Op,C) Rop;
00115 typedef rounding_helper<C> H;
00116 H::set_rounding (H::template translate<Mode>::val);
00117 return Rop::op (x1, x2);
00118 }
00119 };
00120
00121
00122
00123
00124
00125 template<mmx_rounding Mode>
00126 struct rounded_local {
00127 template<typename Op, typename C>
00128 static inline C
00129 rounded_as () {
00130 typedef Return_op(Op,C) Rop;
00131 typedef rounding_helper<C> H;
00132 typename H::rounding_mode old= H::get_rounding ();
00133 H::set_rounding (H::template translate<Mode>::val);
00134 C r= Rop::op ();
00135 H::set_rounding (old);
00136 return r;
00137 }
00138
00139 template<typename Op, typename C, typename X1>
00140 static inline C
00141 rounded_as (const X1& x1) {
00142 typedef Return_op(Op,C) Rop;
00143 typedef rounding_helper<C> H;
00144 typename H::rounding_mode old= H::get_rounding ();
00145 H::set_rounding (H::template translate<Mode>::val);
00146 C r= Rop::op (x1);
00147 H::set_rounding (old);
00148 return r;
00149 }
00150
00151 template<typename Op, typename C, typename X1, typename X2>
00152 static inline C
00153 rounded_as (const X1& x1, const X2& x2) {
00154 typedef Return_op(Op,C) Rop;
00155 typedef rounding_helper<C> H;
00156 typename H::rounding old= H::get_rounding ();
00157 H::set_rounding (H::template translate<Mode>::val);
00158 C r= Rop::op (x1, x2);
00159 H::set_rounding (old);
00160 return r;
00161 }
00162 };
00163
00164
00165
00166
00167
00168 template<typename W, typename V> inline floating<W>
00169 view (const floating<V>& x) {
00170 return *((floating<W>*) ((void*) &x));
00171 }
00172
00173 template<mmx_rounding Mode>
00174 struct rounded_floating {
00175 template<typename Op, typename C>
00176 struct helper {};
00177
00178 template<typename Op, typename V>
00179 struct helper<Op,floating<V> > {
00180 typedef rounding_helper<floating<V> > H;
00181 typedef rnd_floating<V,H::template translate<Mode>::val> RV;
00182 typedef Return_op (Op, floating<RV> ) Rop;
00183
00184 static inline floating<V>
00185 op () {
00186 return view<V> (Rop::op ()); }
00187
00188 template<typename X1>
00189 static inline floating<V>
00190 op (const X1& x1) {
00191 return view<V> (Rop::op (x1)); }
00192
00193 template<typename X1, typename X2>
00194 static inline floating<V>
00195 op (const X1& x1, const X2& x2) {
00196 return view<V> (Rop::op (x1, x2)); }
00197
00198 template<typename W1>
00199 static inline floating<V>
00200 op (const floating<W1>& x1) {
00201 return view<V> (Rop::op (view<RV> (x1))); }
00202
00203 template<typename W1, typename W2>
00204 static inline floating<V>
00205 op (const floating<W1>& x1, const floating<W2>& x2) {
00206 return view<V> (Rop::op (view<RV> (x1), view<RV> (x2))); }
00207 };
00208
00209 template<typename Op, typename C>
00210 static inline C rounded_as () {
00211 return helper<Op,C>::op ();
00212 };
00213
00214 template<typename Op, typename C, typename X1>
00215 static inline C rounded_as (const X1& x1) {
00216 return helper<Op,C>::op (x1);
00217 };
00218
00219 template<typename Op, typename C, typename X1, typename X2>
00220 static inline C rounded_as (const X1& x1, const X2& x2) {
00221 return helper<Op,C>::op (x1, x2);
00222 };
00223 };
00224
00225
00226
00227
00228
00229 template<typename V>
00230 struct rounded_operations: public V {
00231
00232 template<typename Op, typename C>
00233 static inline C
00234 rounded () {
00235 return V::template rounded_as<Op,C> (); }
00236
00237 template<typename Op, typename C>
00238 static inline C
00239 rounded (const C& x1) {
00240 return V::template rounded_as<Op,C> (x1); }
00241
00242 template<typename Op, typename C>
00243 static inline C
00244 rounded (const C& x1, const C& x2) {
00245 return V::template rounded_as<Op,C> (x1, x2); }
00246
00247 TMPL static inline C log2 () {
00248 return rounded<log2_as_op,C> (); }
00249 TMPL static inline C pi () {
00250 return rounded<pi_as_op,C> (); }
00251 TMPL static inline C euler () {
00252 return rounded<euler_as_op,C> (); }
00253 TMPL static inline C catalan () {
00254 return rounded<catalan_as_op,C> (); }
00255 TMPL static inline C smallest () {
00256 return rounded<smallest_as_op,C> (); }
00257 TMPL static inline C largest () {
00258 return rounded<largest_as_op,C> (); }
00259 TMPL static inline C accuracy () {
00260 return rounded<accuracy_as_op,C> (); }
00261
00262 TMPL static inline C add (const C& x1, const C& x2) {
00263 return rounded<add_op> (x1, x2); }
00264 TMPL static inline C sub (const C& x1, const C& x2) {
00265 return rounded<sub_op> (x1, x2); }
00266 TMPL static inline C mul (const C& x1, const C& x2) {
00267 return rounded<mul_op> (x1, x2); }
00268 TMPL static inline C div (const C& x1, const C& x2) {
00269 return rounded<div_op> (x1, x2); }
00270 TMPL static inline C square (const C& x) {
00271 return rounded<square_op> (x); }
00272
00273 TMPL static inline C sqrt (const C& x1) {
00274 return rounded<sqrt_op> (x1); }
00275 TMPL static inline C cbrt (const C& x1) {
00276 return rounded<cbrt_op> (x1); }
00277 TMPL static inline C hypot (const C& x1, const C& x2) {
00278 return rounded<hypot_op> (x1, x2); }
00279 TMPL static inline C exp (const C& x1) {
00280 return rounded<exp_op> (x1); }
00281 TMPL static inline C exp2 (const C& x1) {
00282 return rounded<exp2_op> (x1); }
00283 TMPL static inline C log (const C& x1) {
00284 return rounded<log_op> (x1); }
00285 TMPL static inline C log2 (const C& x1) {
00286 return rounded<log2_op> (x1); }
00287 TMPL static inline C log10 (const C& x1) {
00288 return rounded<log10_op> (x1); }
00289 TMPL static inline C pow (const C& x1, const C& x2) {
00290 return rounded<pow_op> (x1, x2); }
00291 TMPL static inline C cos (const C& x1) {
00292 return rounded<cos_op> (x1); }
00293 TMPL static inline C sin (const C& x1) {
00294 return rounded<sin_op> (x1); }
00295 TMPL static inline C tan (const C& x1) {
00296 return rounded<tan_op> (x1); }
00297 TMPL static inline C acos (const C& x1) {
00298 return rounded<acos_op> (x1); }
00299 TMPL static inline C asin (const C& x1) {
00300 return rounded<asin_op> (x1); }
00301 TMPL static inline C atan (const C& x1) {
00302 return rounded<atan_op> (x1); }
00303 TMPL static inline C atan2 (const C& x1, const C& x2) {
00304 return rounded<atan2_op> (x1, x2); }
00305 TMPL static inline C cosh (const C& x1) {
00306 return rounded<cosh_op> (x1); }
00307 TMPL static inline C sinh (const C& x1) {
00308 return rounded<sinh_op> (x1); }
00309 TMPL static inline C tanh (const C& x1) {
00310 return rounded<tanh_op> (x1); }
00311 TMPL static inline C acosh (const C& x1) {
00312 return rounded<acosh_op> (x1); }
00313 TMPL static inline C asinh (const C& x1) {
00314 return rounded<asinh_op> (x1); }
00315 TMPL static inline C atanh (const C& x1) {
00316 return rounded<atanh_op> (x1); }
00317
00318 TMPL static inline C gamma (const C& x1) {
00319 return rounded<gamma> (x1); }
00320 TMPL static inline C zeta (const C& x1) {
00321 return rounded<zeta> (x1); }
00322 TMPL static inline C erf (const C& x1) {
00323 return rounded<erf> (x1); }
00324
00325 TMPL static inline C floor (const C& x1) {
00326 return rounded<floor_op> (x1); }
00327 TMPL static inline C trunc (const C& x1) {
00328 return rounded<trunc_op> (x1); }
00329 TMPL static inline C ceil (const C& x1) {
00330 return rounded<ceil_op> (x1); }
00331 TMPL static inline C round (const C& x1) {
00332 return rounded<round_op> (x1); }
00333
00334 template<typename C, typename X1>
00335 static inline C
00336 as (const X1& x1) {
00337 return V::template rounded_as<as_op,C> (x1); }
00338
00339 template<typename C, typename X1>
00340 static inline C
00341 abs_as (const X1& x1) {
00342 return V::template rounded_as<abs_as_op,C> (x1); }
00343
00344 template<typename C, typename X1, typename X2>
00345 static inline C
00346 add_as (const X1& x1, const X2& x2) {
00347 return V::template rounded_as<add_op,C> (x1, x2); }
00348
00349 template<typename C, typename X1, typename X2>
00350 static inline C
00351 sub_as (const X1& x1, const X2& x2) {
00352 return V::template rounded_as<sub_op,C> (x1, x2); }
00353
00354 };
00355
00356
00357
00358
00359
00360 template<typename V>
00361 struct rounded_opposite: public V {
00362 TMPL static inline C add (const C& x1, const C& x2) {
00363 C tmp= -x1; tmp -= x2; return -tmp; }
00364 TMPL static inline C sub (const C& x1, const C& x2) {
00365 C tmp= x2; tmp -= x1; return -tmp; }
00366 TMPL static inline C mul (const C& x1, const C& x2) {
00367 C tmp= -x1; tmp *= x2; return -tmp; }
00368 TMPL static inline C div (const C& x1, const C& x2) {
00369 C tmp= -x1; tmp /= x2; return -tmp; }
00370 TMPL static inline C square (const C& x1) {
00371 C tmp= -x1; tmp *= x1; return -tmp; }
00372 };
00373
00374
00375
00376
00377
00378 STMPL
00379 struct rounding_helper<double> {
00380 typedef int rounding_mode;
00381 template<mmx_rounding Mode, typename K= nat> struct translate {
00382 static const rounding_mode val= FE_TONEAREST; };
00383 template<typename K> struct translate<MMX_ROUND_NEAR,K> {
00384 static const rounding_mode val= FE_TONEAREST; };
00385 template<typename K> struct translate<MMX_ROUND_ZERO,K> {
00386 static const rounding_mode val= FE_TOWARDZERO; };
00387 template<typename K> struct translate<MMX_ROUND_UP,K> {
00388 static const rounding_mode val= FE_UPWARD; };
00389 template<typename K> struct translate<MMX_ROUND_DOWN,K> {
00390 static const rounding_mode val= FE_DOWNWARD; };
00391 static inline rounding_mode get_rounding () {
00392 return fegetround (); }
00393 static inline void set_rounding (rounding_mode r) {
00394 fesetround (r); }
00395 static const bool exact_type = false;
00396 typedef rounded_operations<rounded_global<MMX_ROUND_UP> > UV;
00397 typedef rounded_operations<rounded_global<MMX_ROUND_DOWN> > DV;
00398 };
00399
00400 template<typename V>
00401 struct rounding_helper<floating<V> > {
00402 typedef mpfr_rnd_t rounding_mode;
00403 template<mmx_rounding Mode, typename K= nat> struct translate {
00404 static const rounding_mode val= GMP_RNDN; };
00405 template<typename K> struct translate<MMX_ROUND_NEAR,K> {
00406 static const rounding_mode val= GMP_RNDN; };
00407 template<typename K> struct translate<MMX_ROUND_ZERO,K> {
00408 static const rounding_mode val= GMP_RNDZ; };
00409 template<typename K> struct translate<MMX_ROUND_UP,K> {
00410 static const rounding_mode val= GMP_RNDU; };
00411 template<typename K> struct translate<MMX_ROUND_DOWN,K> {
00412 static const rounding_mode val= GMP_RNDD; };
00413 static inline rounding_mode get_rounding () {
00414 return mmx_rounding_mode; }
00415 static inline void set_rounding (rounding_mode r) {
00416 mmx_rounding_mode= r; }
00417 static const bool exact_type = false;
00418 typedef rounded_operations<rounded_floating<MMX_ROUND_UP> > UV;
00419 typedef rounded_operations<rounded_floating<MMX_ROUND_DOWN> > DV;
00420 };
00421
00422 #undef TMPL
00423 }
00424 #endif // __MMX__ROUNDED_HPP