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