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 <stdlib.h> 19 #include <string.h> 20 #include <math.h> 21 #include "JNIHelp.h" 22 #include "JniConstants.h" 23 #include "JniException.h" 24 #include "ScopedUtfChars.h" 25 #include "cbigint.h" 26 27 /* ************************* Defines ************************* */ 28 #if defined(__linux__) || defined(FREEBSD) 29 #define USE_LL 30 #endif 31 32 #define LOW_I32_FROM_VAR(u64) LOW_I32_FROM_LONG64(u64) 33 #define LOW_I32_FROM_PTR(u64ptr) LOW_I32_FROM_LONG64_PTR(u64ptr) 34 #define HIGH_I32_FROM_VAR(u64) HIGH_I32_FROM_LONG64(u64) 35 #define HIGH_I32_FROM_PTR(u64ptr) HIGH_I32_FROM_LONG64_PTR(u64ptr) 36 37 #define MAX_DOUBLE_ACCURACY_WIDTH 17 38 39 #define DEFAULT_DOUBLE_WIDTH MAX_DOUBLE_ACCURACY_WIDTH 40 41 #if defined(USE_LL) 42 #define DOUBLE_INFINITE_LONGBITS (0x7FF0000000000000LL) 43 #else 44 #if defined(USE_L) 45 #define DOUBLE_INFINITE_LONGBITS (0x7FF0000000000000L) 46 #else 47 #define DOUBLE_INFINITE_LONGBITS (0x7FF0000000000000) 48 #endif /* USE_L */ 49 #endif /* USE_LL */ 50 51 #define DOUBLE_MINIMUM_LONGBITS (0x1) 52 53 #if defined(USE_LL) 54 #define DOUBLE_MANTISSA_MASK (0x000FFFFFFFFFFFFFLL) 55 #define DOUBLE_EXPONENT_MASK (0x7FF0000000000000LL) 56 #define DOUBLE_NORMAL_MASK (0x0010000000000000LL) 57 #else 58 #if defined(USE_L) 59 #define DOUBLE_MANTISSA_MASK (0x000FFFFFFFFFFFFFL) 60 #define DOUBLE_EXPONENT_MASK (0x7FF0000000000000L) 61 #define DOUBLE_NORMAL_MASK (0x0010000000000000L) 62 #else 63 #define DOUBLE_MANTISSA_MASK (0x000FFFFFFFFFFFFF) 64 #define DOUBLE_EXPONENT_MASK (0x7FF0000000000000) 65 #define DOUBLE_NORMAL_MASK (0x0010000000000000) 66 #endif /* USE_L */ 67 #endif /* USE_LL */ 68 69 /* Keep a count of the number of times we decrement and increment to 70 * approximate the double, and attempt to detect the case where we 71 * could potentially toggle back and forth between decrementing and 72 * incrementing. It is possible for us to be stuck in the loop when 73 * incrementing by one or decrementing by one may exceed or stay below 74 * the value that we are looking for. In this case, just break out of 75 * the loop if we toggle between incrementing and decrementing for more 76 * than twice. 77 */ 78 #define INCREMENT_DOUBLE(_x, _decCount, _incCount) \ 79 { \ 80 ++DOUBLE_TO_LONGBITS(_x); \ 81 _incCount++; \ 82 if( (_incCount > 2) && (_decCount > 2) ) { \ 83 if( _decCount > _incCount ) { \ 84 DOUBLE_TO_LONGBITS(_x) += _decCount - _incCount; \ 85 } else if( _incCount > _decCount ) { \ 86 DOUBLE_TO_LONGBITS(_x) -= _incCount - _decCount; \ 87 } \ 88 break; \ 89 } \ 90 } 91 #define DECREMENT_DOUBLE(_x, _decCount, _incCount) \ 92 { \ 93 --DOUBLE_TO_LONGBITS(_x); \ 94 _decCount++; \ 95 if( (_incCount > 2) && (_decCount > 2) ) { \ 96 if( _decCount > _incCount ) { \ 97 DOUBLE_TO_LONGBITS(_x) += _decCount - _incCount; \ 98 } else if( _incCount > _decCount ) { \ 99 DOUBLE_TO_LONGBITS(_x) -= _incCount - _decCount; \ 100 } \ 101 break; \ 102 } \ 103 } 104 105 #define allocateU64(x, n) if (!((x) = reinterpret_cast<uint64_t*>(malloc((n) * sizeof(uint64_t))))) goto OutOfMemory; 106 107 /* *********************************************************** */ 108 109 /* ************************ local data ************************ */ 110 static const jdouble double_tens[] = { 111 1.0, 112 1.0e1, 113 1.0e2, 114 1.0e3, 115 1.0e4, 116 1.0e5, 117 1.0e6, 118 1.0e7, 119 1.0e8, 120 1.0e9, 121 1.0e10, 122 1.0e11, 123 1.0e12, 124 1.0e13, 125 1.0e14, 126 1.0e15, 127 1.0e16, 128 1.0e17, 129 1.0e18, 130 1.0e19, 131 1.0e20, 132 1.0e21, 133 1.0e22 134 }; 135 /* *********************************************************** */ 136 137 /* ************** private function declarations ************** */ 138 static jdouble createDouble1 (JNIEnv* env, uint64_t * f, int32_t length, jint e); 139 static jdouble doubleAlgorithm (JNIEnv* env, uint64_t * f, int32_t length, jint e, 140 jdouble z); 141 /* *********************************************************** */ 142 143 #define doubleTenToTheE(e) (*(double_tens + (e))) 144 #define DOUBLE_LOG5_OF_TWO_TO_THE_N 23 145 146 #define sizeOfTenToTheE(e) (((e) / 19) + 1) 147 148 static jdouble createDouble(JNIEnv* env, const char* s, jint e) { 149 /* assumes s is a null terminated string with at least one 150 * character in it */ 151 uint64_t def[DEFAULT_DOUBLE_WIDTH]; 152 uint64_t defBackup[DEFAULT_DOUBLE_WIDTH]; 153 uint64_t* f; 154 uint64_t* fNoOverflow; 155 uint64_t* g; 156 uint64_t* tempBackup; 157 uint32_t overflow; 158 jdouble result; 159 int32_t index = 1; 160 int unprocessedDigits = 0; 161 162 f = def; 163 fNoOverflow = defBackup; 164 *f = 0; 165 tempBackup = g = 0; 166 do 167 { 168 if (*s >= '0' && *s <= '9') 169 { 170 /* Make a back up of f before appending, so that we can 171 * back out of it if there is no more room, i.e. index > 172 * MAX_DOUBLE_ACCURACY_WIDTH. 173 */ 174 memcpy (fNoOverflow, f, sizeof (uint64_t) * index); 175 overflow = 176 simpleAppendDecimalDigitHighPrecision (f, index, *s - '0'); 177 if (overflow) 178 { 179 f[index++] = overflow; 180 /* There is an overflow, but there is no more room 181 * to store the result. We really only need the top 52 182 * bits anyway, so we must back out of the overflow, 183 * and ignore the rest of the string. 184 */ 185 if (index >= MAX_DOUBLE_ACCURACY_WIDTH) 186 { 187 index--; 188 memcpy (f, fNoOverflow, sizeof (uint64_t) * index); 189 break; 190 } 191 if (tempBackup) 192 { 193 fNoOverflow = tempBackup; 194 } 195 } 196 } 197 else 198 index = -1; 199 } 200 while (index > 0 && *(++s) != '\0'); 201 202 /* We've broken out of the parse loop either because we've reached 203 * the end of the string or we've overflowed the maximum accuracy 204 * limit of a double. If we still have unprocessed digits in the 205 * given string, then there are three possible results: 206 * 1. (unprocessed digits + e) == 0, in which case we simply 207 * convert the existing bits that are already parsed 208 * 2. (unprocessed digits + e) < 0, in which case we simply 209 * convert the existing bits that are already parsed along 210 * with the given e 211 * 3. (unprocessed digits + e) > 0 indicates that the value is 212 * simply too big to be stored as a double, so return Infinity 213 */ 214 if ((unprocessedDigits = strlen (s)) > 0) 215 { 216 e += unprocessedDigits; 217 if (index > -1) 218 { 219 if (e == 0) 220 result = toDoubleHighPrecision (f, index); 221 else if (e < 0) 222 result = createDouble1 (env, f, index, e); 223 else 224 { 225 DOUBLE_TO_LONGBITS (result) = DOUBLE_INFINITE_LONGBITS; 226 } 227 } 228 else 229 { 230 LOW_I32_FROM_VAR (result) = -1; 231 HIGH_I32_FROM_VAR (result) = -1; 232 } 233 } 234 else 235 { 236 if (index > -1) 237 { 238 if (e == 0) 239 result = toDoubleHighPrecision (f, index); 240 else 241 result = createDouble1 (env, f, index, e); 242 } 243 else 244 { 245 LOW_I32_FROM_VAR (result) = -1; 246 HIGH_I32_FROM_VAR (result) = -1; 247 } 248 } 249 250 return result; 251 252 } 253 254 static jdouble createDouble1(JNIEnv* env, uint64_t* f, int32_t length, jint e) { 255 int32_t numBits; 256 jdouble result; 257 258 static const jint APPROX_MIN_MAGNITUDE = -309; 259 static const jint APPROX_MAX_MAGNITUDE = 309; 260 261 numBits = highestSetBitHighPrecision (f, length) + 1; 262 numBits -= lowestSetBitHighPrecision (f, length); 263 if (numBits < 54 && e >= 0 && e < DOUBLE_LOG5_OF_TWO_TO_THE_N) 264 { 265 return toDoubleHighPrecision (f, length) * doubleTenToTheE (e); 266 } 267 else if (numBits < 54 && e < 0 && (-e) < DOUBLE_LOG5_OF_TWO_TO_THE_N) 268 { 269 return toDoubleHighPrecision (f, length) / doubleTenToTheE (-e); 270 } 271 else if (e >= 0 && e < APPROX_MAX_MAGNITUDE) 272 { 273 result = toDoubleHighPrecision (f, length) * pow (10.0, e); 274 } 275 else if (e >= APPROX_MAX_MAGNITUDE) 276 { 277 /* Convert the partial result to make sure that the 278 * non-exponential part is not zero. This check fixes the case 279 * where the user enters 0.0e309! */ 280 result = toDoubleHighPrecision (f, length); 281 /* Don't go straight to zero as the fact that x*0 = 0 independent of x might 282 cause the algorithm to produce an incorrect result. Instead try the min value 283 first and let it fall to zero if need be. */ 284 285 if (result == 0.0) 286 { 287 DOUBLE_TO_LONGBITS (result) = DOUBLE_MINIMUM_LONGBITS; 288 } 289 else 290 { 291 DOUBLE_TO_LONGBITS (result) = DOUBLE_INFINITE_LONGBITS; 292 return result; 293 } 294 } 295 else if (e > APPROX_MIN_MAGNITUDE) 296 { 297 result = toDoubleHighPrecision (f, length) / pow (10.0, -e); 298 } 299 300 if (e <= APPROX_MIN_MAGNITUDE) 301 { 302 303 result = toDoubleHighPrecision (f, length) * pow (10.0, e + 52); 304 result = result * pow (10.0, -52); 305 306 } 307 308 /* Don't go straight to zero as the fact that x*0 = 0 independent of x might 309 cause the algorithm to produce an incorrect result. Instead try the min value 310 first and let it fall to zero if need be. */ 311 312 if (result == 0.0) 313 314 DOUBLE_TO_LONGBITS (result) = DOUBLE_MINIMUM_LONGBITS; 315 316 return doubleAlgorithm (env, f, length, e, result); 317 } 318 319 /* The algorithm for the function doubleAlgorithm() below can be found 320 * in: 321 * 322 * "How to Read Floating-Point Numbers Accurately", William D. 323 * Clinger, Proceedings of the ACM SIGPLAN '90 Conference on 324 * Programming Language Design and Implementation, June 20-22, 325 * 1990, pp. 92-101. 326 * 327 * There is a possibility that the function will end up in an endless 328 * loop if the given approximating floating-point number (a very small 329 * floating-point whose value is very close to zero) straddles between 330 * two approximating integer values. We modified the algorithm slightly 331 * to detect the case where it oscillates back and forth between 332 * incrementing and decrementing the floating-point approximation. It 333 * is currently set such that if the oscillation occurs more than twice 334 * then return the original approximation. 335 */ 336 static jdouble doubleAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jdouble z) { 337 uint64_t m; 338 int32_t k, comparison, comparison2; 339 uint64_t* x; 340 uint64_t* y; 341 uint64_t* D; 342 uint64_t* D2; 343 int32_t xLength, yLength, DLength, D2Length, decApproxCount, incApproxCount; 344 345 x = y = D = D2 = 0; 346 xLength = yLength = DLength = D2Length = 0; 347 decApproxCount = incApproxCount = 0; 348 349 do 350 { 351 m = doubleMantissa (z); 352 k = doubleExponent (z); 353 354 if (x && x != f) 355 free(x); 356 357 free(y); 358 free(D); 359 free(D2); 360 361 if (e >= 0 && k >= 0) 362 { 363 xLength = sizeOfTenToTheE (e) + length; 364 allocateU64 (x, xLength); 365 memset (x + length, 0, sizeof (uint64_t) * (xLength - length)); 366 memcpy (x, f, sizeof (uint64_t) * length); 367 timesTenToTheEHighPrecision (x, xLength, e); 368 369 yLength = (k >> 6) + 2; 370 allocateU64 (y, yLength); 371 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1)); 372 *y = m; 373 simpleShiftLeftHighPrecision (y, yLength, k); 374 } 375 else if (e >= 0) 376 { 377 xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1; 378 allocateU64 (x, xLength); 379 memset (x + length, 0, sizeof (uint64_t) * (xLength - length)); 380 memcpy (x, f, sizeof (uint64_t) * length); 381 timesTenToTheEHighPrecision (x, xLength, e); 382 simpleShiftLeftHighPrecision (x, xLength, -k); 383 384 yLength = 1; 385 allocateU64 (y, 1); 386 *y = m; 387 } 388 else if (k >= 0) 389 { 390 xLength = length; 391 x = f; 392 393 yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6); 394 allocateU64 (y, yLength); 395 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1)); 396 *y = m; 397 timesTenToTheEHighPrecision (y, yLength, -e); 398 simpleShiftLeftHighPrecision (y, yLength, k); 399 } 400 else 401 { 402 xLength = length + ((-k) >> 6) + 1; 403 allocateU64 (x, xLength); 404 memset (x + length, 0, sizeof (uint64_t) * (xLength - length)); 405 memcpy (x, f, sizeof (uint64_t) * length); 406 simpleShiftLeftHighPrecision (x, xLength, -k); 407 408 yLength = sizeOfTenToTheE (-e) + 1; 409 allocateU64 (y, yLength); 410 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1)); 411 *y = m; 412 timesTenToTheEHighPrecision (y, yLength, -e); 413 } 414 415 comparison = compareHighPrecision (x, xLength, y, yLength); 416 if (comparison > 0) 417 { /* x > y */ 418 DLength = xLength; 419 allocateU64 (D, DLength); 420 memcpy (D, x, DLength * sizeof (uint64_t)); 421 subtractHighPrecision (D, DLength, y, yLength); 422 } 423 else if (comparison) 424 { /* y > x */ 425 DLength = yLength; 426 allocateU64 (D, DLength); 427 memcpy (D, y, DLength * sizeof (uint64_t)); 428 subtractHighPrecision (D, DLength, x, xLength); 429 } 430 else 431 { /* y == x */ 432 DLength = 1; 433 allocateU64 (D, 1); 434 *D = 0; 435 } 436 437 D2Length = DLength + 1; 438 allocateU64 (D2, D2Length); 439 m <<= 1; 440 multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length); 441 m >>= 1; 442 443 comparison2 = compareHighPrecision (D2, D2Length, y, yLength); 444 if (comparison2 < 0) 445 { 446 if (comparison < 0 && m == DOUBLE_NORMAL_MASK) 447 { 448 simpleShiftLeftHighPrecision (D2, D2Length, 1); 449 if (compareHighPrecision (D2, D2Length, y, yLength) > 0) 450 { 451 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount); 452 } 453 else 454 { 455 break; 456 } 457 } 458 else 459 { 460 break; 461 } 462 } 463 else if (comparison2 == 0) 464 { 465 if ((LOW_U32_FROM_VAR (m) & 1) == 0) 466 { 467 if (comparison < 0 && m == DOUBLE_NORMAL_MASK) 468 { 469 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount); 470 } 471 else 472 { 473 break; 474 } 475 } 476 else if (comparison < 0) 477 { 478 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount); 479 break; 480 } 481 else 482 { 483 INCREMENT_DOUBLE (z, decApproxCount, incApproxCount); 484 break; 485 } 486 } 487 else if (comparison < 0) 488 { 489 DECREMENT_DOUBLE (z, decApproxCount, incApproxCount); 490 } 491 else 492 { 493 if (DOUBLE_TO_LONGBITS (z) == DOUBLE_INFINITE_LONGBITS) 494 break; 495 INCREMENT_DOUBLE (z, decApproxCount, incApproxCount); 496 } 497 } 498 while (1); 499 500 if (x && x != f) 501 free(x); 502 free(y); 503 free(D); 504 free(D2); 505 return z; 506 507 OutOfMemory: 508 if (x && x != f) 509 free(x); 510 free(y); 511 free(D); 512 free(D2); 513 jniThrowOutOfMemoryError(env, NULL); 514 return z; 515 } 516 517 518 519 #define MAX_FLOAT_ACCURACY_WIDTH 8 520 521 #define DEFAULT_FLOAT_WIDTH MAX_FLOAT_ACCURACY_WIDTH 522 523 static jfloat createFloat1(JNIEnv* env, uint64_t* f, int32_t length, jint e); 524 static jfloat floatAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jfloat z); 525 526 static const uint32_t float_tens[] = { 527 0x3f800000, 528 0x41200000, 529 0x42c80000, 530 0x447a0000, 531 0x461c4000, 532 0x47c35000, 533 0x49742400, 534 0x4b189680, 535 0x4cbebc20, 536 0x4e6e6b28, 537 0x501502f9 /* 10 ^ 10 in float */ 538 }; 539 540 #define floatTenToTheE(e) (*reinterpret_cast<const jfloat*>(float_tens + (e))) 541 #define FLOAT_LOG5_OF_TWO_TO_THE_N 11 542 543 #define FLOAT_INFINITE_INTBITS (0x7F800000) 544 #define FLOAT_MINIMUM_INTBITS (1) 545 546 #define FLOAT_MANTISSA_MASK (0x007FFFFF) 547 #define FLOAT_EXPONENT_MASK (0x7F800000) 548 #define FLOAT_NORMAL_MASK (0x00800000) 549 550 /* Keep a count of the number of times we decrement and increment to 551 * approximate the double, and attempt to detect the case where we 552 * could potentially toggle back and forth between decrementing and 553 * incrementing. It is possible for us to be stuck in the loop when 554 * incrementing by one or decrementing by one may exceed or stay below 555 * the value that we are looking for. In this case, just break out of 556 * the loop if we toggle between incrementing and decrementing for more 557 * than twice. 558 */ 559 #define INCREMENT_FLOAT(_x, _decCount, _incCount) \ 560 { \ 561 ++FLOAT_TO_INTBITS(_x); \ 562 _incCount++; \ 563 if( (_incCount > 2) && (_decCount > 2) ) { \ 564 if( _decCount > _incCount ) { \ 565 FLOAT_TO_INTBITS(_x) += _decCount - _incCount; \ 566 } else if( _incCount > _decCount ) { \ 567 FLOAT_TO_INTBITS(_x) -= _incCount - _decCount; \ 568 } \ 569 break; \ 570 } \ 571 } 572 #define DECREMENT_FLOAT(_x, _decCount, _incCount) \ 573 { \ 574 --FLOAT_TO_INTBITS(_x); \ 575 _decCount++; \ 576 if( (_incCount > 2) && (_decCount > 2) ) { \ 577 if( _decCount > _incCount ) { \ 578 FLOAT_TO_INTBITS(_x) += _decCount - _incCount; \ 579 } else if( _incCount > _decCount ) { \ 580 FLOAT_TO_INTBITS(_x) -= _incCount - _decCount; \ 581 } \ 582 break; \ 583 } \ 584 } 585 586 static jfloat createFloat(JNIEnv* env, const char* s, jint e) { 587 /* assumes s is a null terminated string with at least one 588 * character in it */ 589 uint64_t def[DEFAULT_FLOAT_WIDTH]; 590 uint64_t defBackup[DEFAULT_FLOAT_WIDTH]; 591 uint64_t* f; 592 uint64_t* fNoOverflow; 593 uint64_t* g; 594 uint64_t* tempBackup; 595 uint32_t overflow; 596 jfloat result; 597 int32_t index = 1; 598 int unprocessedDigits = 0; 599 600 f = def; 601 fNoOverflow = defBackup; 602 *f = 0; 603 tempBackup = g = 0; 604 do 605 { 606 if (*s >= '0' && *s <= '9') 607 { 608 /* Make a back up of f before appending, so that we can 609 * back out of it if there is no more room, i.e. index > 610 * MAX_FLOAT_ACCURACY_WIDTH. 611 */ 612 memcpy (fNoOverflow, f, sizeof (uint64_t) * index); 613 overflow = 614 simpleAppendDecimalDigitHighPrecision (f, index, *s - '0'); 615 if (overflow) 616 { 617 618 f[index++] = overflow; 619 /* There is an overflow, but there is no more room 620 * to store the result. We really only need the top 52 621 * bits anyway, so we must back out of the overflow, 622 * and ignore the rest of the string. 623 */ 624 if (index >= MAX_FLOAT_ACCURACY_WIDTH) 625 { 626 index--; 627 memcpy (f, fNoOverflow, sizeof (uint64_t) * index); 628 break; 629 } 630 if (tempBackup) 631 { 632 fNoOverflow = tempBackup; 633 } 634 } 635 } 636 else 637 index = -1; 638 } 639 while (index > 0 && *(++s) != '\0'); 640 641 /* We've broken out of the parse loop either because we've reached 642 * the end of the string or we've overflowed the maximum accuracy 643 * limit of a double. If we still have unprocessed digits in the 644 * given string, then there are three possible results: 645 * 1. (unprocessed digits + e) == 0, in which case we simply 646 * convert the existing bits that are already parsed 647 * 2. (unprocessed digits + e) < 0, in which case we simply 648 * convert the existing bits that are already parsed along 649 * with the given e 650 * 3. (unprocessed digits + e) > 0 indicates that the value is 651 * simply too big to be stored as a double, so return Infinity 652 */ 653 if ((unprocessedDigits = strlen (s)) > 0) 654 { 655 e += unprocessedDigits; 656 if (index > -1) 657 { 658 if (e <= 0) 659 { 660 result = createFloat1 (env, f, index, e); 661 } 662 else 663 { 664 FLOAT_TO_INTBITS (result) = FLOAT_INFINITE_INTBITS; 665 } 666 } 667 else 668 { 669 result = INTBITS_TO_FLOAT(index); 670 } 671 } 672 else 673 { 674 if (index > -1) 675 { 676 result = createFloat1 (env, f, index, e); 677 } 678 else 679 { 680 result = INTBITS_TO_FLOAT(index); 681 } 682 } 683 684 return result; 685 686 } 687 688 static jfloat createFloat1 (JNIEnv* env, uint64_t* f, int32_t length, jint e) { 689 int32_t numBits; 690 jdouble dresult; 691 jfloat result; 692 693 numBits = highestSetBitHighPrecision (f, length) + 1; 694 if (numBits < 25 && e >= 0 && e < FLOAT_LOG5_OF_TWO_TO_THE_N) 695 { 696 return ((jfloat) LOW_I32_FROM_PTR (f)) * floatTenToTheE (e); 697 } 698 else if (numBits < 25 && e < 0 && (-e) < FLOAT_LOG5_OF_TWO_TO_THE_N) 699 { 700 return ((jfloat) LOW_I32_FROM_PTR (f)) / floatTenToTheE (-e); 701 } 702 else if (e >= 0 && e < 39) 703 { 704 result = (jfloat) (toDoubleHighPrecision (f, length) * pow (10.0, (double) e)); 705 } 706 else if (e >= 39) 707 { 708 /* Convert the partial result to make sure that the 709 * non-exponential part is not zero. This check fixes the case 710 * where the user enters 0.0e309! */ 711 result = (jfloat) toDoubleHighPrecision (f, length); 712 713 if (result == 0.0) 714 715 FLOAT_TO_INTBITS (result) = FLOAT_MINIMUM_INTBITS; 716 else 717 FLOAT_TO_INTBITS (result) = FLOAT_INFINITE_INTBITS; 718 } 719 else if (e > -309) 720 { 721 int dexp; 722 uint32_t fmant, fovfl; 723 uint64_t dmant; 724 dresult = toDoubleHighPrecision (f, length) / pow (10.0, (double) -e); 725 if (IS_DENORMAL_DBL (dresult)) 726 { 727 FLOAT_TO_INTBITS (result) = 0; 728 return result; 729 } 730 dexp = doubleExponent (dresult) + 51; 731 dmant = doubleMantissa (dresult); 732 /* Is it too small to be represented by a single-precision 733 * float? */ 734 if (dexp <= -155) 735 { 736 FLOAT_TO_INTBITS (result) = 0; 737 return result; 738 } 739 /* Is it a denormalized single-precision float? */ 740 if ((dexp <= -127) && (dexp > -155)) 741 { 742 /* Only interested in 24 msb bits of the 53-bit double mantissa */ 743 fmant = (uint32_t) (dmant >> 29); 744 fovfl = ((uint32_t) (dmant & 0x1FFFFFFF)) << 3; 745 while ((dexp < -127) && ((fmant | fovfl) != 0)) 746 { 747 if ((fmant & 1) != 0) 748 { 749 fovfl |= 0x80000000; 750 } 751 fovfl >>= 1; 752 fmant >>= 1; 753 dexp++; 754 } 755 if ((fovfl & 0x80000000) != 0) 756 { 757 if ((fovfl & 0x7FFFFFFC) != 0) 758 { 759 fmant++; 760 } 761 else if ((fmant & 1) != 0) 762 { 763 fmant++; 764 } 765 } 766 else if ((fovfl & 0x40000000) != 0) 767 { 768 if ((fovfl & 0x3FFFFFFC) != 0) 769 { 770 fmant++; 771 } 772 } 773 FLOAT_TO_INTBITS (result) = fmant; 774 } 775 else 776 { 777 result = (jfloat) dresult; 778 } 779 } 780 781 /* Don't go straight to zero as the fact that x*0 = 0 independent 782 * of x might cause the algorithm to produce an incorrect result. 783 * Instead try the min value first and let it fall to zero if need 784 * be. 785 */ 786 if (e <= -309 || FLOAT_TO_INTBITS (result) == 0) 787 FLOAT_TO_INTBITS (result) = FLOAT_MINIMUM_INTBITS; 788 789 return floatAlgorithm (env, f, length, e, (jfloat) result); 790 } 791 792 /* The algorithm for the function floatAlgorithm() below can be found 793 * in: 794 * 795 * "How to Read Floating-Point Numbers Accurately", William D. 796 * Clinger, Proceedings of the ACM SIGPLAN '90 Conference on 797 * Programming Language Design and Implementation, June 20-22, 798 * 1990, pp. 92-101. 799 * 800 * There is a possibility that the function will end up in an endless 801 * loop if the given approximating floating-point number (a very small 802 * floating-point whose value is very close to zero) straddles between 803 * two approximating integer values. We modified the algorithm slightly 804 * to detect the case where it oscillates back and forth between 805 * incrementing and decrementing the floating-point approximation. It 806 * is currently set such that if the oscillation occurs more than twice 807 * then return the original approximation. 808 */ 809 static jfloat floatAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jfloat z) { 810 uint64_t m; 811 int32_t k, comparison, comparison2; 812 uint64_t* x; 813 uint64_t* y; 814 uint64_t* D; 815 uint64_t* D2; 816 int32_t xLength, yLength, DLength, D2Length; 817 int32_t decApproxCount, incApproxCount; 818 819 x = y = D = D2 = 0; 820 xLength = yLength = DLength = D2Length = 0; 821 decApproxCount = incApproxCount = 0; 822 823 do 824 { 825 m = floatMantissa (z); 826 k = floatExponent (z); 827 828 if (x && x != f) 829 free(x); 830 831 free(y); 832 free(D); 833 free(D2); 834 835 if (e >= 0 && k >= 0) 836 { 837 xLength = sizeOfTenToTheE (e) + length; 838 allocateU64 (x, xLength); 839 memset (x + length, 0, sizeof (uint64_t) * (xLength - length)); 840 memcpy (x, f, sizeof (uint64_t) * length); 841 timesTenToTheEHighPrecision (x, xLength, e); 842 843 yLength = (k >> 6) + 2; 844 allocateU64 (y, yLength); 845 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1)); 846 *y = m; 847 simpleShiftLeftHighPrecision (y, yLength, k); 848 } 849 else if (e >= 0) 850 { 851 xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1; 852 allocateU64 (x, xLength); 853 memset (x + length, 0, sizeof (uint64_t) * (xLength - length)); 854 memcpy (x, f, sizeof (uint64_t) * length); 855 timesTenToTheEHighPrecision (x, xLength, e); 856 simpleShiftLeftHighPrecision (x, xLength, -k); 857 858 yLength = 1; 859 allocateU64 (y, 1); 860 *y = m; 861 } 862 else if (k >= 0) 863 { 864 xLength = length; 865 x = f; 866 867 yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6); 868 allocateU64 (y, yLength); 869 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1)); 870 *y = m; 871 timesTenToTheEHighPrecision (y, yLength, -e); 872 simpleShiftLeftHighPrecision (y, yLength, k); 873 } 874 else 875 { 876 xLength = length + ((-k) >> 6) + 1; 877 allocateU64 (x, xLength); 878 memset (x + length, 0, sizeof (uint64_t) * (xLength - length)); 879 memcpy (x, f, sizeof (uint64_t) * length); 880 simpleShiftLeftHighPrecision (x, xLength, -k); 881 882 yLength = sizeOfTenToTheE (-e) + 1; 883 allocateU64 (y, yLength); 884 memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1)); 885 *y = m; 886 timesTenToTheEHighPrecision (y, yLength, -e); 887 } 888 889 comparison = compareHighPrecision (x, xLength, y, yLength); 890 if (comparison > 0) 891 { /* x > y */ 892 DLength = xLength; 893 allocateU64 (D, DLength); 894 memcpy (D, x, DLength * sizeof (uint64_t)); 895 subtractHighPrecision (D, DLength, y, yLength); 896 } 897 else if (comparison) 898 { /* y > x */ 899 DLength = yLength; 900 allocateU64 (D, DLength); 901 memcpy (D, y, DLength * sizeof (uint64_t)); 902 subtractHighPrecision (D, DLength, x, xLength); 903 } 904 else 905 { /* y == x */ 906 DLength = 1; 907 allocateU64 (D, 1); 908 *D = 0; 909 } 910 911 D2Length = DLength + 1; 912 allocateU64 (D2, D2Length); 913 m <<= 1; 914 multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length); 915 m >>= 1; 916 917 comparison2 = compareHighPrecision (D2, D2Length, y, yLength); 918 if (comparison2 < 0) 919 { 920 if (comparison < 0 && m == FLOAT_NORMAL_MASK) 921 { 922 simpleShiftLeftHighPrecision (D2, D2Length, 1); 923 if (compareHighPrecision (D2, D2Length, y, yLength) > 0) 924 { 925 DECREMENT_FLOAT (z, decApproxCount, incApproxCount); 926 } 927 else 928 { 929 break; 930 } 931 } 932 else 933 { 934 break; 935 } 936 } 937 else if (comparison2 == 0) 938 { 939 if ((m & 1) == 0) 940 { 941 if (comparison < 0 && m == FLOAT_NORMAL_MASK) 942 { 943 DECREMENT_FLOAT (z, decApproxCount, incApproxCount); 944 } 945 else 946 { 947 break; 948 } 949 } 950 else if (comparison < 0) 951 { 952 DECREMENT_FLOAT (z, decApproxCount, incApproxCount); 953 break; 954 } 955 else 956 { 957 INCREMENT_FLOAT (z, decApproxCount, incApproxCount); 958 break; 959 } 960 } 961 else if (comparison < 0) 962 { 963 DECREMENT_FLOAT (z, decApproxCount, incApproxCount); 964 } 965 else 966 { 967 if (FLOAT_TO_INTBITS (z) == FLOAT_EXPONENT_MASK) 968 break; 969 INCREMENT_FLOAT (z, decApproxCount, incApproxCount); 970 } 971 } 972 while (1); 973 974 if (x && x != f) 975 free(x); 976 free(y); 977 free(D); 978 free(D2); 979 return z; 980 981 OutOfMemory: 982 if (x && x != f) 983 free(x); 984 free(y); 985 free(D); 986 free(D2); 987 jniThrowOutOfMemoryError(env, NULL); 988 return z; 989 } 990 991 static jfloat StringToReal_parseFltImpl(JNIEnv* env, jclass, jstring s, jint e) { 992 ScopedUtfChars str(env, s); 993 if (str.c_str() == NULL) { 994 return 0.0; 995 } 996 return createFloat(env, str.c_str(), e); 997 } 998 999 static jdouble StringToReal_parseDblImpl(JNIEnv* env, jclass, jstring s, jint e) { 1000 ScopedUtfChars str(env, s); 1001 if (str.c_str() == NULL) { 1002 return 0.0; 1003 } 1004 return createDouble(env, str.c_str(), e); 1005 } 1006 1007 static JNINativeMethod gMethods[] = { 1008 NATIVE_METHOD(StringToReal, parseFltImpl, "(Ljava/lang/String;I)F"), 1009 NATIVE_METHOD(StringToReal, parseDblImpl, "(Ljava/lang/String;I)D"), 1010 }; 1011 void register_java_lang_StringToReal(JNIEnv* env) { 1012 jniRegisterNativeMethods(env, "java/lang/StringToReal", gMethods, NELEM(gMethods)); 1013 } 1014
