Skip to content

Commit

Permalink
replace double with fp_t, replace int64_t with fp_component, refactor…
Browse files Browse the repository at this point in the history 
… code to use the new typedefs
  • Loading branch information
@jrahlf
jrahlf committed Feb 6, 2022
1 parent 36402dc commit b7e1f52
Showing 1 changed file with 100 additions and 77 deletions.
177 changes: 100 additions & 77 deletions src/printf/printf.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand Down Expand Up @@ -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;
Expand All @@ -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) )

Expand Down Expand Up @@ -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;
Expand All @@ -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) {
Expand All @@ -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;
}
Expand All @@ -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) {
Expand All @@ -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;
}
Expand Down Expand Up @@ -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 <m.jasperse@gmail.com>
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
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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;
Expand All @@ -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.
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down

0 comments on commit b7e1f52

Please sign in to comment.