00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #ifndef __MMX_BALL_ROUNDED_HPP
00014 #define __MMX_BALL_ROUNDED_HPP
00015 #include <numerix/rounded.hpp>
00016 namespace mmx {
00017 #define TMPL template<typename C,typename R,typename V>
00018 #define Ball ball<C,R,V>
00019
00020 TMPL class ball;
00021 TMPL inline C center (const Ball& z);
00022 TMPL inline R radius (const Ball& z);
00023 TMPL inline C& center (Ball& z);
00024 TMPL inline R& radius (Ball& z);
00025
00026
00027
00028
00029
00030 struct ball_rounded {};
00031
00032 template<typename C>
00033 struct ball_variant_helper {
00034 typedef ball_rounded BV;
00035 };
00036
00037
00038
00039
00040
00041 struct ball_rounding {};
00042
00043 template<typename W>
00044 struct implementation<ball_rounding,W,ball_rounded> {
00045
00046 TMPL static inline void
00047 add_rounding_error (Ball& d) {
00048 typedef Round_up(R) Up;
00049 R err= Up::template as<R> (rounding_error (center (d)));
00050 radius (d)= Up::add (radius (d), err);
00051 }
00052
00053 TMPL static inline void
00054 add_additive_error (Ball& d) {
00055 typedef Round_up(R) Up;
00056 R err= Up::template as<R> (additive_error (center (d)));
00057 radius (d)= Up::add (radius (d), err);
00058 }
00059
00060 TMPL static inline void
00061 add_multiplicative_error (Ball& d) {
00062 typedef Round_up(R) Up;
00063 R err= Up::template as<R> (multiplicative_error (center (d)));
00064 radius (d)= Up::add (radius (d), err);
00065 }
00066
00067 TMPL static inline void
00068 add_elementary_error (Ball& d) {
00069 typedef Round_up(R) Up;
00070 R err= Up::template as<R> (elementary_error (center (d)));
00071 radius (d)= Up::add (radius (d), err);
00072 }
00073
00074 TMPL static inline C
00075 lower (const Ball& z) {
00076 typedef Round_down(C) Down;
00077 return Down::template sub_as<C> (center (z), radius (z));
00078 }
00079
00080 TMPL static inline C
00081 upper (const Ball& z) {
00082 typedef Round_up(C) Up;
00083 return Up::template add_as<C> (center (z), radius (z));
00084 }
00085
00086 TMPL static inline void
00087 make_interval (Ball& d, const C& l, const C& r) {
00088 typedef Round_up(C) Up;
00089 C c= decexp2 (l + r, 1);
00090 center (d)= c;
00091 radius (d)= Up::template as<R> (max (Up::sub (r, c), Up::sub (c, l)));
00092 }
00093
00094 TMPL static inline R
00095 abs_down (const Ball& z) {
00096 typedef Round_down(R) Down;
00097 R a= Down::template abs_as<R> (center (z));
00098 if (a <= radius (z)) return R(0);
00099 else return Down::sub (a, radius (z));
00100 }
00101
00102 TMPL static inline R
00103 abs_up (const Ball& z) {
00104 typedef Round_up(R) Up;
00105 return Up::add (Up::template abs_as<R> (center(z)), radius(z));
00106 }
00107
00108 TMPL static inline R
00109 bnd_down (const Ball& z) {
00110 typedef Round_down(R) Down;
00111 return Down::sub (Down::template as<R> (center(z)), radius(z));
00112 }
00113
00114 TMPL static inline R
00115 bnd_up (const Ball& z) {
00116 typedef Round_up(R) Up;
00117 return Up::add (Up::template as<R> (center(z)), radius(z));
00118 }
00119
00120 };
00121
00122
00123
00124
00125
00126 struct ball_additive {};
00127
00128 template<typename W>
00129 struct implementation<ball_additive,W,ball_rounded> {
00130 typedef implementation<ball_rounding,W> Rnd;
00131
00132 TMPL static inline void
00133 neg (Ball& d, const Ball& z) {
00134 center (d)= -center (z);
00135 radius (d)= radius (z);
00136 }
00137
00138 TMPL static inline void
00139 add (Ball& d, const Ball& z1, const Ball& z2) {
00140 typedef Round_up(R) Up;
00141 center (d)= center (z1) + center (z2);
00142 radius (d)= Up::add (radius (z1), radius (z2));
00143 Rnd::add_additive_error (d);
00144 }
00145
00146 TMPL static inline void
00147 add (Ball& d, const C& z1, const Ball& z2) {
00148 center (d)= z1 + center (z2);
00149 radius (d)= radius (z2);
00150 Rnd::add_additive_error (d);
00151 }
00152
00153 TMPL static inline void
00154 add (Ball& d, const Ball& z1, const C& z2) {
00155 center (d)= center (z1) + z2;
00156 radius (d)= radius (z1);
00157 Rnd::add_additive_error (d);
00158 }
00159
00160 TMPL static inline void
00161 sub (Ball& d, const Ball& z1, const Ball& z2) {
00162 typedef Round_up(R) Up;
00163 center (d)= center (z1) - center (z2);
00164 radius (d)= Up::add (radius (z1), radius (z2));
00165 Rnd::add_additive_error (d);
00166 }
00167
00168 TMPL static inline void
00169 sub (Ball& d, const C& z1, const Ball& z2) {
00170 center (d)= z1 - center (z2);
00171 radius (d)= radius (z2);
00172 Rnd::add_additive_error (d);
00173 }
00174
00175 TMPL static inline void
00176 sub (Ball& d, const Ball& z1, const C& z2) {
00177 center (d)= center (z1) - z2;
00178 radius (d)= radius (z1);
00179 Rnd::add_additive_error (d);
00180 }
00181
00182 };
00183
00184
00185
00186
00187
00188 struct ball_multiplicative {};
00189
00190 template<typename W>
00191 struct implementation<ball_multiplicative,W,ball_rounded> {
00192 typedef implementation<ball_rounding,W> Rnd;
00193
00194 TMPL static inline void
00195 mul (Ball& d, const Ball& z1, const Ball& z2) {
00196 typedef Round_up(R) Up;
00197 R a1= Rnd::abs_up (z1);
00198 R a2= Up::template abs_as<R> (center (z2));
00199 center (d)= center (z1) * center (z2);
00200 radius (d)= Up::add (Up::mul (a1, radius (z2)), Up::mul (radius (z1), a2));
00201 Rnd::add_multiplicative_error (d);
00202 }
00203
00204 TMPL static inline void
00205 mul (Ball& d, const C& z1, const Ball& z2) {
00206 typedef Round_up(R) Up;
00207 R a1= Up::template abs_as<R> (z1);
00208 center (d)= z1 * center (z2);
00209 radius (d)= Up::mul (a1, radius (z2));
00210 Rnd::add_multiplicative_error (d);
00211 }
00212
00213 TMPL static inline void
00214 mul (Ball& d, const Ball& z1, const C& z2) {
00215 typedef Round_up(R) Up;
00216 R a2= Up::template abs_as<R> (z2);
00217 center (d)= center (z1) * z2;
00218 radius (d)= Up::mul (radius (z1), a2);
00219 Rnd::add_multiplicative_error (d);
00220 }
00221
00222 TMPL static inline void
00223 square (Ball& d, const Ball& z) {
00224 typedef Round_up(R) Up;
00225 R a= Up::template abs_as<R> (center (z));
00226 center (d)= square_op::op (center (z));
00227 radius (d)= Up::mul (Up::add (a + a, radius (z)), radius (z));
00228 Rnd::add_multiplicative_error (d);
00229 }
00230
00231 TMPL static inline void
00232 invert (Ball& d, const Ball& z) {
00233 typedef Round_up(R) Up;
00234 typedef Round_down(R) Down;
00235 center (d)= invert_op::op (center (z));
00236 radius (d)= Up::div (radius (z), Down::square (Rnd::abs_down (z)));
00237 Rnd::add_multiplicative_error (d);
00238 }
00239
00240 TMPL static inline void
00241 invert (Ball& d, const C& z) {
00242 center (d)= invert_op::op (z);
00243 radius (d)= multiplicative_error (center (d));
00244 }
00245
00246 TMPL static inline void
00247 div (Ball& d, const Ball& z1, const Ball& z2) {
00248 d= z1 * invert_op::op (z2);
00249 }
00250
00251 TMPL static inline void
00252 div (Ball& d, const Ball& z1, const C& z2) {
00253 d= z1 * invert_op::op (z2);
00254 }
00255
00256 TMPL static inline void
00257 div (Ball& d, const C& z1, const Ball& z2) {
00258 d= z1 * invert_op::op (z2);
00259 }
00260
00261 };
00262
00263
00264
00265
00266
00267 struct ball_root {};
00268
00269 template<typename W>
00270 struct implementation<ball_root,W,ball_rounded> {
00271 typedef implementation<ball_rounding,W> Rnd;
00272
00273 TMPL static void
00274 sqrt (Ball& d, const Ball& z) {
00275 typedef Round_up(R) Up;
00276 typedef Round_down(R) Down;
00277 if (is_negative_or_zero (z)) d= Nan (Ball);
00278 else {
00279 R aux= Down::sqrt (Rnd::abs_down (z));
00280 center (d)= sqrt_op::op (center (z));
00281 radius (d)= Up::div (radius (z), aux + aux);
00282 Rnd::add_multiplicative_error (d);
00283 }
00284 }
00285
00286 TMPL static void
00287 hypot (Ball& d, const Ball& x, const Ball& y) {
00288 typedef Round_up(R) Up;
00289 center (d)= hypot_op::op (center (x), center (y));
00290 radius (d)= Up::hypot (radius (x), radius (y));
00291 Rnd::add_multiplicative_error (d);
00292 }
00293
00294 };
00295
00296
00297
00298
00299
00300 struct ball_elementary {};
00301
00302 template<typename W>
00303 struct implementation<ball_elementary,W,ball_rounded> {
00304 typedef implementation<ball_rounding,W> Rnd;
00305
00306 TMPL static void
00307 exp (Ball& d, const Ball& z) {
00308 typedef Round_up(R) Up;
00309 center (d)= exp_op::op (center (z));
00310 radius (d)= Up::mul (Up::exp (Rnd::bnd_up (z)), radius (z));
00311 Rnd::add_elementary_error (d);
00312 }
00313
00314 TMPL static void
00315 log (Ball& d, const Ball& z) {
00316 typedef Round_up(R) Up;
00317 if (is_negative_or_zero (z)) d= Nan (Ball);
00318 else {
00319 center (d)= log_op::op (center (z));
00320 radius (d)= Up::div (radius (z), Rnd::bnd_down (z));
00321 Rnd::add_elementary_error (d);
00322 }
00323 }
00324
00325 TMPL static inline void
00326 pow (Ball& d, const Ball& z, const Ball& y) {
00327 d= exp_op::op (y * log_op::op (z));
00328 }
00329
00330 TMPL static void
00331 cos (Ball& d, const Ball& z) {
00332 typedef Round_up(R) Up;
00333 typedef Round_down(R) Down;
00334 if (radius (z) >= Down::template pi<R> ()) {
00335 center (d)= promote (0, center (z));
00336 radius (d)= promote (1, radius (z));
00337 }
00338 else {
00339 center (d)= cos_op::op (center (z));
00340 radius (d)= radius (z);
00341 Rnd::add_elementary_error (d);
00342 if (radius (z) < abs (center (d))) {
00343 R u1= Up::sin (Rnd::bnd_down (z));
00344 R u2= Up::sin (Rnd::bnd_up (z));
00345 radius (d)= Up::mul (max (abs (u1), abs (u2)), radius (z));
00346 Rnd::add_elementary_error (d);
00347 }
00348 }
00349 }
00350
00351 TMPL static void
00352 sin (Ball& d, const Ball& z) {
00353 typedef Round_up(R) Up;
00354 typedef Round_down(R) Down;
00355 if (radius (z) >= Down::template pi<R> ()) {
00356 center (d)= promote (0, center (z));
00357 radius (d)= promote (1, radius (z));
00358 }
00359 else {
00360 center (d)= sin_op::op (center (z));
00361 radius (d)= radius (z);
00362 Rnd::add_elementary_error (d);
00363 if (radius (z) < abs (center (d))) {
00364 R u1= Up::cos (Rnd::bnd_down (z));
00365 R u2= Up::cos (Rnd::bnd_up (z));
00366 radius (d)= Up::mul (max (abs (u1), abs (u2)), radius (z));
00367 Rnd::add_elementary_error (d);
00368 }
00369 }
00370 }
00371
00372 TMPL static inline void
00373 tan (Ball& d, const Ball& z) {
00374 d= sin_op::op (z) / cos_op::op (z);
00375 }
00376
00377 TMPL static inline void
00378 cosh (Ball& d, const Ball& z) {
00379 typedef Round_up(R) Up;
00380 center (d)= cosh_op::op (center (z));
00381 radius (d)= Up::mul (Up::sinh (Rnd::bnd_up (abs (z))), radius (z));
00382 Rnd::add_elementary_error (d);
00383 }
00384
00385 TMPL static inline void
00386 sinh (Ball& d, const Ball& z) {
00387 typedef Round_up(R) Up;
00388 center (d)= sinh_op::op (center (z));
00389 radius (d)= Up::mul (Up::cosh (Rnd::bnd_up (abs (z))), radius (z));
00390 Rnd::add_elementary_error (d);
00391 }
00392
00393 TMPL static inline void
00394 tanh (Ball& d, const Ball& z) {
00395 d= sinh_op::op (z) / cosh_op::op (z);
00396 }
00397
00398 TMPL static void
00399 acos (Ball& d, const Ball& z) {
00400 typedef Round_up(R) Up;
00401 typedef Round_down(R) Down;
00402 Ball u= promote (1, z) - square_op::op (z);
00403 R v= Rnd::bnd_down (u);
00404 if (v <= 0) d= Nan (Ball);
00405 else {
00406 center (d)= acos_op::op (center (z));
00407 radius (d)= Up::div (radius (z), Down::sqrt (v));
00408 Rnd::add_elementary_error (d);
00409 }
00410 }
00411
00412 TMPL static void
00413 asin (Ball& d, const Ball& z) {
00414 typedef Round_up(R) Up;
00415 typedef Round_down(R) Down;
00416 Ball u= promote (1, z) - square_op::op (z);
00417 R v= Rnd::bnd_down (u);
00418 if (v <= 0) d= Nan (Ball);
00419 else {
00420 center (d)= asin_op::op (center (z));
00421 radius (d)= Up::div (radius (z), Down::sqrt (v));
00422 Rnd::add_elementary_error (d);
00423 }
00424 }
00425
00426 TMPL static void
00427 atan (Ball& d, const Ball& z) {
00428 typedef Round_up(R) Up;
00429 typedef Round_down(R) Down;
00430 center (d)= atan_op::op (center (z));
00431 R u= Rnd::bnd_down (z);
00432 R v= Down::add (promote (1, radius (z)), Down::square (u));
00433 if (u <= 0) radius (d)= radius (z);
00434 else radius (d)= Up::div (radius (z), v);
00435 Rnd::add_elementary_error (d);
00436 }
00437
00438 TMPL static void
00439 atan2 (Ball& d, const Ball& y, const Ball& x) {
00440 typedef Round_up(R) Up;
00441 Ball z= square_op::op (x) + square_op::op (y);
00442 if (z <= 0) d= Nan (Ball);
00443 else {
00444 R b = Up::div (promote (1, radius (y)), Rnd::bnd_down (z));
00445 R dx= Up::mul (Rnd::bnd_up (abs (y)), b);
00446 R dy= Up::mul (Rnd::bnd_up (abs (x)), b);
00447 center (d)= atan2_op::op (center (y), center (x));
00448 radius (d)= Up::add (Up::mul (dx, radius (x)), Up::mul (dy, radius (y)));
00449 Rnd::add_elementary_error (d);
00450 }
00451 }
00452
00453 };
00454
00455
00456
00457
00458
00459 struct ball_ordered {};
00460
00461 template<typename W>
00462 struct implementation<ball_ordered,W,ball_rounded> {
00463 typedef implementation<ball_rounding,W> Rnd;
00464
00465 TMPL static void
00466 min (Ball& d, const Ball& z1, const Ball& z2) {
00467 if (z1 < z2) d= z1;
00468 else if (z2 < z1) d= z2;
00469 else {
00470 center (d)= min_op::op (center (z1), center (z2));
00471 radius (d)= max_op::op (radius (z1), radius (z2));
00472 }
00473 }
00474
00475 TMPL static void
00476 max (Ball& d, const Ball& z1, const Ball& z2) {
00477 if (z1 > z2) d= z1;
00478 else if (z2 > z1) d= z2;
00479 else {
00480 center (d)= max_op::op (center (z1), center (z2));
00481 radius (d)= max_op::op (radius (z1), radius (z2));
00482 }
00483 }
00484
00485 TMPL static void
00486 floor (Ball& d, const Ball& z) {
00487 typedef Round_up(C) Up;
00488 typedef Round_down(C) Down;
00489 C l= Down::floor (lower (z));
00490 C r= Up ::floor (upper (z));
00491 Rnd::make_interval (d, l, r);
00492 }
00493
00494 TMPL static void
00495 trunc (Ball& d, const Ball& z) {
00496 typedef Round_up(C) Up;
00497 typedef Round_down(C) Down;
00498 C l= Down::trunc (lower (z));
00499 C r= Up ::trunc (upper (z));
00500 Rnd::make_interval (d, l, r);
00501 }
00502
00503 TMPL static void
00504 ceil (Ball& d, const Ball& z) {
00505 typedef Round_up(C) Up;
00506 typedef Round_down(C) Down;
00507 C l= Down::ceil (lower (z));
00508 C r= Up ::ceil (upper (z));
00509 Rnd::make_interval (d, l, r);
00510 }
00511
00512 TMPL static void
00513 round (Ball& d, const Ball& z) {
00514 typedef Round_up(C) Up;
00515 typedef Round_down(C) Down;
00516 C l= Down::round (lower (z));
00517 C r= Up ::round (upper (z));
00518 Rnd::make_interval (d, l, r);
00519 }
00520
00521 };
00522
00523
00524
00525
00526
00527 struct ball_abs {};
00528
00529 template<typename W>
00530 struct implementation<ball_abs,W,ball_rounded> {
00531 typedef implementation<ball_rounding,W> Rnd;
00532
00533 template<typename C ,typename R ,typename V,
00534 typename C2,typename V2> static inline void
00535 abs (ball<C2,R,V2>& d, const Ball& z) {
00536 center (d)= abs_op::op (center (z));
00537 radius (d)= radius (z);
00538 Rnd::add_multiplicative_error (d);
00539 }
00540
00541 TMPL static inline void
00542 abs (Ball& d, const Ball& z) {
00543 center (d)= abs_op::op (center (z));
00544 radius (d)= radius (z);
00545 }
00546
00547 };
00548
00549
00550
00551
00552
00553 struct ball_shift {};
00554
00555 template<typename W>
00556 struct implementation<ball_shift,W,ball_rounded> {
00557 typedef implementation<ball_rounding,W> Rnd;
00558
00559 TMPL static inline void
00560 shiftl (Ball& d, const xint& shift) {
00561 center (d)= incexp2 (center (d), shift);
00562 radius (d)= incexp2 (radius (d), shift);
00563 Rnd::add_multiplicative_error (d);
00564 }
00565
00566 TMPL static inline void
00567 shiftl (Ball& d, const Ball& z, const xint& shift) {
00568 center (d)= incexp2 (center (z), shift);
00569 radius (d)= incexp2 (radius (z), shift);
00570 Rnd::add_multiplicative_error (d);
00571 }
00572
00573 TMPL static inline void
00574 shiftr (Ball& d, const xint& shift) {
00575 center (d)= decexp2 (center (d), shift);
00576 radius (d)= decexp2 (radius (d), shift);
00577 Rnd::add_multiplicative_error (d);
00578 }
00579
00580 TMPL static inline void
00581 shiftr (Ball& d, const Ball& z, const xint& shift) {
00582 center (d)= decexp2 (center (z), shift);
00583 radius (d)= decexp2 (radius (z), shift);
00584 Rnd::add_multiplicative_error (d);
00585 }
00586
00587 };
00588
00589 #undef TMPL
00590 #undef Ball
00591 }
00592 #endif // __MMX_BALL_ROUNDED_HPP