Function to print a double - exactly

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.
// Max fractional decimal digits in a double.

/*
* 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 {
size_t flength;
size_t ilength;

/*
* Initialize 2 digits: 0.0
*/
a->ilength = 1;
a->flength = 1;
}

/*
*  a = a*FLT_RADIX + carry
*/
size_t msd = AD_F + a->ilength;
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;
}

/*
*/
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);
}
if (a->ilength > 1 && a->digit[msd - 1] == 0) {
a->ilength--;
}
if (carry) {
carry = carry * 10u;
a->flength++;
assert(carry == 0);
}
}

/*
*/
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]);
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.

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


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'


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++;
}
}


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.

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.
• 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) {. – chux Feb 1 at 13:57
• 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. – chux Feb 1 at 14:04
• Last bit about commenting #define AD_F is spot on. – chux Feb 1 at 14:06
• 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. – Toby Speight Feb 1 at 14:09
• 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. – chux Feb 1 at 14:15