1 /* Copyright (C) 2007-2013 Free Software Foundation, Inc. 2 3 This file is part of GCC. 4 5 GCC is free software; you can redistribute it and/or modify it under 6 the terms of the GNU General Public License as published by the Free 7 Software Foundation; either version 3, or (at your option) any later 8 version. 9 10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY 11 WARRANTY; without even the implied warranty of MERCHANTABILITY or 12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License 13 for more details. 14 15 Under Section 7 of GPL version 3, you are granted additional 16 permissions described in the GCC Runtime Library Exception, version 17 3.1, as published by the Free Software Foundation. 18 19 You should have received a copy of the GNU General Public License and 20 a copy of the GCC Runtime Library Exception along with this program; 21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see 22 <http://www.gnu.org/licenses/>. */ 23 24 /***************************************************************************** 25 * 26 * BID64 encoding: 27 * **************************************** 28 * 63 62 53 52 0 29 * |---|------------------|--------------| 30 * | S | Biased Exp (E) | Coeff (c) | 31 * |---|------------------|--------------| 32 * 33 * bias = 398 34 * number = (-1)^s * 10^(E-398) * c 35 * coefficient range - 0 to (2^53)-1 36 * COEFF_MAX = 2^53-1 = 9007199254740991 37 * 38 ***************************************************************************** 39 * 40 * BID128 encoding: 41 * 1-bit sign 42 * 14-bit biased exponent in [0x21, 0x3020] = [33, 12320] 43 * unbiased exponent in [-6176, 6111]; exponent bias = 6176 44 * 113-bit unsigned binary integer coefficient (49-bit high + 64-bit low) 45 * Note: 10^34-1 ~ 2^112.945555... < 2^113 => coefficient fits in 113 bits 46 * 47 * Note: assume invalid encodings are not passed to this function 48 * 49 * Round a number C with q decimal digits, represented as a binary integer 50 * to q - x digits. Six different routines are provided for different values 51 * of q. The maximum value of q used in the library is q = 3 * P - 1 where 52 * P = 16 or P = 34 (so q <= 111 decimal digits). 53 * The partitioning is based on the following, where Kx is the scaled 54 * integer representing the value of 10^(-x) rounded up to a number of bits 55 * sufficient to ensure correct rounding: 56 * 57 * -------------------------------------------------------------------------- 58 * q x max. value of a max number min. number 59 * of bits in C of bits in Kx 60 * -------------------------------------------------------------------------- 61 * 62 * GROUP 1: 64 bits 63 * round64_2_18 () 64 * 65 * 2 [1,1] 10^1 - 1 < 2^3.33 4 4 66 * ... ... ... ... ... 67 * 18 [1,17] 10^18 - 1 < 2^59.80 60 61 68 * 69 * GROUP 2: 128 bits 70 * round128_19_38 () 71 * 72 * 19 [1,18] 10^19 - 1 < 2^63.11 64 65 73 * 20 [1,19] 10^20 - 1 < 2^66.44 67 68 74 * ... ... ... ... ... 75 * 38 [1,37] 10^38 - 1 < 2^126.24 127 128 76 * 77 * GROUP 3: 192 bits 78 * round192_39_57 () 79 * 80 * 39 [1,38] 10^39 - 1 < 2^129.56 130 131 81 * ... ... ... ... ... 82 * 57 [1,56] 10^57 - 1 < 2^189.35 190 191 83 * 84 * GROUP 4: 256 bits 85 * round256_58_76 () 86 * 87 * 58 [1,57] 10^58 - 1 < 2^192.68 193 194 88 * ... ... ... ... ... 89 * 76 [1,75] 10^76 - 1 < 2^252.47 253 254 90 * 91 * GROUP 5: 320 bits 92 * round320_77_96 () 93 * 94 * 77 [1,76] 10^77 - 1 < 2^255.79 256 257 95 * 78 [1,77] 10^78 - 1 < 2^259.12 260 261 96 * ... ... ... ... ... 97 * 96 [1,95] 10^96 - 1 < 2^318.91 319 320 98 * 99 * GROUP 6: 384 bits 100 * round384_97_115 () 101 * 102 * 97 [1,96] 10^97 - 1 < 2^322.23 323 324 103 * ... ... ... ... ... 104 * 115 [1,114] 10^115 - 1 < 2^382.03 383 384 105 * 106 ****************************************************************************/ 107 108 #include "bid_internal.h" 109 110 void 111 round64_2_18 (int q, 112 int x, 113 UINT64 C, 114 UINT64 * ptr_Cstar, 115 int *incr_exp, 116 int *ptr_is_midpoint_lt_even, 117 int *ptr_is_midpoint_gt_even, 118 int *ptr_is_inexact_lt_midpoint, 119 int *ptr_is_inexact_gt_midpoint) { 120 121 UINT128 P128; 122 UINT128 fstar; 123 UINT64 Cstar; 124 UINT64 tmp64; 125 int shift; 126 int ind; 127 128 // Note: 129 // In round128_2_18() positive numbers with 2 <= q <= 18 will be 130 // rounded to nearest only for 1 <= x <= 3: 131 // x = 1 or x = 2 when q = 17 132 // x = 2 or x = 3 when q = 18 133 // However, for generality and possible uses outside the frame of IEEE 754R 134 // this implementation works for 1 <= x <= q - 1 135 136 // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even, 137 // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are 138 // initialized to 0 by the caller 139 140 // round a number C with q decimal digits, 2 <= q <= 18 141 // to q - x digits, 1 <= x <= 17 142 // C = C + 1/2 * 10^x where the result C fits in 64 bits 143 // (because the largest value is 999999999999999999 + 50000000000000000 = 144 // 0x0e92596fd628ffff, which fits in 60 bits) 145 ind = x - 1; // 0 <= ind <= 16 146 C = C + midpoint64[ind]; 147 // kx ~= 10^(-x), kx = Kx64[ind] * 2^(-Ex), 0 <= ind <= 16 148 // P128 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx 149 // the approximation kx of 10^(-x) was rounded up to 64 bits 150 __mul_64x64_to_128MACH (P128, C, Kx64[ind]); 151 // calculate C* = floor (P128) and f* 152 // Cstar = P128 >> Ex 153 // fstar = low Ex bits of P128 154 shift = Ex64m64[ind]; // in [3, 56] 155 Cstar = P128.w[1] >> shift; 156 fstar.w[1] = P128.w[1] & mask64[ind]; 157 fstar.w[0] = P128.w[0]; 158 // the top Ex bits of 10^(-x) are T* = ten2mxtrunc64[ind], e.g. 159 // if x=1, T*=ten2mxtrunc64[0]=0xcccccccccccccccc 160 // if (0 < f* < 10^(-x)) then the result is a midpoint 161 // if floor(C*) is even then C* = floor(C*) - logical right 162 // shift; C* has q - x decimal digits, correct by Prop. 1) 163 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 164 // shift; C* has q - x decimal digits, correct by Pr. 1) 165 // else 166 // C* = floor(C*) (logical right shift; C has q - x decimal digits, 167 // correct by Property 1) 168 // in the caling function n = C* * 10^(e+x) 169 170 // determine inexactness of the rounding of C* 171 // if (0 < f* - 1/2 < 10^(-x)) then 172 // the result is exact 173 // else // if (f* - 1/2 > T*) then 174 // the result is inexact 175 if (fstar.w[1] > half64[ind] || 176 (fstar.w[1] == half64[ind] && fstar.w[0])) { 177 // f* > 1/2 and the result may be exact 178 // Calculate f* - 1/2 179 tmp64 = fstar.w[1] - half64[ind]; 180 if (tmp64 || fstar.w[0] > ten2mxtrunc64[ind]) { // f* - 1/2 > 10^(-x) 181 *ptr_is_inexact_lt_midpoint = 1; 182 } // else the result is exact 183 } else { // the result is inexact; f2* <= 1/2 184 *ptr_is_inexact_gt_midpoint = 1; 185 } 186 // check for midpoints (could do this before determining inexactness) 187 if (fstar.w[1] == 0 && fstar.w[0] <= ten2mxtrunc64[ind]) { 188 // the result is a midpoint 189 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD] 190 // if floor(C*) is odd C = floor(C*) - 1; the result may be 0 191 Cstar--; // Cstar is now even 192 *ptr_is_midpoint_gt_even = 1; 193 *ptr_is_inexact_lt_midpoint = 0; 194 *ptr_is_inexact_gt_midpoint = 0; 195 } else { // else MP in [ODD, EVEN] 196 *ptr_is_midpoint_lt_even = 1; 197 *ptr_is_inexact_lt_midpoint = 0; 198 *ptr_is_inexact_gt_midpoint = 0; 199 } 200 } 201 // check for rounding overflow, which occurs if Cstar = 10^(q-x) 202 ind = q - x; // 1 <= ind <= q - 1 203 if (Cstar == ten2k64[ind]) { // if Cstar = 10^(q-x) 204 Cstar = ten2k64[ind - 1]; // Cstar = 10^(q-x-1) 205 *incr_exp = 1; 206 } else { // 10^33 <= Cstar <= 10^34 - 1 207 *incr_exp = 0; 208 } 209 *ptr_Cstar = Cstar; 210 } 211 212 213 void 214 round128_19_38 (int q, 215 int x, 216 UINT128 C, 217 UINT128 * ptr_Cstar, 218 int *incr_exp, 219 int *ptr_is_midpoint_lt_even, 220 int *ptr_is_midpoint_gt_even, 221 int *ptr_is_inexact_lt_midpoint, 222 int *ptr_is_inexact_gt_midpoint) { 223 224 UINT256 P256; 225 UINT256 fstar; 226 UINT128 Cstar; 227 UINT64 tmp64; 228 int shift; 229 int ind; 230 231 // Note: 232 // In round128_19_38() positive numbers with 19 <= q <= 38 will be 233 // rounded to nearest only for 1 <= x <= 23: 234 // x = 3 or x = 4 when q = 19 235 // x = 4 or x = 5 when q = 20 236 // ... 237 // x = 18 or x = 19 when q = 34 238 // x = 1 or x = 2 or x = 19 or x = 20 when q = 35 239 // x = 2 or x = 3 or x = 20 or x = 21 when q = 36 240 // x = 3 or x = 4 or x = 21 or x = 22 when q = 37 241 // x = 4 or x = 5 or x = 22 or x = 23 when q = 38 242 // However, for generality and possible uses outside the frame of IEEE 754R 243 // this implementation works for 1 <= x <= q - 1 244 245 // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even, 246 // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are 247 // initialized to 0 by the caller 248 249 // round a number C with q decimal digits, 19 <= q <= 38 250 // to q - x digits, 1 <= x <= 37 251 // C = C + 1/2 * 10^x where the result C fits in 128 bits 252 // (because the largest value is 99999999999999999999999999999999999999 + 253 // 5000000000000000000000000000000000000 = 254 // 0x4efe43b0c573e7e68a043d8fffffffff, which fits is 127 bits) 255 256 ind = x - 1; // 0 <= ind <= 36 257 if (ind <= 18) { // if 0 <= ind <= 18 258 tmp64 = C.w[0]; 259 C.w[0] = C.w[0] + midpoint64[ind]; 260 if (C.w[0] < tmp64) 261 C.w[1]++; 262 } else { // if 19 <= ind <= 37 263 tmp64 = C.w[0]; 264 C.w[0] = C.w[0] + midpoint128[ind - 19].w[0]; 265 if (C.w[0] < tmp64) { 266 C.w[1]++; 267 } 268 C.w[1] = C.w[1] + midpoint128[ind - 19].w[1]; 269 } 270 // kx ~= 10^(-x), kx = Kx128[ind] * 2^(-Ex), 0 <= ind <= 36 271 // P256 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx 272 // the approximation kx of 10^(-x) was rounded up to 128 bits 273 __mul_128x128_to_256 (P256, C, Kx128[ind]); 274 // calculate C* = floor (P256) and f* 275 // Cstar = P256 >> Ex 276 // fstar = low Ex bits of P256 277 shift = Ex128m128[ind]; // in [2, 63] but have to consider two cases 278 if (ind <= 18) { // if 0 <= ind <= 18 279 Cstar.w[0] = (P256.w[2] >> shift) | (P256.w[3] << (64 - shift)); 280 Cstar.w[1] = (P256.w[3] >> shift); 281 fstar.w[0] = P256.w[0]; 282 fstar.w[1] = P256.w[1]; 283 fstar.w[2] = P256.w[2] & mask128[ind]; 284 fstar.w[3] = 0x0ULL; 285 } else { // if 19 <= ind <= 37 286 Cstar.w[0] = P256.w[3] >> shift; 287 Cstar.w[1] = 0x0ULL; 288 fstar.w[0] = P256.w[0]; 289 fstar.w[1] = P256.w[1]; 290 fstar.w[2] = P256.w[2]; 291 fstar.w[3] = P256.w[3] & mask128[ind]; 292 } 293 // the top Ex bits of 10^(-x) are T* = ten2mxtrunc64[ind], e.g. 294 // if x=1, T*=ten2mxtrunc128[0]=0xcccccccccccccccccccccccccccccccc 295 // if (0 < f* < 10^(-x)) then the result is a midpoint 296 // if floor(C*) is even then C* = floor(C*) - logical right 297 // shift; C* has q - x decimal digits, correct by Prop. 1) 298 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 299 // shift; C* has q - x decimal digits, correct by Pr. 1) 300 // else 301 // C* = floor(C*) (logical right shift; C has q - x decimal digits, 302 // correct by Property 1) 303 // in the caling function n = C* * 10^(e+x) 304 305 // determine inexactness of the rounding of C* 306 // if (0 < f* - 1/2 < 10^(-x)) then 307 // the result is exact 308 // else // if (f* - 1/2 > T*) then 309 // the result is inexact 310 if (ind <= 18) { // if 0 <= ind <= 18 311 if (fstar.w[2] > half128[ind] || 312 (fstar.w[2] == half128[ind] && (fstar.w[1] || fstar.w[0]))) { 313 // f* > 1/2 and the result may be exact 314 // Calculate f* - 1/2 315 tmp64 = fstar.w[2] - half128[ind]; 316 if (tmp64 || fstar.w[1] > ten2mxtrunc128[ind].w[1] || (fstar.w[1] == ten2mxtrunc128[ind].w[1] && fstar.w[0] > ten2mxtrunc128[ind].w[0])) { // f* - 1/2 > 10^(-x) 317 *ptr_is_inexact_lt_midpoint = 1; 318 } // else the result is exact 319 } else { // the result is inexact; f2* <= 1/2 320 *ptr_is_inexact_gt_midpoint = 1; 321 } 322 } else { // if 19 <= ind <= 37 323 if (fstar.w[3] > half128[ind] || (fstar.w[3] == half128[ind] && 324 (fstar.w[2] || fstar.w[1] 325 || fstar.w[0]))) { 326 // f* > 1/2 and the result may be exact 327 // Calculate f* - 1/2 328 tmp64 = fstar.w[3] - half128[ind]; 329 if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mxtrunc128[ind].w[1] || (fstar.w[1] == ten2mxtrunc128[ind].w[1] && fstar.w[0] > ten2mxtrunc128[ind].w[0])) { // f* - 1/2 > 10^(-x) 330 *ptr_is_inexact_lt_midpoint = 1; 331 } // else the result is exact 332 } else { // the result is inexact; f2* <= 1/2 333 *ptr_is_inexact_gt_midpoint = 1; 334 } 335 } 336 // check for midpoints (could do this before determining inexactness) 337 if (fstar.w[3] == 0 && fstar.w[2] == 0 && 338 (fstar.w[1] < ten2mxtrunc128[ind].w[1] || 339 (fstar.w[1] == ten2mxtrunc128[ind].w[1] && 340 fstar.w[0] <= ten2mxtrunc128[ind].w[0]))) { 341 // the result is a midpoint 342 if (Cstar.w[0] & 0x01) { // Cstar is odd; MP in [EVEN, ODD] 343 // if floor(C*) is odd C = floor(C*) - 1; the result may be 0 344 Cstar.w[0]--; // Cstar is now even 345 if (Cstar.w[0] == 0xffffffffffffffffULL) { 346 Cstar.w[1]--; 347 } 348 *ptr_is_midpoint_gt_even = 1; 349 *ptr_is_inexact_lt_midpoint = 0; 350 *ptr_is_inexact_gt_midpoint = 0; 351 } else { // else MP in [ODD, EVEN] 352 *ptr_is_midpoint_lt_even = 1; 353 *ptr_is_inexact_lt_midpoint = 0; 354 *ptr_is_inexact_gt_midpoint = 0; 355 } 356 } 357 // check for rounding overflow, which occurs if Cstar = 10^(q-x) 358 ind = q - x; // 1 <= ind <= q - 1 359 if (ind <= 19) { 360 if (Cstar.w[1] == 0x0ULL && Cstar.w[0] == ten2k64[ind]) { 361 // if Cstar = 10^(q-x) 362 Cstar.w[0] = ten2k64[ind - 1]; // Cstar = 10^(q-x-1) 363 *incr_exp = 1; 364 } else { 365 *incr_exp = 0; 366 } 367 } else if (ind == 20) { 368 // if ind = 20 369 if (Cstar.w[1] == ten2k128[0].w[1] 370 && Cstar.w[0] == ten2k128[0].w[0]) { 371 // if Cstar = 10^(q-x) 372 Cstar.w[0] = ten2k64[19]; // Cstar = 10^(q-x-1) 373 Cstar.w[1] = 0x0ULL; 374 *incr_exp = 1; 375 } else { 376 *incr_exp = 0; 377 } 378 } else { // if 21 <= ind <= 37 379 if (Cstar.w[1] == ten2k128[ind - 20].w[1] && 380 Cstar.w[0] == ten2k128[ind - 20].w[0]) { 381 // if Cstar = 10^(q-x) 382 Cstar.w[0] = ten2k128[ind - 21].w[0]; // Cstar = 10^(q-x-1) 383 Cstar.w[1] = ten2k128[ind - 21].w[1]; 384 *incr_exp = 1; 385 } else { 386 *incr_exp = 0; 387 } 388 } 389 ptr_Cstar->w[1] = Cstar.w[1]; 390 ptr_Cstar->w[0] = Cstar.w[0]; 391 } 392 393 394 void 395 round192_39_57 (int q, 396 int x, 397 UINT192 C, 398 UINT192 * ptr_Cstar, 399 int *incr_exp, 400 int *ptr_is_midpoint_lt_even, 401 int *ptr_is_midpoint_gt_even, 402 int *ptr_is_inexact_lt_midpoint, 403 int *ptr_is_inexact_gt_midpoint) { 404 405 UINT384 P384; 406 UINT384 fstar; 407 UINT192 Cstar; 408 UINT64 tmp64; 409 int shift; 410 int ind; 411 412 // Note: 413 // In round192_39_57() positive numbers with 39 <= q <= 57 will be 414 // rounded to nearest only for 5 <= x <= 42: 415 // x = 23 or x = 24 or x = 5 or x = 6 when q = 39 416 // x = 24 or x = 25 or x = 6 or x = 7 when q = 40 417 // ... 418 // x = 41 or x = 42 or x = 23 or x = 24 when q = 57 419 // However, for generality and possible uses outside the frame of IEEE 754R 420 // this implementation works for 1 <= x <= q - 1 421 422 // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even, 423 // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are 424 // initialized to 0 by the caller 425 426 // round a number C with q decimal digits, 39 <= q <= 57 427 // to q - x digits, 1 <= x <= 56 428 // C = C + 1/2 * 10^x where the result C fits in 192 bits 429 // (because the largest value is 430 // 999999999999999999999999999999999999999999999999999999999 + 431 // 50000000000000000000000000000000000000000000000000000000 = 432 // 0x2ad282f212a1da846afdaf18c034ff09da7fffffffffffff, which fits in 190 bits) 433 ind = x - 1; // 0 <= ind <= 55 434 if (ind <= 18) { // if 0 <= ind <= 18 435 tmp64 = C.w[0]; 436 C.w[0] = C.w[0] + midpoint64[ind]; 437 if (C.w[0] < tmp64) { 438 C.w[1]++; 439 if (C.w[1] == 0x0) { 440 C.w[2]++; 441 } 442 } 443 } else if (ind <= 37) { // if 19 <= ind <= 37 444 tmp64 = C.w[0]; 445 C.w[0] = C.w[0] + midpoint128[ind - 19].w[0]; 446 if (C.w[0] < tmp64) { 447 C.w[1]++; 448 if (C.w[1] == 0x0) { 449 C.w[2]++; 450 } 451 } 452 tmp64 = C.w[1]; 453 C.w[1] = C.w[1] + midpoint128[ind - 19].w[1]; 454 if (C.w[1] < tmp64) { 455 C.w[2]++; 456 } 457 } else { // if 38 <= ind <= 57 (actually ind <= 55) 458 tmp64 = C.w[0]; 459 C.w[0] = C.w[0] + midpoint192[ind - 38].w[0]; 460 if (C.w[0] < tmp64) { 461 C.w[1]++; 462 if (C.w[1] == 0x0ull) { 463 C.w[2]++; 464 } 465 } 466 tmp64 = C.w[1]; 467 C.w[1] = C.w[1] + midpoint192[ind - 38].w[1]; 468 if (C.w[1] < tmp64) { 469 C.w[2]++; 470 } 471 C.w[2] = C.w[2] + midpoint192[ind - 38].w[2]; 472 } 473 // kx ~= 10^(-x), kx = Kx192[ind] * 2^(-Ex), 0 <= ind <= 55 474 // P384 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx 475 // the approximation kx of 10^(-x) was rounded up to 192 bits 476 __mul_192x192_to_384 (P384, C, Kx192[ind]); 477 // calculate C* = floor (P384) and f* 478 // Cstar = P384 >> Ex 479 // fstar = low Ex bits of P384 480 shift = Ex192m192[ind]; // in [1, 63] but have to consider three cases 481 if (ind <= 18) { // if 0 <= ind <= 18 482 Cstar.w[2] = (P384.w[5] >> shift); 483 Cstar.w[1] = (P384.w[5] << (64 - shift)) | (P384.w[4] >> shift); 484 Cstar.w[0] = (P384.w[4] << (64 - shift)) | (P384.w[3] >> shift); 485 fstar.w[5] = 0x0ULL; 486 fstar.w[4] = 0x0ULL; 487 fstar.w[3] = P384.w[3] & mask192[ind]; 488 fstar.w[2] = P384.w[2]; 489 fstar.w[1] = P384.w[1]; 490 fstar.w[0] = P384.w[0]; 491 } else if (ind <= 37) { // if 19 <= ind <= 37 492 Cstar.w[2] = 0x0ULL; 493 Cstar.w[1] = P384.w[5] >> shift; 494 Cstar.w[0] = (P384.w[5] << (64 - shift)) | (P384.w[4] >> shift); 495 fstar.w[5] = 0x0ULL; 496 fstar.w[4] = P384.w[4] & mask192[ind]; 497 fstar.w[3] = P384.w[3]; 498 fstar.w[2] = P384.w[2]; 499 fstar.w[1] = P384.w[1]; 500 fstar.w[0] = P384.w[0]; 501 } else { // if 38 <= ind <= 57 502 Cstar.w[2] = 0x0ULL; 503 Cstar.w[1] = 0x0ULL; 504 Cstar.w[0] = P384.w[5] >> shift; 505 fstar.w[5] = P384.w[5] & mask192[ind]; 506 fstar.w[4] = P384.w[4]; 507 fstar.w[3] = P384.w[3]; 508 fstar.w[2] = P384.w[2]; 509 fstar.w[1] = P384.w[1]; 510 fstar.w[0] = P384.w[0]; 511 } 512 513 // the top Ex bits of 10^(-x) are T* = ten2mxtrunc192[ind], e.g. if x=1, 514 // T*=ten2mxtrunc192[0]=0xcccccccccccccccccccccccccccccccccccccccccccccccc 515 // if (0 < f* < 10^(-x)) then the result is a midpoint 516 // if floor(C*) is even then C* = floor(C*) - logical right 517 // shift; C* has q - x decimal digits, correct by Prop. 1) 518 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 519 // shift; C* has q - x decimal digits, correct by Pr. 1) 520 // else 521 // C* = floor(C*) (logical right shift; C has q - x decimal digits, 522 // correct by Property 1) 523 // in the caling function n = C* * 10^(e+x) 524 525 // determine inexactness of the rounding of C* 526 // if (0 < f* - 1/2 < 10^(-x)) then 527 // the result is exact 528 // else // if (f* - 1/2 > T*) then 529 // the result is inexact 530 if (ind <= 18) { // if 0 <= ind <= 18 531 if (fstar.w[3] > half192[ind] || (fstar.w[3] == half192[ind] && 532 (fstar.w[2] || fstar.w[1] 533 || fstar.w[0]))) { 534 // f* > 1/2 and the result may be exact 535 // Calculate f* - 1/2 536 tmp64 = fstar.w[3] - half192[ind]; 537 if (tmp64 || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) { // f* - 1/2 > 10^(-x) 538 *ptr_is_inexact_lt_midpoint = 1; 539 } // else the result is exact 540 } else { // the result is inexact; f2* <= 1/2 541 *ptr_is_inexact_gt_midpoint = 1; 542 } 543 } else if (ind <= 37) { // if 19 <= ind <= 37 544 if (fstar.w[4] > half192[ind] || (fstar.w[4] == half192[ind] && 545 (fstar.w[3] || fstar.w[2] 546 || fstar.w[1] || fstar.w[0]))) { 547 // f* > 1/2 and the result may be exact 548 // Calculate f* - 1/2 549 tmp64 = fstar.w[4] - half192[ind]; 550 if (tmp64 || fstar.w[3] || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) { // f* - 1/2 > 10^(-x) 551 *ptr_is_inexact_lt_midpoint = 1; 552 } // else the result is exact 553 } else { // the result is inexact; f2* <= 1/2 554 *ptr_is_inexact_gt_midpoint = 1; 555 } 556 } else { // if 38 <= ind <= 55 557 if (fstar.w[5] > half192[ind] || (fstar.w[5] == half192[ind] && 558 (fstar.w[4] || fstar.w[3] 559 || fstar.w[2] || fstar.w[1] 560 || fstar.w[0]))) { 561 // f* > 1/2 and the result may be exact 562 // Calculate f* - 1/2 563 tmp64 = fstar.w[5] - half192[ind]; 564 if (tmp64 || fstar.w[4] || fstar.w[3] || fstar.w[2] > ten2mxtrunc192[ind].w[2] || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] > ten2mxtrunc192[ind].w[1]) || (fstar.w[2] == ten2mxtrunc192[ind].w[2] && fstar.w[1] == ten2mxtrunc192[ind].w[1] && fstar.w[0] > ten2mxtrunc192[ind].w[0])) { // f* - 1/2 > 10^(-x) 565 *ptr_is_inexact_lt_midpoint = 1; 566 } // else the result is exact 567 } else { // the result is inexact; f2* <= 1/2 568 *ptr_is_inexact_gt_midpoint = 1; 569 } 570 } 571 // check for midpoints (could do this before determining inexactness) 572 if (fstar.w[5] == 0 && fstar.w[4] == 0 && fstar.w[3] == 0 && 573 (fstar.w[2] < ten2mxtrunc192[ind].w[2] || 574 (fstar.w[2] == ten2mxtrunc192[ind].w[2] && 575 fstar.w[1] < ten2mxtrunc192[ind].w[1]) || 576 (fstar.w[2] == ten2mxtrunc192[ind].w[2] && 577 fstar.w[1] == ten2mxtrunc192[ind].w[1] && 578 fstar.w[0] <= ten2mxtrunc192[ind].w[0]))) { 579 // the result is a midpoint 580 if (Cstar.w[0] & 0x01) { // Cstar is odd; MP in [EVEN, ODD] 581 // if floor(C*) is odd C = floor(C*) - 1; the result may be 0 582 Cstar.w[0]--; // Cstar is now even 583 if (Cstar.w[0] == 0xffffffffffffffffULL) { 584 Cstar.w[1]--; 585 if (Cstar.w[1] == 0xffffffffffffffffULL) { 586 Cstar.w[2]--; 587 } 588 } 589 *ptr_is_midpoint_gt_even = 1; 590 *ptr_is_inexact_lt_midpoint = 0; 591 *ptr_is_inexact_gt_midpoint = 0; 592 } else { // else MP in [ODD, EVEN] 593 *ptr_is_midpoint_lt_even = 1; 594 *ptr_is_inexact_lt_midpoint = 0; 595 *ptr_is_inexact_gt_midpoint = 0; 596 } 597 } 598 // check for rounding overflow, which occurs if Cstar = 10^(q-x) 599 ind = q - x; // 1 <= ind <= q - 1 600 if (ind <= 19) { 601 if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == 0x0ULL && 602 Cstar.w[0] == ten2k64[ind]) { 603 // if Cstar = 10^(q-x) 604 Cstar.w[0] = ten2k64[ind - 1]; // Cstar = 10^(q-x-1) 605 *incr_exp = 1; 606 } else { 607 *incr_exp = 0; 608 } 609 } else if (ind == 20) { 610 // if ind = 20 611 if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == ten2k128[0].w[1] && 612 Cstar.w[0] == ten2k128[0].w[0]) { 613 // if Cstar = 10^(q-x) 614 Cstar.w[0] = ten2k64[19]; // Cstar = 10^(q-x-1) 615 Cstar.w[1] = 0x0ULL; 616 *incr_exp = 1; 617 } else { 618 *incr_exp = 0; 619 } 620 } else if (ind <= 38) { // if 21 <= ind <= 38 621 if (Cstar.w[2] == 0x0ULL && Cstar.w[1] == ten2k128[ind - 20].w[1] && 622 Cstar.w[0] == ten2k128[ind - 20].w[0]) { 623 // if Cstar = 10^(q-x) 624 Cstar.w[0] = ten2k128[ind - 21].w[0]; // Cstar = 10^(q-x-1) 625 Cstar.w[1] = ten2k128[ind - 21].w[1]; 626 *incr_exp = 1; 627 } else { 628 *incr_exp = 0; 629 } 630 } else if (ind == 39) { 631 if (Cstar.w[2] == ten2k256[0].w[2] && Cstar.w[1] == ten2k256[0].w[1] 632 && Cstar.w[0] == ten2k256[0].w[0]) { 633 // if Cstar = 10^(q-x) 634 Cstar.w[0] = ten2k128[18].w[0]; // Cstar = 10^(q-x-1) 635 Cstar.w[1] = ten2k128[18].w[1]; 636 Cstar.w[2] = 0x0ULL; 637 *incr_exp = 1; 638 } else { 639 *incr_exp = 0; 640 } 641 } else { // if 40 <= ind <= 56 642 if (Cstar.w[2] == ten2k256[ind - 39].w[2] && 643 Cstar.w[1] == ten2k256[ind - 39].w[1] && 644 Cstar.w[0] == ten2k256[ind - 39].w[0]) { 645 // if Cstar = 10^(q-x) 646 Cstar.w[0] = ten2k256[ind - 40].w[0]; // Cstar = 10^(q-x-1) 647 Cstar.w[1] = ten2k256[ind - 40].w[1]; 648 Cstar.w[2] = ten2k256[ind - 40].w[2]; 649 *incr_exp = 1; 650 } else { 651 *incr_exp = 0; 652 } 653 } 654 ptr_Cstar->w[2] = Cstar.w[2]; 655 ptr_Cstar->w[1] = Cstar.w[1]; 656 ptr_Cstar->w[0] = Cstar.w[0]; 657 } 658 659 660 void 661 round256_58_76 (int q, 662 int x, 663 UINT256 C, 664 UINT256 * ptr_Cstar, 665 int *incr_exp, 666 int *ptr_is_midpoint_lt_even, 667 int *ptr_is_midpoint_gt_even, 668 int *ptr_is_inexact_lt_midpoint, 669 int *ptr_is_inexact_gt_midpoint) { 670 671 UINT512 P512; 672 UINT512 fstar; 673 UINT256 Cstar; 674 UINT64 tmp64; 675 int shift; 676 int ind; 677 678 // Note: 679 // In round256_58_76() positive numbers with 58 <= q <= 76 will be 680 // rounded to nearest only for 24 <= x <= 61: 681 // x = 42 or x = 43 or x = 24 or x = 25 when q = 58 682 // x = 43 or x = 44 or x = 25 or x = 26 when q = 59 683 // ... 684 // x = 60 or x = 61 or x = 42 or x = 43 when q = 76 685 // However, for generality and possible uses outside the frame of IEEE 754R 686 // this implementation works for 1 <= x <= q - 1 687 688 // assume *ptr_is_midpoint_lt_even, *ptr_is_midpoint_gt_even, 689 // *ptr_is_inexact_lt_midpoint, and *ptr_is_inexact_gt_midpoint are 690 // initialized to 0 by the caller 691 692 // round a number C with q decimal digits, 58 <= q <= 76 693 // to q - x digits, 1 <= x <= 75 694 // C = C + 1/2 * 10^x where the result C fits in 256 bits 695 // (because the largest value is 9999999999999999999999999999999999999999 696 // 999999999999999999999999999999999999 + 500000000000000000000000000 697 // 000000000000000000000000000000000000000000000000 = 698 // 0x1736ca15d27a56cae15cf0e7b403d1f2bd6ebb0a50dc83ffffffffffffffffff, 699 // which fits in 253 bits) 700 ind = x - 1; // 0 <= ind <= 74 701 if (ind <= 18) { // if 0 <= ind <= 18 702 tmp64 = C.w[0]; 703 C.w[0] = C.w[0] + midpoint64[ind]; 704 if (C.w[0] < tmp64) { 705 C.w[1]++; 706 if (C.w[1] == 0x0) { 707 C.w[2]++; 708 if (C.w[2] == 0x0) { 709 C.w[3]++; 710 } 711 } 712 } 713 } else if (ind <= 37) { // if 19 <= ind <= 37 714 tmp64 = C.w[0]; 715 C.w[0] = C.w[0] + midpoint128[ind - 19].w[0]; 716 if (C.w[0] < tmp64) { 717 C.w[1]++; 718 if (C.w[1] == 0x0) { 719 C.w[2]++; 720 if (C.w[2] == 0x0) { 721 C.w[3]++; 722 } 723 } 724 } 725 tmp64 = C.w[1]; 726 C.w[1] = C.w[1] + midpoint128[ind - 19].w[1]; 727 if (C.w[1] < tmp64) { 728 C.w[2]++; 729 if (C.w[2] == 0x0) { 730 C.w[3]++; 731 } 732 } 733 } else if (ind <= 57) { // if 38 <= ind <= 57 734 tmp64 = C.w[0]; 735 C.w[0] = C.w[0] + midpoint192[ind - 38].w[0]; 736 if (C.w[0] < tmp64) { 737 C.w[1]++; 738 if (C.w[1] == 0x0ull) { 739 C.w[2]++; 740 if (C.w[2] == 0x0) { 741 C.w[3]++; 742 } 743 } 744 } 745 tmp64 = C.w[1]; 746 C.w[1] = C.w[1] + midpoint192[ind - 38].w[1]; 747 if (C.w[1] < tmp64) { 748 C.w[2]++; 749 if (C.w[2] == 0x0) { 750 C.w[3]++; 751 } 752 } 753 tmp64 = C.w[2]; 754 C.w[2] = C.w[2] + midpoint192[ind - 38].w[2]; 755 if (C.w[2] < tmp64) { 756 C.w[3]++; 757 } 758 } else { // if 58 <= ind <= 76 (actually 58 <= ind <= 74) 759 tmp64 = C.w[0]; 760 C.w[0] = C.w[0] + midpoint256[ind - 58].w[0]; 761 if (C.w[0] < tmp64) { 762 C.w[1]++; 763 if (C.w[1] == 0x0ull) { 764 C.w[2]++; 765 if (C.w[2] == 0x0) { 766 C.w[3]++; 767 } 768 } 769 } 770 tmp64 = C.w[1]; 771 C.w[1] = C.w[1] + midpoint256[ind - 58].w[1]; 772 if (C.w[1] < tmp64) { 773 C.w[2]++; 774 if (C.w[2] == 0x0) { 775 C.w[3]++; 776 } 777 } 778 tmp64 = C.w[2]; 779 C.w[2] = C.w[2] + midpoint256[ind - 58].w[2]; 780 if (C.w[2] < tmp64) { 781 C.w[3]++; 782 } 783 C.w[3] = C.w[3] + midpoint256[ind - 58].w[3]; 784 } 785 // kx ~= 10^(-x), kx = Kx256[ind] * 2^(-Ex), 0 <= ind <= 74 786 // P512 = (C + 1/2 * 10^x) * kx * 2^Ex = (C + 1/2 * 10^x) * Kx 787 // the approximation kx of 10^(-x) was rounded up to 192 bits 788 __mul_256x256_to_512 (P512, C, Kx256[ind]); 789 // calculate C* = floor (P512) and f* 790 // Cstar = P512 >> Ex 791 // fstar = low Ex bits of P512 792 shift = Ex256m256[ind]; // in [0, 63] but have to consider four cases 793 if (ind <= 18) { // if 0 <= ind <= 18 794 Cstar.w[3] = (P512.w[7] >> shift); 795 Cstar.w[2] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift); 796 Cstar.w[1] = (P512.w[6] << (64 - shift)) | (P512.w[5] >> shift); 797 Cstar.w[0] = (P512.w[5] << (64 - shift)) | (P512.w[4] >> shift); 798 fstar.w[7] = 0x0ULL; 799 fstar.w[6] = 0x0ULL; 800 fstar.w[5] = 0x0ULL; 801 fstar.w[4] = P512.w[4] & mask256[ind]; 802 fstar.w[3] = P512.w[3]; 803 fstar.w[2] = P512.w[2]; 804 fstar.w[1] = P512.w[1]; 805 fstar.w[0] = P512.w[0]; 806 } else if (ind <= 37) { // if 19 <= ind <= 37 807 Cstar.w[3] = 0x0ULL; 808 Cstar.w[2] = P512.w[7] >> shift; 809 Cstar.w[1] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift); 810 Cstar.w[0] = (P512.w[6] << (64 - shift)) | (P512.w[5] >> shift); 811 fstar.w[7] = 0x0ULL; 812 fstar.w[6] = 0x0ULL; 813 fstar.w[5] = P512.w[5] & mask256[ind]; 814 fstar.w[4] = P512.w[4]; 815 fstar.w[3] = P512.w[3]; 816 fstar.w[2] = P512.w[2]; 817 fstar.w[1] = P512.w[1]; 818 fstar.w[0] = P512.w[0]; 819 } else if (ind <= 56) { // if 38 <= ind <= 56 820 Cstar.w[3] = 0x0ULL; 821 Cstar.w[2] = 0x0ULL; 822 Cstar.w[1] = P512.w[7] >> shift; 823 Cstar.w[0] = (P512.w[7] << (64 - shift)) | (P512.w[6] >> shift); 824 fstar.w[7] = 0x0ULL; 825 fstar.w[6] = P512.w[6] & mask256[ind]; 826 fstar.w[5] = P512.w[5]; 827 fstar.w[4] = P512.w[4]; 828 fstar.w[3] = P512.w[3]; 829 fstar.w[2] = P512.w[2]; 830 fstar.w[1] = P512.w[1]; 831 fstar.w[0] = P512.w[0]; 832 } else if (ind == 57) { 833 Cstar.w[3] = 0x0ULL; 834 Cstar.w[2] = 0x0ULL; 835 Cstar.w[1] = 0x0ULL; 836 Cstar.w[0] = P512.w[7]; 837 fstar.w[7] = 0x0ULL; 838 fstar.w[6] = P512.w[6]; 839 fstar.w[5] = P512.w[5]; 840 fstar.w[4] = P512.w[4]; 841 fstar.w[3] = P512.w[3]; 842 fstar.w[2] = P512.w[2]; 843 fstar.w[1] = P512.w[1]; 844 fstar.w[0] = P512.w[0]; 845 } else { // if 58 <= ind <= 74 846 Cstar.w[3] = 0x0ULL; 847 Cstar.w[2] = 0x0ULL; 848 Cstar.w[1] = 0x0ULL; 849 Cstar.w[0] = P512.w[7] >> shift; 850 fstar.w[7] = P512.w[7] & mask256[ind]; 851 fstar.w[6] = P512.w[6]; 852 fstar.w[5] = P512.w[5]; 853 fstar.w[4] = P512.w[4]; 854 fstar.w[3] = P512.w[3]; 855 fstar.w[2] = P512.w[2]; 856 fstar.w[1] = P512.w[1]; 857 fstar.w[0] = P512.w[0]; 858 } 859 860 // the top Ex bits of 10^(-x) are T* = ten2mxtrunc256[ind], e.g. if x=1, 861 // T*=ten2mxtrunc256[0]= 862 // 0xcccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 863 // if (0 < f* < 10^(-x)) then the result is a midpoint 864 // if floor(C*) is even then C* = floor(C*) - logical right 865 // shift; C* has q - x decimal digits, correct by Prop. 1) 866 // else if floor(C*) is odd C* = floor(C*)-1 (logical right 867 // shift; C* has q - x decimal digits, correct by Pr. 1) 868 // else 869 // C* = floor(C*) (logical right shift; C has q - x decimal digits, 870 // correct by Property 1) 871 // in the caling function n = C* * 10^(e+x) 872 873 // determine inexactness of the rounding of C* 874 // if (0 < f* - 1/2 < 10^(-x)) then 875 // the result is exact 876 // else // if (f* - 1/2 > T*) then 877 // the result is inexact 878 if (ind <= 18) { // if 0 <= ind <= 18 879 if (fstar.w[4] > half256[ind] || (fstar.w[4] == half256[ind] && 880 (fstar.w[3] || fstar.w[2] 881 || fstar.w[1] || fstar.w[0]))) { 882 // f* > 1/2 and the result may be exact 883 // Calculate f* - 1/2 884 tmp64 = fstar.w[4] - half256[ind]; 885 if (tmp64 || fstar.w[3] > ten2mxtrunc256[ind].w[2] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) { // f* - 1/2 > 10^(-x) 886 *ptr_is_inexact_lt_midpoint = 1; 887 } // else the result is exact 888 } else { // the result is inexact; f2* <= 1/2 889 *ptr_is_inexact_gt_midpoint = 1; 890 } 891 } else if (ind <= 37) { // if 19 <= ind <= 37 892 if (fstar.w[5] > half256[ind] || (fstar.w[5] == half256[ind] && 893 (fstar.w[4] || fstar.w[3] 894 || fstar.w[2] || fstar.w[1] 895 || fstar.w[0]))) { 896 // f* > 1/2 and the result may be exact 897 // Calculate f* - 1/2 898 tmp64 = fstar.w[5] - half256[ind]; 899 if (tmp64 || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) { // f* - 1/2 > 10^(-x) 900 *ptr_is_inexact_lt_midpoint = 1; 901 } // else the result is exact 902 } else { // the result is inexact; f2* <= 1/2 903 *ptr_is_inexact_gt_midpoint = 1; 904 } 905 } else if (ind <= 57) { // if 38 <= ind <= 57 906 if (fstar.w[6] > half256[ind] || (fstar.w[6] == half256[ind] && 907 (fstar.w[5] || fstar.w[4] 908 || fstar.w[3] || fstar.w[2] 909 || fstar.w[1] || fstar.w[0]))) { 910 // f* > 1/2 and the result may be exact 911 // Calculate f* - 1/2 912 tmp64 = fstar.w[6] - half256[ind]; 913 if (tmp64 || fstar.w[5] || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) { // f* - 1/2 > 10^(-x) 914 *ptr_is_inexact_lt_midpoint = 1; 915 } // else the result is exact 916 } else { // the result is inexact; f2* <= 1/2 917 *ptr_is_inexact_gt_midpoint = 1; 918 } 919 } else { // if 58 <= ind <= 74 920 if (fstar.w[7] > half256[ind] || (fstar.w[7] == half256[ind] && 921 (fstar.w[6] || fstar.w[5] 922 || fstar.w[4] || fstar.w[3] 923 || fstar.w[2] || fstar.w[1] 924 || fstar.w[0]))) { 925 // f* > 1/2 and the result may be exact 926 // Calculate f* - 1/2 927 tmp64 = fstar.w[7] - half256[ind]; 928 if (tmp64 || fstar.w[6] || fstar.w[5] || fstar.w[4] || fstar.w[3] > ten2mxtrunc256[ind].w[3] || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] > ten2mxtrunc256[ind].w[2]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] > ten2mxtrunc256[ind].w[1]) || (fstar.w[3] == ten2mxtrunc256[ind].w[3] && fstar.w[2] == ten2mxtrunc256[ind].w[2] && fstar.w[1] == ten2mxtrunc256[ind].w[1] && fstar.w[0] > ten2mxtrunc256[ind].w[0])) { // f* - 1/2 > 10^(-x) 929 *ptr_is_inexact_lt_midpoint = 1; 930 } // else the result is exact 931 } else { // the result is inexact; f2* <= 1/2 932 *ptr_is_inexact_gt_midpoint = 1; 933 } 934 } 935 // check for midpoints (could do this before determining inexactness) 936 if (fstar.w[7] == 0 && fstar.w[6] == 0 && 937 fstar.w[5] == 0 && fstar.w[4] == 0 && 938 (fstar.w[3] < ten2mxtrunc256[ind].w[3] || 939 (fstar.w[3] == ten2mxtrunc256[ind].w[3] && 940 fstar.w[2] < ten2mxtrunc256[ind].w[2]) || 941 (fstar.w[3] == ten2mxtrunc256[ind].w[3] && 942 fstar.w[2] == ten2mxtrunc256[ind].w[2] && 943 fstar.w[1] < ten2mxtrunc256[ind].w[1]) || 944 (fstar.w[3] == ten2mxtrunc256[ind].w[3] && 945 fstar.w[2] == ten2mxtrunc256[ind].w[2] && 946 fstar.w[1] == ten2mxtrunc256[ind].w[1] && 947 fstar.w[0] <= ten2mxtrunc256[ind].w[0]))) { 948 // the result is a midpoint 949 if (Cstar.w[0] & 0x01) { // Cstar is odd; MP in [EVEN, ODD] 950 // if floor(C*) is odd C = floor(C*) - 1; the result may be 0 951 Cstar.w[0]--; // Cstar is now even 952 if (Cstar.w[0] == 0xffffffffffffffffULL) { 953 Cstar.w[1]--; 954 if (Cstar.w[1] == 0xffffffffffffffffULL) { 955 Cstar.w[2]--; 956 if (Cstar.w[2] == 0xffffffffffffffffULL) { 957 Cstar.w[3]--; 958 } 959 } 960 } 961 *ptr_is_midpoint_gt_even = 1; 962 *ptr_is_inexact_lt_midpoint = 0; 963 *ptr_is_inexact_gt_midpoint = 0; 964 } else { // else MP in [ODD, EVEN] 965 *ptr_is_midpoint_lt_even = 1; 966 *ptr_is_inexact_lt_midpoint = 0; 967 *ptr_is_inexact_gt_midpoint = 0; 968 } 969 } 970 // check for rounding overflow, which occurs if Cstar = 10^(q-x) 971 ind = q - x; // 1 <= ind <= q - 1 972 if (ind <= 19) { 973 if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL && 974 Cstar.w[1] == 0x0ULL && Cstar.w[0] == ten2k64[ind]) { 975 // if Cstar = 10^(q-x) 976 Cstar.w[0] = ten2k64[ind - 1]; // Cstar = 10^(q-x-1) 977 *incr_exp = 1; 978 } else { 979 *incr_exp = 0; 980 } 981 } else if (ind == 20) { 982 // if ind = 20 983 if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL && 984 Cstar.w[1] == ten2k128[0].w[1] 985 && Cstar.w[0] == ten2k128[0].w[0]) { 986 // if Cstar = 10^(q-x) 987 Cstar.w[0] = ten2k64[19]; // Cstar = 10^(q-x-1) 988 Cstar.w[1] = 0x0ULL; 989 *incr_exp = 1; 990 } else { 991 *incr_exp = 0; 992 } 993 } else if (ind <= 38) { // if 21 <= ind <= 38 994 if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == 0x0ULL && 995 Cstar.w[1] == ten2k128[ind - 20].w[1] && 996 Cstar.w[0] == ten2k128[ind - 20].w[0]) { 997 // if Cstar = 10^(q-x) 998 Cstar.w[0] = ten2k128[ind - 21].w[0]; // Cstar = 10^(q-x-1) 999 Cstar.w[1] = ten2k128[ind - 21].w[1]; 1000 *incr_exp = 1; 1001 } else { 1002 *incr_exp = 0; 1003 } 1004 } else if (ind == 39) { 1005 if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == ten2k256[0].w[2] && 1006 Cstar.w[1] == ten2k256[0].w[1] 1007 && Cstar.w[0] == ten2k256[0].w[0]) { 1008 // if Cstar = 10^(q-x) 1009 Cstar.w[0] = ten2k128[18].w[0]; // Cstar = 10^(q-x-1) 1010 Cstar.w[1] = ten2k128[18].w[1]; 1011 Cstar.w[2] = 0x0ULL; 1012 *incr_exp = 1; 1013 } else { 1014 *incr_exp = 0; 1015 } 1016 } else if (ind <= 57) { // if 40 <= ind <= 57 1017 if (Cstar.w[3] == 0x0ULL && Cstar.w[2] == ten2k256[ind - 39].w[2] && 1018 Cstar.w[1] == ten2k256[ind - 39].w[1] && 1019 Cstar.w[0] == ten2k256[ind - 39].w[0]) { 1020 // if Cstar = 10^(q-x) 1021 Cstar.w[0] = ten2k256[ind - 40].w[0]; // Cstar = 10^(q-x-1) 1022 Cstar.w[1] = ten2k256[ind - 40].w[1]; 1023 Cstar.w[2] = ten2k256[ind - 40].w[2]; 1024 *incr_exp = 1; 1025 } else { 1026 *incr_exp = 0; 1027 } 1028 // else if (ind == 58) is not needed becauae we do not have ten2k192[] yet 1029 } else { // if 58 <= ind <= 77 (actually 58 <= ind <= 74) 1030 if (Cstar.w[3] == ten2k256[ind - 39].w[3] && 1031 Cstar.w[2] == ten2k256[ind - 39].w[2] && 1032 Cstar.w[1] == ten2k256[ind - 39].w[1] && 1033 Cstar.w[0] == ten2k256[ind - 39].w[0]) { 1034 // if Cstar = 10^(q-x) 1035 Cstar.w[0] = ten2k256[ind - 40].w[0]; // Cstar = 10^(q-x-1) 1036 Cstar.w[1] = ten2k256[ind - 40].w[1]; 1037 Cstar.w[2] = ten2k256[ind - 40].w[2]; 1038 Cstar.w[3] = ten2k256[ind - 40].w[3]; 1039 *incr_exp = 1; 1040 } else { 1041 *incr_exp = 0; 1042 } 1043 } 1044 ptr_Cstar->w[3] = Cstar.w[3]; 1045 ptr_Cstar->w[2] = Cstar.w[2]; 1046 ptr_Cstar->w[1] = Cstar.w[1]; 1047 ptr_Cstar->w[0] = Cstar.w[0]; 1048 1049 } 1050