Home | History | Annotate | Download | only in native
      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(__APPLE__)
     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