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