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 #define HIGH_IN_U64(u64) ((u64) >> 32) 22 #define LOW_IN_U64(u64) ((u64) & 0x00000000FFFFFFFFLL) 23 24 #define TEN_E1 (0xALL) 25 #define TEN_E2 (0x64LL) 26 #define TEN_E3 (0x3E8LL) 27 #define TEN_E4 (0x2710LL) 28 #define TEN_E5 (0x186A0LL) 29 #define TEN_E6 (0xF4240LL) 30 #define TEN_E7 (0x989680LL) 31 #define TEN_E8 (0x5F5E100LL) 32 #define TEN_E9 (0x3B9ACA00LL) 33 #define TEN_E19 (0x8AC7230489E80000LL) 34 35 #define TIMES_TEN(x) (((x) << 3) + ((x) << 1)) 36 #define bitSection(x, mask, shift) (((x) & (mask)) >> (shift)) 37 #define CREATE_DOUBLE_BITS(normalizedM, e) (((normalizedM) & MANTISSA_MASK) | ((static_cast<uint64_t>((e) + E_OFFSET)) << 52)) 38 39 #define MANTISSA_MASK (0x000FFFFFFFFFFFFFLL) 40 #define EXPONENT_MASK (0x7FF0000000000000LL) 41 #define NORMAL_MASK (0x0010000000000000LL) 42 #define SIGN_MASK (0x8000000000000000LL) 43 44 #define E_OFFSET (1075) 45 46 #define FLOAT_MANTISSA_MASK (0x007FFFFF) 47 #define FLOAT_EXPONENT_MASK (0x7F800000) 48 #define FLOAT_NORMAL_MASK (0x00800000) 49 #define FLOAT_E_OFFSET (150) 50 51 int32_t 52 simpleAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2) 53 { 54 /* assumes length > 0 */ 55 int32_t index = 1; 56 57 *arg1 += arg2; 58 if (arg2 <= *arg1) 59 return 0; 60 else if (length == 1) 61 return 1; 62 63 while (++arg1[index] == 0 && ++index < length) { 64 } 65 66 return index == length; 67 } 68 69 int32_t 70 addHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2) 71 { 72 /* addition is limited by length of arg1 as it this function is 73 * storing the result in arg1 */ 74 /* fix for cc (GCC) 3.2 20020903 (Red Hat Linux 8.0 3.2-7): code generated does not 75 * do the temp1 + temp2 + carry addition correct. carry is 64 bit because gcc has 76 * subtle issues when you mix 64 / 32 bit maths. */ 77 uint64_t temp1, temp2, temp3; /* temporary variables to help the SH-4, and gcc */ 78 uint64_t carry; 79 int32_t index; 80 81 if (length1 == 0 || length2 == 0) 82 { 83 return 0; 84 } 85 else if (length1 < length2) 86 { 87 length2 = length1; 88 } 89 90 carry = 0; 91 index = 0; 92 do 93 { 94 temp1 = arg1[index]; 95 temp2 = arg2[index]; 96 temp3 = temp1 + temp2; 97 arg1[index] = temp3 + carry; 98 if (arg2[index] < arg1[index]) 99 carry = 0; 100 else if (arg2[index] != arg1[index]) 101 carry = 1; 102 } 103 while (++index < length2); 104 if (!carry) 105 return 0; 106 else if (index == length1) 107 return 1; 108 109 while (++arg1[index] == 0 && ++index < length1) { 110 } 111 112 return index == length1; 113 } 114 115 void 116 subtractHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2) 117 { 118 /* assumes arg1 > arg2 */ 119 int32_t index; 120 for (index = 0; index < length1; ++index) 121 arg1[index] = ~arg1[index]; 122 simpleAddHighPrecision (arg1, length1, 1); 123 124 while (length2 > 0 && arg2[length2 - 1] == 0) 125 --length2; 126 127 addHighPrecision (arg1, length1, arg2, length2); 128 129 for (index = 0; index < length1; ++index) 130 arg1[index] = ~arg1[index]; 131 simpleAddHighPrecision (arg1, length1, 1); 132 } 133 134 static uint32_t simpleMultiplyHighPrecision(uint64_t* arg1, int32_t length, uint64_t arg2) { 135 /* assumes arg2 only holds 32 bits of information */ 136 uint64_t product; 137 int32_t index; 138 139 index = 0; 140 product = 0; 141 142 do 143 { 144 product = 145 HIGH_IN_U64 (product) + arg2 * LOW_U32_FROM_PTR (arg1 + index); 146 LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product); 147 product = 148 HIGH_IN_U64 (product) + arg2 * HIGH_U32_FROM_PTR (arg1 + index); 149 HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (product); 150 } 151 while (++index < length); 152 153 return HIGH_U32_FROM_VAR (product); 154 } 155 156 static void 157 simpleMultiplyAddHighPrecision (uint64_t * arg1, int32_t length, uint64_t arg2, 158 uint32_t * result) 159 { 160 /* Assumes result can hold the product and arg2 only holds 32 bits 161 of information */ 162 uint64_t product; 163 int32_t index, resultIndex; 164 165 index = resultIndex = 0; 166 product = 0; 167 168 do 169 { 170 product = 171 HIGH_IN_U64 (product) + result[resultIndex] + 172 arg2 * LOW_U32_FROM_PTR (arg1 + index); 173 result[resultIndex] = LOW_U32_FROM_VAR (product); 174 ++resultIndex; 175 product = 176 HIGH_IN_U64 (product) + result[resultIndex] + 177 arg2 * HIGH_U32_FROM_PTR (arg1 + index); 178 result[resultIndex] = LOW_U32_FROM_VAR (product); 179 ++resultIndex; 180 } 181 while (++index < length); 182 183 result[resultIndex] += HIGH_U32_FROM_VAR (product); 184 if (result[resultIndex] < HIGH_U32_FROM_VAR (product)) 185 { 186 /* must be careful with ++ operator and macro expansion */ 187 ++resultIndex; 188 while (++result[resultIndex] == 0) 189 ++resultIndex; 190 } 191 } 192 193 void 194 multiplyHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2, 195 uint64_t * result, int32_t length) 196 { 197 /* assumes result is large enough to hold product */ 198 uint64_t* temp; 199 uint32_t* resultIn32; 200 int32_t count, index; 201 202 if (length1 < length2) 203 { 204 temp = arg1; 205 arg1 = arg2; 206 arg2 = temp; 207 count = length1; 208 length1 = length2; 209 length2 = count; 210 } 211 212 memset (result, 0, sizeof (uint64_t) * length); 213 214 /* length1 > length2 */ 215 resultIn32 = reinterpret_cast<uint32_t*>(result); 216 index = -1; 217 for (count = 0; count < length2; ++count) 218 { 219 simpleMultiplyAddHighPrecision (arg1, length1, LOW_IN_U64 (arg2[count]), 220 resultIn32 + (++index)); 221 simpleMultiplyAddHighPrecision(arg1, length1, HIGH_IN_U64(arg2[count]), resultIn32 + (++index)); 222 } 223 } 224 225 uint32_t 226 simpleAppendDecimalDigitHighPrecision (uint64_t * arg1, int32_t length, uint64_t digit) 227 { 228 /* assumes digit is less than 32 bits */ 229 uint64_t arg; 230 int32_t index = 0; 231 232 digit <<= 32; 233 do 234 { 235 arg = LOW_IN_U64 (arg1[index]); 236 digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg); 237 LOW_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit); 238 239 arg = HIGH_IN_U64 (arg1[index]); 240 digit = HIGH_IN_U64 (digit) + TIMES_TEN (arg); 241 HIGH_U32_FROM_PTR (arg1 + index) = LOW_U32_FROM_VAR (digit); 242 } 243 while (++index < length); 244 245 return HIGH_U32_FROM_VAR (digit); 246 } 247 248 void 249 simpleShiftLeftHighPrecision (uint64_t * arg1, int32_t length, int32_t arg2) 250 { 251 /* assumes length > 0 */ 252 int32_t index, offset; 253 if (arg2 >= 64) 254 { 255 offset = arg2 >> 6; 256 index = length; 257 258 while (--index - offset >= 0) 259 arg1[index] = arg1[index - offset]; 260 do 261 { 262 arg1[index] = 0; 263 } 264 while (--index >= 0); 265 266 arg2 &= 0x3F; 267 } 268 269 if (arg2 == 0) 270 return; 271 while (--length > 0) 272 { 273 arg1[length] = arg1[length] << arg2 | arg1[length - 1] >> (64 - arg2); 274 } 275 *arg1 <<= arg2; 276 } 277 278 int32_t 279 highestSetBit (uint64_t * y) 280 { 281 uint32_t x; 282 int32_t result; 283 284 if (*y == 0) 285 return 0; 286 287 if (*y & 0xFFFFFFFF00000000LL) 288 { 289 x = HIGH_U32_FROM_PTR (y); 290 result = 32; 291 } 292 else 293 { 294 x = LOW_U32_FROM_PTR (y); 295 result = 0; 296 } 297 298 if (x & 0xFFFF0000) 299 { 300 x = bitSection (x, 0xFFFF0000, 16); 301 result += 16; 302 } 303 if (x & 0xFF00) 304 { 305 x = bitSection (x, 0xFF00, 8); 306 result += 8; 307 } 308 if (x & 0xF0) 309 { 310 x = bitSection (x, 0xF0, 4); 311 result += 4; 312 } 313 if (x > 0x7) 314 return result + 4; 315 else if (x > 0x3) 316 return result + 3; 317 else if (x > 0x1) 318 return result + 2; 319 else 320 return result + 1; 321 } 322 323 int32_t 324 lowestSetBit (uint64_t * y) 325 { 326 uint32_t x; 327 int32_t result; 328 329 if (*y == 0) 330 return 0; 331 332 if (*y & 0x00000000FFFFFFFFLL) 333 { 334 x = LOW_U32_FROM_PTR (y); 335 result = 0; 336 } 337 else 338 { 339 x = HIGH_U32_FROM_PTR (y); 340 result = 32; 341 } 342 343 if (!(x & 0xFFFF)) 344 { 345 x = bitSection (x, 0xFFFF0000, 16); 346 result += 16; 347 } 348 if (!(x & 0xFF)) 349 { 350 x = bitSection (x, 0xFF00, 8); 351 result += 8; 352 } 353 if (!(x & 0xF)) 354 { 355 x = bitSection (x, 0xF0, 4); 356 result += 4; 357 } 358 359 if (x & 0x1) 360 return result + 1; 361 else if (x & 0x2) 362 return result + 2; 363 else if (x & 0x4) 364 return result + 3; 365 else 366 return result + 4; 367 } 368 369 int32_t 370 highestSetBitHighPrecision (uint64_t * arg, int32_t length) 371 { 372 int32_t highBit; 373 374 while (--length >= 0) 375 { 376 highBit = highestSetBit (arg + length); 377 if (highBit) 378 return highBit + 64 * length; 379 } 380 381 return 0; 382 } 383 384 int32_t 385 lowestSetBitHighPrecision (uint64_t * arg, int32_t length) 386 { 387 int32_t lowBit, index = -1; 388 389 while (++index < length) 390 { 391 lowBit = lowestSetBit (arg + index); 392 if (lowBit) 393 return lowBit + 64 * index; 394 } 395 396 return 0; 397 } 398 399 int32_t 400 compareHighPrecision (uint64_t * arg1, int32_t length1, uint64_t * arg2, int32_t length2) 401 { 402 while (--length1 >= 0 && arg1[length1] == 0) { 403 } 404 while (--length2 >= 0 && arg2[length2] == 0) { 405 } 406 407 if (length1 > length2) 408 return 1; 409 else if (length1 < length2) 410 return -1; 411 else if (length1 > -1) 412 { 413 do 414 { 415 if (arg1[length1] > arg2[length1]) 416 return 1; 417 else if (arg1[length1] < arg2[length1]) 418 return -1; 419 } 420 while (--length1 >= 0); 421 } 422 423 return 0; 424 } 425 426 jdouble 427 toDoubleHighPrecision (uint64_t * arg, int32_t length) 428 { 429 int32_t highBit; 430 uint64_t mantissa, test64; 431 uint32_t test; 432 jdouble result; 433 434 while (length > 0 && arg[length - 1] == 0) 435 --length; 436 437 if (length == 0) 438 result = 0.0; 439 else if (length > 16) 440 { 441 DOUBLE_TO_LONGBITS (result) = EXPONENT_MASK; 442 } 443 else if (length == 1) 444 { 445 highBit = highestSetBit (arg); 446 if (highBit <= 53) 447 { 448 highBit = 53 - highBit; 449 mantissa = *arg << highBit; 450 DOUBLE_TO_LONGBITS (result) = 451 CREATE_DOUBLE_BITS (mantissa, -highBit); 452 } 453 else 454 { 455 highBit -= 53; 456 mantissa = *arg >> highBit; 457 DOUBLE_TO_LONGBITS (result) = 458 CREATE_DOUBLE_BITS (mantissa, highBit); 459 460 /* perform rounding, round to even in case of tie */ 461 test = (LOW_U32_FROM_PTR (arg) << (11 - highBit)) & 0x7FF; 462 if (test > 0x400 || ((test == 0x400) && (mantissa & 1))) 463 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1; 464 } 465 } 466 else 467 { 468 highBit = highestSetBit (arg + (--length)); 469 if (highBit <= 53) 470 { 471 highBit = 53 - highBit; 472 if (highBit > 0) 473 { 474 mantissa = 475 (arg[length] << highBit) | (arg[length - 1] >> 476 (64 - highBit)); 477 } 478 else 479 { 480 mantissa = arg[length]; 481 } 482 DOUBLE_TO_LONGBITS (result) = 483 CREATE_DOUBLE_BITS (mantissa, length * 64 - highBit); 484 485 /* perform rounding, round to even in case of tie */ 486 test64 = arg[--length] << highBit; 487 if (test64 > SIGN_MASK || ((test64 == SIGN_MASK) && (mantissa & 1))) 488 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1; 489 else if (test64 == SIGN_MASK) 490 { 491 while (--length >= 0) 492 { 493 if (arg[length] != 0) 494 { 495 DOUBLE_TO_LONGBITS (result) = 496 DOUBLE_TO_LONGBITS (result) + 1; 497 break; 498 } 499 } 500 } 501 } 502 else 503 { 504 highBit -= 53; 505 mantissa = arg[length] >> highBit; 506 DOUBLE_TO_LONGBITS (result) = 507 CREATE_DOUBLE_BITS (mantissa, length * 64 + highBit); 508 509 /* perform rounding, round to even in case of tie */ 510 test = (LOW_U32_FROM_PTR (arg + length) << (11 - highBit)) & 0x7FF; 511 if (test > 0x400 || ((test == 0x400) && (mantissa & 1))) 512 DOUBLE_TO_LONGBITS (result) = DOUBLE_TO_LONGBITS (result) + 1; 513 else if (test == 0x400) 514 { 515 do 516 { 517 if (arg[--length] != 0) 518 { 519 DOUBLE_TO_LONGBITS (result) = 520 DOUBLE_TO_LONGBITS (result) + 1; 521 break; 522 } 523 } 524 while (length > 0); 525 } 526 } 527 } 528 529 return result; 530 } 531 532 static uint64_t simpleMultiplyHighPrecision64(uint64_t* arg1, int32_t length, uint64_t arg2); 533 534 int32_t 535 timesTenToTheEHighPrecision (uint64_t * result, int32_t length, jint e) 536 { 537 /* assumes result can hold value */ 538 uint64_t overflow; 539 int exp10 = e; 540 541 if (e == 0) 542 return length; 543 544 /* bad O(n) way of doing it, but simple */ 545 /* 546 do { 547 overflow = simpleAppendDecimalDigitHighPrecision(result, length, 0); 548 if (overflow) 549 result[length++] = overflow; 550 } while (--e); 551 */ 552 /* Replace the current implementaion which performs a 553 * "multiplication" by 10 e number of times with an actual 554 * multiplication. 10e19 is the largest exponent to the power of ten 555 * that will fit in a 64-bit integer, and 10e9 is the largest exponent to 556 * the power of ten that will fit in a 64-bit integer. Not sure where the 557 * break-even point is between an actual multiplication and a 558 * simpleAappendDecimalDigit() so just pick 10e3 as that point for 559 * now. 560 */ 561 while (exp10 >= 19) 562 { 563 overflow = simpleMultiplyHighPrecision64 (result, length, TEN_E19); 564 if (overflow) 565 result[length++] = overflow; 566 exp10 -= 19; 567 } 568 while (exp10 >= 9) 569 { 570 overflow = simpleMultiplyHighPrecision (result, length, TEN_E9); 571 if (overflow) 572 result[length++] = overflow; 573 exp10 -= 9; 574 } 575 if (exp10 == 0) 576 return length; 577 else if (exp10 == 1) 578 { 579 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0); 580 if (overflow) 581 result[length++] = overflow; 582 } 583 else if (exp10 == 2) 584 { 585 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0); 586 if (overflow) 587 result[length++] = overflow; 588 overflow = simpleAppendDecimalDigitHighPrecision (result, length, 0); 589 if (overflow) 590 result[length++] = overflow; 591 } 592 else if (exp10 == 3) 593 { 594 overflow = simpleMultiplyHighPrecision (result, length, TEN_E3); 595 if (overflow) 596 result[length++] = overflow; 597 } 598 else if (exp10 == 4) 599 { 600 overflow = simpleMultiplyHighPrecision (result, length, TEN_E4); 601 if (overflow) 602 result[length++] = overflow; 603 } 604 else if (exp10 == 5) 605 { 606 overflow = simpleMultiplyHighPrecision (result, length, TEN_E5); 607 if (overflow) 608 result[length++] = overflow; 609 } 610 else if (exp10 == 6) 611 { 612 overflow = simpleMultiplyHighPrecision (result, length, TEN_E6); 613 if (overflow) 614 result[length++] = overflow; 615 } 616 else if (exp10 == 7) 617 { 618 overflow = simpleMultiplyHighPrecision (result, length, TEN_E7); 619 if (overflow) 620 result[length++] = overflow; 621 } 622 else if (exp10 == 8) 623 { 624 overflow = simpleMultiplyHighPrecision (result, length, TEN_E8); 625 if (overflow) 626 result[length++] = overflow; 627 } 628 return length; 629 } 630 631 uint64_t 632 doubleMantissa (jdouble z) 633 { 634 uint64_t m = DOUBLE_TO_LONGBITS (z); 635 636 if ((m & EXPONENT_MASK) != 0) 637 m = (m & MANTISSA_MASK) | NORMAL_MASK; 638 else 639 m = (m & MANTISSA_MASK); 640 641 return m; 642 } 643 644 int32_t 645 doubleExponent (jdouble z) 646 { 647 /* assumes positive double */ 648 int32_t k = HIGH_U32_FROM_VAR (z) >> 20; 649 650 if (k) 651 k -= E_OFFSET; 652 else 653 k = 1 - E_OFFSET; 654 655 return k; 656 } 657 658 uint32_t floatMantissa(jfloat z) { 659 uint32_t m = FLOAT_TO_INTBITS (z); 660 661 if ((m & FLOAT_EXPONENT_MASK) != 0) 662 m = (m & FLOAT_MANTISSA_MASK) | FLOAT_NORMAL_MASK; 663 else 664 m = (m & FLOAT_MANTISSA_MASK); 665 666 return m; 667 } 668 669 int32_t 670 floatExponent (jfloat z) 671 { 672 /* assumes positive float */ 673 int32_t k = FLOAT_TO_INTBITS (z) >> 23; 674 if (k) 675 k -= FLOAT_E_OFFSET; 676 else 677 k = 1 - FLOAT_E_OFFSET; 678 679 return k; 680 } 681 682 /* Allow a 64-bit value in arg2 */ 683 uint64_t 684 simpleMultiplyHighPrecision64 (uint64_t * arg1, int32_t length, uint64_t arg2) 685 { 686 uint64_t intermediate, carry1, carry2, prod1, prod2, sum; 687 uint64_t* pArg1; 688 int32_t index; 689 uint32_t buf32; 690 691 index = 0; 692 intermediate = 0; 693 pArg1 = arg1 + index; 694 carry1 = carry2 = 0; 695 696 do 697 { 698 if ((*pArg1 != 0) || (intermediate != 0)) 699 { 700 prod1 = 701 static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1)); 702 sum = intermediate + prod1; 703 if ((sum < prod1) || (sum < intermediate)) 704 { 705 carry1 = 1; 706 } 707 else 708 { 709 carry1 = 0; 710 } 711 prod1 = 712 static_cast<uint64_t>(LOW_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(HIGH_U32_FROM_PTR (pArg1)); 713 prod2 = 714 static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(LOW_U32_FROM_PTR (pArg1)); 715 intermediate = carry2 + HIGH_IN_U64 (sum) + prod1 + prod2; 716 if ((intermediate < prod1) || (intermediate < prod2)) 717 { 718 carry2 = 1; 719 } 720 else 721 { 722 carry2 = 0; 723 } 724 LOW_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (sum); 725 buf32 = HIGH_U32_FROM_PTR (pArg1); 726 HIGH_U32_FROM_PTR (pArg1) = LOW_U32_FROM_VAR (intermediate); 727 intermediate = carry1 + HIGH_IN_U64 (intermediate) 728 + static_cast<uint64_t>(HIGH_U32_FROM_VAR (arg2)) * static_cast<uint64_t>(buf32); 729 } 730 pArg1++; 731 } 732 while (++index < length); 733 return intermediate; 734 } 735
