Summary
On Pelles C 14.50 (x64), the conversion of a double to decimal text produces
digits that are correct to about fifteen significant figures. Digits
requested beyond that are produced, but they are wrong, not merely
imprecise.
The practical consequence is that the round-trip guarantee implied by
DECIMAL_DIG does not hold:
-------------------------------- ------------------------------------------
DECIMAL_DIG declared by 17
<float.h>
DBL_DECIMAL_DIG declared by 17
<float.h>
finite normal doubles printed 48,191 of 100,000 -- 48.19%
with %.17g that do not read back
as the same value
-------------------------------- ------------------------------------------
The same test program, built with another toolchain's runtime on the same
machine and run on the same values, reports zero failures.
The input side is not affected: strtod reads a correctly rounded
seventeen-digit decimal back to the exact bit pattern. Hexadecimal output,
%a, is not affected either. The fault appears to lie solely in the
binary-to-decimal conversion used by the printf family.
A single self-contained test program is attached. It needs no data files, no
arguments and no other libraries, and it runs in under a second.
Environment
------------------------------------- -------------------------------------
Compiler Pelles ISO C Compiler, Version
14.50.0, x64
Switches -Tx64-coff -std:C17 -Ze -Zx
Link crt64.lib kernel32.lib, console
subsystem
Optimisation irrelevant: identical output
unoptimised and at -Ot, -Os and -Ox
FLT_EVAL_METHOD 0
------------------------------------- -------------------------------------
The optimiser was eliminated first, by rebuilding at every optimisation
level and comparing the output byte for byte. All four builds agree with one
another and all four disagree with the correctly rounded digits, so this is
a runtime library matter and not a code generation one.
How it was found
The work in progress was unrelated: making a C++ image processing library
usable from Pelles C by compiling its C++ once into an object file and
presenting a pure C interface, so that a Pelles C program can call it
without a DLL. In that arrangement the C++ code is compiled by another
compiler, but every C function it calls -- malloc, memcpy, fopen, printf,
cos, sin -- is resolved at the final link by Pelles C's own crt64.lib. That
is the point of the design: one C runtime, one heap, one set of file
handles.
To check that the C interface altered nothing, thirteen image processing
computations were written twice: once in C, calling through the interface,
and once in C++, calling the library directly. The two sets of output images
were compared byte for byte. Twelve of the thirteen were identical. One
differed in a single pixel, by one grey level.
Chasing that one pixel is what turned this up. The chain of elimination was:
1. The optimiser was ruled out, as described above.
2. Both computations were instrumented to compare their floating point in
memory, stage by stage, before anything was cast to a byte. All nine
stages were bit-identical, so the interface was not altering anything.
3. The small amount of arithmetic in the caller was compiled by both
compilers from the same source and compared: 51,062 floats,
bit-identical.
4. What remained was the mathematical library. Because the design routes
every C call to crt64.lib, the transform in question was using Pelles
C's cos and sin where the C++ comparison used the other toolchain's.
Measured over 4,096 angles: 2,061 of 8,192 values differ, by at most one
unit in the last place. That is not a fault in anything -- neither
library claims correctly rounded transcendentals -- and it fully
explains one pixel in 65,536 landing on the boundary between grey 31 and
grey 32.
The present report arises from step 4. While printing the cos values for
comparison, the two runtimes printed different decimal text for identical
bit patterns. That is not a question about mathematics; it is a question
about printf, and it turned out to be a real one.
The minimal test
decimal_dig_test.c, attached, does two things.
Part one takes three doubles as exact IEEE-754 bit patterns, so that no
decimal parsing happens before the test, and compares printf's output
against the correctly rounded digits. Those expected strings were obtained
independently, by arbitrary-precision arithmetic, and are written into the
program as literals.
Part two generates 100,000 finite normal doubles from a four-line xorshift
generator -- so the sample is identical on every machine and every run --
prints each with %.*g at DECIMAL_DIG, reads it back with strtod, and
compares the bits.
strtod is used as the reader because it was verified separately, against
exact bit patterns, to be correct on this implementation. If it were not,
part one would fail in the same direction and the two parts would disagree.
Results
Part one, on Pelles C 14.50
bits 3feffff9606a6fbc
%.17g expected 0.99999684172825409 got 0.999996841728254 <-- DIFFERS
%.20f expected 0.99999684172825409334 got 0.99999684172825400310 <-- DIFFERS
bits 3fefe9cdad01883a
%.17g expected 0.99729045667869021 got 0.99729045667869013 <-- DIFFERS
%.20f expected 0.99729045667869020697 got 0.99729045667869013413 <-- DIFFERS
bits 400921fb54442d18
%.17g expected 3.1415926535897931 got 3.1415926535897932 <-- DIFFERS
%.20f expected 3.14159265358979311600 got 3.14159265358979318989 <-- DIFFERS
Two features are worth separating.
The digits beyond the fifteenth are wrong. In the first value the correct
digits are ...254093 34, and what is produced is ...254003 10. This is not a
rounding difference at the last digit; the sixteenth digit onwards is a
different number.
Fewer digits are produced than were requested. %.17g on the first value
returns 0.999996841728254, which is fifteen significant figures, not
seventeen. In the 100,000-value sample this happens in 9.84% of cases.
Part two, on Pelles C 14.50
values tested 100000
failed to read back as the same double 48191 (48.19%)
printed with fewer digits than requested 9836 (9.84%)
The same program, another runtime, same machine, same values
values tested 100000
failed to read back as the same double 0 (0.00%)
How it varies with the precision requested
Over 4,000 finite normal doubles spanning the whole exponent range, counting
strings that differ from the correctly rounded conversion:
conversion Pelles C 14.50 control
------------ ---------------- ---------
%.15g 3.33% 0.00%
%.16g 29.93% 0.00%
%.17g 85.47% 0.00%
%.18g 97.75% 0.00%
And the round-trip failures over the same sample:
conversion Pelles C 14.50 control
------------ ---------------- ---------
%.16g 61.88% 44.77%
%.17g 48.10% 0.00%
%.18g 47.45% 0.00%
The %.16g row is expected to fail on both, since sixteen digits cannot name
every double; it is included as a control on the method. The %.17g row is
the one that matters, and the two columns should both be zero.
What is not affected
These were tested and are correct, which may help to localise the fault.
strtod is correct. Five exact seventeen-digit decimals were parsed and each
produced the expected bit pattern.
%a is correct. The hexadecimal form was checked across every binary exponent
from 2^-1070 to 2^1023, 2,094 values, and separately on DBL_MAX, DBL_MIN,
the largest subnormal and the smallest subnormal. All were exact. %a is
therefore a sound workaround for anyone who needs to write a double as text
without losing it.
The optimiser is not involved. Unoptimised, -Ot, -Os and -Ox all produce
identical output.
Why it matters
<float.h> on this implementation declares DECIMAL_DIG = 17 and
DBL_DECIMAL_DIG = 17. The C standard describes DECIMAL_DIG as the number of
decimal digits such that any floating value may be converted to decimal and
back with no change of value (C11, 5.2.4.2.2, paragraph 12). Seventeen is
the correct value for IEEE-754 double. The declaration is therefore right;
what is missing is the conversion that would make it true.
The consequence is quiet, which is what makes it serious. Any program that
writes double data as text -- a configuration file, a comma-separated
export, a log, a checkpoint to be resumed from, a scientific data file, a
coordinate in a header -- and that writes more than fifteen significant
figures in the belief that it is preserving the value, is altering the value
at the moment of writing, and nothing in the program will say so. The loss
is at the level of one or two units in the last place, which is exactly the
level at which it survives casual inspection and then accumulates.
A false lead, recorded so that nobody repeats it
An earlier version of the measuring program appeared to show %a producing a
grossly wrong result for DBL_MAX: 0x8.181818181818p-255 in place of
0x1.fffffffffffffp+1023. That was a fault in the measuring program, not in
Pelles C. The same program had earlier printed DBL_MAX with %.20f into a
128-byte buffer, and that conversion needs more than 330 characters. The
buffer overflow corrupted what followed.
It is recorded here because it was nearly reported as a defect, and because
the only reason it was not is that the result was checked a second time by a
different route.
How to reproduce
pocc -Tx64-coff -std:C17 -Ze -Zx decimal_dig_test.c
polink -subsystem:console -machine:x64 decimal_dig_test.obj ^
crt64.lib kernel32.lib -out:decimal_dig_test.exe
decimal_dig_test
The program returns 0 if both parts pass and 1 otherwise. For a comparison,
the identical source builds unchanged with any other C compiler.
Suggested area to look at
The pattern -- correct to about fifteen significant figures, then wrong
rather than absent, together with %g sometimes yielding fewer significant
digits than requested -- is consistent with the conversion being carried out
through an intermediate of limited precision, with the remaining digits
generated from a value that has already lost them, rather than by an exact
arbitrary-precision method such as Steele and White's Dragon4 or Grisu with
an exact fallback.
strtod on the same implementation is exact, so the necessary
arbitrary-precision machinery may already be present on the input side.
File attached 'decimal_dig_test.c':
/* ===========================================================================
* decimal_dig_test.c
* ---------------------------------------------------------------------------
* A self-contained test of one thing only: whether printf can write a double
* with DECIMAL_DIG significant digits such that strtod reads back the same
* value.
*
* <float.h> declares
*
* DECIMAL_DIG the number of decimal digits needed to hold any
* floating value of the widest supported type so that
* converting it to decimal and back changes nothing
* DBL_DECIMAL_DIG the same, for double
*
* The C standard describes DECIMAL_DIG (C11 5.2.4.2.2, paragraph 12) as the
* number of digits for which that round trip is guaranteed. This program
* simply takes the implementation at its word and checks.
*
* IT NEEDS NOTHING. No data files, no libraries beyond the standard ones,
* no command line arguments. It runs in under a second.
*
* HOW IT AVOIDS TESTING ITSELF
*
* Part one sets its values from EXACT BIT PATTERNS, so no decimal
* parsing takes place before the test. Each expected string was
* obtained independently, from arbitrary-precision arithmetic, and is
* written into this file as a literal.
*
* Part two uses strtod to read the printed text back. That is safe
* here: strtod was checked separately against exact bit patterns and is
* correct on this implementation. If it were not, part one would fail
* in the same way and the two parts would disagree.
*
* BUILD, Pelles C:
* pocc -Tx64-coff -std:C17 -Ze -Zx decimal_dig_test.c
* polink -subsystem:console -machine:x64 decimal_dig_test.obj \
* crt64.lib kernel32.lib -out:decimal_dig_test.exe
*
* BUILD, any other C compiler, for comparison:
* cc -std=c11 -O1 decimal_dig_test.c -o decimal_dig_test
* ======================================================================== */