1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 18 #include <string.h> 19 #include "cbigint.h" 20 21 #if defined(__linux__) || defined(FREEBSD) || defined(ZOS) || defined(MACOSX) || defined(AIX) 22 #define USE_LL 23 #endif 24 25 #if __BYTE_ORDER == __LITTLE_ENDIAN 26 #define at(i) (i) 27 #else 28 #define at(i) ((i)^1) 29 /* the sequence for halfAt is -1, 2, 1, 4, 3, 6, 5, 8... */ 30 /* and it should correspond to 0, 1, 2, 3, 4, 5, 6, 7... */ 31 #define halfAt(i) (-((-(i)) ^ 1)) 32 #endif 33 34 #define HIGH_IN_U64(u64) ((u64) >> 32) 35 #if defined(USE_LL) 36 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFLL) 37 #else 38 #if defined(USE_L) 39 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFL) 40 #else 41 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFF) 42 #endif /* USE_L */ 43 #endif /* USE_LL */ 44 45 #if defined(USE_LL) 46 #define TEN_E1 (0xALL) 47 #define TEN_E2 (0x64LL) 48 #define TEN_E3 (0x3E8LL) 49 #define TEN_E4 (0x2710LL) 50 #define TEN_E5 (0x186A0LL) 51 #define TEN_E6 (0xF4240LL) 52 #define TEN_E7 (0x989680LL) 53 #define TEN_E8 (0x5F5E100LL) 54 #define TEN_E9 (0x3B9ACA00LL) 55 #define TEN_E19 (0x8AC7230489E80000LL) 56 #else 57 #if defined(USE_L) 58 #define TEN_E1 (0xAL) 59 #define TEN_E2 (0x64L) 60 #define TEN_E3 (0x3E8L) 61 #define TEN_E4 (0x2710L) 62 #define TEN_E5 (0x186A0L) 63 #define TEN_E6 (0xF4240L) 64 #define TEN_E7 (0x989680L) 65 #define TEN_E8 (0x5F5E100L) 66 #define TEN_E9 (0x3B9ACA00L) 67 #define TEN_E19 (0x8AC7230489E80000L) 68 #else 69 #define TEN_E1 (0xA) 70 #define TEN_E2 (0x64) 71 #define TEN_E3 (0x3E8) 72 #define TEN_E4 (0x2710) 73 #define TEN_E5 (0x186A0) 74 #define TEN_E6 (0xF4240) 75 #define TEN_E7 (0x989680) 76 #define TEN_E8 (0x5F5E100) 77 #define TEN_E9 (0x3B9ACA00) 78 #define TEN_E19 (0x8AC7230489E80000) 79 #endif /* USE_L */ 80 #endif /* USE_LL */ 81 82 #define TIMES_TEN(x) (((x) << 3) + ((x) << 1)) 83 #define bitSection(x, mask, shift) (((x) & (mask)) >> (shift)) 84 #define CREATE_DOUBLE_BITS(normalizedM, e) (((normalizedM) & MANTISSA_MASK) | ((static_cast<uint64_t>((e) + E_OFFSET)) << 52)) 85 86 #if defined(USE_LL) 87 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFLL) 88 #define EXPONENT_MASK (0x7FF0000000000000LL) 89 #define NORMAL_MASK (0x0010000000000000LL) 90 #define SIGN_MASK (0x8000000000000000LL) 91 #else 92 #if defined(USE_L) 93 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFL) 94 #define EXPONENT_MASK (0x7FF0000000000000L) 95 #define NORMAL_MASK (0x0010000000000000L) 96 #define SIGN_MASK (0x8000000000000000L) 97 #else 98 #define MANTISSA_MASK (0x000FFFFFFFFFFFFF) 99 #define EXPONENT_MASK (0x7FF0000000000000) 100 #define NORMAL_MASK (0x0010000000000000) 101 #define SIGN_MASK (0x8000000000000000) 102 #endif /* USE_L */ 103 #endif /* USE_LL */ 104 105 #define E_OFFSET (1075) 106 107 #define FLOAT_MANTISSA_MASK (0x007FFFFF) 108 #define FLOAT_EXPONENT_MASK (0x7F800000) 109 #define FLOAT_NORMAL_MASK (0x00800000) 110 #define FLOAT_E_OFFSET (150) 111 112 int32_t 113 simpleAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2) 114 { 115 /* assumes length > 0 */ 116 int32_t index = 1; 117 118 *arg1 += arg2; 119 if (arg2 <= *arg1) 120 return 0; 121 else if (length == 1) 122 return 1; 123 124 while (++arg1[index] == 0 && ++index < length); 125 126 return index == length; 127 } 128 129 int32_t 130 addHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2) 131 { 132 /* addition is limited by length of arg1 as it this function is 133 * storing the result in arg1 */ 134 /* fix for cc (GCC) 3.2 20020903 (Red Hat Linux 8.0 3.2-7): code generated does not 135 * do the temp1 + temp2 + carry addition correct. carry is 64 bit because gcc has 136 * subtle issues when you mix 64 / 32 bit maths. */ 137 uint64_t temp1, temp2, temp3; /* temporary variables to help the SH-4, and gcc */ 138 uint64_t carry; 139 int32_t index; 140 141 if (length1 == 0 || length2 == 0) 142 { 143 return 0; 144 } 145 else if (length1 < length2) 146 { 147 length2 = length1; 148 } 149 150 carry = 0; 151 index = 0; 152 do 153 { 154 temp1 = arg1[index]; 155 temp2 = arg2[index]; 156 temp3 = temp1 + temp2; 157 arg1[index] = temp3 + carry; 158 if (arg2[index] < arg1[index]) 159 carry = 0; 160 else if (arg2[index] != arg1[index]) 161 carry = 1; 162 } 163 while (++index < length2); 164 if (!carry) 165 return 0; 166 else if (index == length1) 167 return 1; 168 169 while (++arg1[index] == 0 && ++index < length1); 170 171 return index == length1; 172 } 173 174 void 175 subtractHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2) 176 { 177 /* assumes arg1 > arg2 */ 178 int32_t index; 179 for (index = 0; index < length1; ++index) 180 arg1[index] = ~arg1[index]; 181 simpleAddHighPrecision (arg1, length1, 1); 182 183 while (length2 > 0 && arg2[length2 - 1] == 0) 184 --length2; 185 186 addHighPrecision (arg1, length1, arg2, length2); 187 188 for (index = 0; index < length1; ++index) 189 arg1[index] = ~arg1[index]; 190 simpleAddHighPrecision (arg1, length1, 1); 191 } 192 193 static uint32_t simpleMultiplyHighPrecision(uint64_t* arg1, int32_t length, uint64_t arg2) { 194 /* assumes arg2 only holds 32 bits of information */ 195 uint64_t product; 196 int32_t index; 197 198 index = 0; 199 product = 0; 200 201 do 202 { 203 product = 204 HIGH_IN_U64 (product) + arg2 * LOW_U32_FROM_PTR (arg1 + index); 205 LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product); 206 product = 207 HIGH_IN_U64 (product) + arg2 * HIGH_U32_FROM_PTR (arg1 + index); 208 HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product); 209 } 210 while (++index < length); 211 212 return HIGH_U32_FROM_VAR (product); 213 } 214 215 static void 216 simpleMultiplyAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2, 217 uint32_t * result) 218 { 219 /* Assumes result can hold the product and arg2 only holds 32 bits 220 of information */ 221 uint64_t product; 222 int32_t index, resultIndex; 223 224 index = resultIndex = 0; 225 product = 0; 226 227 do 228 { 229 product = 230 HIGH_IN_U64 (product) + result[at (resultIndex)] + 231 arg2 * LOW_U32_FROM_PTR (arg1 + index); 232 result[at (resultIndex)] = LOW_U32_FROM_VAR (product); 233 ++resultIndex; 234 product = 235 HIGH_IN_U64 (product) + result[at (resultIndex)] + 236 arg2 * HIGH_U32_FROM_PTR (arg1 + index); 237 result[at (resultIndex)] = LOW_U32_FROM_VAR (product); 238 ++resultIndex; 239 } 240 while (++index < length); 241 242 result[at (resultIndex)] += HIGH_U32_FROM_VAR (product); 243 if (result[at (resultIndex)] < HIGH_U32_FROM_VAR (product)) 244 { 245 /* must be careful with ++ operator and macro expansion */ 246 ++resultIndex; 247 while (++result[at (resultIndex)] == 0) 248 ++resultIndex; 249 } 250 } 251 252 #if __BYTE_ORDER != __LITTLE_ENDIAN 253 void simpleMultiplyAddHighPrecisionBigEndianFix(uint64_t* arg1, int32_t length, uint64_t arg2, uint32_t* result) { 254 /* Assumes result can hold the product and arg2 only holds 32 bits 255 of information */ 256 uint64_t product; 257 int32_t index, resultIndex; 258 259 index = resultIndex = 0; 260 product = 0; 261 262 do { 263 product = HIGH_IN_U64(product) + result[halfAt(resultIndex)] + arg2 * LOW_U32_FROM_PTR(arg1 + index); 264 result[halfAt(resultIndex)] = LOW_U32_FROM_VAR(product); 265 ++resultIndex; 266 product = HIGH_IN_U64(product) + result[halfAt(resultIndex)] + arg2 * HIGH_U32_FROM_PTR(arg1 + index); 267 result[halfAt(resultIndex)] = LOW_U32_FROM_VAR(product); 268 ++resultIndex; 269 } while (++index < length); 270 271 result[halfAt(resultIndex)] += HIGH_U32_FROM_VAR(product); 272 if (result[halfAt(resultIndex)] < HIGH_U32_FROM_VAR(product)) { 273 /* must be careful with ++ operator and macro expansion */ 274 ++resultIndex; 275 while (++result[halfAt(resultIndex)] == 0) ++resultIndex; 276 } 277 } 278 #endif 279 280 void 281 multiplyHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2, 282 uint64_t * result, int32_t length) 283 { 284 /* assumes result is large enough to hold product */ 285 uint64_t* temp; 286 uint32_t* resultIn32; 287 int32_t count, index; 288 289 if (length1 < length2) 290 { 291 temp = arg1; 292 arg1 = arg2; 293 arg2 = temp; 294 count = length1; 295 length1 = length2; 296 length2 = count; 297 } 298 299 memset (result, 0, sizeof (uint64_t) * length); 300 301 /* length1 > length2 */ 302 resultIn32 = reinterpret_cast<uint32_t*>(result); 303 index = -1; 304 for (count = 0; count < length2; ++count) 305 { 306 simpleMultiplyAddHighPrecision (arg1, length1, LOW_IN_U64 (arg2[count]), 307 resultIn32 + (++index)); 308 #if __BYTE_ORDER == __LITTLE_ENDIAN 309 simpleMultiplyAddHighPrecision(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index)); 310 #else 311 simpleMultiplyAddHighPrecisionBigEndianFix(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index)); 312 #endif 313 } 314 } 315 316 uint32_t 317 simpleAppendDecimalDigitHighPrecision (uint64_t * arg1, int32_t length, uint64_t digit) 318 { 319 /* assumes digit is less than 32 bits */ 320 uint64_t arg; 321 int32_t index = 0; 322 323 digit <<= 32; 324 do 325 { 326 arg = LOW_IN_U64 (arg1[index]); 327 digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg); 328 LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit); 329 330 arg = HIGH_IN_U64 (arg1[index]); 331 digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg); 332 HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit); 333 } 334 while (++index < length); 335 336 return HIGH_U32_FROM_VAR (digit); 337 } 338 339 void 340 simpleShiftLeftHighPrecision (uint64_t * arg1, int32_t length, int32_t arg2) 341 { 342 /* assumes length > 0 */ 343 int32_t index, offset; 344 if (arg2 >= 64) 345 { 346 offset = arg2 >> 6; 347 index = length; 348 349 while (--index - offset >= 0) 350 arg1[index] = arg1[index - offset]; 351 do 352 { 353 arg1[index] = 0; 354 } 355 while (--index >= 0); 356 357 arg2 &= 0x3F; 358 } 359 360 if (arg2 == 0) 361 return; 362 while (--length > 0) 363 { 364 arg1[length] = arg1[length] << arg2 | arg1[length - 1] >> (64 - arg2); 365 } 366 *arg1 <<= arg2; 367 } 368 369 int32_t 370 highestSetBit (uint64_t * y) 371 { 372 uint32_t x; 373 int32_t result; 374 375 if (*y == 0) 376 return 0; 377 378 #if defined(USE_LL) 379 if (*y & 0xFFFFFFFF00000000LL) 380 { 381 x = HIGH_U32_FROM_PTR (y); 382 result = 32; 383 } 384 else 385 { 386 x = LOW_U32_FROM_PTR (y); 387 result = 0; 388 } 389 #else 390 #if defined(USE_L) 391 if (*y & 0xFFFFFFFF00000000L) 392 { 393 x = HIGH_U32_FROM_PTR (y); 394 result = 32; 395 } 396 else 397 { 398 x = LOW_U32_FROM_PTR (y); 399 result = 0; 400 } 401 #else 402 if (*y & 0xFFFFFFFF00000000) 403 { 404 x = HIGH_U32_FROM_PTR (y); 405 result = 32; 406 } 407 else 408 { 409 x = LOW_U32_FROM_PTR (y); 410 result = 0; 411 } 412 #endif /* USE_L */ 413 #endif /* USE_LL */ 414 415 if (x & 0xFFFF0000) 416 { 417 x = bitSection (x, 0xFFFF0000, 16); 418 result += 16; 419 } 420 if (x & 0xFF00) 421 { 422 x = bitSection (x, 0xFF00, 8); 423 result += 8; 424 } 425 if (x & 0xF0) 426 { 427 x = bitSection (x, 0xF0, 4); 428 result += 4; 429 } 430 if (x > 0x7) 431 return result + 4; 432 else if (x > 0x3) 433 return result + 3; 434 else if (x > 0x1) 435 return result + 2; 436 else 437 return result + 1; 438 } 439 440 int32_t 441 lowestSetBit (uint64_t * y) 442 { 443 uint32_t x; 444 int32_t result; 445 446 if (*y == 0) 447 return 0; 448 449 #if defined(USE_LL) 450 if (*y & 0x00000000FFFFFFFFLL) 451 { 452 x = LOW_U32_FROM_PTR (y); 453 result = 0; 454 } 455 else 456 { 457 x = HIGH_U32_FROM_PTR (y); 458 result = 32; 459 } 460 #else 461 #if defined(USE_L) 462 if (*y & 0x00000000FFFFFFFFL) 463 { 464 x = LOW_U32_FROM_PTR (y); 465 result = 0; 466 } 467 else 468 { 469 x = HIGH_U32_FROM_PTR (y); 470 result = 32; 471 } 472 #else 473 if (*y & 0x00000000FFFFFFFF) 474 { 475 x = LOW_U32_FROM_PTR (y); 476 result = 0; 477 } 478 else 479 { 480 x = HIGH_U32_FROM_PTR (y); 481 result = 32; 482 } 483 #endif /* USE_L */ 484 #endif /* USE_LL */ 485 486 if (!(x & 0xFFFF)) 487 { 488 x = bitSection (x, 0xFFFF0000, 16); 489 result += 16; 490 } 491 if (!(x & 0xFF)) 492 { 493 x = bitSection (x, 0xFF00, 8); 494 result += 8; 495 } 496 if (!(x & 0xF)) 497 { 498 x = bitSection (x, 0xF0, 4); 499 result += 4; 500 } 501 502 if (x & 0x1) 503 return result + 1; 504 else if (x & 0x2) 505 return result + 2; 506 else if (x & 0x4) 507 return result + 3; 508 else 509 return result + 4; 510 } 511 512 int32_t 513 highestSetBitHighPrecision (uint64_t * arg, int32_t length) 514 { 515 int32_t highBit; 516 517 while (--length >= 0) 518 { 519 highBit = highestSetBit (arg + length); 520 if (highBit) 521 return highBit + 64 * length; 522 } 523 524 return 0; 525 } 526 527 int32_t 528 lowestSetBitHighPrecision (uint64_t * arg, int32_t length) 529 { 530 int32_t lowBit, index = -1; 531 532 while (++index < length) 533 { 534 lowBit = lowestSetBit (arg + index); 535 if (lowBit) 536 return lowBit + 64 * index; 537 } 538 539 return 0; 540 } 541 542 int32_t 543 compareHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2) 544 { 545 while (--length1 >= 0 && arg1[length1] == 0); 546 while (--length2 >= 0 && arg2[length2] == 0); 547 548 if (length1 > length2) 549 return 1; 550 else if (length1 < length2) 551 return -1; 552 else if (length1 > -1) 553 { 554 do 555 { 556 if (arg1[length1] > arg2[length1]) 557 return 1; 558 else if (arg1[length1] < arg2[length1]) 559 return -1; 560 } 561 while (--length1 >= 0); 562 } 563 564 return 0; 565 } 566 567 jdouble 568 toDoubleHighPrecision (uint64_t * arg, int32_t length) 569 { 570 int32_t highBit; 571 uint64_t mantissa, test64; 572 uint32_t test; 573 jdouble result; 574 575 while (length > 0 && arg[length - 1] == 0) 576 --length; 577 578 if (length == 0) 579 result = 0.0; 580 else if (length > 16) 581 { 582 DOUBLE_TO_LONGBITS (result) = EXPONENT_MASK; 583 } 584 else if (length == 1) 585 { 586 highBit = highestSetBit (arg); 587 if (highBit <= 53) 588 { 589 highBit = 53 - highBit; 590 mantissa = *arg << highBit; 591 DOUBLE_TO_LONGBITS (result) = 592 CREATE_DOUBLE_BITS (mantissa, -highBit); 593 } 594 else 595 { 596 highBit -= 53; 597 mantissa = *arg >> highBit; 598 DOUBLE_TO_LONGBITS (result) = 599 CREATE_DOUBLE_BITS (mantissa, highBit); 600 601 /* perform rounding, round to even in case of tie */ 602 test = (LOW_U32_FROM_PTR (arg) << (11 - highBit)) & 0x7FF; 603 if (test > 0x400 || ((test == 0x400) && (mantissa & 1))) 604 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1; 605 } 606 } 607 else 608 { 609 highBit = highestSetBit (arg + (--length)); 610 if (highBit <= 53) 611 { 612 highBit = 53 - highBit; 613 if (highBit > 0) 614 { 615 mantissa = 616 (arg[length] << highBit) | (arg[length - 1] >> 617 (64 - highBit)); 618 } 619 else 620 { 621 mantissa = arg[length]; 622 } 623 DOUBLE_TO_LONGBITS (result) = 624 CREATE_DOUBLE_BITS (mantissa, length * 64 - highBit); 625 626 /* perform rounding, round to even in case of tie */ 627 test64 = arg[--length] << highBit; 628 if (test64 > SIGN_MASK || ((test64 == SIGN_MASK) && (mantissa & 1))) 629 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1; 630 else if (test64 == SIGN_MASK) 631 { 632 while (--length >= 0) 633 { 634 if (arg[length] != 0) 635 { 636 DOUBLE_TO_LONGBITS (result) = 637 DOUBLE_TO_LONGBITS (result) + 1; 638 break; 639 } 640 } 641 } 642 } 643 else 644 { 645 highBit -= 53; 646 mantissa = arg[length] >> highBit; 647 DOUBLE_TO_LONGBITS (result) = 648 CREATE_DOUBLE_BITS (mantissa, length * 64 + highBit); 649 650 /* perform rounding, round to even in case of tie */ 651 test = (LOW_U32_FROM_PTR (arg + length) << (11 - highBit)) & 0x7FF; 652 if (test > 0x400 || ((test == 0x400) && (mantissa & 1))) 653 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1; 654 else if (test == 0x400) 655 { 656 do 657 { 658 if (arg[--length] != 0) 659 { 660 DOUBLE_TO_LONGBITS (result) = 661 DOUBLE_TO_LONGBITS (result) + 1; 662 break; 663 } 664 } 665 while (length > 0); 666 } 667 } 668 } 669 670 return result; 671 } 672 673 static uint64_t simpleMultiplyHighPrecision64(uint64_t* arg1, int32_t length, uint64_t arg2); 674 675 int32_t 676 timesTenToTheEHighPrecision (uint64_t * result, int32_t length, jint e) 677 { 678 /* assumes result can hold value */ 679 uint64_t overflow; 680 int exp10 = e; 681 682 if (e == 0) 683 return length; 684 685 /* bad O(n) way of doing it, but simple */ 686 /* 687 do { 688 overflow = simpleAppendDecimalDigitHighPrecision(result, length, 0); 689 if (overflow) 690 result[length++] = overflow; 691 } while (--e); 692 */ 693 /* Replace the current implementaion which performs a 694 * "multiplication" by 10 e number of times with an actual 695 * multiplication. 10e19 is the largest exponent to the power of ten 696 * that will fit in a 64-bit integer, and 10e9 is the largest exponent to 697 * the power of ten that will fit in a 64-bit integer. Not sure where the 698 * break-even point is between an actual multiplication and a 699 * simpleAappendDecimalDigit() so just pick 10e3 as that point for 700 * now. 701 */ 702 while (exp10 >= 19) 703 { 704 overflow = simpleMultiplyHighPrecision64 (result, length, TEN_E19); 705 if (overflow) 706 result[length++] = overflow; 707 exp10 -= 19; 708 } 709 while (exp10 >= 9) 710 { 711 overflow = simpleMultiplyHighPrecision (result, length, TEN_E9); 712 if (overflow) 713 result[length++] = overflow; 714 exp10 -= 9; 715 } 716 if (exp10 == 0) 717 return length; 718 else if (exp10 == 1) 719 { 720 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0); 721 if (overflow) 722 result[length++] = overflow; 723 } 724 else if (exp10 == 2) 725 { 726 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0); 727 if (overflow) 728 result[length++] = overflow; 729 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0); 730 if (overflow) 731 result[length++] = overflow; 732 } 733 else if (exp10 == 3) 734 { 735 overflow = simpleMultiplyHighPrecision (result, length, TEN_E3); 736 if (overflow) 737 result[length++] = overflow; 738 } 739 else if (exp10 == 4) 740 { 741 overflow = simpleMultiplyHighPrecision (result, length, TEN_E4); 742 if (overflow) 743 result[length++] = overflow; 744 } 745 else if (exp10 == 5) 746 { 747 overflow = simpleMultiplyHighPrecision (result, length, TEN_E5); 748 if (overflow) 749 result[length++] = overflow; 750 } 751 else if (exp10 == 6) 752 { 753 overflow = simpleMultiplyHighPrecision (result, length, TEN_E6); 754 if (overflow) 755 result[length++] = overflow; 756 } 757 else if (exp10 == 7) 758 { 759 overflow = simpleMultiplyHighPrecision (result, length, TEN_E7); 760 if (overflow) 761 result[length++] = overflow; 762 } 763 else if (exp10 == 8) 764 { 765 overflow = simpleMultiplyHighPrecision (result, length, TEN_E8); 766 if (overflow) 767 result[length++] = overflow; 768 } 769 return length; 770 } 771 772 uint64_t 773 doubleMantissa (jdouble z) 774 { 775 uint64_t m = DOUBLE_TO_LONGBITS (z); 776 777 if ((m & EXPONENT_MASK) != 0) 778 m = (m & MANTISSA_MASK) | NORMAL_MASK; 779 else 780 m = (m & MANTISSA_MASK); 781 782 return m; 783 } 784 785 int32_t 786 doubleExponent (jdouble z) 787 { 788 /* assumes positive double */ 789 int32_t k = HIGH_U32_FROM_VAR (z) >> 20; 790 791 if (k) 792 k -= E_OFFSET; 793 else 794 k = 1 - E_OFFSET; 795 796 return k; 797 } 798 799 uint32_t floatMantissa(jfloat z) { 800 uint32_t m = FLOAT_TO_INTBITS (z); 801 802 if ((m & FLOAT_EXPONENT_MASK) != 0) 803 m = (m & FLOAT_MANTISSA_MASK) | FLOAT_NORMAL_MASK; 804 else 805 m = (m & FLOAT_MANTISSA_MASK); 806 807 return m; 808 } 809 810 int32_t 811 floatExponent (jfloat z) 812 { 813 /* assumes positive float */ 814 int32_t k = FLOAT_TO_INTBITS (z) >> 23; 815 if (k) 816 k -= FLOAT_E_OFFSET; 817 else 818 k = 1 - FLOAT_E_OFFSET; 819 820 return k; 821 } 822 823 /* Allow a 64-bit value in arg2 */ 824 uint64_t 825 simpleMultiplyHighPrecision64 (uint64_t * arg1, int32_t length, uint64_t arg2) 826 { 827 uint64_t intermediate, carry1, carry2, prod1, prod2, sum; 828 uint64_t* pArg1; 829 int32_t index; 830 uint32_t buf32; 831 832 index = 0; 833 intermediate = 0; 834 pArg1 = arg1 + index; 835 carry1 = carry2 = 0; 836 837 do 838 { 839 if ((*pArg1 != 0) || (intermediate != 0)) 840 { 841 prod1 = 842 static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1)); 843 sum = intermediate + prod1; 844 if ((sum < prod1) || (sum < intermediate)) 845 { 846 carry1 = 1; 847 } 848 else 849 { 850 carry1 = 0; 851 } 852 prod1 = 853 static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(HIGH_U32_FROM_PTR (pArg1)); 854 prod2 = 855 static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1)); 856 intermediate = carry2 + HIGH_IN_U64 (sum) + prod1 + prod2; 857 if ((intermediate < prod1) || (intermediate < prod2)) 858 { 859 carry2 = 1; 860 } 861 else 862 { 863 carry2 = 0; 864 } 865 LOW_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (sum); 866 buf32 = HIGH_U32_FROM_PTR (pArg1); 867 HIGH_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (intermediate); 868 intermediate = carry1 + HIGH_IN_U64 (intermediate) 869 + static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(buf32); 870 } 871 pArg1++; 872 } 873 while (++index < length); 874 return intermediate; 875 } 876
