path: root/luni/src/main/native/java_lang_StringToReal.cpp

/*
 *  Licensed to the Apache Software Foundation (ASF) under one or more
 *  contributor license agreements.  See the NOTICE file distributed with
 *  this work for additional information regarding copyright ownership.
 *  The ASF licenses this file to You under the Apache License, Version 2.0
 *  (the "License"); you may not use this file except in compliance with
 *  the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 */

#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "JNIHelp.h"
#include "JniConstants.h"
#include "JniException.h"
#include "ScopedUtfChars.h"
#include "cbigint.h"

/* ************************* Defines ************************* */

#define LOW_I32_FROM_VAR(u64)     LOW_I32_FROM_LONG64(u64)
#define LOW_I32_FROM_PTR(u64ptr)  LOW_I32_FROM_LONG64_PTR(u64ptr)
#define HIGH_I32_FROM_VAR(u64)    HIGH_I32_FROM_LONG64(u64)
#define HIGH_I32_FROM_PTR(u64ptr) HIGH_I32_FROM_LONG64_PTR(u64ptr)

#define MAX_DOUBLE_ACCURACY_WIDTH 17

#define DEFAULT_DOUBLE_WIDTH MAX_DOUBLE_ACCURACY_WIDTH

#define DOUBLE_INFINITE_LONGBITS (0x7FF0000000000000LL)

#define DOUBLE_MINIMUM_LONGBITS (0x1)

#define DOUBLE_MANTISSA_MASK (0x000FFFFFFFFFFFFFLL)
#define DOUBLE_EXPONENT_MASK (0x7FF0000000000000LL)
#define DOUBLE_NORMAL_MASK   (0x0010000000000000LL)

#define allocateU64(x, n) if (!((x) = reinterpret_cast<uint64_t*>(malloc((n) * sizeof(uint64_t))))) goto OutOfMemory;

/* *********************************************************** */

/* ************************ local data ************************ */
static const jdouble double_tens[] = {
  1.0,
  1.0e1,
  1.0e2,
  1.0e3,
  1.0e4,
  1.0e5,
  1.0e6,
  1.0e7,
  1.0e8,
  1.0e9,
  1.0e10,
  1.0e11,
  1.0e12,
  1.0e13,
  1.0e14,
  1.0e15,
  1.0e16,
  1.0e17,
  1.0e18,
  1.0e19,
  1.0e20,
  1.0e21,
  1.0e22
};
/* *********************************************************** */

/* ************** private function declarations ************** */
static jdouble createDouble1   (JNIEnv* env, uint64_t * f, int32_t length, jint e);
static jdouble doubleAlgorithm (JNIEnv* env, uint64_t * f, int32_t length, jint e,
                                jdouble z);
/* *********************************************************** */

#define doubleTenToTheE(e) (*(double_tens + (e)))
#define DOUBLE_LOG5_OF_TWO_TO_THE_N 23

#define sizeOfTenToTheE(e) (((e) / 19) + 1)

static jdouble createDouble(JNIEnv* env, const char* s, jint e) {
  /* assumes s is a null terminated string with at least one
   * character in it */
  uint64_t def[DEFAULT_DOUBLE_WIDTH];
  uint64_t defBackup[DEFAULT_DOUBLE_WIDTH];
  uint64_t* f;
  uint64_t* fNoOverflow;
  uint64_t* g;
  uint64_t* tempBackup;
  uint32_t overflow;
  jdouble result;
  int32_t index = 1;
  int unprocessedDigits = 0;

  f = def;
  fNoOverflow = defBackup;
  *f = 0;
  tempBackup = g = 0;
  do
    {
      if (*s >= '0' && *s <= '9')
        {
          /* Make a back up of f before appending, so that we can
           * back out of it if there is no more room, i.e. index >
           * MAX_DOUBLE_ACCURACY_WIDTH.
           */
          memcpy (fNoOverflow, f, sizeof (uint64_t) * index);
          overflow =
            simpleAppendDecimalDigitHighPrecision (f, index, *s - '0');
          if (overflow)
            {
              f[index++] = overflow;
              /* There is an overflow, but there is no more room
               * to store the result. We really only need the top 52
               * bits anyway, so we must back out of the overflow,
               * and ignore the rest of the string.
               */
              if (index >= MAX_DOUBLE_ACCURACY_WIDTH)
                {
                  index--;
                  memcpy (f, fNoOverflow, sizeof (uint64_t) * index);
                  break;
                }
              if (tempBackup)
                {
                  fNoOverflow = tempBackup;
                }
            }
        }
      else
        index = -1;
    }
  while (index > 0 && *(++s) != '\0');

  /* We've broken out of the parse loop either because we've reached
   * the end of the string or we've overflowed the maximum accuracy
   * limit of a double. If we still have unprocessed digits in the
   * given string, then there are three possible results:
   *   1. (unprocessed digits + e) == 0, in which case we simply
   *      convert the existing bits that are already parsed
   *   2. (unprocessed digits + e) < 0, in which case we simply
   *      convert the existing bits that are already parsed along
   *      with the given e
   *   3. (unprocessed digits + e) > 0 indicates that the value is
   *      simply too big to be stored as a double, so return Infinity
   */
  if ((unprocessedDigits = strlen (s)) > 0)
    {
      e += unprocessedDigits;
      if (index > -1)
        {
          if (e == 0)
            result = toDoubleHighPrecision (f, index);
          else if (e < 0)
            result = createDouble1 (env, f, index, e);
          else
            {
              DOUBLE_TO_LONGBITS (result) = DOUBLE_INFINITE_LONGBITS;
            }
        }
      else
        {
          LOW_I32_FROM_VAR  (result) = -1;
          HIGH_I32_FROM_VAR (result) = -1;
        }
    }
  else
    {
      if (index > -1)
        {
          if (e == 0)
            result = toDoubleHighPrecision (f, index);
          else
            result = createDouble1 (env, f, index, e);
        }
      else
        {
          LOW_I32_FROM_VAR  (result) = -1;
          HIGH_I32_FROM_VAR (result) = -1;
        }
    }

  return result;
}

static jdouble createDouble1(JNIEnv* env, uint64_t* f, int32_t length, jint e) {
  int32_t numBits;
  jdouble result;

  static const jint APPROX_MIN_MAGNITUDE = -309;
  static const jint APPROX_MAX_MAGNITUDE = 309;

  numBits = highestSetBitHighPrecision (f, length) + 1;
  numBits -= lowestSetBitHighPrecision (f, length);
  if (numBits < 54 && e >= 0 && e < DOUBLE_LOG5_OF_TWO_TO_THE_N)
    {
      return toDoubleHighPrecision (f, length) * doubleTenToTheE (e);
    }
  else if (numBits < 54 && e < 0 && (-e) < DOUBLE_LOG5_OF_TWO_TO_THE_N)
    {
      return toDoubleHighPrecision (f, length) / doubleTenToTheE (-e);
    }
  else if (e >= 0 && e < APPROX_MAX_MAGNITUDE)
    {
      result = toDoubleHighPrecision (f, length) * pow (10.0, e);
    }
  else if (e >= APPROX_MAX_MAGNITUDE)
    {
      /* Convert the partial result to make sure that the
       * non-exponential part is not zero. This check fixes the case
       * where the user enters 0.0e309! */
      result = toDoubleHighPrecision (f, length);
      /* Don't go straight to zero as the fact that x*0 = 0 independent of x might
         cause the algorithm to produce an incorrect result.  Instead try the min value
         first and let it fall to zero if need be. */

      if (result == 0.0)
        {
          DOUBLE_TO_LONGBITS (result) = DOUBLE_MINIMUM_LONGBITS;
        }
      else
        {
          DOUBLE_TO_LONGBITS (result) = DOUBLE_INFINITE_LONGBITS;
          return result;
        }
    }
  else if (e > APPROX_MIN_MAGNITUDE)
    {
      result = toDoubleHighPrecision (f, length) / pow (10.0, -e);
    }

  if (e <= APPROX_MIN_MAGNITUDE)
    {

      result = toDoubleHighPrecision (f, length) * pow (10.0, e + 52);
      result = result * pow (10.0, -52);

    }

  /* Don't go straight to zero as the fact that x*0 = 0 independent of x might
     cause the algorithm to produce an incorrect result.  Instead try the min value
     first and let it fall to zero if need be. */

  if (result == 0.0)
    DOUBLE_TO_LONGBITS (result) = DOUBLE_MINIMUM_LONGBITS;

  return doubleAlgorithm (env, f, length, e, result);
}

/* The algorithm for the function doubleAlgorithm() below can be found
 * in:
 *
 *      "How to Read Floating-Point Numbers Accurately", William D.
 *      Clinger, Proceedings of the ACM SIGPLAN '90 Conference on
 *      Programming Language Design and Implementation, June 20-22,
 *      1990, pp. 92-101.
 */
static jdouble doubleAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jdouble z) {
  uint64_t m;
  int32_t k, comparison, comparison2;
  uint64_t* x;
  uint64_t* y;
  uint64_t* D;
  uint64_t* D2;
  int32_t xLength, yLength, DLength, D2Length;

  x = y = D = D2 = 0;
  xLength = yLength = DLength = D2Length = 0;

  do
    {
      m = doubleMantissa (z);
      k = doubleExponent (z);

      if (x && x != f)
          free(x);

      free(y);
      free(D);
      free(D2);

      if (e >= 0 && k >= 0)
        {
          xLength = sizeOfTenToTheE (e) + length;
          allocateU64 (x, xLength);
          memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
          memcpy (x, f, sizeof (uint64_t) * length);
          timesTenToTheEHighPrecision (x, xLength, e);

          yLength = (k >> 6) + 2;
          allocateU64 (y, yLength);
          memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
          *y = m;
          simpleShiftLeftHighPrecision (y, yLength, k);
        }
      else if (e >= 0)
        {
          xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1;
          allocateU64 (x, xLength);
          memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
          memcpy (x, f, sizeof (uint64_t) * length);
          timesTenToTheEHighPrecision (x, xLength, e);
          simpleShiftLeftHighPrecision (x, xLength, -k);

          yLength = 1;
          allocateU64 (y, 1);
          *y = m;
        }
      else if (k >= 0)
        {
          xLength = length;
          x = f;

          yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6);
          allocateU64 (y, yLength);
          memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
          *y = m;
          timesTenToTheEHighPrecision (y, yLength, -e);
          simpleShiftLeftHighPrecision (y, yLength, k);
        }
      else
        {
          xLength = length + ((-k) >> 6) + 1;
          allocateU64 (x, xLength);
          memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
          memcpy (x, f, sizeof (uint64_t) * length);
          simpleShiftLeftHighPrecision (x, xLength, -k);

          yLength = sizeOfTenToTheE (-e) + 1;
          allocateU64 (y, yLength);
          memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
          *y = m;
          timesTenToTheEHighPrecision (y, yLength, -e);
        }

      comparison = compareHighPrecision (x, xLength, y, yLength);
      if (comparison > 0)
        {                       /* x > y */
          DLength = xLength;
          allocateU64 (D, DLength);
          memcpy (D, x, DLength * sizeof (uint64_t));
          subtractHighPrecision (D, DLength, y, yLength);
        }
      else if (comparison)
        {                       /* y > x */
          DLength = yLength;
          allocateU64 (D, DLength);
          memcpy (D, y, DLength * sizeof (uint64_t));
          subtractHighPrecision (D, DLength, x, xLength);
        }
      else
        {                       /* y == x */
          DLength = 1;
          allocateU64 (D, 1);
          *D = 0;
        }

      D2Length = DLength + 1;
      allocateU64 (D2, D2Length);
      m <<= 1;
      multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length);
      m >>= 1;

      comparison2 = compareHighPrecision (D2, D2Length, y, yLength);
      if (comparison2 < 0)
        {
          if (comparison < 0 && m == DOUBLE_NORMAL_MASK
              && DOUBLE_TO_LONGBITS(z) != DOUBLE_NORMAL_MASK)
            {
              simpleShiftLeftHighPrecision (D2, D2Length, 1);
              if (compareHighPrecision (D2, D2Length, y, yLength) > 0)
                {
                  --DOUBLE_TO_LONGBITS (z);
                }
              else
                {
                  break;
                }
            }
          else
            {
              break;
            }
        }
      else if (comparison2 == 0)
        {
          if ((LOW_U32_FROM_VAR (m) & 1) == 0)
            {
              if (comparison < 0 && m == DOUBLE_NORMAL_MASK)
                {
                  --DOUBLE_TO_LONGBITS (z);
                }
              else
                {
                  break;
                }
            }
          else if (comparison < 0)
            {
              --DOUBLE_TO_LONGBITS (z);
              break;
            }
          else
            {
              ++DOUBLE_TO_LONGBITS (z);
              break;
            }
        }
      else if (comparison < 0)
        {
          --DOUBLE_TO_LONGBITS (z);
        }
      else
        {
          if (DOUBLE_TO_LONGBITS (z) == DOUBLE_INFINITE_LONGBITS)
            break;
          ++DOUBLE_TO_LONGBITS (z);
        }
    }
  while (1);

  if (x && x != f)
     free(x);
  free(y);
  free(D);
  free(D2);
  return z;

OutOfMemory:
  if (x && x != f)
      free(x);
  free(y);
  free(D);
  free(D2);
  jniThrowOutOfMemoryError(env, NULL);
  return z;
}



#define MAX_FLOAT_ACCURACY_WIDTH 8

#define DEFAULT_FLOAT_WIDTH MAX_FLOAT_ACCURACY_WIDTH

static jfloat createFloat1(JNIEnv* env, uint64_t* f, int32_t length, jint e);
static jfloat floatAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jfloat z);

static const uint32_t float_tens[] = {
  0x3f800000,
  0x41200000,
  0x42c80000,
  0x447a0000,
  0x461c4000,
  0x47c35000,
  0x49742400,
  0x4b189680,
  0x4cbebc20,
  0x4e6e6b28,
  0x501502f9                    /* 10 ^ 10 in float */
};

#define floatTenToTheE(e) (*reinterpret_cast<const jfloat*>(float_tens + (e)))
#define FLOAT_LOG5_OF_TWO_TO_THE_N 11

#define FLOAT_INFINITE_INTBITS (0x7F800000)
#define FLOAT_MINIMUM_INTBITS (1)

#define FLOAT_MANTISSA_MASK (0x007FFFFF)
#define FLOAT_EXPONENT_MASK (0x7F800000)
#define FLOAT_NORMAL_MASK   (0x00800000)

static jfloat createFloat(JNIEnv* env, const char* s, jint e) {
  /* assumes s is a null terminated string with at least one
   * character in it */
  uint64_t def[DEFAULT_FLOAT_WIDTH];
  uint64_t defBackup[DEFAULT_FLOAT_WIDTH];
  uint64_t* f;
  uint64_t* fNoOverflow;
  uint64_t* g;
  uint64_t* tempBackup;
  uint32_t overflow;
  jfloat result;
  int32_t index = 1;
  int unprocessedDigits = 0;

  f = def;
  fNoOverflow = defBackup;
  *f = 0;
  tempBackup = g = 0;
  do
    {
      if (*s >= '0' && *s <= '9')
        {
          /* Make a back up of f before appending, so that we can
           * back out of it if there is no more room, i.e. index >
           * MAX_FLOAT_ACCURACY_WIDTH.
           */
          memcpy (fNoOverflow, f, sizeof (uint64_t) * index);
          overflow =
            simpleAppendDecimalDigitHighPrecision (f, index, *s - '0');
          if (overflow)
            {

              f[index++] = overflow;
              /* There is an overflow, but there is no more room
               * to store the result. We really only need the top 52
               * bits anyway, so we must back out of the overflow,
               * and ignore the rest of the string.
               */
              if (index >= MAX_FLOAT_ACCURACY_WIDTH)
                {
                  index--;
                  memcpy (f, fNoOverflow, sizeof (uint64_t) * index);
                  break;
                }
              if (tempBackup)
                {
                  fNoOverflow = tempBackup;
                }
            }
        }
      else
        index = -1;
    }
  while (index > 0 && *(++s) != '\0');

  /* We've broken out of the parse loop either because we've reached
   * the end of the string or we've overflowed the maximum accuracy
   * limit of a double. If we still have unprocessed digits in the
   * given string, then there are three possible results:
   *   1. (unprocessed digits + e) == 0, in which case we simply
   *      convert the existing bits that are already parsed
   *   2. (unprocessed digits + e) < 0, in which case we simply
   *      convert the existing bits that are already parsed along
   *      with the given e
   *   3. (unprocessed digits + e) > 0 indicates that the value is
   *      simply too big to be stored as a double, so return Infinity
   */
  if ((unprocessedDigits = strlen (s)) > 0)
    {
      e += unprocessedDigits;
      if (index > -1)
        {
          if (e <= 0)
            {
              result = createFloat1 (env, f, index, e);
            }
          else
            {
              FLOAT_TO_INTBITS (result) = FLOAT_INFINITE_INTBITS;
            }
        }
      else
        {
          result = INTBITS_TO_FLOAT(index);
        }
    }
  else
    {
      if (index > -1)
        {
          result = createFloat1 (env, f, index, e);
        }
      else
        {
          result = INTBITS_TO_FLOAT(index);
        }
    }

  return result;
}

static jfloat createFloat1 (JNIEnv* env, uint64_t* f, int32_t length, jint e) {
  int32_t numBits;
  jdouble dresult;
  jfloat result;

  numBits = highestSetBitHighPrecision (f, length) + 1;
  if (numBits < 25 && e >= 0 && e < FLOAT_LOG5_OF_TWO_TO_THE_N)
    {
      return ((jfloat) LOW_I32_FROM_PTR (f)) * floatTenToTheE (e);
    }
  else if (numBits < 25 && e < 0 && (-e) < FLOAT_LOG5_OF_TWO_TO_THE_N)
    {
      return ((jfloat) LOW_I32_FROM_PTR (f)) / floatTenToTheE (-e);
    }
  else if (e >= 0 && e < 39)
    {
      result = (jfloat) (toDoubleHighPrecision (f, length) * pow (10.0, (double) e));
    }
  else if (e >= 39)
    {
      /* Convert the partial result to make sure that the
       * non-exponential part is not zero. This check fixes the case
       * where the user enters 0.0e309! */
      result = (jfloat) toDoubleHighPrecision (f, length);

      if (result == 0.0)

        FLOAT_TO_INTBITS (result) = FLOAT_MINIMUM_INTBITS;
      else
        FLOAT_TO_INTBITS (result) = FLOAT_INFINITE_INTBITS;
    }
  else if (e > -309)
    {
      int dexp;
      uint32_t fmant, fovfl;
      uint64_t dmant;
      dresult = toDoubleHighPrecision (f, length) / pow (10.0, (double) -e);
      if (IS_DENORMAL_DBL (dresult))
        {
          FLOAT_TO_INTBITS (result) = 0;
          return result;
        }
      dexp = doubleExponent (dresult) + 51;
      dmant = doubleMantissa (dresult);
      /* Is it too small to be represented by a single-precision
       * float? */
      if (dexp <= -155)
        {
          FLOAT_TO_INTBITS (result) = 0;
          return result;
        }
      /* Is it a denormalized single-precision float? */
      if ((dexp <= -127) && (dexp > -155))
        {
          /* Only interested in 24 msb bits of the 53-bit double mantissa */
          fmant = (uint32_t) (dmant >> 29);
          fovfl = ((uint32_t) (dmant & 0x1FFFFFFF)) << 3;
          while ((dexp < -127) && ((fmant | fovfl) != 0))
            {
              if ((fmant & 1) != 0)
                {
                  fovfl |= 0x80000000;
                }
              fovfl >>= 1;
              fmant >>= 1;
              dexp++;
            }
          if ((fovfl & 0x80000000) != 0)
            {
              if ((fovfl & 0x7FFFFFFC) != 0)
                {
                  fmant++;
                }
              else if ((fmant & 1) != 0)
                {
                  fmant++;
                }
            }
          else if ((fovfl & 0x40000000) != 0)
            {
              if ((fovfl & 0x3FFFFFFC) != 0)
                {
                  fmant++;
                }
            }
          FLOAT_TO_INTBITS (result) = fmant;
        }
      else
        {
          result = (jfloat) dresult;
        }
    }

  /* Don't go straight to zero as the fact that x*0 = 0 independent
   * of x might cause the algorithm to produce an incorrect result.
   * Instead try the min  value first and let it fall to zero if need
   * be.
   */
  if (e <= -309 || FLOAT_TO_INTBITS (result) == 0)
    FLOAT_TO_INTBITS (result) = FLOAT_MINIMUM_INTBITS;

  return floatAlgorithm (env, f, length, e, (jfloat) result);
}

/* The algorithm for the function floatAlgorithm() below can be found
 * in:
 *
 *      "How to Read Floating-Point Numbers Accurately", William D.
 *      Clinger, Proceedings of the ACM SIGPLAN '90 Conference on
 *      Programming Language Design and Implementation, June 20-22,
 *      1990, pp. 92-101.
 */
static jfloat floatAlgorithm(JNIEnv* env, uint64_t* f, int32_t length, jint e, jfloat z) {
  uint64_t m;
  int32_t k, comparison, comparison2;
  uint64_t* x;
  uint64_t* y;
  uint64_t* D;
  uint64_t* D2;
  int32_t xLength, yLength, DLength, D2Length;

  x = y = D = D2 = 0;
  xLength = yLength = DLength = D2Length = 0;

  do
    {
      m = floatMantissa (z);
      k = floatExponent (z);

      if (x && x != f)
          free(x);

      free(y);
      free(D);
      free(D2);

      if (e >= 0 && k >= 0)
        {
          xLength = sizeOfTenToTheE (e) + length;
          allocateU64 (x, xLength);
          memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
          memcpy (x, f, sizeof (uint64_t) * length);
          timesTenToTheEHighPrecision (x, xLength, e);

          yLength = (k >> 6) + 2;
          allocateU64 (y, yLength);
          memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
          *y = m;
          simpleShiftLeftHighPrecision (y, yLength, k);
        }
      else if (e >= 0)
        {
          xLength = sizeOfTenToTheE (e) + length + ((-k) >> 6) + 1;
          allocateU64 (x, xLength);
          memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
          memcpy (x, f, sizeof (uint64_t) * length);
          timesTenToTheEHighPrecision (x, xLength, e);
          simpleShiftLeftHighPrecision (x, xLength, -k);

          yLength = 1;
          allocateU64 (y, 1);
          *y = m;
        }
      else if (k >= 0)
        {
          xLength = length;
          x = f;

          yLength = sizeOfTenToTheE (-e) + 2 + (k >> 6);
          allocateU64 (y, yLength);
          memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
          *y = m;
          timesTenToTheEHighPrecision (y, yLength, -e);
          simpleShiftLeftHighPrecision (y, yLength, k);
        }
      else
        {
          xLength = length + ((-k) >> 6) + 1;
          allocateU64 (x, xLength);
          memset (x + length, 0, sizeof (uint64_t) * (xLength - length));
          memcpy (x, f, sizeof (uint64_t) * length);
          simpleShiftLeftHighPrecision (x, xLength, -k);

          yLength = sizeOfTenToTheE (-e) + 1;
          allocateU64 (y, yLength);
          memset (y + 1, 0, sizeof (uint64_t) * (yLength - 1));
          *y = m;
          timesTenToTheEHighPrecision (y, yLength, -e);
        }

      comparison = compareHighPrecision (x, xLength, y, yLength);
      if (comparison > 0)
        {                       /* x > y */
          DLength = xLength;
          allocateU64 (D, DLength);
          memcpy (D, x, DLength * sizeof (uint64_t));
          subtractHighPrecision (D, DLength, y, yLength);
        }
      else if (comparison)
        {                       /* y > x */
          DLength = yLength;
          allocateU64 (D, DLength);
          memcpy (D, y, DLength * sizeof (uint64_t));
          subtractHighPrecision (D, DLength, x, xLength);
        }
      else
        {                       /* y == x */
          DLength = 1;
          allocateU64 (D, 1);
          *D = 0;
        }

      D2Length = DLength + 1;
      allocateU64 (D2, D2Length);
      m <<= 1;
      multiplyHighPrecision (D, DLength, &m, 1, D2, D2Length);
      m >>= 1;

      comparison2 = compareHighPrecision (D2, D2Length, y, yLength);
      if (comparison2 < 0)
        {
          if (comparison < 0 && m == FLOAT_NORMAL_MASK
              && FLOAT_TO_INTBITS(z) != FLOAT_NORMAL_MASK)
            {
              simpleShiftLeftHighPrecision (D2, D2Length, 1);
              if (compareHighPrecision (D2, D2Length, y, yLength) > 0)
                {
                  --FLOAT_TO_INTBITS (z);
                }
              else
                {
                  break;
                }
            }
          else
            {
              break;
            }
        }
      else if (comparison2 == 0)
        {
          if ((m & 1) == 0)
            {
              if (comparison < 0 && m == FLOAT_NORMAL_MASK)
                {
                  --FLOAT_TO_INTBITS (z);
                }
              else
                {
                  break;
                }
            }
          else if (comparison < 0)
            {
              --FLOAT_TO_INTBITS (z);
              break;
            }
          else
            {
              ++FLOAT_TO_INTBITS (z);
              break;
            }
        }
      else if (comparison < 0)
        {
          --FLOAT_TO_INTBITS (z);
        }
      else
        {
          if (FLOAT_TO_INTBITS (z) == FLOAT_EXPONENT_MASK)
            break;
          ++FLOAT_TO_INTBITS (z);
        }
    }
  while (1);

  if (x && x != f)
      free(x);
  free(y);
  free(D);
  free(D2);
  return z;

OutOfMemory:
  if (x && x != f)
      free(x);
  free(y);
  free(D);
  free(D2);
  jniThrowOutOfMemoryError(env, NULL);
  return z;
}

static jfloat StringToReal_parseFltImpl(JNIEnv* env, jclass, jstring s, jint e) {
    ScopedUtfChars str(env, s);
    if (str.c_str() == NULL) {
        return 0.0;
    }
    return createFloat(env, str.c_str(), e);
}

static jdouble StringToReal_parseDblImpl(JNIEnv* env, jclass, jstring s, jint e) {
    ScopedUtfChars str(env, s);
    if (str.c_str() == NULL) {
        return 0.0;
    }
    return createDouble(env, str.c_str(), e);
}

static JNINativeMethod gMethods[] = {
    NATIVE_METHOD(StringToReal, parseFltImpl, "(Ljava/lang/String;I)F"),
    NATIVE_METHOD(StringToReal, parseDblImpl, "(Ljava/lang/String;I)D"),
};
void register_java_lang_StringToReal(JNIEnv* env) {
    jniRegisterNativeMethods(env, "java/lang/StringToReal", gMethods, NELEM(gMethods));
}