Function to print a double, in decimal - exactly

Question

Rarely is the exact decimal 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'

Edit: Added output for FLT_MAX, FLT_MIN, -FLT_TRUE_MIN for reference.

 3.4028234663852886e+38   '340282346638528859811704183484516925440.0'
 1.1754943508222875e-38   '0.000000000000000000000000000000000000011754943508222875079687365372222456778186655567720875215087517062784172594547271728515625'
-1.4012984643248171e-45   '-0.00000000000000000000000000000000000000000000140129846432481707092372958328991613128026194187651577175706828388979108268586060148663818836212158203125'

One would imagine that FLT_RADIX == 10 should be easy to support, too... — Toby Speight
– Toby Speight, Commented Aug 31, 2021 at 7:26
@TobySpeight True, but hard to test. IIRC, the latest platform with FLT_RADIX == 10 is over 15 years ago. C23 may introduce a new set of FP: decimal types, distinct from float, double, etc. — chux
– chux, Commented Aug 31, 2021 at 13:20

Toby Speight · Accepted Answer · 2019-02-01 14:18:21Z

6

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.

edited Feb 1, 2019 at 14:18

answered Feb 1, 2019 at 13:49

Toby Speight

88.7k14 gold badges104 silver badges327 bronze badges

\$\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
– chux

2019-02-01 13:57:31 +00:00
Commented Feb 1, 2019 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
– chux

2019-02-01 14:04:53 +00:00
Commented Feb 1, 2019 at 14:04
\$\begingroup\$ Last bit about commenting #define AD_F is spot on. \$\endgroup\$

chux
– chux

2019-02-01 14:06:42 +00:00
Commented Feb 1, 2019 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
– Toby Speight

2019-02-01 14:09:15 +00:00
Commented Feb 1, 2019 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
– chux

2019-02-01 14:15:28 +00:00
Commented Feb 1, 2019 at 14:15

Add a comment |

Stack Exchange Network

Function to print a double, in decimal - exactly

1 Answer 1

You must log in to answer this question.

Linked

Hot Network Questions

Function to print a double, in decimal - exactly

1 Answer 1

You must log in to answer this question.

Linked

Related

Hot Network Questions