6
\$\begingroup\$

Rarely is the exact value of a double needed to be printed and only its leading significant digits, after rounding, are needed.

It is a curiosity to see the exact value of a double as all finite double are exact. (Even if the math used on them generates mathematical approximation.)

An exact decimal output can be used to evaluate rounded outputs. Example


Below I seek a general review of print_double() and its support functions with emphasis on portability, assuming FLT_RADIX == 2, though not necessarily binary64.

Efficiency is a lesser concern.

/*
 * print_double.c
 * chux 2019
 */

#include <assert.h>
#include <float.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

_Static_assert(FLT_RADIX == 2, "TBD code for non-binary FP");

// Max whole number decimal digits in a `double`.
#define AD_I (DBL_MAX_10_EXP + 1)
// Max fractional decimal digits in a `double`.
#define AD_F (DBL_MANT_DIG - DBL_MIN_EXP)

/*
 * Managed array of decimal digits
 * Code first uses the middle of the array and works its way out to the ends
 * Most significant digits in highest indexed elements.
 */
typedef struct {
  unsigned char digit[AD_F + AD_I];
  size_t flength;
  size_t ilength;
} ad;

/*
 * Initialize 2 digits: 0.0
 */
static void ad_zero(ad *a) {
  a->ilength = 1;
  a->flength = 1;
  a->digit[AD_F] = 0;
  a->digit[AD_F-1] = 0;
}

/*
 *  a = a*FLT_RADIX + carry
 */
static int ad_mul(ad *a, unsigned carry) {
  size_t msd = AD_F + a->ilength;
  size_t lsd = AD_F;
  for (size_t i = lsd; i < msd; i++) {
    carry += (unsigned char) (a->digit[i] * FLT_RADIX);
    a->digit[i] = (unsigned char) (carry % 10);
    carry /= 10;
  }
  if (carry) {
   //printf("xx %zu\n", a->ilength);
    a->digit[AD_F + a->ilength++] = (unsigned char) carry;
    assert(carry / 10 == 0);
  }
  return 0;
}

/*
 * Divide by FLT_RADIX,  a /= FLT_RADIX
 */
static void ad_div(ad *a) {
  size_t msd = AD_F + a->ilength;
  size_t lsd = AD_F - a->flength;
  unsigned carry = 0;
  for (size_t i = msd; i > lsd;) {
    i--;
    carry = carry * 10u + a->digit[i];
    a->digit[i] = (unsigned char) (carry / FLT_RADIX);
    carry %= FLT_RADIX;
  }
  if (a->ilength > 1 && a->digit[msd - 1] == 0) {
    a->ilength--;
  }
  if (carry) {
    carry = carry * 10u;
    a->flength++;
    a->digit[AD_F - a->flength] = (unsigned char) (carry / FLT_RADIX);
    carry %= FLT_RADIX;
    assert(carry == 0);
  }
}

/*
 * Print ad
 */
static void ad_print(const ad *a) {
  size_t msd = AD_F + a->ilength;
  size_t lsd = AD_F - a->flength;
  for (size_t i = msd; i > lsd;) {
    printf("%d", a->digit[--i]);
    if (i == AD_F) {
      putchar('.');
    }
  }
}

/*
 * Print the exact value of double
 */
void print_double(double d) {
  if (!isfinite(d)) {
    printf("%f", d);
    return;
  }
  if (signbit(d)) {
    d = -d;
    putchar('-');
  }

  // Array to hold all the digits.    
  ad a;
  ad_zero(&a);

  int expo;
  d = frexp(d, &expo);
  while (d > 0) {
    expo--;
    double ipart;
    d = modf(d, &ipart) * FLT_RADIX;
    ad_mul(&a, (unsigned char) ipart);
  }
  expo++;
  while (expo > 0) {
    expo--;
    ad_mul(&a, 0);
  }
  if (expo < 0) {
    while (expo < 0) {
      expo++;
      ad_div(&a);
    }
  }
  ad_print(&a);
}

Sample usage

#include <float.h>
#include <stdio.h>

// Usually I'd put this in a .h file.
void print_double(double d);  

void print_double_wrap(double d) {
  // Print via printf()
  printf("% -*.*e '", DBL_DECIMAL_DIG + 8, DBL_DECIMAL_DIG - 1, d);

  print_double(d);

  puts("'");
}

int main(void) {
  print_double_wrap(0.0 / 0.0);
  print_double_wrap(1.0 / 0.0);
  print_double_wrap(-1.0 / 0.0);
  print_double_wrap(0);
  print_double_wrap(-0.0);
  print_double_wrap(1);
  print_double_wrap(123);
  print_double_wrap(1234);
  print_double_wrap(1e10);
  print_double_wrap(-1e20);
  print_double_wrap(1e30);
  print_double_wrap(1e300);
  print_double_wrap(123.5);
  print_double_wrap(0.01);
  print_double_wrap(DBL_MAX);
  print_double_wrap(DBL_MIN);
  print_double_wrap(-DBL_TRUE_MIN);
}

Output

-nan                      '-nan'
 inf                      'inf'
-inf                      '-inf'
 0.0000000000000000e+00   '0.0'
-0.0000000000000000e+00   '-0.0'
 1.0000000000000000e+00   '1.0'
 1.2300000000000000e+02   '123.0'
 1.2340000000000000e+03   '1234.0'
 1.0000000000000000e+10   '10000000000.0'
-1.0000000000000000e+20   '-100000000000000000000.0'
 1.0000000000000000e+30   '1000000000000000019884624838656.0'
 1.0000000000000001e+300  '1000000000000000052504760255204420248704468581108159154915854115511802457988908195786371375080447864043704443832883878176942523235360430575644792184786706982848387200926575803737830233794788090059368953234970799945081119038967640880074652742780142494579258788820056842838115669472196386865459400540160.0'
 1.2350000000000000e+02   '123.5'
 1.0000000000000000e-02   '0.01000000000000000020816681711721685132943093776702880859375'
 1.7976931348623157e+308  '179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.0'
 2.2250738585072014e-308  '0.00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000002225073858507201383090232717332404064219215980462331830553327416887204434813918195854283159012511020564067339731035811005152434161553460108856012385377718821130777993532002330479610147442583636071921565046942503734208375250806650616658158948720491179968591639648500635908770118304874799780887753749949451580451605050915399856582470818645113537935804992115981085766051992433352114352390148795699609591288891602992641511063466313393663477586513029371762047325631781485664350872122828637642044846811407613911477062801689853244110024161447421618567166150540154285084716752901903161322778896729707373123334086988983175067838846926092773977972858659654941091369095406136467568702398678315290680984617210924625396728515625'
-4.9406564584124654e-324  '-0.000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000004940656458412465441765687928682213723650598026143247644255856825006755072702087518652998363616359923797965646954457177309266567103559397963987747960107818781263007131903114045278458171678489821036887186360569987307230500063874091535649843873124733972731696151400317153853980741262385655911710266585566867681870395603106249319452715914924553293054565444011274801297099995419319894090804165633245247571478690147267801593552386115501348035264934720193790268107107491703332226844753335720832431936092382893458368060106011506169809753078342277318329247904982524730776375927247874656084778203734469699533647017972677717585125660551199131504891101451037862738167250955837389733598993664809941164205702637090279242767544565229087538682506419718265533447265625'
\$\endgroup\$
1
\$\begingroup\$

An interesting curiosity, as you say. Because 10 is an exact multiple of 2, all binary fractions have an exact, non-repeating decimal representation (but not vice versa). That means we're guaranteed to terminate.

The code compiles cleanly with an aggressive set of warnings enabled, and Valgrind is completely happy with the test program. But I expect you already know that.


There's very little I would change, but I did notice a redundancy here:

  if (expo < 0) {
    while (expo < 0) {
      expo++;
      ad_div(&a);
    }
  }

The outer if is pointless (and a good compiler will ignore it), so it's just useless clutter.


On the style side, I'd move the mutation of i from the body of these for loops into the loop-control part (inside the ( )), so that the control variable is constant within the body of each iteration:

  for (size_t i = msd; i > lsd;) {
    i--;
    //...
  }
  for (size_t i = msd; i > lsd;) {
    printf("%d", a->digit[--i]);
    //...
  }

to become, respectively:

  for (size_t i = msd; i-- > lsd;) {
    //...
  }
  for (size_t i = msd; i-- > lsd;) {
    printf("%d", a->digit[i]);
    //...
  }

That makes no difference to the functionality, but is less surprising and makes it easier to reason about the code. (If you want to be cute, you can write those operators together without a space, aka the infamous "goes to" operator, -->.)


One line that might need more comments is this one:

// Max fractional decimal digits in a `double`.
#define AD_F (DBL_MANT_DIG - DBL_MIN_EXP)

As we all know, DBL_MANT_DIG and DBL_MIN_EXP are defined in terms of FLT_RADIX, so it looks like an error to use these to infer the number of decimal digits required. A small amount of mathematical reasoning shows that each added bit requires one more decimal digit for the representation; I recommend that you summarise that in the comment to show that it's not a mistake.

\$\endgroup\$
  • \$\begingroup\$ Yes if (expo < 0) { while (expo < 0) { simplification good idea. Code came from former code that had a if (expo < 0) { putchar('.'); /* or other do something interesting here */ while (expo < 0) {. \$\endgroup\$ – chux Feb 1 at 13:57
  • \$\begingroup\$ Second style idea: Solitary i-- comes from a noted tendency to reduce use of post-fix notation in line (as in suggested i-- > lsd). It seems the newer crop of coders tend to have trouble with the latter. But I did miss the fine goto operator opportunity. Could have coded it that way just to encourage feedback. \$\endgroup\$ – chux Feb 1 at 14:04
  • \$\begingroup\$ Last bit about commenting #define AD_F is spot on. \$\endgroup\$ – chux Feb 1 at 14:06
  • \$\begingroup\$ As I said, it's a style/opinion point about i-- in the body - my preference is not to modify within the body unless it's natural for the step to be part-way through (and then consider commenting it). But I'm open to other opinions on that. \$\endgroup\$ – Toby Speight Feb 1 at 14:09
  • \$\begingroup\$ The other common reason for for (i = msd; i > lsd;) { i--; vs. for (i = msd; i-- > lsd;) { is when i defined and used outside the loop and the decrement must not be taken on loop exit. \$\endgroup\$ – chux Feb 1 at 14:15

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.