Fixes #114, #115: Better handling of powers near the extrema of the r…

…epresentation range.

* Special-casing the power-of-10 function for 10^-308.
* Altered `get_normalized_components()` to not strive for better precision for near-extremal powers (at both ends) and settle for the "rougher" way to obtain the integral and decimal parts.
* Added a test-case of a near-extremal-power which fails due to the re-normalization even if the powering is correct.

src/printf/printf.c

@@ -223,13 +223,17 @@ typedef unsigned int printf_flags_t;
 typedef uint32_t double_uint_t;
 #define DOUBLE_EXPONENT_MASK 0xFFU
 #define DOUBLE_BASE_EXPONENT 127
+#define DOUBLE_MAX_SUBNORMAL_EXPONENT_OF_10 -38
+#define DOUBLE_MAX_SUBNORMAL_POWER_OF_10 1e-38
 
 #elif DBL_MANT_DIG == 53
 
 #define DOUBLE_SIZE_IN_BITS 64
 typedef uint64_t double_uint_t;
 #define DOUBLE_EXPONENT_MASK 0x7FFU
 #define DOUBLE_BASE_EXPONENT 1023
+#define DOUBLE_MAX_SUBNORMAL_EXPONENT_OF_10 -308
+#define DOUBLE_MAX_SUBNORMAL_POWER_OF_10 1e-308
 
 #else
 #error "Unsupported double type configuration"
@@ -650,11 +654,20 @@ static struct scaling_factor update_normalization(struct scaling_factor sf, doub
   return result;
 }
 
-static struct double_components get_normalized_components(bool negative, printf_size_t precision, double non_normalized, struct scaling_factor normalization)
+static struct double_components get_normalized_components(bool negative, printf_size_t precision, double non_normalized, struct scaling_factor normalization, int floored_exp10)
 {
   struct double_components components;
   components.is_negative = negative;
-  components.integral = (int_fast64_t) apply_scaling(non_normalized, normalization);
+  double scaled = apply_scaling(non_normalized, normalization);
+
+  bool close_to_representation_extremum = ( (-floored_exp10 + (int) precision) >= DBL_MAX_10_EXP - 1 );
+  if (close_to_representation_extremum) {
+    // We can't have a normalization factor which also accounts for the precision, i.e. moves
+    // some decimal digits into the mantissa, since it's unrepresentable, or nearly unrepresentable.
+    // So, we'll give up early on getting extra precision...
+    return get_components(negative ? -scaled : scaled, precision);
+  }
+  components.integral = (int_fast64_t) scaled;
   double remainder = non_normalized - unapply_scaling((double) components.integral, normalization);
   double prec_power_of_10 = powers_of_10[precision];
   struct scaling_factor account_for_precision = update_normalization(normalization, prec_power_of_10);
@@ -819,6 +832,10 @@ static double log10_of_positive(double positive_number)
 
 static double pow10_of_int(int floored_exp10)
 {
+  // A crude hack for avoiding undesired behavior with barely-normal or slightly-subnormal values.
+  if (floored_exp10 == DOUBLE_MAX_SUBNORMAL_EXPONENT_OF_10) {
+    return DOUBLE_MAX_SUBNORMAL_POWER_OF_10;
+  }
   // Compute 10^(floored_exp10) but (try to) make sure that doesn't overflow
   double_with_bit_access dwba;
   int exp2 = bastardized_floor(floored_exp10 * 3.321928094887362 + 0.5);
@@ -887,7 +904,7 @@ static void print_exponential_number(output_gadget_t* output, double number, pri
   struct double_components decimal_part_components =
     should_skip_normalization ?
     get_components(negative ? -abs_number : abs_number, precision) :
-    get_normalized_components(negative, precision, abs_number, normalization);
+    get_normalized_components(negative, precision, abs_number, normalization, floored_exp10);
 
   // Account for roll-over, e.g. rounding from 9.99 to 100.0 - which effects
   // the exponent and may require additional tweaking of the parts

test/test_suite.cpp

@@ -1002,6 +1002,7 @@ TEST_CASE("floating-point specifiers, precision and flags", "[]" ) {
   PRINTING_CHECK("0.0012",           ==, sprintf_, buffer, "%.2G", 0.001234);
   PRINTING_CHECK(" +0.001234",       ==, sprintf_, buffer, "%+10.4G", 0.001234);
   PRINTING_CHECK("+001.234e-05",     ==, sprintf_, buffer, "%+012.4g", 0.00001234);
+  // Note: The following two values are _barely_ normal; make their mantissa 1.1 and they lose their normality.
   PRINTING_CHECK("-1.23e-308",       ==, sprintf_, buffer, "%.3g", -1.2345e-308);
   PRINTING_CHECK("+1.230E+308",      ==, sprintf_, buffer, "%+.3E", 1.23e+308);
   PRINTING_CHECK("1.000e+01",        ==, sprintf_, buffer, "%.3e", 9.9996);
@@ -1012,6 +1013,8 @@ TEST_CASE("floating-point specifiers, precision and flags", "[]" ) {
   PRINTING_CHECK("-4e+04",           ==, sprintf_, buffer, "%.1g", -40661.5);
   PRINTING_CHECK("-4.e+04",          ==, sprintf_, buffer, "%#.1g", -40661.5);
   PRINTING_CHECK("100.",             ==, sprintf_, buffer, "%#.3g", 99.998580932617187500);
+  // Note: The following value is _barely_ normal; make the mantissa 1.1 and it loses its normality.
+  PRINTING_CHECK("1.2345678901e-308", ==, sprintf_, buffer, "%.10e", 1.2345678901e-308);
   // Rounding-focused checks
   PRINTING_CHECK("4.895512e+04",     ==, sprintf_, buffer, "%e", 48955.125);
   PRINTING_CHECK("9.2524e+04",       ==, sprintf_, buffer, "%.4e", 92523.5);