replace double with fp_t, replace int64_t with fp_component, refactor code to use the new typedefs

Author: jrahlf

src/printf/printf.c

@@ -103,6 +103,13 @@ extern "C" {
 #define PRINTF_DEFAULT_FLOAT_PRECISION  6
 #endif
 
+// Enable/disable high precission double type for internal processing of floating point values (i.e. va_arg(double))
+// If 0, uses  float which has less precision, but may be faster and reduce code size on devices without double FPU
+// default is 1 as this corresponds to the C standard implementation
+#ifndef PRINTF_USE_DOUBLE_FLOATING_POINT_PRECISION
+#define PRINTF_USE_DOUBLE_FLOATING_POINT_PRECISION  1
+#endif
+
 // According to the C languages standard, printf() and related functions must be able to print any
 // integral number in floating-point notation, regardless of length, when using the %f specifier -
 // possibly hundreds of characters, potentially overflowing your buffers. In this implementation,
@@ -204,24 +211,40 @@ typedef unsigned int printf_flags_t;
 #error "Non-binary-radix floating-point types are unsupported."
 #endif
 
-#if DBL_MANT_DIG == 24
+#if PRINTF_USE_DOUBLE_FLOATING_POINT_PRECISION == 1
+typedef double fp_t;
+#else
+typedef float fp_t;
+#endif
+
+#if (DBL_MANT_DIG == 24) || (PRINTF_USE_DOUBLE_FLOATING_POINT_PRECISION == 0)
 
-#define DOUBLE_SIZE_IN_BITS 32
-typedef uint32_t double_uint_t;
-#define DOUBLE_EXPONENT_MASK 0xFFU
-#define DOUBLE_BASE_EXPONENT 127
+#define FLOATING_POINT_SIZE_IN_BITS 32
+typedef uint32_t floating_point_uint_t;
+#define FLOATING_POINT_EXPONENT_MASK 0xFFU
+#define FLOATING_POINT_BASE_EXPONENT 127
+#define FLOATING_POINT_MANT_DIG 24
+#define FLOATING_POINT_MAX FLT_MAX
 
 #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 FLOATING_POINT_SIZE_IN_BITS 64
+typedef uint64_t floating_point_uint_t;
+#define FLOATING_POINT_EXPONENT_MASK 0x7FFU
+#define FLOATING_POINT_BASE_EXPONENT 1023
+#define FLOATING_POINT_MANT_DIG 53
+#define FLOATING_POINT_MAX DBL_MAX
 
 #else
 #error "Unsupported double type configuration"
 #endif
-#define DOUBLE_STORED_MANTISSA_BITS (DBL_MANT_DIG - 1)
+#define FLOATING_POINT_STORED_MANTISSA_BITS (FLOATING_POINT_MANT_DIG - 1)
+
+#if PRINTF_MAX_INTEGRAL_DIGITS_FOR_DECIMAL <= 9
+typedef int_fast32_t fp_component;
+#else
+typedef int_fast64_t fp_component;
+#endif
 
 #if PRINTF_SUPPORT_LONG_LONG
 typedef unsigned long long printf_unsigned_value_t;
@@ -245,35 +268,35 @@ typedef unsigned int printf_size_t;
   // trailing '\0'.
 
 typedef union {
-  double_uint_t U;
-  double        F;
-} double_with_bit_access;
+  floating_point_uint_t U;
+  fp_t                  F;
+} floating_point_with_bit_access;
 
 // This is unnecessary in C99, since compound initializers can be used,
 // but:
 // 1. Some compilers are finicky about this;
 // 2. Some people may want to convert this to C89;
 // 3. If you try to use it as C++, only C++20 supports compound literals
-static inline double_with_bit_access get_bit_access(double x)
+static inline floating_point_with_bit_access get_bit_access(fp_t x)
 {
-  double_with_bit_access dwba;
+  floating_point_with_bit_access dwba;
   dwba.F = x;
   return dwba;
 }
 
-static inline int get_sign_bit(double x)
+static inline int get_sign_bit(fp_t x)
 {
   // The sign is stored in the highest bit
-  return (int) (get_bit_access(x).U >> (DOUBLE_SIZE_IN_BITS - 1));
+  return (int) (get_bit_access(x).U >> (FLOATING_POINT_SIZE_IN_BITS - 1));
 }
 
-static inline int get_exp2(double_with_bit_access x)
+static inline int get_exp2(floating_point_with_bit_access x)
 {
   // The exponent in an IEEE-754 floating-point number occupies a contiguous
   // sequence of bits (e.g. 52..62 for 64-bit doubles), but with a non-trivial representation: An
   // unsigned offset from some negative value (with the extremal offset values reserved for
   // special use).
-  return (int)((x.U >> DOUBLE_STORED_MANTISSA_BITS ) & DOUBLE_EXPONENT_MASK) - DOUBLE_BASE_EXPONENT;
+  return (int)((x.U >> FLOATING_POINT_STORED_MANTISSA_BITS ) & FLOATING_POINT_EXPONENT_MASK) - FLOATING_POINT_BASE_EXPONENT;
 }
 #define PRINTF_ABS(_x) ( (_x) > 0 ? (_x) : -(_x) )
 
@@ -539,56 +562,56 @@ static void print_integer(output_gadget_t* output, printf_unsigned_value_t value
 
 #if (PRINTF_SUPPORT_DECIMAL_SPECIFIERS || PRINTF_SUPPORT_EXPONENTIAL_SPECIFIERS)
 
-// Stores a fixed-precision representation of a double relative
+// Stores a fixed-precision representation of a float/double relative
 // to a fixed precision (which cannot be determined by examining this structure)
-struct double_components {
-  int_fast64_t integral;
-  int_fast64_t fractional;
-    // ... truncation of the actual fractional part of the double value, scaled
+struct floating_point_components {
+  fp_component integral;
+  fp_component fractional;
+    // ... truncation of the actual fractional part of the float/double value, scaled
     // by the precision value
   bool is_negative;
 };
 
 #define NUM_DECIMAL_DIGITS_IN_INT64_T 18
 #define PRINTF_MAX_PRECOMPUTED_POWER_OF_10  NUM_DECIMAL_DIGITS_IN_INT64_T
-static const double powers_of_10[NUM_DECIMAL_DIGITS_IN_INT64_T] = {
-  1e00, 1e01, 1e02, 1e03, 1e04, 1e05, 1e06, 1e07, 1e08,
-  1e09, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17
+static const fp_t powers_of_10[NUM_DECIMAL_DIGITS_IN_INT64_T] = {
+  (fp_t) 1e00, (fp_t) 1e01, (fp_t) 1e02, (fp_t) 1e03, (fp_t) 1e04, (fp_t) 1e05, (fp_t) 1e06, (fp_t) 1e07, (fp_t) 1e08,
+  (fp_t) 1e09, (fp_t) 1e10, (fp_t) 1e11, (fp_t) 1e12, (fp_t) 1e13, (fp_t) 1e14, (fp_t) 1e15, (fp_t) 1e16, (fp_t) 1e17
 };
 
 #define PRINTF_MAX_SUPPORTED_PRECISION NUM_DECIMAL_DIGITS_IN_INT64_T - 1
 
 
-// Break up a double number - which is known to be a finite non-negative number -
+// Break up a float/double number - which is known to be a finite non-negative number -
 // into its base-10 parts: integral - before the decimal point, and fractional - after it.
 // Taken the precision into account, but does not change it even internally.
-static struct double_components get_components(double number, printf_size_t precision)
+static struct floating_point_components get_components(fp_t number, printf_size_t precision)
 {
-  struct double_components number_;
+  struct floating_point_components number_;
   number_.is_negative = get_sign_bit(number);
-  double abs_number = (number_.is_negative) ? -number : number;
-  number_.integral = (int_fast64_t)abs_number;
-  double remainder = (abs_number - (double) number_.integral) * powers_of_10[precision];
-  number_.fractional = (int_fast64_t)remainder;
+  fp_t abs_number = (number_.is_negative) ? -number : number;
+  number_.integral = (fp_component)abs_number;
+  fp_t remainder = (abs_number - (fp_t) number_.integral) * powers_of_10[precision];
+  number_.fractional = (fp_component)remainder;
 
-  remainder -= (double) number_.fractional;
+  remainder -= (fp_t) number_.fractional;
 
-  if (remainder > 0.5) {
+  if (remainder > (fp_t) 0.5) {
     ++number_.fractional;
     // handle rollover, e.g. case 0.99 with precision 1 is 1.0
-    if ((double) number_.fractional >= powers_of_10[precision]) {
+    if ((fp_t) number_.fractional >= powers_of_10[precision]) {
       number_.fractional = 0;
       ++number_.integral;
     }
   }
-  else if ((remainder == 0.5) && ((number_.fractional == 0U) || (number_.fractional & 1U))) {
+  else if ((remainder == (fp_t) 0.5) && ((number_.fractional == 0U) || (number_.fractional & 1U))) {
     // if halfway, round up if odd OR if last digit is 0
     ++number_.fractional;
   }
 
   if (precision == 0U) {
-    remainder = abs_number - (double) number_.integral;
-    if ((!(remainder < 0.5) || (remainder > 0.5)) && (number_.integral & 1)) {
+    remainder = abs_number - (fp_t) number_.integral;
+    if ((!(remainder < (fp_t) 0.5) || (remainder > (fp_t) 0.5)) && (number_.integral & 1)) {
       // exactly 0.5 and ODD, then round up
       // 1.5 -> 2, but 2.5 -> 2
       ++number_.integral;
@@ -599,21 +622,21 @@ static struct double_components get_components(double number, printf_size_t prec
 
 #if PRINTF_SUPPORT_EXPONENTIAL_SPECIFIERS
 struct scaling_factor {
-  double raw_factor;
+  fp_t raw_factor;
   bool multiply; // if true, need to multiply by raw_factor; otherwise need to divide by it
 };
 
-static double apply_scaling(double num, struct scaling_factor normalization)
+static fp_t apply_scaling(fp_t num, struct scaling_factor normalization)
 {
   return normalization.multiply ? num * normalization.raw_factor : num / normalization.raw_factor;
 }
 
-static double unapply_scaling(double normalized, struct scaling_factor normalization)
+static fp_t unapply_scaling(fp_t normalized, struct scaling_factor normalization)
 {
   return normalization.multiply ? normalized / normalization.raw_factor : normalized * normalization.raw_factor;
 }
 
-static struct scaling_factor update_normalization(struct scaling_factor sf, double extra_multiplicative_factor)
+static struct scaling_factor update_normalization(struct scaling_factor sf, fp_t extra_multiplicative_factor)
 {
   struct scaling_factor result;
   if (sf.multiply) {
@@ -637,37 +660,37 @@ 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 floating_point_components get_normalized_components(bool negative, printf_size_t precision, fp_t non_normalized, struct scaling_factor normalization)
 {
-  struct double_components components;
+  struct floating_point_components components;
   components.is_negative = negative;
-  components.integral = (int_fast64_t) apply_scaling(non_normalized, normalization);
-  double remainder = non_normalized - unapply_scaling((double) components.integral, normalization);
-  double prec_power_of_10 = powers_of_10[precision];
+  components.integral = (fp_component) apply_scaling(non_normalized, normalization);
+  fp_t remainder = non_normalized - unapply_scaling((fp_t) components.integral, normalization);
+  fp_t prec_power_of_10 = powers_of_10[precision];
   struct scaling_factor account_for_precision = update_normalization(normalization, prec_power_of_10);
-  double scaled_remainder = apply_scaling(remainder, account_for_precision);
-  double rounding_threshold = 0.5;
+  fp_t scaled_remainder = apply_scaling(remainder, account_for_precision);
+  fp_t rounding_threshold = 0.5;
 
   if (precision == 0U) {
     components.fractional = 0;
     components.integral += (scaled_remainder >= rounding_threshold);
     if (scaled_remainder == rounding_threshold) {
       // banker's rounding: Round towards the even number (making the mean error 0)
-      components.integral &= ~((int_fast64_t) 0x1);
+      components.integral &= ~((fp_component) 0x1);
     }
   }
   else {
-    components.fractional = (int_fast64_t) scaled_remainder;
-    scaled_remainder -= (double) components.fractional;
+    components.fractional = (fp_component) scaled_remainder;
+    scaled_remainder -= (fp_t) components.fractional;
 
     components.fractional += (scaled_remainder >= rounding_threshold);
     if (scaled_remainder == rounding_threshold) {
       // banker's rounding: Round towards the even number (making the mean error 0)
-      components.fractional &= ~((int_fast64_t) 0x1);
+      components.fractional &= ~((fp_component) 0x1);
     }
     // handle rollover, e.g. the case of 0.99 with precision 1 becoming (0,100),
     // and must then be corrected into (1, 0).
-    if ((double) components.fractional >= prec_power_of_10) {
+    if ((fp_t) components.fractional >= prec_power_of_10) {
       components.fractional = 0;
       ++components.integral;
     }
@@ -677,7 +700,7 @@ static struct double_components get_normalized_components(bool negative, printf_
 #endif // PRINTF_SUPPORT_EXPONENTIAL_SPECIFIERS
 
 static void print_broken_up_decimal(
-  struct double_components number_, output_gadget_t* output, printf_size_t precision,
+  struct floating_point_components number_, output_gadget_t* output, printf_size_t precision,
   printf_size_t width, printf_flags_t flags, char *buf, printf_size_t len)
 {
   if (precision != 0U) {
@@ -688,7 +711,7 @@ static void print_broken_up_decimal(
     // %g/%G mandates we skip the trailing 0 digits...
     if ((flags & FLAGS_ADAPT_EXP) && !(flags & FLAGS_HASH) && (number_.fractional > 0)) {
       while(true) {
-        int_fast64_t digit = number_.fractional % 10U;
+        fp_component digit = number_.fractional % 10U;
         if (digit != 0) {
           break;
         }
@@ -759,44 +782,44 @@ static void print_broken_up_decimal(
 }
 
       // internal ftoa for fixed decimal floating point
-static void print_decimal_number(output_gadget_t* output, double number, printf_size_t precision, printf_size_t width, printf_flags_t flags, char* buf, printf_size_t len)
+static void print_decimal_number(output_gadget_t* output, fp_t number, printf_size_t precision, printf_size_t width, printf_flags_t flags, char* buf, printf_size_t len)
 {
-  struct double_components value_ = get_components(number, precision);
+  struct floating_point_components value_ = get_components(number, precision);
   print_broken_up_decimal(value_, output, precision, width, flags, buf, len);
 }
 
 #if PRINTF_SUPPORT_EXPONENTIAL_SPECIFIERS
 // internal ftoa variant for exponential floating-point type, contributed by Martijn Jasperse <[email protected]>
-static void print_exponential_number(output_gadget_t* output, double number, printf_size_t precision, printf_size_t width, printf_flags_t flags, char* buf, printf_size_t len)
+static void print_exponential_number(output_gadget_t* output, fp_t number, printf_size_t precision, printf_size_t width, printf_flags_t flags, char* buf, printf_size_t len)
 {
   const bool negative = get_sign_bit(number);
   // This number will decrease gradually (by factors of 10) as we "extract" the exponent out of it
-  double abs_number =  negative ? -number : number;
+  fp_t abs_number =  negative ? -number : number;
 
   int exp10;
   bool abs_exp10_covered_by_powers_table;
   struct scaling_factor normalization;
 
 
   // Determine the decimal exponent
-  if (abs_number == 0.0) {
+  if (abs_number == (fp_t) 0.0) {
     // TODO: This is a special-case for 0.0 (and -0.0); but proper handling is required for denormals more generally.
     exp10 = 0; // ... and no need to set a normalization factor or check the powers table
   }
   else  {
-    double_with_bit_access conv = get_bit_access(abs_number);
+    floating_point_with_bit_access conv = get_bit_access(abs_number);
     {
       // based on the algorithm by David Gay (https://www.ampl.com/netlib/fp/dtoa.c)
       int exp2 = get_exp2(conv);
 	  // drop the exponent, so conv.F comes into the range [1,2)
-      conv.U = (conv.U & (( (double_uint_t)(1) << DOUBLE_STORED_MANTISSA_BITS) - 1U)) | ((double_uint_t) DOUBLE_BASE_EXPONENT << DOUBLE_STORED_MANTISSA_BITS);
+      conv.U = (conv.U & (( (floating_point_uint_t)(1) << FLOATING_POINT_STORED_MANTISSA_BITS) - 1U)) | ((floating_point_uint_t) FLOATING_POINT_BASE_EXPONENT << FLOATING_POINT_STORED_MANTISSA_BITS);
       // now approximate log10 from the log2 integer part and an expansion of ln around 1.5
-      exp10 = (int)(0.1760912590558 + exp2 * 0.301029995663981 + (conv.F - 1.5) * 0.289529654602168);
+      exp10 = (int)((fp_t) 0.1760912590558 + ((fp_t) exp2) * ((fp_t) 0.301029995663981) + (conv.F - (fp_t) 1.5) * ((fp_t) 0.289529654602168));
       // now we want to compute 10^exp10 but we want to be sure it won't overflow
-      exp2 = (int)(exp10 * 3.321928094887362 + 0.5);
-      const double z  = exp10 * 2.302585092994046 - exp2 * 0.6931471805599453;
-      const double z2 = z * z;
-      conv.U = ((double_uint_t)(exp2) + DOUBLE_BASE_EXPONENT) << DOUBLE_STORED_MANTISSA_BITS;
+      exp2 = (int)(((fp_t) exp10) * ((fp_t) 3.321928094887362) + (fp_t) 0.5);
+      const fp_t z  = ((fp_t) exp10) * ((fp_t) 2.302585092994046) - ((fp_t) exp2) * ((fp_t) 0.6931471805599453);
+      const fp_t z2 = z * z;
+      conv.U = ((floating_point_uint_t)(exp2) + FLOATING_POINT_BASE_EXPONENT) << FLOATING_POINT_STORED_MANTISSA_BITS;
       // compute exp(z) using continued fractions, see https://en.wikipedia.org/wiki/Exponential_function#Continued_fractions_for_ex
       conv.F *= 1 + 2 * z / (2 - z + (z2 / (6 + (z2 / (10 + z2 / 14)))));
       // correct for rounding errors
@@ -833,7 +856,7 @@ static void print_exponential_number(output_gadget_t* output, double number, pri
 
   normalization.multiply = (exp10 < 0 && abs_exp10_covered_by_powers_table);
   bool should_skip_normalization = (fall_back_to_decimal_only_mode || exp10 == 0);
-  struct double_components decimal_part_components =
+  struct floating_point_components decimal_part_components =
     should_skip_normalization ?
     get_components(negative ? -abs_number : abs_number, precision) :
     get_normalized_components(negative, precision, abs_number, normalization);
@@ -894,7 +917,7 @@ static void print_exponential_number(output_gadget_t* output, double number, pri
 }
 #endif  // PRINTF_SUPPORT_EXPONENTIAL_SPECIFIERS
 
-static void print_floating_point(output_gadget_t* output, double value, printf_size_t precision, printf_size_t width, printf_flags_t flags, bool prefer_exponential)
+static void print_floating_point(output_gadget_t* output, fp_t value, printf_size_t precision, printf_size_t width, printf_flags_t flags, bool prefer_exponential)
 {
   char buf[PRINTF_FTOA_BUFFER_SIZE];
   printf_size_t len = 0U;
@@ -904,17 +927,17 @@ static void print_floating_point(output_gadget_t* output, double value, printf_s
     out_rev_(output, "nan", 3, width, flags);
     return;
   }
-  if (value < -DBL_MAX) {
+  if (value < -FLOATING_POINT_MAX) {
     out_rev_(output, "fni-", 4, width, flags);
     return;
   }
-  if (value > DBL_MAX) {
+  if (value > FLOATING_POINT_MAX) {
     out_rev_(output, (flags & FLAGS_PLUS) ? "fni+" : "fni", (flags & FLAGS_PLUS) ? 4U : 3U, width, flags);
     return;
   }
 
   if (!prefer_exponential &&
-      ((value > PRINTF_FLOAT_NOTATION_THRESHOLD) || (value < -PRINTF_FLOAT_NOTATION_THRESHOLD))) {
+      ((value > (fp_t) PRINTF_FLOAT_NOTATION_THRESHOLD) || (value < - (fp_t) PRINTF_FLOAT_NOTATION_THRESHOLD))) {
     // The required behavior of standard printf is to print _every_ integral-part digit -- which could mean
     // printing hundreds of characters, overflowing any fixed internal buffer and necessitating a more complicated
     // implementation.
@@ -1171,7 +1194,7 @@ static int _vsnprintf(output_gadget_t* output, const char* format, va_list args)
       case 'f' :
       case 'F' :
         if (*format == 'F') flags |= FLAGS_UPPERCASE;
-        print_floating_point(output, va_arg(args, double), precision, width, flags, PRINTF_PREFER_DECIMAL);
+        print_floating_point(output, (fp_t) va_arg(args, double), precision, width, flags, PRINTF_PREFER_DECIMAL);
         format++;
         break;
 #endif
@@ -1182,7 +1205,7 @@ static int _vsnprintf(output_gadget_t* output, const char* format, va_list args)
       case 'G':
         if ((*format == 'g')||(*format == 'G')) flags |= FLAGS_ADAPT_EXP;
         if ((*format == 'E')||(*format == 'G')) flags |= FLAGS_UPPERCASE;
-        print_floating_point(output, va_arg(args, double), precision, width, flags, PRINTF_PREFER_EXPONENTIAL);
+        print_floating_point(output, (fp_t) va_arg(args, double), precision, width, flags, PRINTF_PREFER_EXPONENTIAL);
         format++;
         break;
 #endif  // PRINTF_SUPPORT_EXPONENTIAL_SPECIFIERS