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

Question

I'm new to C and learning out of a book / off the internet. I'm trying to write a function that I can pass any double to and get returned an int to be used in a printf("%.*lf" ... statement such that the int returned will neither diminish precision nor produce trailing zeros.

I have a working function but it is pretty big since it is written for readability and all commented up.

To summarize the function, I count how many divisions by 10 it takes to get the double in the range 10 > d >= 0, take only the fractional part and put it into a string with n decimal places where n = 15 - number_of_digits_left_of_decimal (I read that type double can only keep track of 15 digits), check the string from right to left for trailing zeroes and keep count, and, finally, return an int that represents the number of non-zero digits right of the decimal.

Is there an easier way? Thanks.

int get_number_of_digits_after_decimal(double d)
{
  int i = 0;      /* sometimes you need an int */
  int pl = 0;     /* precision left = 15 - sigfigs */
  int sigfigs = 1; /* the number of digits in d */
  char line[20];  /* used to find last non-zero digit right of the decimal place */
  double temp;    /* a copy of d used for destructive calculations */

  /* find digits to right of decimal */
  temp = d;
  while(sigfigs < 15)
  {
    if(temp < 0)
      temp *= -1;
    if(temp < 10)
      break;
    temp /= 10;
    ++sigfigs;
  }
  /* at this point 10 > temp >= 0
  * decrement temp unitl 1 > temp >=0 */
  while(temp > 1)
  {
    --temp;
  }
  if(temp == 0)
    return(0);
  pl = 15 - sigfigs;   /* if n digits left of decimal, 15-n to right */
  switch(pl)
  {
  case 14:
    sprintf(line, "%.14lf", d);
    break;
  case 13:
    sprintf(line, "%.13lf", d);
    break;
  case 12:
    sprintf(line, "%.12lf", d);
    break;
  case 11:
    sprintf(line, "%.11lf", d);
    break;
  case 10:
    sprintf(line, "%.10lf", d);
    break;
  case 9:
    sprintf(line, "%.9f", d);
    break;
  case 8:
    sprintf(line, "%.8lf", d);
    break;
  case 7:
    sprintf(line, "%.7lf", d);
    break;
  case 6:
    sprintf(line, "%.6lf", d);
    break;
  case 5:
    sprintf(line, "%.5lf", d);
    break;
  case 4:
    sprintf(line, "%.4lf", d);
    break;
  case 3:
    sprintf(line, "%.3lf", d);
    break;
  case 2:
    sprintf(line, "%.2lf", d);
    break;
  case 1:
    sprintf(line, "%.1lf", d);
    break;
  case 0:
    return(0);
    break;
  }
  i = (strlen(line) - 1); /* last meaningful digit char */
  while(1) /* start at end of string, move left checking for first non-zero */
  {
    if(line[i] == '0') /* if 0 at end */
    {
      --i;
      --pl;
    }
    else
    {
      break;
    }
  }
  return(pl);
}

Answer 1

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⁺⁰ = 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).

Answer 2

The first thing I notice is that you're dividing temp by 10 and that is causing a loss of precision.

Not to shut you down or discourage you from trying again, but a correct implementation of this is considerably more involved than what you've shown.

Guy L. Steele and Jon L. White wrote a paper called "How to print floating-point numbers accurately" that details some of the pitfalls and presents a working algorithm for printing floating-point numbers. It's a good read.