1 /* 2 * triple_double.h 3 * 4 * This file contains useful tools and data for triple double data representation. 5 * 6 */ 7 8 #ifndef TRIPLE_DOUBLE_H 9 #define TRIPLE_DOUBLE_H 1 10 11 #include "scs_lib/scs.h" 12 #include "scs_lib/scs_private.h" 13 14 /* undef all the variables that might have been defined in 15 scs_lib/scs_private.h */ 16 #undef VERSION 17 #undef PACKAGE 18 #undef HAVE_GMP_H 19 #undef HAVE_MPFR_H 20 #undef HAVE_MATHLIB_H 21 /* then include the proper definitions */ 22 #include "crlibm_config.h" 23 24 #ifdef HAVE_INTTYPES_H 25 #include <inttypes.h> 26 #endif 27 28 /* Set to O for larger but faster functions. 29 As it only impacts the second step, smaller is preferred */ 30 #define TRIPLEDOUBLE_AS_FUNCTIONS 0 31 32 33 #define Renormalize3(resh, resm, resl, ah, am, al) DoRenormalize3(resh, resm, resl, ah, am, al) 34 /*extern void Renormalize3(double* resh, double* resm, double* resl, double ah, double am, double al) ;*/ 35 36 #if TRIPLEDOUBLE_AS_FUNCTIONS 37 extern void Mul23(double* resh, double* resm, double* resl, double ah, double al, double bh, double bl); 38 extern void Mul233(double* resh, double* resm, double* resl, double ah, double al, double bh, double bm, double bl); 39 extern void Mul33(double* resh, double* resm, double* resl, double ah, double am, double al, double bh, double bm, double bl); 40 extern void Mul133(double* resh, double* resm, double* resl, double a, double bh, double bm, double bl); 41 extern void Mul123(double* resh, double* resm, double* resl, double a, double bh, double bl); 42 extern void Sqrt13(double* resh, double* resm, double* resl, double x); 43 extern void Recpr33(double* resh, double* resm, double* resl, double dh, double dm, double dl); 44 #else 45 #define Mul23(resh, resm, resl, ah, al, bh, bl) DoMul23(resh, resm, resl, ah, al, bh, bl) 46 #define Mul233(resh, resm, resl, ah, al, bh, bm, bl) DoMul233(resh, resm, resl, ah, al, bh, bm, bl) 47 #define Mul33(resh, resm, resl, ah, am, al, bh, bm, bl) DoMul33(resh, resm, resl, ah, am, al, bh, bm, bl) 48 #define Mul133(resh, resm, resl, a, bh, bm, bl) DoMul133(resh, resm, resl, a, bh, bm, bl) 49 #define Mul123(resh, resm, resl, a, bh, bl) DoMul123(resh, resm, resl, a, bh, bl) 50 #define Sqrt13(resh, resm, resl , x) DoSqrt13(resh, resm, resl , x) 51 #define Recpr33(resh, resm, resl, dh, dm, dl) DoRecpr33(resh, resm, resl, dh, dm, dl) 52 #endif 53 54 55 56 /* Renormalize3 57 58 Procedure for renormalizing a triple double number, i.e. 59 computing exactly an equivalent sum of three non-overlapping 60 double numbers 61 62 63 Arguments: a triple double number ah, am, al 64 65 Results: a triple double number resh, resm, resl 66 67 Preconditions: abs(ah) > abs(am) > abs(al) 68 ah and am are overlapping not more than 51 bits 69 am and al are overlapping not more than 51 bits 70 71 Guarantees: abs(resh) > abs(resm) > abs(resl) 72 resh and resm are non-overlapping 73 resm and resl are non-overlapping 74 resm = round-to-nearest(resm + resl) 75 76 Details: resh, resm and resl are considered to be pointers 77 78 */ 79 #define DoRenormalize3(resh, resm, resl, ah, am, al) \ 80 { \ 81 double _t1h, _t1l, _t2l; \ 82 \ 83 Add12(_t1h, _t1l, (am), (al)); \ 84 Add12((*(resh)), _t2l, (ah), (_t1h)); \ 85 Add12((*(resm)), (*(resl)), _t2l, _t1l); \ 86 } 87 88 89 /* Mul23 90 91 Procedure for multiplying two double double numbers resulting 92 in a triple double number 93 94 95 Arguments: two double double numbers: 96 ah, al and 97 bh, bl 98 99 Results: a triple double number resh, resm, resl 100 101 Preconditions: abs(ah) > abs(al) 102 ah and al do not overlap 103 ah = round-to-nearest(ah + al) 104 abs(bh) > abs(bl) 105 bh and bl do not overlap 106 bh = round-to-nearest(bh + bl) 107 108 Guarantees: resm and resl are non-overlapping 109 resm = round-to-nearest(resm + resl) 110 abs(resm) <= 2^(-49) * abs(resh) 111 resh+resm+resl = (ah+al) * (bh+bl) * (1 + eps) 112 where 113 abs(eps) <= 2^(-149) 114 115 Details: resh, resm and resl are considered to be pointers 116 */ 117 #define DoMul23(resh, resm, resl, ah, al, bh, bl) \ 118 { \ 119 double _t1, _t2, _t3, _t4, _t5, _t6, _t7, _t8, _t9, _t10; \ 120 \ 121 Mul12((resh),&_t1,(ah),(bh)); \ 122 Mul12(&_t2,&_t3,(ah),(bl)); \ 123 Mul12(&_t4,&_t5,(al),(bh)); \ 124 _t6 = (al) * (bl); \ 125 Add22Cond(&_t7,&_t8,_t2,_t3,_t4,_t5); \ 126 Add12(_t9,_t10,_t1,_t6); \ 127 Add22Cond((resm),(resl),_t7,_t8,_t9,_t10); \ 128 } 129 130 131 132 /* Mul233 133 134 Procedure for multiplying a double double number by 135 a triple double number resulting in a triple double number 136 137 138 Arguments: a double double number ah, al 139 a triple double number bh, bm, bl 140 141 Results: a triple double number resh, resm, resl 142 143 Preconditions: abs(ah) > abs(al) 144 ah and al do not overlap 145 ah = round-to-nearest(ah + al) 146 abs(bm) <= 2^(-b_o) * abs(bh) 147 abs(bl) <= 2^(-b_u) * abs(bm) 148 where 149 b_o >= 2 150 b_u >= 1 151 152 Guarantees: resm and resl are non-overlapping 153 resm = round-to-nearest(resm + resl) 154 abs(resm) <= 2^(\gamma) * abs(resh) 155 where 156 \gamma >= min(48,b_o-4,b_o+b_u-4) 157 resh+resm+resl=(ah+al) * (bh+bm+bl) * (1+eps) 158 where 159 abs(eps) <= 160 (2^(-99-b_o) + 2^(-99-b_o-b_u) + 2^(-152)) / 161 (1 - 2^(-53) - 2^(-b_o+1) - 2^(-b_o-b_u+1)) 162 163 Details: resh, resm and resl are considered to be pointers 164 */ 165 #define DoMul233(resh, resm, resl, ah, al, bh, bm, bl) \ 166 { \ 167 double _t1, _t2, _t3, _t4, _t5, _t6, _t7, _t8, _t9, _t10; \ 168 double _t11, _t12, _t13, _t14, _t15, _t16, _t17, _t18; \ 169 \ 170 Mul12((resh),&_t1,(ah),(bh)); \ 171 Mul12(&_t2,&_t3,(ah),(bm)); \ 172 Mul12(&_t4,&_t5,(ah),(bl)); \ 173 Mul12(&_t6,&_t7,(al),(bh)); \ 174 Mul12(&_t8,&_t9,(al),(bm)); \ 175 _t10 = (al) * (bl); \ 176 Add22Cond(&_t11,&_t12,_t2,_t3,_t4,_t5); \ 177 Add22Cond(&_t13,&_t14,_t6,_t7,_t8,_t9); \ 178 Add22Cond(&_t15,&_t16,_t11,_t12,_t13,_t14); \ 179 Add12Cond(_t17,_t18,_t1,_t10); \ 180 Add22Cond((resm),(resl),_t17,_t18,_t15,_t16); \ 181 } 182 183 184 185 186 /* Add33 187 188 Procedure for adding two triple double numbers resulting 189 in a triple double number 190 191 192 Arguments: two triple double numbers: 193 ah, am, al and 194 bh, bm, bl 195 196 Results: a triple double number resh, resm, resl 197 198 Preconditions: abs(bh) <= 0.75 * abs(ah) OR ( sign(bh) = sign(ah) AND abs(bh) <= abs(ah)) (i) 199 abs(am) <= 2^(-a_o) * abs(ah) 200 abs(al) <= 2^(-a_u) * abs(am) 201 abs(bm) <= 2^(-b_o) * abs(bh) 202 abs(bl) <= 2^(-b_u) * abs(bm) 203 where 204 b_o >= a_o >= 4 205 b_u >= a_u >= 4 206 207 Condition (i) may not be respected if 208 one can assume in this case that ah=am=al=0 209 210 Guarantees: resm and resl are non-overlapping 211 resm = round-to-nearest(resm + resl) 212 abs(resm) <= 2^(-min(a_o,b_o) + 5) * abs(resh) 213 resh+resm+resl = (ah+am+al + bh+bm+bl) * (1+eps) 214 where 215 abs(eps) <= 2^(-min(a_o+a_u,b_o+b_u)-47) + 2^(-min(a_o,a_u)-98) 216 217 Details: resh, resm and resl are considered to be pointers 218 */ 219 220 #define Add33(resh, resm, resl, ah, am, al, bh, bm, bl) \ 221 { \ 222 double _t1, _t2, _t3, _t4, _t5, _t6, _t7, _t8; \ 223 \ 224 Add12((*(resh)),_t1,(ah),(bh)); \ 225 Add12Cond(_t2,_t3,(am),(bm)); \ 226 _t6 = (al) + (bl); \ 227 Add12Cond(_t7,_t4,_t1,_t2); \ 228 _t5 = _t3 + _t4; \ 229 _t8 = _t5 + _t6; \ 230 Add12Cond((*(resm)),(*(resl)),_t7,_t8); \ 231 } 232 233 /* Add33Cond 234 235 Procedure for adding two triple double numbers resulting 236 in a triple double number 237 238 239 Arguments: two triple double numbers: 240 ah, am, al and 241 bh, bm, bl 242 243 Results: a triple double number resh, resm, resl 244 245 Preconditions: abs(am) <= 2^(-a_o) * abs(ah) 246 abs(al) <= 2^(-a_u) * abs(am) 247 abs(bm) <= 2^(-b_o) * abs(bh) 248 abs(bl) <= 2^(-b_u) * abs(bm) 249 where 250 b_o >= a_o >= 4 251 b_u >= a_u >= 4 252 253 Condition (i) may not be respected if 254 one can assume in this case that ah=am=al=0 255 256 Guarantees: TODO 257 Details: resh, resm and resl are considered to be pointers 258 */ 259 260 #define Add33Cond(resh, resm, resl, ah, am, al, bh, bm, bl) \ 261 { \ 262 double _t1, _t2, _t3, _t4, _t5, _t6, _t7, _t8; \ 263 \ 264 Add12Cond((*(resh)),_t1,(ah),(bh)); \ 265 Add12Cond(_t2,_t3,(am),(bm)); \ 266 _t6 = (al) + (bl); \ 267 Add12Cond(_t7,_t4,_t1,_t2); \ 268 _t5 = _t3 + _t4; \ 269 _t8 = _t5 + _t6; \ 270 Add12Cond((*(resm)),(*(resl)),_t7,_t8); \ 271 } 272 273 274 275 /* Add233 276 277 Procedure for adding a double double number to a triple 278 double number resulting in a triple double number 279 280 281 Arguments: a double double number ah, al 282 a triple double number bh, bm, bl 283 284 Results: a triple double number resh, resm, resl 285 286 Preconditions: abs(ah) > abs(al) 287 ah and al do not overlap 288 ah = round-to-nearest(ah + al) 289 abs(bh) <= 2^(-2) * abs(ah) 290 abs(bm) <= 2^(-b_o) * abs(bh) 291 abs(bl) <= 2^(-b_u) * abs(bm) 292 where 293 b_o >= 2 294 b_u >= 1 295 296 Guarantees: resm and resl are non-overlapping 297 resm = round-to-nearest(resm + resl) 298 abs(resm) <= 2^(\gamma) * abs(resh) 299 where 300 \gamma >= min(45,b_o-4,b_o+b_u-2) 301 resh+resm+resl=((ah+al) + (bh+bm+bl)) * (1+eps) 302 where 303 abs(eps) <= 304 <= 2^(-b_o-b_u-52) + 2^(-b_o-104) + 2^(-153) 305 306 Details: resh, resm and resl are considered to be pointers 307 */ 308 #define Add233(resh, resm, resl, ah, al, bh, bm, bl) \ 309 { \ 310 double _t1, _t2, _t3, _t4, _t5, _t6, _t7; \ 311 \ 312 Add12((*(resh)),_t1,(ah),(bh)); \ 313 Add12Cond(_t2,_t3,(al),(bm)); \ 314 Add12Cond(_t4,_t5,_t1,_t2); \ 315 _t6 = _t3 + (bl); \ 316 _t7 = _t6 + _t5; \ 317 Add12Cond((*(resm)),(*(resl)),_t4,_t7); \ 318 } 319 320 /* Add123 321 322 Procedure for adding a double number to a double 323 double number resulting in a triple double number 324 325 326 Arguments: a double number a 327 a double double number bh, bl 328 329 Results: a triple double number resh, resm, resl 330 331 Preconditions: abs(bh) <= 2^(-2) * abs(a) 332 abs(bl) <= 2^(-53) * abs(bh) 333 334 Guarantees: resm and resl are non-overlapping 335 resm = round-to-nearest(resm + resl) 336 abs(resm) <= 2^(-\gamma) * abs(resh) 337 where 338 339 \gamma >= 52 340 341 resh+resm+resl=(a + (bh+bm+bl)) exactly 342 343 344 Details: resh, resm and resl are considered to be pointers 345 */ 346 #define Add123(resh, resm, resl, a, bh, bl) \ 347 { \ 348 double _t1; \ 349 \ 350 Add12((*(resh)),_t1,(a),(bh)); \ 351 Add12((*(resm)),(*(resl)),_t1,(bl)); \ 352 } 353 354 /* Add213 355 356 Procedure for adding a double double number to a double 357 number resulting in a triple double number 358 359 360 Arguments: a double double number ah, al 361 a double number b 362 363 Results: a triple double number resh, resm, resl 364 365 Preconditions: abs(b) <= 2^(-2) * abs(ah) 366 abs(al) <= 2^(-53) * abs(ah) 367 368 Guarantees: resm and resl are non-overlapping 369 resm = round-to-nearest(resm + resl) 370 abs(resm) <= 2^(-\gamma) * abs(resh) 371 where 372 373 \gamma >= 52 374 375 resh+resm+resl=(a + (bh+bm+bl)) exactly 376 377 378 Details: resh, resm and resl are considered to be pointers 379 */ 380 #define Add213(resh, resm, resl, ah, al, b) \ 381 { \ 382 double _t1; \ 383 \ 384 Add12((*(resh)),_t1,(ah),(b)); \ 385 Add12Cond((*(resm)),(*(resl)),(al),(b)); \ 386 } 387 388 389 390 /* Add23 391 392 Procedure for adding a double-double number to a double-double 393 number resulting in a triple double number 394 395 396 Arguments: a double double number ah, al 397 a double double number bh, bl 398 399 Results: a triple double number resh, resm, resl 400 401 Preconditions: abs(bh) <= 2^(-2) * abs(ah) 402 abs(al) <= 2^(-53) * abs(ah) 403 abs(bl) <= 2^(-53) * abs(bh) 404 405 Guarantees: TO DO 406 407 408 Details: resh, resm and resl are considered to be pointers 409 */ 410 #define Add23(resh, resm, resl, ah, al, bh, bl) \ 411 { \ 412 double _t1, _t2, _t3, _t4, _t5, _t6; \ 413 \ 414 Add12((*(resh)),_t1,(ah),(bh)); \ 415 Add12Cond(_t2,_t3,(al),(bl)); \ 416 Add12Cond(_t4,_t5,_t1,_t2); \ 417 _t6 = _t3 + _t5; \ 418 Add12Cond((*(resm)),(*(resl)),_t4,_t6); \ 419 } 420 421 422 423 424 /* Add133 425 426 Procedure for adding a double number to a triple 427 double number resulting in a triple double number 428 429 430 Arguments: a double number a 431 a triple double number bh, bm, bl 432 433 Results: a triple double number resh, resm, resl 434 435 Preconditions: abs(bh) <= 2^(-2) * abs(a) 436 abs(bm) <= 2^(-b_o) * abs(bh) 437 abs(bl) <= 2^(-b_u) * abs(bm) 438 where 439 b_o >= 2 440 b_u >= 1 441 442 Guarantees: resm and resl are non-overlapping 443 resm = round-to-nearest(resm + resl) 444 abs(resm) <= 2^(\gamma) * abs(resh) 445 where 446 \gamma >= min(47,2-b_o,1-b_o-b_u) 447 resh+resm+resl=(a + (bh+bm+bl)) * (1+eps) 448 where 449 abs(eps) <= 450 <= 2^(-52-b_o-b_u) + 2^(-154) 451 452 453 Details: resh, resm and resl are considered to be pointers 454 */ 455 #define Add133(resh, resm, resl, a, bh, bm, bl) \ 456 { \ 457 double _t1, _t2, _t3, _t4; \ 458 \ 459 Add12((*(resh)),_t1,(a),(bh)); \ 460 Add12Cond(_t2,_t3,_t1,(bm)); \ 461 _t4 = _t3 + (bl); \ 462 Add12Cond((*(resm)),(*(resl)),_t2,_t4); \ 463 } 464 465 /* Add133Cond 466 467 Procedure for adding a double number to a triple 468 double number resulting in a triple double number 469 470 471 Arguments: a double number a 472 a triple double number bh, bm, bl 473 474 Results: a triple double number resh, resm, resl 475 476 Preconditions: abs(bm) <= 2^(-b_o) * abs(bh) 477 abs(bl) <= 2^(-b_u) * abs(bm) 478 where 479 b_o >= 2 480 b_u >= 1 481 482 Guarantees: resm and resl are non-overlapping 483 resm = round-to-nearest(resm + resl) 484 abs(resm) <= 2^(\gamma) * abs(resh) 485 where 486 487 TODO 488 489 resh+resm+resl=(a + (bh+bm+bl)) * (1+eps) 490 where 491 abs(eps) <= 492 493 TODO 494 495 496 Details: resh, resm and resl are considered to be pointers 497 */ 498 #define Add133Cond(resh, resm, resl, a, bh, bm, bl) \ 499 { \ 500 double _t1, _t2, _t3, _t4; \ 501 \ 502 Add12Cond((*(resh)),_t1,(a),(bh)); \ 503 Add12Cond(_t2,_t3,_t1,(bm)); \ 504 _t4 = _t3 + (bl); \ 505 Add12Cond((*(resm)),(*(resl)),_t2,_t4); \ 506 } 507 508 509 510 /* Add233Cond 511 512 Procedure for adding a double double number to a triple 513 double number resulting in a triple double number 514 515 516 Arguments: a double double number ah, al 517 a triple double number bh, bm, bl 518 519 Results: a triple double number resh, resm, resl 520 521 Preconditions: abs(ah) > abs(al) 522 ah and al do not overlap 523 ah = round-to-nearest(ah + al) 524 abs(bm) <= 2^(-b_o) * abs(bh) 525 abs(bl) <= 2^(-b_u) * abs(bm) 526 where 527 b_o >= 2 528 b_u >= 1 529 530 Guarantees: resm and resl are non-overlapping 531 resm = round-to-nearest(resm + resl) 532 abs(resm) <= 2^(\gamma) * abs(resh) 533 where 534 \gamma >= ???? 535 resh+resm+resl=((ah+al) + (bh+bm+bl)) * (1+eps) 536 where 537 abs(eps) <= 538 <= ???? 539 540 Details: resh, resm and resl are considered to be pointers 541 */ 542 #define Add233Cond(resh, resm, resl, ah, al, bh, bm, bl) \ 543 { \ 544 double _t1, _t2, _t3, _t4, _t5, _t6, _t7; \ 545 \ 546 Add12Cond((*(resh)),_t1,(ah),(bh)); \ 547 Add12Cond(_t2,_t3,(al),(bm)); \ 548 Add12Cond(_t4,_t5,_t1,_t2); \ 549 _t6 = _t3 + (bl); \ 550 _t7 = _t6 + _t5; \ 551 Add12Cond((*(resm)),(*(resl)),_t4,_t7); \ 552 } 553 554 555 556 557 /* Mul33 558 559 Procedure for multiplying two triple double numbers resulting 560 in a triple double number 561 562 563 Arguments: two triple double numbers: 564 ah, am, al and 565 bh, bm, bl 566 567 Results: a triple double number resh, resm, resl 568 569 Preconditions: abs(am) <= 2^(-a_o) * abs(ah) 570 abs(al) <= 2^(-a_u) * abs(am) 571 abs(bm) <= 2^(-b_o) * abs(bh) 572 abs(bl) <= 2^(-b_u) * abs(bm) 573 where 574 b_o, a_o >= 5 575 b_u, a_u >= 5 576 577 578 Guarantees: resm and resl are non-overlapping 579 resm = round-to-nearest(resm + resl) 580 abs(resm) <= 2^(-g_o) * abs(resh) 581 with 582 g_o > min(48,-4+a_o,-4+b_o,-4+a_o-b_o) 583 resh+resm+resl = (ah+am+al) * (bh+bm+bl) * (1+eps) 584 where 585 abs(eps) <= 2^-151 + 2^-99-a_o + 2^-99-b_o + 586 + 2^-49-a_o-a_u + 2^-49-b_o-b_u + 2^50-a_o-b_o-b_u + 587 + 2^50-a_o-b_o-b_u + 2^-101-a_o-b_o + 2^-52-a_o-a_u-b_o-b_u 588 589 Details: resh, resm and resl are considered to be pointers 590 */ 591 592 #define DoMul33(resh, resm, resl, ah, am, al, bh, bm, bl) \ 593 { \ 594 double _t1, _t2, _t3, _t4, _t5, _t6, _t7, _t8, _t9; \ 595 double _t10, _t11, _t12, _t13, _t14, _t15, _t16, _t17; \ 596 double _t18, _t19, _t20, _t21, _t22; \ 597 \ 598 Mul12((resh),&_t1,(ah),(bh)); \ 599 Mul12(&_t2,&_t3,(ah),(bm)); \ 600 Mul12(&_t4,&_t5,(am),(bh)); \ 601 Mul12(&_t6,&_t7,(am),(bm)); \ 602 _t8 = (ah) * (bl); \ 603 _t9 = (al) * (bh); \ 604 _t10 = (am) * (bl); \ 605 _t11 = (al) * (bm); \ 606 _t12 = _t8 + _t9; \ 607 _t13 = _t10 + _t11; \ 608 Add12Cond(_t14,_t15,_t1,_t6); \ 609 _t16 = _t7 + _t15; \ 610 _t17 = _t12 + _t13; \ 611 _t18 = _t16 + _t17; \ 612 Add12Cond(_t19,_t20,_t14,_t18); \ 613 Add22Cond(&_t21,&_t22,_t2,_t3,_t4,_t5); \ 614 Add22Cond((resm),(resl),_t21,_t22,_t19,_t20); \ 615 } 616 617 618 /* Mul133 619 620 Procedure for multiplying double by a triple double number resulting 621 in a triple double number 622 623 624 Arguments: a double a 625 a triple double bh, bm, bl 626 627 Results: a triple double number resh, resm, resl 628 629 Preconditions: abs(bm) <= 2^(-b_o) * abs(bh) 630 abs(bl) <= 2^(-b_u) * abs(bm) 631 where 632 b_o >= 2 633 b_u >= 2 634 635 636 Guarantees: resm and resl are non-overlapping 637 resm = round-to-nearest(resm + resl) 638 abs(resm) <= 2^(-g_o) * abs(resh) 639 with 640 g_o > min(47,-5-b_o,-5+b_o+b_u) 641 resh+resm+resl = a * (bh+bm+bl) * (1+eps) 642 where 643 abs(eps) <= 2^-49-b_o-b_u + 2^-101-b_o + 2^-156 644 645 Details: resh, resm and resl are considered to be pointers 646 */ 647 #define DoMul133(resh, resm, resl, a, bh, bm, bl) \ 648 { \ 649 double _t2, _t3, _t4, _t5, _t7, _t8, _t9, _t10; \ 650 \ 651 Mul12((resh),&_t2,(a),(bh)); \ 652 Mul12(&_t3,&_t4,(a),(bm)); \ 653 _t5 = (a) * (bl); \ 654 Add12Cond(_t9,_t7,_t2,_t3); \ 655 _t8 = _t4 + _t5; \ 656 _t10 = _t7 + _t8; \ 657 Add12Cond((*(resm)),(*(resl)),_t9,_t10); \ 658 } 659 660 /* Mul123 661 662 Procedure for multiplying double by a double double number resulting 663 in a triple double number 664 665 666 Arguments: a double a 667 a double double bh, bl 668 669 Results: a triple double number resh, resm, resl 670 671 Guarantees: resm and resl are non-overlapping 672 resm = round-to-nearest(resm + resl) 673 abs(resm) <= 2^(-g_o) * abs(resh) 674 with 675 g_o > 47 676 resh+resm+resl = a * (bh+bm) * (1+eps) 677 where 678 abs(eps) <= 2^-154 679 680 Details: resh, resm and resl are considered to be pointers 681 */ 682 #define DoMul123(resh, resm, resl, a, bh, bl) \ 683 { \ 684 double _t1, _t2, _t3, _t4, _t5, _t6; \ 685 \ 686 Mul12((resh),&_t1,(a),(bh)); \ 687 Mul12(&_t2,&_t3,(a),(bl)); \ 688 Add12Cond(_t5,_t4,_t1,_t2); \ 689 _t6 = _t3 + _t4; \ 690 Add12Cond((*(resm)),(*(resl)),_t5,_t6); \ 691 } 692 693 694 695 /* ReturnRoundToNearest3 696 697 Procedure for rounding a triple to a double number 698 in round-to-nearest-ties-to-even mode. 699 700 701 Arguments: a triple double number xh, xm, xl 702 703 Results: a double number xprime 704 returned by a return-statement 705 706 Preconditions: xh, xm and xl are non-overlapping 707 xm = RN(xm +math xl) 708 xh != 0, xm != 0 709 xl = 0 iff xm != +/- 0.5 * ulp(xh) (0.25 if xh = 2^e) 710 711 Guarantees: xprime = RN(xh + xm + xl) 712 713 Sideeffects: returns, i.e. leaves the function 714 715 */ 716 #define ReturnRoundToNearest3(xh,xm,xl) \ 717 { \ 718 double _t1, _t2, _t3, _t4, _t5, _t6; \ 719 db_number _xp, _xn; \ 720 \ 721 _xp.d = (xh); \ 722 _xn.i[HI] = _xp.i[HI]; \ 723 _xn.i[LO] = _xp.i[LO]; \ 724 _xn.l--; \ 725 _t1 = _xn.d; \ 726 _xp.l++; \ 727 _t4 = _xp.d; \ 728 _t2 = (xh) - _t1; \ 729 _t3 = _t2 * -0.5; \ 730 _t5 = _t4 - (xh); \ 731 _t6 = _t5 * 0.5; \ 732 if (((xm) != _t3) && ((xm) != _t6)) return ((xh) + (xm)); \ 733 if ((xm) * (xl) > 0.0) { \ 734 if ((xh) * (xl) > 0.0) \ 735 return _t4; \ 736 else \ 737 return _t1; \ 738 } else return (xh); \ 739 } 740 741 /* ReturnRoundToNearest3Other 742 743 ATTENTION: THIS CURRENTLY UNPROVEN CODE !!! 744 745 Procedure for rounding a triple to a double number 746 in round-to-nearest-ties-to-even mode. 747 748 749 Arguments: a triple double number xh, xm, xl 750 751 Results: a double number xprime 752 returned by a return-statement 753 754 Preconditions: |xm + xl| <= 2^(-5) * |xh| 755 756 Guarantees: xprime = RN(xh + xm + xl) 757 758 Sideeffects: returns, i.e. leaves the function 759 760 */ 761 #define ReturnRoundToNearest3Other(xh,xm,xl) \ 762 { \ 763 double _t3, _t4; \ 764 db_number _t3db; \ 765 \ 766 Add12(_t3,_t4,(xm),(xl)); \ 767 if (_t4 != 0.0) { \ 768 _t3db.d = _t3; \ 769 if (!(_t3db.i[LO] & 0x00000001)) { \ 770 if ((_t4 > 0.0) ^ ((_t3db.i[HI] & 0x80000000) != 0)) \ 771 _t3db.l++; \ 772 else \ 773 _t3db.l--; \ 774 _t3 = _t3db.d; \ 775 } \ 776 } \ 777 return (xh) + _t3; \ 778 } 779 780 781 782 /* ReturnRoundUpwards3 783 784 Procedure for rounding a triple to a double number 785 in round-upwards mode. 786 787 788 Arguments: a triple double number xh, xm, xl 789 790 Results: a double number xprime 791 returned by a return-statement 792 793 Preconditions: xh, xm and xl are non-overlapping 794 xm = RN(xm +math xl) 795 xh != 0, xm != 0 796 797 Exact algebraic images have already 798 been filtered out. 799 800 Guarantees: xprime = RU(xh + xm + xl) 801 802 Sideeffects: returns, i.e. leaves the function 803 804 */ 805 #define ReturnRoundUpwards3(xh,xm,xl) \ 806 { \ 807 double _t1, _t2, _t3; \ 808 db_number _tdb; \ 809 \ 810 Add12(_t1,_t2,(xh),(xm)); \ 811 _t3 = _t2 + (xl); \ 812 if (_t3 > 0.0) { \ 813 if (_t1 > 0.0) { \ 814 _tdb.d = _t1; \ 815 _tdb.l++; \ 816 return _tdb.d; \ 817 } else { \ 818 _tdb.d = _t1; \ 819 _tdb.l--; \ 820 return _tdb.d; \ 821 } \ 822 } else return _t1; \ 823 } 824 825 826 /* ReturnRoundDownwards3 827 828 Procedure for rounding a triple to a double number 829 in round-downwards mode. 830 831 832 Arguments: a triple double number xh, xm, xl 833 834 Results: a double number xprime 835 returned by a return-statement 836 837 Preconditions: xh, xm and xl are non-overlapping 838 xm = RN(xm +math xl) 839 xh != 0, xm != 0 840 841 Exact algebraic images have already 842 been filtered out. 843 844 Guarantees: xprime = RD(xh + xm + xl) 845 846 Sideeffects: returns, i.e. leaves the function 847 848 */ 849 #define ReturnRoundDownwards3(xh,xm,xl) \ 850 { \ 851 double _t1, _t2, _t3; \ 852 db_number _tdb; \ 853 \ 854 Add12(_t1,_t2,(xh),(xm)); \ 855 _t3 = _t2 + (xl); \ 856 if (_t3 < 0.0) { \ 857 if (_t1 > 0.0) { \ 858 _tdb.d = _t1; \ 859 _tdb.l--; \ 860 return _tdb.d; \ 861 } else { \ 862 _tdb.d = _t1; \ 863 _tdb.l++; \ 864 return _tdb.d; \ 865 } \ 866 } else return _t1; \ 867 } 868 869 870 /* ReturnRoundTowardsZero3 871 872 Procedure for rounding a triple to a double number 873 in round-towards-zero mode. 874 875 876 Arguments: a triple double number xh, xm, xl 877 878 Results: a double number xprime 879 returned by a return-statement 880 881 Preconditions: xh, xm and xl are non-overlapping 882 xm = RN(xm +math xl) 883 xh != 0, xm != 0 884 885 Exact algebraic images have already 886 been filtered out. 887 888 Guarantees: xprime = RZ(xh + xm + xl) 889 890 Sideeffects: returns, i.e. leaves the function 891 892 */ 893 #define ReturnRoundTowardsZero3(xh,xm,xl) \ 894 { \ 895 double _t1, _t2, _t3; \ 896 db_number _tdb; \ 897 \ 898 Add12(_t1,_t2,(xh),(xm)); \ 899 _t3 = _t2 + (xl); \ 900 if (_t1 > 0.0) { \ 901 if (_t3 < 0.0) { \ 902 _tdb.d = _t1; \ 903 _tdb.l--; \ 904 return _tdb.d; \ 905 } else return _t1; \ 906 } else { \ 907 if (_t3 > 0.0) { \ 908 _tdb.d = _t1; \ 909 _tdb.l--; \ 910 return _tdb.d; \ 911 } else return _t1; \ 912 } \ 913 } 914 915 916 /* ReturnRoundUpwards3Unfiltered 917 918 Procedure for rounding a triple to a double number 919 in round-upwards mode. 920 921 922 Arguments: a triple double number xh, xm, xl 923 a double constant wca representing 2^k 924 where 2^-k is Lefevre's worst case accuracy 925 926 Results: a double number xprime 927 returned by a return-statement 928 929 Preconditions: xh, xm and xl are non-overlapping 930 xm = RN(xm +math xl) 931 xh != 0, xm != 0 932 933 Guarantees: xprime = RU(xh + xm + xl) 934 935 Sideeffects: returns, i.e. leaves the function 936 937 */ 938 #define ReturnRoundUpwards3Unfiltered(xh,xm,xl,wca) \ 939 { \ 940 double _t1, _t2, _t3; \ 941 db_number _tdb, _tdb2; \ 942 \ 943 Add12(_t1,_t2,(xh),(xm)); \ 944 _t3 = _t2 + (xl); \ 945 if (_t3 > 0.0) { \ 946 _tdb2.d = wca * _t3; \ 947 _tdb.d = _t1; \ 948 if ((_tdb2.i[HI] & 0x7ff00000) < (_tdb.i[HI] & 0x7ff00000)) \ 949 return _t1; \ 950 if (_t1 > 0.0) { \ 951 _tdb.l++; \ 952 return _tdb.d; \ 953 } else { \ 954 _tdb.l--; \ 955 return _tdb.d; \ 956 } \ 957 } else return _t1; \ 958 } 959 960 961 962 /* ReturnRoundDownwards3Unfiltered 963 964 Procedure for rounding a triple to a double number 965 in round-downwards mode. 966 967 968 Arguments: a triple double number xh, xm, xl 969 a double constant wca representing 2^k 970 where 2^-k is Lefevre's worst case accuracy 971 972 Results: a double number xprime 973 returned by a return-statement 974 975 Preconditions: xh, xm and xl are non-overlapping 976 xm = RN(xm +math xl) 977 xh != 0, xm != 0 978 979 Guarantees: xprime = RD(xh + xm + xl) 980 981 Sideeffects: returns, i.e. leaves the function 982 983 */ 984 #define ReturnRoundDownwards3Unfiltered(xh,xm,xl,wca) \ 985 { \ 986 double _t1, _t2, _t3; \ 987 db_number _tdb, _tdb2; \ 988 \ 989 Add12(_t1,_t2,(xh),(xm)); \ 990 _t3 = _t2 + (xl); \ 991 if (_t3 < 0.0) { \ 992 _tdb2.d = wca * _t3; \ 993 _tdb.d = _t1; \ 994 if ((_tdb2.i[HI] & 0x7ff00000) < (_tdb.i[HI] & 0x7ff00000)) \ 995 return _t1; \ 996 if (_t1 > 0.0) { \ 997 _tdb.l--; \ 998 return _tdb.d; \ 999 } else { \ 1000 _tdb.l++; \ 1001 return _tdb.d; \ 1002 } \ 1003 } else return _t1; \ 1004 } 1005 1006 /* ReturnRoundTowardsZero3Unfiltered 1007 1008 Procedure for rounding a triple to a double number 1009 in round-towards-zero mode. 1010 1011 1012 Arguments: a triple double number xh, xm, xl 1013 a double constant wca representing 2^k 1014 where 2^-k is Lefevre's worst case accuracy 1015 1016 Results: a double number xprime 1017 returned by a return-statement 1018 1019 Preconditions: xh, xm and xl are non-overlapping 1020 xm = RN(xm +math xl) 1021 xh != 0, xm != 0 1022 1023 Guarantees: xprime = RZ(xh + xm + xl) 1024 1025 Sideeffects: returns, i.e. leaves the function 1026 1027 */ 1028 #define ReturnRoundTowardsZero3Unfiltered(xh,xm,xl,wca) \ 1029 { \ 1030 if ((xh) > 0) \ 1031 ReturnRoundDownwards3Unfiltered((xh),(xm),(xl),(wca)) \ 1032 else \ 1033 ReturnRoundUpwards3Unfiltered((xh),(xm),(xl),(wca)) \ 1034 } 1035 1036 /* RoundToNearest3 1037 1038 Procedure for rounding a triple to a double number 1039 in round-to-nearest-ties-to-even mode. 1040 1041 1042 Arguments: a triple double number xh, xm, xl 1043 1044 Results: a double number xprime 1045 returned by a return-statement 1046 1047 Preconditions: xh, xm and xl are non-overlapping 1048 xm = RN(xm +math xl) 1049 xh != 0, xm != 0 1050 xl = 0 iff xm != +/- 0.5 * ulp(xh) (0.25 if xh = 2^e) 1051 1052 Guarantees: xprime = RN(xh + xm + xl) 1053 1054 Details: res is considered to be a pointer 1055 1056 */ 1057 #define RoundToNearest3(res,xh,xm,xl) \ 1058 { \ 1059 double _t1, _t2, _t3, _t4, _t5, _t6; \ 1060 db_number _xp, _xn; \ 1061 \ 1062 _xp.d = (xh); \ 1063 _xn.i[HI] = _xp.i[HI]; \ 1064 _xn.i[LO] = _xp.i[LO]; \ 1065 _xn.l--; \ 1066 _t1 = _xn.d; \ 1067 _xp.l++; \ 1068 _t4 = _xp.d; \ 1069 _t2 = (xh) - _t1; \ 1070 _t3 = _t2 * -0.5; \ 1071 _t5 = _t4 - (xh); \ 1072 _t6 = _t5 * 0.5; \ 1073 if (((xm) != _t3) && ((xm) != _t6)) \ 1074 (*(res)) = ((xh) + (xm)); \ 1075 else { \ 1076 if ((xm) * (xl) > 0.0) { \ 1077 if ((xh) * (xl) > 0.0) \ 1078 (*(res)) = _t4; \ 1079 else \ 1080 (*(res)) = _t1; \ 1081 } else (*(res)) = (xh); \ 1082 } \ 1083 } 1084 1085 /* RoundUpwards3 1086 1087 Procedure for rounding a triple to a double number 1088 in round-upwards mode. 1089 1090 1091 Arguments: a triple double number xh, xm, xl 1092 1093 Results: a double number xprime 1094 returned by a return-statement 1095 1096 Preconditions: xh, xm and xl are non-overlapping 1097 xm = RN(xm +math xl) 1098 xh != 0, xm != 0 1099 1100 Exact algebraic images have already 1101 been filtered out. 1102 1103 Guarantees: xprime = RU(xh + xm + xl) 1104 1105 Details: res is considered to be a pointer 1106 1107 */ 1108 #define RoundUpwards3(res,xh,xm,xl) \ 1109 { \ 1110 double _t1, _t2, _t3; \ 1111 db_number _tdb; \ 1112 \ 1113 Add12(_t1,_t2,(xh),(xm)); \ 1114 _t3 = _t2 + (xl); \ 1115 if (_t3 > 0.0) { \ 1116 if (_t1 > 0.0) { \ 1117 _tdb.d = _t1; \ 1118 _tdb.l++; \ 1119 (*(res)) = _tdb.d; \ 1120 } else { \ 1121 _tdb.d = _t1; \ 1122 _tdb.l--; \ 1123 (*(res)) = _tdb.d; \ 1124 } \ 1125 } else (*(res)) = _t1; \ 1126 } 1127 1128 1129 /* RoundDownwards3 1130 1131 Procedure for rounding a triple to a double number 1132 in round-downwards mode. 1133 1134 1135 Arguments: a triple double number xh, xm, xl 1136 1137 Results: a double number xprime 1138 returned by a return-statement 1139 1140 Preconditions: xh, xm and xl are non-overlapping 1141 xm = RN(xm +math xl) 1142 xh != 0, xm != 0 1143 1144 Exact algebraic images have already 1145 been filtered out. 1146 1147 Guarantees: xprime = RD(xh + xm + xl) 1148 1149 Details: res is considered to be a pointer 1150 1151 */ 1152 #define RoundDownwards3(res,xh,xm,xl) \ 1153 { \ 1154 double _t1, _t2, _t3; \ 1155 db_number _tdb; \ 1156 \ 1157 Add12(_t1,_t2,(xh),(xm)); \ 1158 _t3 = _t2 + (xl); \ 1159 if (_t3 < 0.0) { \ 1160 if (_t1 > 0.0) { \ 1161 _tdb.d = _t1; \ 1162 _tdb.l--; \ 1163 (*(res)) = _tdb.d; \ 1164 } else { \ 1165 _tdb.d = _t1; \ 1166 _tdb.l++; \ 1167 (*(res)) = _tdb.d; \ 1168 } \ 1169 } else (*(res)) = _t1; \ 1170 } 1171 1172 1173 /* RoundTowardsZero3 1174 1175 Procedure for rounding a triple to a double number 1176 in round-towards-zero mode. 1177 1178 1179 Arguments: a triple double number xh, xm, xl 1180 1181 Results: a double number xprime 1182 returned by a return-statement 1183 1184 Preconditions: xh, xm and xl are non-overlapping 1185 xm = RN(xm +math xl) 1186 xh != 0, xm != 0 1187 1188 Exact algebraic images have already 1189 been filtered out. 1190 1191 Guarantees: xprime = RZ(xh + xm + xl) 1192 1193 Details: res is considered to be a pointer 1194 1195 */ 1196 #define RoundTowardsZero3(res,xh,xm,xl) \ 1197 { \ 1198 double _t1, _t2, _t3; \ 1199 db_number _tdb; \ 1200 \ 1201 Add12(_t1,_t2,(xh),(xm)); \ 1202 _t3 = _t2 + (xl); \ 1203 if (_t1 > 0.0) { \ 1204 if (_t3 < 0.0) { \ 1205 _tdb.d = _t1; \ 1206 _tdb.l--; \ 1207 (*(res)) = _tdb.d; \ 1208 } else (*(res)) = _t1; \ 1209 } else { \ 1210 if (_t3 > 0.0) { \ 1211 _tdb.d = _t1; \ 1212 _tdb.l--; \ 1213 (*(res)) = _tdb.d; \ 1214 } else (*(res)) = _t1; \ 1215 } \ 1216 } 1217 1218 /* sqrt13 1219 1220 Computes a triple-double approximation of sqrt(x) 1221 1222 Should be provable to be exact to at least 140 bits. 1223 1224 Only handles the following special cases: 1225 - x == 0 1226 - subnormal x 1227 The following cases are not handled: 1228 - x < 0 1229 - x = +/-Infty, NaN 1230 1231 */ 1232 1233 1234 #define DoSqrt13(resh, resm, resl , x) \ 1235 { \ 1236 db_number _xdb; \ 1237 int _E; \ 1238 double _m, _r0, _r1, _r2, _r3h, _r3l, _r4h, _r4l; \ 1239 double _r5h, _r5m, _r5l, _srtmh, _srtml, _srtmm; \ 1240 double _r2PHr2h, _r2PHr2l, _r2Sqh, _r2Sql; \ 1241 double _mMr2h, _mMr2l, _mMr2Ch, _mMr2Cl; \ 1242 double _MHmMr2Ch, _MHmMr2Cl; \ 1243 double _r3Sqh, _r3Sql, _mMr3Sqh, _mMr3Sql; \ 1244 double _srtmhover,_srtmmover,_srtmlover; \ 1245 double _HmMr4Sqm,_HmMr4Sql, _mMr4Sqhover, _mMr4Sqmover, _mMr4Sqlover; \ 1246 double _mMr4Sqh, _mMr4Sqm, _mMr4Sql, _r4Sqh, _r4Sqm, _r4Sql; \ 1247 \ 1248 /* Special case x = 0 */ \ 1249 if ((x) == 0) { \ 1250 (*(resh)) = (x); \ 1251 (*(resm)) = 0; \ 1252 (*(resl)) = 0; \ 1253 } else { \ 1254 \ 1255 _E = 0; \ 1256 \ 1257 /* Convert to integer format */ \ 1258 _xdb.d = (x); \ 1259 \ 1260 /* Handle subnormal case */ \ 1261 if (_xdb.i[HI] < 0x00100000) { \ 1262 _E = -52; \ 1263 _xdb.d *= ((db_number) ((double) SQRTTWO52)).d; \ 1264 /* make x a normal number */ \ 1265 } \ 1266 \ 1267 /* Extract exponent E and mantissa m */ \ 1268 _E += (_xdb.i[HI]>>20)-1023; \ 1269 _xdb.i[HI] = (_xdb.i[HI] & 0x000fffff) | 0x3ff00000; \ 1270 _m = _xdb.d; \ 1271 \ 1272 /* Make exponent even */ \ 1273 if (_E & 0x00000001) { \ 1274 _E++; \ 1275 _m *= 0.5; /* Suppose now 1/2 <= m <= 2 */ \ 1276 } \ 1277 \ 1278 /* Construct sqrt(2^E) = 2^(E/2) */ \ 1279 _xdb.i[HI] = (_E/2 + 1023) << 20; \ 1280 _xdb.i[LO] = 0; \ 1281 \ 1282 /* Compute initial approximation to r = 1/sqrt(m) */ \ 1283 \ 1284 _r0 = SQRTPOLYC0 + \ 1285 _m * (SQRTPOLYC1 + _m * (SQRTPOLYC2 + _m * (SQRTPOLYC3 + _m * SQRTPOLYC4))); \ 1286 \ 1287 /* Iterate two times on double precision */ \ 1288 \ 1289 _r1 = 0.5 * _r0 * (3 - _m * (_r0 * _r0)); \ 1290 _r2 = 0.5 * _r1 * (3 - _m * (_r1 * _r1)); \ 1291 \ 1292 /* Iterate two times on double-double precision */ \ 1293 \ 1294 Mul12(&_r2Sqh, &_r2Sql, _r2, _r2); \ 1295 Add12(_r2PHr2h, _r2PHr2l, _r2, (0.5 * _r2)); \ 1296 Mul12(&_mMr2h, &_mMr2l, _m, _r2); \ 1297 Mul22(&_mMr2Ch, &_mMr2Cl, _mMr2h, _mMr2l, _r2Sqh, _r2Sql); \ 1298 \ 1299 _MHmMr2Ch = -0.5 * _mMr2Ch; \ 1300 _MHmMr2Cl = -0.5 * _mMr2Cl; \ 1301 \ 1302 Add22(&_r3h, &_r3l, _r2PHr2h, _r2PHr2l, _MHmMr2Ch, _MHmMr2Cl); \ 1303 \ 1304 Mul22(&_r3Sqh, &_r3Sql, _r3h, _r3l, _r3h, _r3l); \ 1305 Mul22(&_mMr3Sqh, &_mMr3Sql, _m, 0.0, _r3Sqh, _r3Sql); \ 1306 /* To prove: mMr3Sqh = 1.0 in each case */ \ 1307 \ 1308 Mul22(&_r4h, &_r4l, _r3h, _r3l, 1.0, (-0.5 * _mMr3Sql)); \ 1309 \ 1310 /* Iterate once on triple-double precision */ \ 1311 \ 1312 Mul23(&_r4Sqh, &_r4Sqm, &_r4Sql, _r4h, _r4l, _r4h, _r4l); \ 1313 Mul133(&_mMr4Sqhover, &_mMr4Sqmover, &_mMr4Sqlover, _m, _r4Sqh, _r4Sqm, _r4Sql); \ 1314 Renormalize3(&_mMr4Sqh, &_mMr4Sqm, &_mMr4Sql, _mMr4Sqhover, _mMr4Sqmover, _mMr4Sqlover); \ 1315 /* To prove: mMr4Sqh = 1.0 in each case */ \ 1316 \ 1317 _HmMr4Sqm = -0.5 * _mMr4Sqm; \ 1318 _HmMr4Sql = -0.5 * _mMr4Sql; \ 1319 \ 1320 Mul233(&_r5h,&_r5m,&_r5l,_r4h,_r4l,1.0,_HmMr4Sqm,_HmMr4Sql); \ 1321 \ 1322 /* Multiply obtained reciprocal square root by m */ \ 1323 \ 1324 Mul133(&_srtmhover, &_srtmmover, &_srtmlover,_m,_r5h,_r5m,_r5l); \ 1325 \ 1326 Renormalize3(&_srtmh,&_srtmm,&_srtml,_srtmhover,_srtmmover,_srtmlover); \ 1327 \ 1328 /* Multiply componentwise by sqrt(2^E) */ \ 1329 /* which is an integer power of 2 that may not produce a subnormal */ \ 1330 \ 1331 (*(resh)) = _xdb.d * _srtmh; \ 1332 (*(resm)) = _xdb.d * _srtmm; \ 1333 (*(resl)) = _xdb.d * _srtml; \ 1334 \ 1335 } /* End: special case 0 */ \ 1336 } 1337 1338 1339 /* recpr33() 1340 1341 Computes a triple-double reciprocal of a triple-double 1342 1343 Should be provable to be exact to at least 140 bits 1344 1345 No special case handling is done 1346 1347 dh + dm + dl must be renormalized 1348 1349 The result is renormalized 1350 1351 */ 1352 1353 1354 #define DoRecpr33(resh, resm, resl, dh, dm, dl) \ 1355 { \ 1356 double _rec_r1, _rec_t1, _rec_t2, _rec_t3, _rec_t4, _rec_t5, _rec_t6, _rec_t7, _rec_t8, _rec_t9, _rec_t10, _rec_t11, _rec_t12, _rec_t13, _rec_t14; \ 1357 double _rec_r2h, _rec_r2l, _rec_t15, _rec_t16, _rec_t17, _rec_t18, _rec_t19, _rec_t20, _rec_t21, _rec_t22, _rec_t23; \ 1358 \ 1359 _rec_r1 = 1.0 / (dh); \ 1360 Mul12(&_rec_t1,&_rec_t2,_rec_r1,(dh)); \ 1361 _rec_t3 = _rec_t1 - 1.0; \ 1362 Add12Cond(_rec_t4,_rec_t5,_rec_t3,_rec_t2); \ 1363 Mul12(&_rec_t6,&_rec_t7,_rec_r1,(dm)); \ 1364 Add12(_rec_t8,_rec_t9,-1.0,_rec_t6); \ 1365 _rec_t10 = _rec_t9 + _rec_t7; \ 1366 Add12(_rec_t11,_rec_t12,_rec_t8,_rec_t10); \ 1367 _rec_r1 = -_rec_r1; \ 1368 Add22Cond(&_rec_t13,&_rec_t14,_rec_t4,_rec_t5,_rec_t11,_rec_t12); \ 1369 Mul122(&_rec_r2h,&_rec_r2l,_rec_r1,_rec_t13,_rec_t14); \ 1370 Mul233(&_rec_t15,&_rec_t16,&_rec_t17,_rec_r2h,_rec_r2l,(dh),(dm),(dl)); \ 1371 Renormalize3(&_rec_t18,&_rec_t19,&_rec_t20,_rec_t15,_rec_t16,_rec_t17); \ 1372 _rec_t18 = -1.0; \ 1373 Mul233(&_rec_t21,&_rec_t22,&_rec_t23,_rec_r2h,_rec_r2l,_rec_t18,_rec_t19,_rec_t20); \ 1374 _rec_t21 = -_rec_t21; _rec_t22 = -_rec_t22; _rec_t23 = -_rec_t23; \ 1375 Renormalize3((resh),(resm),(resl),_rec_t21,_rec_t22,_rec_t23); \ 1376 } 1377 1378 1379 1380 #endif /*TRIPLE_rec_DOUBLE_rec_H*/ 1381