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