C dynamically printf double, no loss of precision and no trailing zeroes

by

There’s probably no easier way. It’s a quite involved problem.

Your code isn’t solving it right for several reasons:

  • Most practical implementations of floating-point arithmetic aren’t decimal, they are binary. So, when you multiply a floating-point number by 10 or divide it by 10, you may lose precision (this depends on the number).
  • Even though the standard 64-bit IEEE-754 floating-point format reserves 53 bits for the mantissa, which is equivalent to floor(log10(2 ^ 53)) = 15 decimal digits, a valid number in this format may need up to some 1080 decimal digits in the fractional part when printed exactly, which is what you appear to be asking about.

One way of solving this is to use the %a format type specifier in snprintf(), which is going to print the floating-point value using hexadecimal digits for the mantissa and the C standard from 1999 guarantees that this will print all significant digits if the floating-point format is radix-2 (AKA base-2 or simply binary). So, with this you can obtain all the binary digits of the mantissa of the number. And from here you will be able to figure out how many decimal digits are in the fractional part.

Now, observe that:

1.00000 = 2+0 = 1.00000 (binary)
0.50000 = 2-1 = 0.10000
0.25000 = 2-2 = 0.01000
0.12500 = 2-3 = 0.00100
0.06250 = 2-4 = 0.00010
0.03125 = 2-5 = 0.00001

and so on.

You can clearly see here that a binary digit at i-th position to the right of the point in the binary representation produces the last non-zero decimal digit also in i-th position to the right of the point in the decimal representation.

So, if you know where the least significant non-zero bit is in a binary floating-point number, you can figure out how many decimal digits are needed to print the fractional part of the number exactly.

And this is what my program is doing.

Code:

// file: PrintFullFraction.c
//
// compile with gcc 4.6.2 or better:
//   gcc -Wall -Wextra -std=c99 -O2 PrintFullFraction.c -o PrintFullFraction.exe
#include <limits.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include <assert.h>

#if FLT_RADIX != 2
#error currently supported only FLT_RADIX = 2
#endif

int FractionalDigits(double d)
{
  char buf[
           1 + // sign, '-' or '+'
           (sizeof(d) * CHAR_BIT + 3) / 4 + // mantissa hex digits max
           1 + // decimal point, '.'
           1 + // mantissa-exponent separator, 'p'
           1 + // mantissa sign, '-' or '+'
           (sizeof(d) * CHAR_BIT + 2) / 3 + // exponent decimal digits max
           1 // string terminator, '\0'
          ];
  int n;
  char *pp, *p;
  int e, lsbFound, lsbPos;

  // convert d into "+/- 0x h.hhhh p +/- ddd" representation and check for errors
  if ((n = snprintf(buf, sizeof(buf), "%+a", d)) < 0 ||
      (unsigned)n >= sizeof(buf))
    return -1;

//printf("{%s}", buf);

  // make sure the conversion didn't produce something like "nan" or "inf"
  // instead of "+/- 0x h.hhhh p +/- ddd"
  if (strstr(buf, "0x") != buf + 1 ||
      (pp = strchr(buf, 'p')) == NULL)
    return 0;

  // extract the base-2 exponent manually, checking for overflows
  e = 0;
  p = pp + 1 + (pp[1] == '-' || pp[1] == '+'); // skip the exponent sign at first
  for (; *p != '\0'; p++)
  {
    if (e > INT_MAX / 10)
      return -2;
    e *= 10;
    if (e > INT_MAX - (*p - '0'))
      return -2;
    e += *p - '0';
  }
  if (pp[1] == '-') // apply the sign to the exponent
    e = -e;

//printf("[%s|%d]", buf, e);

  // find the position of the least significant non-zero bit
  lsbFound = lsbPos = 0;
  for (p = pp - 1; *p != 'x'; p--)
  {
    if (*p == '.')
      continue;
    if (!lsbFound)
    {
      int hdigit = (*p >= 'a') ? (*p - 'a' + 10) : (*p - '0'); // assuming ASCII chars
      if (hdigit)
      {
        static const int lsbPosInNibble[16] = { 0,4,3,4,  2,4,3,4, 1,4,3,4, 2,4,3,4 };
        lsbFound = 1;
        lsbPos = -lsbPosInNibble[hdigit];
      }
    }
    else
    {
      lsbPos -= 4;
    }
  }
  lsbPos += 4;

  if (!lsbFound)
    return 0; // d is 0 (integer)

  // adjust the least significant non-zero bit position
  // by the base-2 exponent (just add them), checking
  // for overflows

  if (lsbPos >= 0 && e >= 0)
    return 0; // lsbPos + e >= 0, d is integer

  if (lsbPos < 0 && e < 0)
    if (lsbPos < INT_MIN - e)
      return -2; // d isn't integer and needs too many fractional digits

  if ((lsbPos += e) >= 0)
    return 0; // d is integer

  if (lsbPos == INT_MIN && -INT_MAX != INT_MIN)
    return -2; // d isn't integer and needs too many fractional digits

  return -lsbPos;
}

const double testData[] =
{
  0,
  1, // 2 ^ 0
  0.5, // 2 ^ -1
  0.25, // 2 ^ -2
  0.125,
  0.0625, // ...
  0.03125,
  0.015625,
  0.0078125, // 2 ^ -7
  1.0/256, // 2 ^ -8
  1.0/256/256, // 2 ^ -16
  1.0/256/256/256, // 2 ^ -24
  1.0/256/256/256/256, // 2 ^ -32
  1.0/256/256/256/256/256/256/256/256, // 2 ^ -64
  3.14159265358979323846264338327950288419716939937510582097494459,
  0.1,
  INFINITY,
#ifdef NAN
  NAN,
#endif
  DBL_MIN
};

int main(void)
{
  unsigned i;
  for (i = 0; i < sizeof(testData) / sizeof(testData[0]); i++)
  {
    int digits = FractionalDigits(testData[i]);
    assert(digits >= 0);
    printf("%f %e %.*f\n", testData[i], testData[i], digits, testData[i]);
  }
  return 0;
}

Output (ideone):

0.000000 0.000000e+00 0
1.000000 1.000000e+00 1
0.500000 5.000000e-01 0.5
0.250000 2.500000e-01 0.25
0.125000 1.250000e-01 0.125
0.062500 6.250000e-02 0.0625
0.031250 3.125000e-02 0.03125
0.015625 1.562500e-02 0.015625
0.007812 7.812500e-03 0.0078125
0.003906 3.906250e-03 0.00390625
0.000015 1.525879e-05 0.0000152587890625
0.000000 5.960464e-08 0.000000059604644775390625
0.000000 2.328306e-10 0.00000000023283064365386962890625
0.000000 5.421011e-20 0.0000000000000000000542101086242752217003726400434970855712890625
3.141593 3.141593e+00 3.141592653589793115997963468544185161590576171875
0.100000 1.000000e-01 0.1000000000000000055511151231257827021181583404541015625
inf inf inf
nan nan nan
0.000000 2.225074e-308 0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002225073858507201383090232717332404064219215980462331830553327416887204434813918195854283159012511020564067339731035811005152434161553460108856012385377718821130777993532002330479610147442583636071921565046942503734208375250806650616658158948720491179968591639648500635908770118304874799780887753749949451580451605050915399856582470818645113537935804992115981085766051992433352114352390148795699609591288891602992641511063466313393663477586513029371762047325631781485664350872122828637642044846811407613911477062801689853244110024161447421618567166150540154285084716752901903161322778896729707373123334086988983175067838846926092773977972858659654941091369095406136467568702398678315290680984617210924625396728515625

You can see that π and 0.1 are only true up to 15 decimal digits and the rest of the digits show what the numbers got really rounded to since these numbers cannot be represented exactly in a binary floating-point format.

You can also see that DBL_MIN, the smallest positive normalized double value, has 1022 digits in the fractional part and of those there are 715 significant digits.

Possible issues with this solution:

  • Your compiler’s printf() functions do not support %a or do not correctly print all digits requested by the precision (this is quite possible).
  • Your computer uses non-binary floating-point formats (this is extremely rare).

More Related Contents:

  1. Printf width specifier to maintain precision of floating-point value
  2. Correct format specifier for double in printf
  3. What happens to a float variable when %d is used in a printf?
  4. Why is printf round floating point numbers up?
  5. How does printf and co differentiate between float and double
  6. What precisely does the %g printf specifier mean?
  7. C printf using %d and %f
  8. What does %d mean? Why d and not another letter? [closed]
  9. Why does >= operator work as > operator in C when compared with 0.9 [duplicate]
  10. printf(“%p”) and casting to (void *)
  11. How to get the sign, mantissa and exponent of a floating point number
  12. Printf variable number of decimals in float
  13. How to format a number using comma as thousands separator in C?
  14. What does “%.*s” mean in printf?
  15. printf specify integer format string for float
  16. Simd matmul program gives different numerical results
  17. How to control the number of exponent digits after ‘e’ in C printf %e?
  18. C: Casting minimum 32-bit integer (-2147483648) to float gives positive number (2147483648.0)
  19. Are floating point operations in C associative?
  20. Unsigned values in C
  21. How to combine two 32-bit integers into one 64-bit integer?
  22. Centering strings with printf()
  23. What is the conversion specifier for printf that formats a long?
  24. Using printf with a non-null terminated string
  25. printf string, variable length item
  26. Does C have a string type? [closed]
  27. Variadic function (va_arg) doesn’t work with float?
  28. printf rounding behavior for doubles
  29. sin, cos, tan and rounding error
  30. What is the difference between %f and %lf in C?

Leave a Comment