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