vendor/bela/src/bela/numbers.cc
inline ExpDigitsA SplitToSixA(const double value) {
  ExpDigitsA exp_dig;
  int exp = 5;
  double d = value;
  // First step: calculate a close approximation of the output, where the
  // value d will be between 100,000 and 999,999, representing the digits
  // in the output ASCII array, and exp is the base-10 exponent.  It would be
  // faster to use a table here, and to look up the base-2 exponent of value,
  // however value is an IEEE-754 64-bit number, so the table would have 2,000
  // entries, which is not cache-friendly.
  if (d >= 999999.5) {
    if (d >= 1e+261) {
      exp += 256;
      d *= 1e-256;
    }
    if (d >= 1e+133) {
      exp += 128;
      d *= 1e-128;
    }
    if (d >= 1e+69) {
      exp += 64;
      d *= 1e-64;
    }
    if (d >= 1e+37) {
      exp += 32;
      d *= 1e-32;
    }
    if (d >= 1e+21) {
      exp += 16;
      d *= 1e-16;
    }
    if (d >= 1e+13) {
      exp += 8;
      d *= 1e-8;
    }
    if (d >= 1e+9) {
      exp += 4;
      d *= 1e-4;
    }
    if (d >= 1e+7) {
      exp += 2;
      d *= 1e-2;
    }
    if (d >= 1e+6) {
      exp += 1;
      d *= 1e-1;
    }
  } else {
    if (d < 1e-250) {
      exp -= 256;
      d *= 1e256;
    }
    if (d < 1e-122) {
      exp -= 128;
      d *= 1e128;
    }
    if (d < 1e-58) {
      exp -= 64;
      d *= 1e64;
    }
    if (d < 1e-26) {
      exp -= 32;
      d *= 1e32;
    }
    if (d < 1e-10) {
      exp -= 16;
      d *= 1e16;
    }
    if (d < 1e-2) {
      exp -= 8;
      d *= 1e8;
    }
    if (d < 1e+2) {
      exp -= 4;
      d *= 1e4;
    }
    if (d < 1e+4) {
      exp -= 2;
      d *= 1e2;
    }
    if (d < 1e+5) {
      exp -= 1;
      d *= 1e1;
    }
  }
  // At this point, d is in the range [99999.5..999999.5) and exp is in the
  // range [-324..308]. Since we need to round d up, we want to add a half
  // and truncate.
  // However, the technique above may have lost some precision, due to its
  // repeated multiplication by constants that each may be off by half a bit
  // of precision.  This only matters if we're close to the edge though.
  // Since we'd like to know if the fractional part of d is close to a half,
  // we multiply it by 65536 and see if the fractional part is close to 32768.
  // (The number doesn't have to be a power of two,but powers of two are faster)
  uint64_t d64k = static_cast<uint64_t>(d * 65536);
  int dddddd; // A 6-digit decimal integer.
  if ((d64k % 65536) == 32767 || (d64k % 65536) == 32768) {
    // OK, it's fairly likely that precision was lost above, which is
    // not a surprise given only 52 mantissa bits are available.  Therefore
    // redo the calculation using 128-bit numbers.  (64 bits are not enough).

    // Start out with digits rounded down; maybe add one below.
    dddddd = static_cast<int>(d64k / 65536);

    // mantissa is a 64-bit integer representing M.mmm... * 2^63.  The actual
    // value we're representing, of course, is M.mmm... * 2^exp2.
    int exp2;
    double m = std::frexp(value, &exp2);
    uint64_t mantissa = static_cast<uint64_t>(m * (32768.0 * 65536.0 * 65536.0 * 65536.0));
    // std::frexp returns an m value in the range [0.5, 1.0), however we
    // can't multiply it by 2^64 and convert to an integer because some FPUs
    // throw an exception when converting an number higher than 2^63 into an
    // integer - even an unsigned 64-bit integer!  Fortunately it doesn't matter
    // since m only has 52 significant bits anyway.
    mantissa <<= 1;
    exp2 -= 64; // not needed, but nice for debugging

    // OK, we are here to compare:
    //     (dddddd + 0.5) * 10^(exp-5)  vs.  mantissa * 2^exp2
    // so we can round up dddddd if appropriate.  Those values span the full
    // range of 600 orders of magnitude of IEE 64-bit floating-point.
    // Fortunately, we already know they are very close, so we don't need to
    // track the base-2 exponent of both sides.  This greatly simplifies the
    // the math since the 2^exp2 calculation is unnecessary and the power-of-10
    // calculation can become a power-of-5 instead.

    std::pair<uint64_t, uint64_t> edge, val;
    if (exp >= 6) {
      // Compare (dddddd + 0.5) * 5 ^ (exp - 5) to mantissa
      // Since we're tossing powers of two, 2 * dddddd + 1 is the
      // same as dddddd + 0.5
      edge = PowFive(2 * dddddd + 1, exp - 5);

      val.first = mantissa;
      val.second = 0;
    } else {
      // We can't compare (dddddd + 0.5) * 5 ^ (exp - 5) to mantissa as we did
      // above because (exp - 5) is negative.  So we compare (dddddd + 0.5) to
      // mantissa * 5 ^ (5 - exp)
      edge = PowFive(2 * dddddd + 1, 0);

      val = PowFive(mantissa, 5 - exp);
    }
    // printf("exp=%d %016lx %016lx vs %016lx %016lx\n", exp, val.first,
    //        val.second, edge.first, edge.second);
    if (val > edge) {
      dddddd++;
    } else if (val == edge) {
      dddddd += (dddddd & 1);
    }
  } else {
    // Here, we are not close to the edge.
    dddddd = static_cast<int>((d64k + 32768) / 65536);
  }
  if (dddddd == 1000000) {
    dddddd = 100000;
    exp += 1;
  }
  exp_dig.exponent = exp;

  int two_digits = dddddd / 10000;
  dddddd -= two_digits * 10000;
  PutTwoDigits(two_digits, &exp_dig.digits[0]);

  two_digits = dddddd / 100;
  dddddd -= two_digits * 100;
  PutTwoDigits(two_digits, &exp_dig.digits[2]);

  PutTwoDigits(dddddd, &exp_dig.digits[4]);
  return exp_dig;
}
// SplitToSix converts value, a positive double-precision floating-point number,
// into a base-10 exponent and 6 ASCII digits, where the first digit is never
// zero.  For example, SplitToSix(1) returns an exponent of zero and a digits
// array of {'1', '0', '0', '0', '0', '0'}.  If value is exactly halfway between
// two possible representations, e.g. value = 100000.5, then "round to even" is
// performed.
SplitToSix converts value, a positive double-precision floating-point number,
into a base-10 exponent and 6 ASCII digits, where the first digit is never
zero.
SplitToSix converts value a positive double - precision floating - point number into a base - 10 exponent and 6 ASCII digits where the first digit is never zero .
inline ExpDigitsA SplitToSixA(const double value) {
  ExpDigitsA exp_dig;
  int exp = 5;
  double d = value;
  if (d >= 999999.5) {
    if (d >= 1e+261) {
      exp += 256;
      d *= 1e-256;
    }
    if (d >= 1e+133) {
      exp += 128;
      d *= 1e-128;
    }
    if (d >= 1e+69) {
      exp += 64;
      d *= 1e-64;
    }
    if (d >= 1e+37) {
      exp += 32;
      d *= 1e-32;
    }
    if (d >= 1e+21) {
      exp += 16;
      d *= 1e-16;
    }
    if (d >= 1e+13) {
      exp += 8;
      d *= 1e-8;
    }
    if (d >= 1e+9) {
      exp += 4;
      d *= 1e-4;
    }
    if (d >= 1e+7) {
      exp += 2;
      d *= 1e-2;
    }
    if (d >= 1e+6) {
      exp += 1;
      d *= 1e-1;
    }
  } else {
    if (d < 1e-250) {
      exp -= 256;
      d *= 1e256;
    }
    if (d < 1e-122) {
      exp -= 128;
      d *= 1e128;
    }
    if (d < 1e-58) {
      exp -= 64;
      d *= 1e64;
    }
    if (d < 1e-26) {
      exp -= 32;
      d *= 1e32;
    }
    if (d < 1e-10) {
      exp -= 16;
      d *= 1e16;
    }
    if (d < 1e-2) {
      exp -= 8;
      d *= 1e8;
    }
    if (d < 1e+2) {
      exp -= 4;
      d *= 1e4;
    }
    if (d < 1e+4) {
      exp -= 2;
      d *= 1e2;
    }
    if (d < 1e+5) {
      exp -= 1;
      d *= 1e1;
    }
  }
  uint64_t d64k = static_cast<uint64_t>(d * 65536);
  int dddddd; 
  if ((d64k % 65536) == 32767 || (d64k % 65536) == 32768) {
    dddddd = static_cast<int>(d64k / 65536);
    int exp2;
    double m = std::frexp(value, &exp2);
    uint64_t mantissa = static_cast<uint64_t>(m * (32768.0 * 65536.0 * 65536.0 * 65536.0));
    mantissa <<= 1;
    exp2 -= 64; 
    std::pair<uint64_t, uint64_t> edge, val;
    if (exp >= 6) {
      edge = PowFive(2 * dddddd + 1, exp - 5);
      val.first = mantissa;
      val.second = 0;
    } else {
      edge = PowFive(2 * dddddd + 1, 0);
      val = PowFive(mantissa, 5 - exp);
    }
    if (val > edge) {
      dddddd++;
    } else if (val == edge) {
      dddddd += (dddddd & 1);
    }
  } else {
    dddddd = static_cast<int>((d64k + 32768) / 65536);
  }
  if (dddddd == 1000000) {
    dddddd = 100000;
    exp += 1;
  }
  exp_dig.exponent = exp;
  int two_digits = dddddd / 10000;
  dddddd -= two_digits * 10000;
  PutTwoDigits(two_digits, &exp_dig.digits[0]);
  two_digits = dddddd / 100;
  dddddd -= two_digits * 100;
  PutTwoDigits(two_digits, &exp_dig.digits[2]);
  PutTwoDigits(dddddd, &exp_dig.digits[4]);
  return exp_dig;
}
inline ExpDigitsA SplitToSixA ( const double value ) { ExpDigitsA exp_dig ; int exp = 5 ; double d = value ; if ( d >= 999999.5 ) { if ( d >= 1e+261 ) { exp += 256 ; d *= 1e-256 ; } if ( d >= 1e+133 ) { exp += 128 ; d *= 1e-128 ; } if ( d >= 1e+69 ) { exp += 64 ; d *= 1e-64 ; } if ( d >= 1e+37 ) { exp += 32 ; d *= 1e-32 ; } if ( d >= 1e+21 ) { exp += 16 ; d *= 1e-16 ; } if ( d >= 1e+13 ) { exp += 8 ; d *= 1e-8 ; } if ( d >= 1e+9 ) { exp += 4 ; d *= 1e-4 ; } if ( d >= 1e+7 ) { exp += 2 ; d *= 1e-2 ; } if ( d >= 1e+6 ) { exp += 1 ; d *= 1e-1 ; } } else { if ( d < 1e-250 ) { exp -= 256 ; d *= 1e256 ; } if ( d < 1e-122 ) { exp -= 128 ; d *= 1e128 ; } if ( d < 1e-58 ) { exp -= 64 ; d *= 1e64 ; } if ( d < 1e-26 ) { exp -= 32 ; d *= 1e32 ; } if ( d < 1e-10 ) { exp -= 16 ; d *= 1e16 ; } if ( d < 1e-2 ) { exp -= 8 ; d *= 1e8 ; } if ( d < 1e+2 ) { exp -= 4 ; d *= 1e4 ; } if ( d < 1e+4 ) { exp -= 2 ; d *= 1e2 ; } if ( d < 1e+5 ) { exp -= 1 ; d *= 1e1 ; } } uint64_t d64k = static_cast < uint64_t > ( d * 65536 ) ; int dddddd ; if ( ( d64k % 65536 ) == 32767 || ( d64k % 65536 ) == 32768 ) { dddddd = static_cast < int > ( d64k / 65536 ) ; int exp2 ; double m = std :: frexp ( value , & exp2 ) ; uint64_t mantissa = static_cast < uint64_t > ( m * ( 32768.0 * 65536.0 * 65536.0 * 65536.0 ) ) ; mantissa <<= 1 ; exp2 -= 64 ; std :: pair < uint64_t , uint64_t > edge , val ; if ( exp >= 6 ) { edge = PowFive ( 2 * dddddd + 1 , exp - 5 ) ; val . first = mantissa ; val . second = 0 ; } else { edge = PowFive ( 2 * dddddd + 1 , 0 ) ; val = PowFive ( mantissa , 5 - exp ) ; } if ( val > edge ) { dddddd ++ ; } else if ( val == edge ) { dddddd += ( dddddd & 1 ) ; } } else { dddddd = static_cast < int > ( ( d64k + 32768 ) / 65536 ) ; } if ( dddddd == 1000000 ) { dddddd = 100000 ; exp += 1 ; } exp_dig . exponent = exp ; int two_digits = dddddd / 10000 ; dddddd -= two_digits * 10000 ; PutTwoDigits ( two_digits , & exp_dig . digits [ 0 ] ) ; two_digits = dddddd / 100 ; dddddd -= two_digits * 100 ; PutTwoDigits ( two_digits , & exp_dig . digits [ 2 ] ) ; PutTwoDigits ( dddddd , & exp_dig . digits [ 4 ] ) ; return exp_dig ; }
SplitToSix converts value, a positive double-precision floating-point number,
into a base-10 exponent and 6 ASCII digits, where the first digit is never
zero.
SplitToSix converts value a positive double - precision floating - point number into a base - 10 exponent and 6 ASCII digits where the first digit is never zero .
// First step: calculate a close approximation of the output, where the // value d will be between 100,000 and 999,999, representing the digits // in the output ASCII array, and exp is the base-10 exponent.  It would be // faster to use a table here, and to look up the base-2 exponent of value, // however value is an IEEE-754 64-bit number, so the table would have 2,000 // entries, which is not cache-friendly. // At this point, d is in the range [99999.5..999999.5) and exp is in the // range [-324..308]. Since we need to round d up, we want to add a half // and truncate. // However, the technique above may have lost some precision, due to its // repeated multiplication by constants that each may be off by half a bit // of precision.  This only matters if we're close to the edge though. // Since we'd like to know if the fractional part of d is close to a half, // we multiply it by 65536 and see if the fractional part is close to 32768. // (The number doesn't have to be a power of two,but powers of two are faster) // A 6-digit decimal integer. // OK, it's fairly likely that precision was lost above, which is // not a surprise given only 52 mantissa bits are available.  Therefore // redo the calculation using 128-bit numbers.  (64 bits are not enough). // Start out with digits rounded down; maybe add one below. // mantissa is a 64-bit integer representing M.mmm... * 2^63.  The actual // value we're representing, of course, is M.mmm... * 2^exp2. // std::frexp returns an m value in the range [0.5, 1.0), however we // can't multiply it by 2^64 and convert to an integer because some FPUs // throw an exception when converting an number higher than 2^63 into an // integer - even an unsigned 64-bit integer!  Fortunately it doesn't matter // since m only has 52 significant bits anyway. // not needed, but nice for debugging // OK, we are here to compare: //     (dddddd + 0.5) * 10^(exp-5)  vs.  mantissa * 2^exp2 // so we can round up dddddd if appropriate.  Those values span the full // range of 600 orders of magnitude of IEE 64-bit floating-point. // Fortunately, we already know they are very close, so we don't need to // track the base-2 exponent of both sides.  This greatly simplifies the // the math since the 2^exp2 calculation is unnecessary and the power-of-10 // calculation can become a power-of-5 instead. // Compare (dddddd + 0.5) * 5 ^ (exp - 5) to mantissa // Since we're tossing powers of two, 2 * dddddd + 1 is the // same as dddddd + 0.5 // We can't compare (dddddd + 0.5) * 5 ^ (exp - 5) to mantissa as we did // above because (exp - 5) is negative.  So we compare (dddddd + 0.5) to // mantissa * 5 ^ (5 - exp) // printf("exp=%d %016lx %016lx vs %016lx %016lx\n", exp, val.first, //        val.second, edge.first, edge.second); // Here, we are not close to the edge.