00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef __MMX_BALL_ROUGH_HPP
00022 #define __MMX_BALL_ROUGH_HPP
00023 #include <numerix/ball_rounded.hpp>
00024 namespace mmx {
00025 #define TMPL template<typename C,typename R,typename V>
00026 #define Ball ball<C,R,V>
00027
00028
00029
00030
00031
00032 struct ball_rough {};
00033
00034 template<typename F, typename V>
00035 struct implementation<F,V,ball_rough>:
00036 public implementation<F,V,ball_rounded> {};
00037
00038 STMPL
00039 struct ball_variant_helper<double> {
00040 typedef ball_rough BV;
00041 };
00042
00043
00044
00045
00046
00047 TMPL inline void
00048 add_rough_additive_error (Ball& d) {
00049 radius (d) += additive_error (center (d));
00050 radius (d) *= (1 + incexp2 (Accuracy (R), 2));
00051 }
00052
00053 template<typename V> inline void
00054 add_rough_additive_error (ball<double,double,V>& d) {
00055 const double alpha= 1.000000000000008882;
00056 radius (d)= alpha * (radius (d) + ldexp (abs (center (d)), -50));
00057 }
00058
00059 template<typename W>
00060 struct implementation<ball_additive,W,ball_rough> {
00061
00062 TMPL static inline void
00063 neg (Ball& d, const Ball& z) {
00064 center (d)= -center (z);
00065 radius (d)= radius (z);
00066 }
00067
00068 TMPL static inline void
00069 add (Ball& d, const Ball& z1, const Ball& z2) {
00070 center (d)= center (z1) + center (z2);
00071 radius (d)= radius (z1) + radius (z2);
00072 add_rough_additive_error (d);
00073 }
00074
00075 TMPL static inline void
00076 add (Ball& d, const C& z1, const Ball& z2) {
00077 center (d)= z1 + center (z2);
00078 radius (d)= radius (z2);
00079 add_rough_additive_error (d);
00080 }
00081
00082 TMPL static inline void
00083 add (Ball& d, const Ball& z1, const C& z2) {
00084 center (d)= center (z1) + z2;
00085 radius (d)= radius (z1);
00086 add_rough_additive_error (d);
00087 }
00088
00089 TMPL static inline void
00090 sub (Ball& d, const Ball& z1, const Ball& z2) {
00091 center (d)= center (z1) - center (z2);
00092 radius (d)= radius (z1) + radius (z2);
00093 add_rough_additive_error (d);
00094 }
00095
00096 TMPL static inline void
00097 sub (Ball& d, const C& z1, const Ball& z2) {
00098 center (d)= z1 - center (z2);
00099 radius (d)= radius (z2);
00100 add_rough_additive_error (d);
00101 }
00102
00103 TMPL static inline void
00104 sub (Ball& d, const Ball& z1, const C& z2) {
00105 center (d)= center (z1) - z2;
00106 radius (d)= radius (z1);
00107 add_rough_additive_error (d);
00108 }
00109
00110 };
00111
00112
00113
00114
00115
00116 TMPL inline void
00117 add_rough_multiplicative_error (Ball& d) {
00118 radius (d) += multiplicative_error (center (d));
00119 radius (d) *= (1 + incexp2 (Accuracy (R), 4));
00120 radius (d) += incexp2 (Smallest (R), 4);
00121 }
00122
00123 template<typename V> inline void
00124 add_rough_multiplicative_error (ball<double,double,V>& d) {
00125 const double alpha= 1.000000000000001777;
00126 const double beta = 3.952525166729972354e-323;
00127 radius (d)= alpha * (radius (d) + ldexp (abs (center (d)), -50)) + beta;
00128 }
00129
00130 template<typename C> inline C
00131 rough_next_above (const C& z) {
00132 return (1 + Accuracy (C)) * z;
00133 }
00134
00135 inline double
00136 rough_next_above (const double& z) {
00137 const double alpha= 1.000000000000000223;
00138 return alpha * z;
00139 }
00140
00141 template<typename W>
00142 struct implementation<ball_multiplicative,W,ball_rough> {
00143
00144 TMPL static inline void
00145 mul (Ball& d, const Ball& z1, const Ball& z2) {
00146 R a1= as<R> (abs (center (z1))) + radius (z1);
00147 R a2= as<R> (abs (center (z2)));
00148 center (d)= center (z1) * center (z2);
00149 radius (d)= a1 * radius (z2) + radius (z1) * a2;
00150 add_rough_multiplicative_error (d);
00151 }
00152
00153 TMPL static inline void
00154 mul (Ball& d, const C& z1, const Ball& z2) {
00155 R a1= as<R> (abs (z1));
00156 center (d)= z1 * center (z2);
00157 radius (d)= a1 * radius (z2);
00158 add_rough_multiplicative_error (d);
00159 }
00160
00161 TMPL static inline void
00162 mul (Ball& d, const Ball& z1, const C& z2) {
00163 R a2= as<R> (abs (z2));
00164 center (d)= center (z1) * z2;
00165 radius (d)= radius (z1) * a2;
00166 add_rough_multiplicative_error (d);
00167 }
00168
00169 TMPL static inline void
00170 square (Ball& d, const Ball& z) {
00171 R a= as<R> (abs (center (z)));
00172 center (d)= square_op::op (center (z));
00173 radius (d)= (a + a + radius (z)) * radius (z);
00174 add_rough_multiplicative_error (d);
00175 }
00176
00177 TMPL static inline void
00178 invert (Ball& d, const Ball& z) {
00179 R a= max (as<R> (abs (center (z))) - rough_next_above (radius (z)), R(0));
00180 center (d)= invert_op::op (center (z));
00181 radius (d)= radius (z) / square_op::op (a);
00182 add_rough_multiplicative_error (d);
00183 }
00184
00185 TMPL static inline void
00186 invert (Ball& d, const C& z) {
00187 center (d)= invert_op::op (z);
00188 radius (d)= multiplicative_error (center (d));
00189 }
00190
00191 TMPL static inline void
00192 div (Ball& d, const Ball& z1, const Ball& z2) {
00193 d= z1 * invert_op::op (z2);
00194 }
00195
00196 TMPL static inline void
00197 div (Ball& d, const Ball& z1, const C& z2) {
00198 d= z1 * invert_op::op (z2);
00199 }
00200
00201 TMPL static inline void
00202 div (Ball& d, const C& z1, const Ball& z2) {
00203 d= z1 * invert_op::op (z2);
00204 }
00205
00206 };
00207
00208
00209
00210
00211
00212 template<typename W>
00213 struct implementation<ball_root,W,ball_rough> {
00214
00215 TMPL static void
00216 sqrt (Ball& d, const Ball& z) {
00217 R l= as<R> (center (z)) - rough_next_above (radius (z));
00218 if (l <= 0) d= Nan (Ball);
00219 else {
00220 center (d)= sqrt_op::op (center (z));
00221 radius (d)= decexp2 (radius (z)) / sqrt_op::op (l);
00222 add_rough_multiplicative_error (d);
00223 }
00224 }
00225
00226 TMPL static void
00227 hypot (Ball& d, const Ball& x, const Ball& y) {
00228 center (d)= hypot_op::op (center (x), center (y));
00229 radius (d)= hypot_op::op (radius (x), radius (y));
00230 add_rough_multiplicative_error (d);
00231 }
00232
00233 };
00234
00235
00236
00237
00238
00239 template<typename W>
00240 struct implementation<ball_shift,W,ball_rough> {
00241
00242 TMPL static inline void
00243 shiftl (Ball& d, const xint& shift) {
00244 center (d)= incexp2 (center (d), shift);
00245 radius (d)= incexp2 (radius (d), shift) + Smallest (R);
00246 }
00247
00248 TMPL static inline void
00249 shiftl (Ball& d, const Ball& z, const xint& shift) {
00250 center (d)= incexp2 (center (z), shift);
00251 radius (d)= incexp2 (radius (z), shift) + Smallest (R);
00252 }
00253
00254 TMPL static inline void
00255 shiftr (Ball& d, const xint& shift) {
00256 center (d)= decexp2 (center (d), shift);
00257 radius (d)= decexp2 (radius (d), shift) + Smallest (R);
00258 }
00259
00260 TMPL static inline void
00261 shiftr (Ball& d, const Ball& z, const xint& shift) {
00262 center (d)= decexp2 (center (z), shift);
00263 radius (d)= decexp2 (radius (z), shift) + Smallest (R);
00264 }
00265
00266 };
00267
00268 #undef TMPL
00269 #undef Ball
00270 }
00271 #endif // __MMX_BALL_ROUGH_HPP