For example, consider the 6-digit decimal number 7.23496. Rounded directly to 3 digits — using round-to-nearest, round half to even rounding — it’s 7.23; rounded first to 5 digits (7.2350) and then to 3 digits it’s 7.24. The value 7.24 is incorrect, reflecting a double rounding error.

In a computer, double rounding occurs in binary floating-point arithmetic; the typical example is a calculated result that’s rounded to fit into an x87 FPU extended precision register and then rounded again to fit into a double-precision variable. But I’ve discovered another context in which double rounding occurs: conversion from a decimal floating-point literal to a single-precision floating-point variable. The double rounding is from full-precision binary to double-precision, and then from double-precision to single-precision.

In this article, I’ll show example conversions in C that are tainted by double rounding errors, and how attaching the ‘f’ suffix to floating-point literals prevents them — in gcc C at least, but not in Visual C++!

Double Rounding Error Example 1

The following C program demonstrates a double rounding error that occurs when a decimal floating-point literal is converted to a float by the compiler. The program converts the same decimal number three times: first to a double d; then to a float f, using the ‘f’ (float) suffix; and then to a float fd, without using the ‘f’ suffix.

#include <stdio.h>

int main (void)
{
 double d;
 float f, fd; 

 d =  0.5000000894069671353303618843710864894092082977294921875;
 f =  0.5000000894069671353303618843710864894092082977294921875f;
 fd = 0.5000000894069671353303618843710864894092082977294921875;

 printf ("d =  %a\n",d);
 printf ("f =  %a\n",f);
 printf ("fd = %a\n",fd);
}

Output on Linux

Here’s the output when the program is compiled and run on Linux using gcc:

d =  0x1.000003p-1
f =  0x1.000002p-1
fd = 0x1.000004p-1

Translated from hexadecimal floating-point constants to binary, the output is

d =  0.1000000000000000000000011
f =  0.100000000000000000000001
fd = 0.10000000000000000000001

To interpret the output, first consider the full-precision binary equivalent of the decimal number in the program (it is a terminating binary fraction):

0.1000000000000000000000010111111111111111111111111111111

Now let’s look at the three variables in turn:

• d is the correctly rounded conversion of the full-precision binary value to double-precision — 53 significant bits. The full-precision binary value is rounded up, since the value of bits 54 and beyond is greater than 1/2 ULP:

  0.1000000000000000000000010111111111111111111111111111111

• f is the correctly rounded conversion of the full-precision binary value to single-precision — 24 significant bits. The full-precision binary value is rounded down, since the value of bits 25 and beyond is less than 1/2 ULP:

  0.1000000000000000000000010111111111111111111111111111111

• fd is an incorrectly rounded conversion to single-precision. It’s 1 ULP above the correctly rounded value, due to a double rounding error. The compiler converted it by first rounding it up to 53 bits, giving 0.1000000000000000000000011, and then rounding that up to 24 bits, giving 0.10000000000000000000001. In other words, the first rounding causes a carry to propagate to bit 25, creating a halfway case that gets resolved by rounding up to nearest even.

Summary

With the ‘f’ suffix on a decimal literal, gcc converts a full-precision binary value directly to single-precision, avoiding double rounding. Without the ‘f’ suffix, gcc converts the full-precision binary value to double-precision and then to single-precision; this indirect conversion from full-precision to single-precision is double rounding, and in this case results in a double rounding error.

Output on Windows

Here’s the output when the program is compiled and run on Windows using Visual C++:

d =  0x1.000003p-1
f =  0x1.000004p-1
fd = 0x1.000004p-1

Like with gcc, fd is doubly rounded; but Visual C++ doubly rounds f as well! It seems Visual C++ ignores the ‘f’ suffix, first rounding to double-precision and then rounding to single-precision. (I am presuming this is the case. It’s also possible that Visual C++ is honoring the ‘f’ suffix, only it’s doing the conversion incorrectly — due to reasons other than double rounding.)

Double Rounding Error Example 2

Here’s a similar example, except that it demonstrates a double rounding error that results in a value that is 1 ULP below the correctly rounded value:

#include <stdio.h>

int main (void)
{
 double d;
 float f, fd; 

 d =  0.5000000298023224154508881156289135105907917022705078125;
 f =  0.5000000298023224154508881156289135105907917022705078125f;
 fd = 0.5000000298023224154508881156289135105907917022705078125;

 printf ("d =  %a\n",d);
 printf ("f =  %a\n",f);
 printf ("fd = %a\n",fd);
}

Output on Linux

Here’s the output when the program is compiled and run on Linux using gcc:

d =  0x1.000001p-1
f =  0x1.000002p-1
fd = 0x1p-1

Translated from hexadecimal floating-point constants to binary, the output is

d =  0.1000000000000000000000001
f =  0.100000000000000000000001
fd = 0.1

The full-precision binary equivalent of the decimal number in the program is

0.1000000000000000000000001000000000000000000000000000001

Again, let’s look at the three variables:

• d is the correctly rounded conversion of the full-precision binary value to 53 significant bits. The full-precision binary value is rounded down, since the value of bits 54 and beyond is less than 1/2 ULP:

  0.1000000000000000000000001000000000000000000000000000001

• f is the correctly rounded conversion of the full-precision binary value to 24 significant bits. The full-precision binary value is rounded up, since the value of bits 25 and beyond is greater than 1/2 ULP:

  0.1000000000000000000000001000000000000000000000000000001

• fd is an incorrectly rounded conversion to single-precision. It’s 1 ULP below the correctly rounded value, due to a double rounding error. The compiler converted it by first rounding it down to 53 bits, giving 0.1000000000000000000000001, and then rounding that down to 24 bits, giving 0.1. In other words, the first rounding creates a halfway case that gets resolved by rounding down to nearest even.

As in example 1, use of the ‘f’ suffix prevents double rounding.

Output on Windows

Here’s the output when the program is compiled and run on Windows using Visual C++:

d =  0x1.000001p-1
f =  0x1.000000p-1
fd = 0x1.000000p-1

As in example 1, Visual C++ ignores the ‘f’ suffix, resulting in double rounding for both f and fd.

Summary

Floating-point literals are subject to double rounding when assigned to single-precision variables, resulting in incorrectly rounded decimal to floating-point conversions. If you’re using the gcc C compiler, you can avoid this by attaching the ‘f’ suffix to your literals.

A Note on the Examples

The examples above are admittedly contrived, but they serve to illustrate the potential for double rounding errors during conversion.

Please Try My Examples

If you are using a compiler other than Visual C++ or gcc, please try out my examples. In particular, I’d like to know if your compiler doubly rounds decimal literals when the ‘f’ suffix is used.

Bug Report

I wrote a bug report against Visual C++; this was Microsoft’s response:

“This is indeed an issue with the compiler but we regret that we cannot fix it for the next release due to its priority.”

References

Here is some background reading on double rounding:

• Handbook of Floating-Point Arithmetic (see pages 75-77).
• When Double Rounding is Odd (see the two pages titled “Double rounding going bad”).
• When Is Double Rounding Innocuous? (see the first two pages).

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Double Rounding Errors in Floating-Point Conversions

For example, consider the decimal number 33.75. It converts to a double with a sign field of 0, an exponent field of 10000000100, and a fraction field of 0000111000000000000000000000000000000000000000000000. The 0 in the sign field means it’s a positive number (1 would mean it’s negative). The value of 10000000100 in the exponent field, which equals 1028 in decimal, means the exponent of the power of two is 5 (the exponent field value is offset, or biased, by 1023). The fraction field represents the normalized binary fraction 1.0000111. Written in binary scientific notation — following the convention that the fraction is written in binary and the power of two is written in decimal — 33.75 equals 1.0000111 x 2⁵.

In this article, I’ll show you the C function I wrote to display a double in binary scientific notation. This function is useful, for example, when verifying that decimal to floating-point conversions are correctly rounded.

The Code

I wrote a function called print_double_binsci() that prints a double-precision floating-point number in binary scientific notation. It is based on a call to my function parse_double(), which isolates the three fields of a double.

I declared and defined this function in files I named binsci.h and binsci.c, respectively.

binsci.h

/***********************************************************/
/* binsci.h: Function to print an IEEE double-precision    */
/*           floating-point number in binary scientific    */
/*           notation                                      */
/***********************************************************/
void print_double_binsci(double d);

binsci.c

/***********************************************************/
/* binsci.c: Function to print an IEEE double-precision    */
/*           floating-point number in binary scientific    */
/*           notation                                      */
/***********************************************************/
#include <stdio.h>
#include "rawdouble.h"
#include "binsci.h"

void print_double_binsci(double d)
{
 int i;
 unsigned char sign_field;
 unsigned short exponent_field;
 unsigned long long fraction_field;

 //Isolate the three fields of the double
 parse_double(d,&sign_field,&exponent_field,&fraction_field);

 //Print a minus sign, if necessary
 if (sign_field == 1)
   printf("-");

 //Print the implicit leading 1 and radix point
 printf("1.");

 //Print the fraction field
 for (i=0; i<=51; i++)
   if ((fraction_field >> (51-i)) & 1)
     printf("1");
   else
     printf("0");

 //Print the power, taking into account exponent bias
 printf(" x 2^%d\n",exponent_field - 1023);
}

Notes

• All 53 bits are printed, including superfluous trailing 0s.
• Numbers that are not raised to a power are printed with the suffix “x 2^0”.
• Not-a-number (NaN) and infinity values are not handled.

Examples of Usage

I wrote a program, called binsciTest.c, that shows some example calls to print_double_binsci():

/***********************************************************/
/* binsciTest.c: Program to test printing of IEEE double   */
/*               precision floating-point numbers in       */
/*               binary scientific notation                */
/***********************************************************/
#include <stdio.h>
#include "binsci.h"

int main (void)
{
 printf("33.75 =\n");
 print_double_binsci(33.75);
 printf("\n");

 printf("0.1 =\n");
 print_double_binsci(0.1);
 printf("\n");

 printf("-0.6 =\n");
 print_double_binsci(-0.6);
 printf("\n");

 printf("3.518437208883201171875e13 =\n");
 print_double_binsci(3.518437208883201171875e13);
 printf("\n");

 printf("9214843084008499.0 =\n");
 print_double_binsci(9214843084008499.0);
 printf("\n");

 printf("30078505129381147446200.0 =\n");
 print_double_binsci(30078505129381147446200.0);
 printf("\n");

 printf("1777820000000000000001.0 =\n");
 print_double_binsci(1777820000000000000001.0);
 printf("\n");

 printf("0.3932922657273 =\n");
 print_double_binsci(0.3932922657273);

 return (0);
}

(Some of these examples were taken from my articles Incorrectly Rounded Conversions in Visual C++ and Incorrectly Rounded Conversions in GCC and GLIBC.)

I compiled and ran it on both Windows and Linux:

• On Windows, I built a project in Visual C++ with files binsci.c, binsci.h, binsciTest.c, rawdouble.c, and rawdouble.h, and compiled and ran it in there.
• On Linux, I compiled with “gcc binsciTest.c binsci.c rawdouble.c -o binsciTest” and then ran it with “./binsciTest”.

Output

33.75 =
1.0000111000000000000000000000000000000000000000000000 x 2^5

0.1 =
1.1001100110011001100110011001100110011001100110011010 x 2^-4

-0.6 =
-1.0011001100110011001100110011001100110011001100110011 x 2^-1

3.518437208883201171875e13 =
1.0000000000000000000000000000000000000000000000000001 x 2^45

9214843084008499.0 =
1.0000010111100110110011101100010101110111011000011001 x 2^53

30078505129381147446200.0 =
1.1001011110100011110001110010011100011011000000100000 x 2^74

1777820000000000000001.0 =
1.1000000110000000110101011011101011010010111000111110 x 2^70

0.3932922657273 =
1.1001001010111011001101010010110001000110001000111001 x 2^-2

Except for 33.75, which is exact, all the other examples are 53 significant bit approximations to the decimal numbers they stand in for (type ‘0.1’ into my decimal to binary converter and see for yourself).

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Displaying IEEE Doubles in Binary Scientific Notation

(Update 7/11/10: dtoa.c is now fixed — see below.)

Discovering the Problem

As a follow up to my article “Visual C++ and GLIBC strtod() Ignore Rounding Mode,” I asked Mark Dickinson — the developer who maintains Python’s copy of dtoa.c — why Python does not enable what I call directed conversions. Out of that discussion came Mark’s conjecture: the directed rounding code was not correct. He soon followed up with an example to prove it.

An Automated Way to Find Examples

I wanted to find more examples, but I wanted a C program to do it for me. I needed a way to generate examples whose correct outcomes were known (in my evaluations of round-to-nearest conversions, David Gay’s strtod() served this purpose). Inspired by Mark’s example, this is the algorithm I devised:

1. Generate a random decimal string that represents a number that is not an integer and is not dyadic.
2. Convert the decimal string to a double (I used David Gay’s strtod(), but the built-in strtod() would have been fine — correct rounding does not matter here).
3. Convert the double to a decimal string, giving the exact decimal representation of a floating-point number (I used David Gay’s g_fmt() function, since I’m running on Windows and printf() in Visual C++ can’t print all the decimal digits of a floating-point number).
4. Change the last digit of the new decimal string — which by the way is always 5 (this article indirectly explains why) — to either 4 or 6, depending on whether you want a number just below or just above the floating-point number.

Except for Mark’s example, the examples below were found by my C program implementation of this algorithm. I ran the program against a copy of dtoa.c most recently changed on “Tue Feb 2 23:05:34 MST 2010,” which I downloaded on May 27, 2010.

Examples: Decimal Numbers Just Above a Floating-Point Number

This section shows two examples that David Gay’s strtod() (hereafter referred to only as strtod()) rounds incorrectly; they are very close to, but just above (in absolute value), a floating-point number.

Example 1

This is Mark’s example, and it shows a positive number rounding to nearest in round toward positive infinity mode:

1.100000000000000088817841970012523233890533447265626

It converts to this binary number, which has a non-terminating binary fraction portion:

1.000110011001100110011001100110011001100110011001101000000000000000
00000000000000000000000000000000000000000000000000000000000000000000 000000000000000000000000000000000001…

Rounded toward positive infinity the result should be

1.0001100110011001100110011001100110011001100110011011

In hexadecimal “%a” notation, this equals 0x1.199999999999bp+0.

strtod() converts it incorrectly to 0x1.199999999999ap+0, which is one ULP below the correctly rounded value. 0x1.199999999999ap+0 is the floating-point number equal to decimal 1.100000000000000088817841970012523233890533447265625, which differs from the input decimal value only in the last decimal digit (5 vs 6).

Example 2

This example shows a negative number rounding to nearest in round toward negative infinity mode:

-0.91276999999999997026378650843980722129344940185546876

It converts to this non-terminating binary number:

-0.111010011010101101001011011100101100010100011001011110000000000000
0000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000001…

Rounded toward negative infinity the result should be -0x1.d35696e58a330p-1; strtod() converts it incorrectly to -0x1.d35696e58a32fp-1.

Examples: Decimal Numbers Just Below a Floating-Point Number

This section shows two examples that strtod() rounds incorrectly; they are very close to, but just below (in absolute value), a floating-point number.

Example 3

This example shows a negative number rounding to nearest in round toward positive infinity mode:

-266.240000000000009094947017729282379150390624

Rounded toward positive infinity it should be -0x1.0a3d70a3d70a3p+8; strtod() rounds it incorrectly to -0x1.0a3d70a3d70a4p+8.

Example 4

This example shows a positive number rounding to nearest in both round toward negative infinity and round toward zero modes:

8.2556288587679180024720432899523381023022507640626854730214745359262
45152950286865234374e-17

Rounded toward negative infinity and rounded toward zero it should be 0x1.7cb9433617c9bp-54; strtod() rounds it incorrectly — in both cases — to 0x1.7cb9433617c9cp-54.

Discussion

In my article “Visual C++ and GLIBC strtod() Ignore Rounding Mode” I said

The percentage of incorrect conversions for both Visual C++ and glibc ranged between 49.00% and 49.06%, suggesting they were not reacting to the rounding mode at all. (Why it wasn’t almost exactly 50% I don’t know.)

It might not have been 50% because, although I didn’t know it at the time, strtod()’s conversions weren’t 100% correct! (I’ll rerun those tests when dtoa.c is fixed.) (Update: the percentages are still near 49%, even with the fixed dtoa.c, so scratch that theory.)

Do You Use Directed Rounding?

If your project uses dtoa.c, could you let me know whether it enables directed conversions? That is, does it define the “Honor_FLT_ROUNDS” flag?

You should be aware that strtod() reacts to the rounding mode, although incorrectly, when Honor_FLT_ROUNDS is not defined.

Bug Report

I reported this to David Gay; I’ll keep you posted.

Update: dtoa.c is Fixed

David Gay updated dtoa.c with a fix, and it appears to work: all four examples above now round correctly, as did the millions of random examples I checked.

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Incorrect Directed Conversions in David Gay’s strtod()

Common practice considers a decimal number correctly rounded when the nearest of the two floating-point numbers is chosen (and when both are equally near, when the one with significant bit number 53 equal to 0 is chosen). This makes sense intuitively, and also reflects the default IEEE 754 rounding mode — round-to-nearest. However, there are three other IEEE 754 rounding modes, which allow for directed rounding: round toward positive infinity, round toward negative infinity, and round toward zero. For a conversion to be considered truly correctly rounded, it must honor all four rounding modes — whichever is currently in effect.

I evaluated the Visual C++ and glibc strtod() functions under the three directed rounding modes, like I did for round-to-nearest mode in my articles “Incorrectly Rounded Conversions in Visual C++” and “Incorrectly Rounded Conversions in GCC and GLIBC.”. What I discovered was this: they only convert correctly about half the time — pure chance! — because they ignore the rounding mode altogether.

Example: Correct Rounding in All Four Rounding Modes

The diagram below shows how correct conversion works in each of the four rounding modes; I used the decimal values 0.1 and -0.1 as examples:

How 0.1 and -0.1 Convert Correctly in Each of the Four Rounding Modes

The conversions — shown in hexadecimal “%a” notation — were done by David Gay’s strtod() function and also verified by hand. In terms of absolute value, there are two different results; here they are in binary:

• 0x1.9999999999999p-4 =

  1.1001100110011001100110011001100110011001100110011001 x 2^-4

• 0x1.999999999999ap-4 =

  1.100110011001100110011001100110011001100110011001101 x 2^-4

Example: Incorrect Rounding in Visual C++ and glibc

Visual C++ and glibc strtod() do four of the eight conversions incorrectly:

How 0.1 and -0.1 Convert Incorrectly in Visual C++ and glibc

0.1, in round toward zero and round toward negative infinity modes, converts incorrectly to 0x1.999999999999ap-4. -0.1, in round toward zero and round toward positive infinity modes, converts incorrectly to -0x1.999999999999ap-4. That is, Visual C++ and glibc are converting 0.1 and -0.1 to a single absolute value in every case: 0x1.999999999999ap-4. That, not coincidentally, is the nearest floating-point number.

Visual C++ and glibc Always Round to Nearest

To see what percentage of conversions Visual C++ and glibc got right, I tested strtod() with random decimal values — both positive and negative — in each of the three directed rounding modes. I compared their results to those computed by David Gay’s strtod(). The percentage of incorrect conversions for both Visual C++ and glibc ranged between 49.00% and 49.06%, suggesting they were not reacting to the rounding mode at all. (Why it wasn’t almost exactly 50% I don’t know.) To confirm this, I generated random decimal values and converted each in all four modes; in every case, the four conversions were the same.

Conclusion: Visual C++ and glibc are using arbitrary-precision arithmetic in their calculations, and their algorithms are hard-coded to perform round-to-nearest rounding.

(Update: I was assuming David Gay’s strtod() got the conversions right in the directed rounding modes, but it turns out there was a bug — now fixed — that caused some conversions to come out incorrectly rounded. In any case, my conclusion is unaffected, since ultimately I didn’t need David Gay’s strtod() to reach it.)

Enabling Rounding Mode Detection in David Gay’s strtod()

David Gay’s strtod() function does not honor rounding mode by default; it must be enabled, by building dtoa.c with the flag “Honor_FLT_ROUNDS”. This forces the inclusion of file fenv.h, which defines the function fegetround(). fegetround() is called to get the current processor rounding mode, which strtod() uses to round its results correctly.

To make this work in Visual C++ — which doesn’t have the file fenv.h — I had to simulate fegetround(). I did this by calling _controlfp_s(&cur_cw,0,0) and translating the rounding mode masks (I am using the x87 FPU).

Conversion Algorithm Uses Floating-Point Arithmetic

David Gay’s strtod() function uses floating-point as well as arbitrary-precision calculations in its conversion, which you can see by reading the source code. You can also see this by building dtoa.c without “Honor_FLT_ROUNDS” and running it in the different rounding modes. Although the conversions won’t be correct (in general), they will differ, meaning they were computed with the help of the processor’s floating-point unit.

Discussion

Compile Time Vs. Run Time Conversions

I’ve been discussing strtod(), which does conversions at run time. But conversions can be done at compile time too — using floating-point literals. For the sake of argument, let’s say both strtod() and compile time conversions honor the current rounding mode (I don’t know if compile time conversions honor the rounding mode, but I expect they don’t). Assume also that both give the same result for the same rounding mode for the same decimal input. The compiled and strtod() conversions can still differ — if the rounding mode is different at run time than at compile time.

GLIBC Bug Report

The glibc bug report “Bugzilla Bug 3479: Incorrect rounding in strtod()” mentions, as an aside, strtod()’s failure to respond to the current rounding mode (see comment 6).

Are There Any Applications That Depend on Directed Conversions?

Are there any Visual C++ or glibc based applications that depend on these “directed” conversions? I hope not — it doesn’t work.

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Visual C++ and GLIBC strtod() Ignore Rounding Mode

(To understand my methodology and terminology, please read my article “Incorrectly Rounded Conversions in Visual C++.”)

My Analysis

I ran my tests on an Intel Core Duo processor running Ubuntu Linux version 10.04, and I used the gcc C compiler version 4.4.3 and the embedded glibc library (eglibc) version 2.11.1. (As far as I can tell, eglibc uses the same conversion code as glibc, so I will refer to eglibc as glibc for the rest of this article.) (Update: I reran my tests on glibc, getting the same results.)

I did automated testing of glibc strtod() string conversions and manual testing of gcc floating-point literal conversions. This is what found in my automated testing of strtod():

• Conversions were done incorrectly only about four ten-thousandths of a percent of the time.
• No conversion was incorrect by more than one unit in the last place (ULP).
• Incorrect conversions included both “hard” and “easy” cases.
• There were incorrect conversions covering all three cases: “halfway”, “greater than halfway”, and “less than halfway.”

My manual testing of gcc was not as comprehensive, of course, but I did manage to find some examples of incorrect conversions.

For the remainder of this article, I will discuss six examples of incorrect conversions — examples uncovered by my automated and manual testing, as well as examples found in Bugzilla bug reports for gcc and glibc.

Incorrect Conversion by Both GCC and GLIBC

Here’s an example that’s converted incorrectly by both gcc and glibc — to the same incorrect value:

Example 1 (Halfway Case)

0.500000000000000166533453693773481063544750213623046875,

which equals 2^-1 + 2^-53 + 2^-54, converts to this binary number:

0.100000000000000000000000000000000000000000000000000011

It has 54 significant bits, with bits 53 and 54 equal to 1. The correctly rounded value — in binary scientific notation — is

1.000000000000000000000000000000000000000000000000001 x 2^-1

gcc and glibc strtod() both compute the converted value as

1.0000000000000000000000000000000000000000000000000001 x 2^-1

which is one ULP below the correctly rounded value.

(I constructed this example by hand.)

Incorrect Conversions by GLIBC Only

Here are three examples that are converted incorrectly by glibc but converted correctly by gcc:

Example 2 (Halfway Case)

3.518437208883201171875e13 converts to this binary number:

1000000000000000000000000000000000000000000000.00000011

It has 54 significant bits, with bits 53 and 54 equal to 1. The correctly rounded value is

1.000000000000000000000000000000000000000000000000001 x 2⁴⁵

glibc strtod() computes the converted value as

1.0000000000000000000000000000000000000000000000000001 x 2⁴⁵

which is one ULP below the correctly rounded value.

(This example was taken from “Bugzilla Bug 3479: Incorrect rounding in strtod().”)

Example 3 (Greater than Halfway Case)

62.5364939768271845828 converts to this binary number, which has a non-terminating binary fraction portion:

111110.1000100101010111101010110101010011111001010001110000000000001…

Bits 53 and 54 are equal to 1, so the correctly rounded value is

1.11110100010010101011110101011010101001111100101001 x 2⁵

glibc strtod() computes the converted value as

1.1111010001001010101111010101101010100111110010100011 x 2⁵

which is one ULP below the correctly rounded value.

(I discovered this example during random testing.)

Example 4 (Less than Halfway Case)

8.10109172351e-10 converts to this non-terminating binary fraction:

0.000000000000000000000000000000110111101010111001011101011101111000
01111110100001100011111111111110…

Significant bit 54 is 0, so it should be rounded down to get the correctly rounded value

1.10111101010111001011101011101111000011111101000011 x 2^-31

glibc strtod() computes the converted value as

1.1011110101011100101110101110111100001111110100001101 x 2^-31

which is one ULP above the correctly rounded value.

(I discovered this example during random testing.)

Incorrect Conversions by GCC Only

Here are two examples that are converted incorrectly by gcc but converted correctly by glibc:

Example 5 (Halfway Case)

The number

1.50000000000000011102230246251565404236316680908203125

which equals 1 + 2^-1 + 2^-53, converts to this binary number:

1.10000000000000000000000000000000000000000000000000001

It rounds correctly to 1.1, since bit 53 is 0.

gcc converts it to

1.1000000000000000000000000000000000000000000000000001

which is one ULP above the correctly rounded value.

(I constructed this example by hand.)

Example 6 (Less than Halfway Case)

9007199254740991.4999999999999999999999999999999995 converts to this binary number, which has a non-terminating binary fraction portion:

11111111111111111111111111111111111111111111111111111.01111111111111
1111111111111111111111111111111111111111111111111111111111111111111
111111111111111111111111111110…

Bit 54 is 0, so it converts correctly to

1.1111111111111111111111111111111111111111111111111111 x 2⁵²

gcc converts it to

1.0 x 2⁵³

which is one ULP above the correctly rounded value.

(This example was taken from “Bugzilla Bug 21718: real.c rounding not perfect.”)

A C Program that Reports the Behavior

Here’s a C program you can use to confirm incorrect rounding on your own Linux system:

// gcc -o examples examples.c
#include <stdio.h>
#include <stdlib.h>

#define NUM_CONVERT 6

int main (void)
{
 unsigned int i;

 char decimal[NUM_CONVERT][60] = {
  "0.500000000000000166533453693773481063544750213623046875",
  "3.518437208883201171875e13",
  "62.5364939768271845828",
  "8.10109172351e-10",
  "1.50000000000000011102230246251565404236316680908203125",
  "9007199254740991.4999999999999999999999999999999995"
 };

 double gcc_conversion[NUM_CONVERT] = {
  0.500000000000000166533453693773481063544750213623046875,
  3.518437208883201171875e13,
  62.5364939768271845828,
  8.10109172351e-10,
  1.50000000000000011102230246251565404236316680908203125,
  9007199254740991.4999999999999999999999999999999995
};

 char correct[NUM_CONVERT][30] = {
  "0x1.0000000000002p-1",
  "0x1.0000000000002p+45",
  "0x1.f44abd5aa7ca4p+5",
  "0x1.bd5cbaef0fd0cp-31",
  "0x1.8p+0",
  "0x1.fffffffffffffp+52"
 };

 for (i=0; i < NUM_CONVERT; i++)
 {
  printf ("%s\n",decimal[i]);
  printf (" Correct = %s\n",correct[i]);
  printf (" gcc =     %a\n",gcc_conversion[i]);
  printf (" strtod =  %a\n\n",strtod(decimal[i],NULL));
 }
}

This is the output on my system (I shaded the incorrect conversions):

0.500000000000000166533453693773481063544750213623046875
 Correct = 0x1.0000000000002p-1
 gcc =     0x1.0000000000001p-1
 strtod =  0x1.0000000000001p-1

3.518437208883201171875e13
 Correct = 0x1.0000000000002p+45
 gcc =     0x1.0000000000002p+45
 strtod =  0x1.0000000000001p+45

62.5364939768271845828
 Correct = 0x1.f44abd5aa7ca4p+5
 gcc =     0x1.f44abd5aa7ca4p+5
 strtod =  0x1.f44abd5aa7ca3p+5

8.10109172351e-10
 Correct = 0x1.bd5cbaef0fd0cp-31
 gcc =     0x1.bd5cbaef0fd0cp-31
 strtod =  0x1.bd5cbaef0fd0dp-31

1.50000000000000011102230246251565404236316680908203125
 Correct = 0x1.8p+0
 gcc =     0x1.8000000000001p+0
 strtod =  0x1.8p+0

9007199254740991.4999999999999999999999999999999995
 Correct = 0x1.fffffffffffffp+52
 gcc =     0x1p+53
 strtod =  0x1.fffffffffffffp+52

The values for the correct conversions were obtained using David Gay’s strtod() function. I ran a separate program to compute them and then hard coded them in this program. (I could also have constructed the correct hex formats by hand, using my analysis in the examples above.)

If you try this on your system, let me know what happens; I expect the same output. Please include the type of Linux, and version numbers if possible.

Bug Reports

There are existing bug reports for gcc and glibc:

• Bugzilla Bug 21718: real.c rounding not perfect.
• Bugzilla Bug 3479: Incorrect rounding in strtod().

It’s not clear whether these bugs will (or should?) be fixed; if they are, I’d imagine David Gay’s strtod() function would be used to replace the current code.

I would guess this problem affects the other GCC compilers, although I haven’t tested.

GCC/GLIBC Vs. Visual C++

Here are the differences I found between gcc/glibc and Visual C++ in my two analyses:

• glibc strtod() converted incorrectly about two orders of magnitude less frequently than Visual C++ strtod().
• I found no integer values that rounded incorrectly in gcc or glibc, unlike in Visual C++.
• gcc and glibc strtod() can convert the same decimal number differently; I found no such example in Visual C++ and its strtod() function.
• Of the five examples Visual C++ converts incorrectly, only one — example 2 — is converted incorrectly by gcc or glibc (in this case, by both). Of the six examples in this article converted incorrectly by gcc or glibc, two are converted incorrectly by Visual C++: examples 1 and 2.
• When gcc or glibc rounds halfway cases incorrectly, it does so either up or down. I found Visual C++ to consistently round incorrectly down.
• I found an incorrect conversion for a “less than halfway case” in gcc/glibc; I found no such example in Visual C++.

Summary

I’ve demonstrated that a decimal floating-point number is not guaranteed to convert to its closest binary floating-point number. Also, I’ve shown that a decimal floating-point number may be converted to different binary floating-point numbers, depending on which compiler or run-time library is used.

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Incorrectly Rounded Conversions in GCC and GLIBC

My Analysis

I did my analysis using Visual C++ Express Edition — both the 2008 and 2010 versions. I converted decimal numbers to double-precision floating-point, via the strtod() library function and floating-point literals converted implicitly by the compiler. I discovered examples of incorrect conversion in three ways: by random testing (using only strtod() so that it could be automated easily), by trial and error, and by direct construction.

My random testing was extensive, but not exhaustive. I generated random decimal strings (positive numbers only) in scientific notation, with a random number of digits (up to 21) and a random exponent (between -308 and +308). I compared the Visual C++ conversion to that given by David Gay’s strtod() function, the de facto standard correct conversion routine. This is what I found:

• Conversions were done incorrectly about three hundredths of a percent of the time.
• No conversion was incorrect by more than one unit in the last place (ULP).
• Incorrect conversions included both “hard” and “easy” cases (defined below).
• All incorrect conversions that I verified by hand were either “halfway” or “greater than halfway” cases (defined below).

For all the examples below, I used my C function fp2bin() to inspect the converted floating-point values. I verified that both strtod() and the compiler converted the same decimal number to the same double-precision floating-point value. I also verified the correctly converted value by hand, using the full-precision binary numbers generated by my decimal to binary converter.

Halfway Cases

A halfway case is when a decimal number converts to a full-precision binary number that is halfway between two binary floating-point values. By convention, in correct rounding, the tie is broken by choosing the floating-point value with its least significant bit equal to 0.

In terms of double-precision floating-point, the halfway case happens when significant bit 54 of the full-precision binary number is equal to 1, and all bits beyond are equal to 0. If bit 53 is equal to 0, then the full-precision binary number is rounded down to 53 bits; if bit 53 is equal to 1, then the full-precision binary number is rounded up to 53 bits. In either case, the correctly rounded result has 0 in its 53rd significant bit.

My testing shows that Visual C++ consistently rounds halfway cases incorrectly; specifically, it always rounds down, even when bit 53 of the full-precision binary number equals 1.

Example 1 (Incorrect)

The integer 9,214,843,084,008,499 converts to this binary number:

100000101111001101100111011000101011101110110000110011

Bit 54 is 1, and with no other bits beyond, is a halfway case. Bit 53 is 1, so the correctly rounded value — in binary scientific notation — is

1.000001011110011011001110110001010111011101100001101 x 2⁵³

Visual C++ computes the converted value as

1.0000010111100110110011101100010101110111011000011001 x 2⁵³

which is one ULP below the correctly rounded value.

Example 2 (Incorrect)

The dyadic fraction

0.500000000000000166533453693773481063544750213623046875,

which equals 2^-1 + 2^-53 + 2^-54, converts to this binary number:

0.100000000000000000000000000000000000000000000000000011

It has 54 significant bits, with bits 53 and 54 equal to 1. The correctly rounded value is

1.000000000000000000000000000000000000000000000000001 x 2^-1

Visual C++ computes the converted value as

1.0000000000000000000000000000000000000000000000000001 x 2^-1

which is one ULP below the correctly rounded value.

Discussion

Examples 1 and 2 would be considered “easy” cases in correct rounding parlance. They require only one extra bit beyond double precision to make the correct rounding decision. These examples support my assertion that Visual C++ universally adopts truncation for halfway cases.

Greater than Halfway Cases

The greater than halfway case occurs when a decimal number converts to a full-precision binary number that is between two binary floating-point values, but closer to the one above. Visual C++ gets some of these wrong and some of these right.

Example 3 (Incorrect)

The integer 30,078,505,129,381,147,446,200 converts to this binary number:

11001011110100011110001110010011100011011000000100000100000000010111
0111000

Bit 54 is 1 and there is a 1 bit beyond bit 54, so the correctly rounded value is

1.1001011110100011110001110010011100011011000000100001 x 2⁷⁴

Visual C++ computes the converted value as

1.10010111101000111100011100100111000110110000001 x 2⁷⁴

which is one ULP below the correctly rounded value.

Example 4 (Correct)

The integer 1,777,820,000,000,000,000,001 converts to this binary number:

11000000110000000110101011011101011010010111000111101100000000000000
001

Just like the prior example, bit 54 is 1 and there is a 1 bit beyond it — but Visual C++ converts this correctly, to

1.100000011000000011010101101110101101001011100011111 x 2⁷⁰

Example 5 (Incorrect)

The dyadic fraction

0.500000000000000166547006220929549868969843373633921146392822265625

which equals 2^-1 + 2^-53 + 2^-54 + 2^-66, converts to this binary number:

0.100000000000000000000000000000000000000000000000000011000000000001

It converts correctly to

1.000000000000000000000000000000000000000000000000001 x 2^-1

But Visual C++ converts it to this, one ULP below the correctly rounded value:

1.0000000000000000000000000000000000000000000000000001 x 2^-1

Example 6 (Correct)

The dyadic fraction

0.50000000000000016656055874808561867439493653364479541778564453125

which equals 2^-1 + 2^-53 + 2^-54 + 2^-65, converts to this binary number:

0.10000000000000000000000000000000000000000000000000001100000000001

Visual C++ converts this correctly, to

1.000000000000000000000000000000000000000000000000001 x 2^-1

Example 7 (Incorrect)

The decimal fraction 0.3932922657273 is non-terminating in binary (it’s not dyadic); here is the relevant portion of it — its first 82 places:

0.011001001010111011001101010010110001000110001000111001100000000000
0000000000000001…

Significant bits 53 and 54 are both 1, as is significant bit 81; it rounds correctly to

1.100100101011101100110101001011000100011000100011101 x 2^-2

Visual C++ computes the value one ULP below:

1.1001001010111011001101010010110001000110001000111001 x 2^-2

Discussion

Visual C++ gets example 4 correct, even though it appears harder than example 3, which it converts incorrectly; example 4 has six more zeros beyond bit 54.

Examples 5 and 6 are dyadic fractions, which theoretically should be easy to convert. Example 5 is 66 digits long, and example 6 is 65 digits long. (Dyadic fractions have the same number of bits as decimal digits.) Example 6 converts correctly, but example 5 does not. Is it because Visual C++ stops looking beyond 65 places?

Example 7 is a hard case, which Visual C++ does get wrong.

My Bug Report

In April 2009, I submitted a bug report to Microsoft, titled “Incorrectly Rounded Decimal to Binary Conversions in Visual C++.” Microsoft closed the bug, but said they would “re-consider [sic] for future versions of VC++”. We’ll see.

You could argue that being only one ULP away from the correct result is inconsequential. Technically, I think it’s ‘legal’ from the perspective of the IEEE floating-point specification. (This seems to be conventional wisdom, but I can’t find the exact wording in the spec — can someone enlighten me?) However, widely used open source code — David Gay’s strtod() function in dtoa.c — gets these conversions right. Should Visual C++ do the same?

(I would imagine this ‘problem’ also affects the rest of Visual Studio, but I haven’t investigated.)

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Incorrectly Rounded Conversions in Visual C++

The quick and dirty program is implemented in native C, so it’s limited to double-precision floating-point arithmetic (although on some systems, extended precision may be used). Higher precision arithmetic — in fact, arbitrary precision arithmetic — is needed to ensure that all decimal inputs are converted correctly. I will demonstrate the need for high precision by analyzing three examples, all taken from Vern Paxson’s paper “A Program for Testing IEEE Decimal–Binary Conversion”.

Correct Rounding

A correctly rounded conversion from a decimal string to a double-precision floating-point number is the 53 significant bit value closest to the decimal input — and in the event of a tie, the value with significant bit number 53 equal to 0. (This is the prevailing definition, being based on the default IEEE 754 rounding mode. But to be truly correctly rounded — theoretically at least — all four IEEE 754 rounding modes must be honored.) The closest 53-bit number is obtained by rounding what would be the full-precision binary number to 53 bits. This can be done correctly only if the correct value of bit 54 — and sometimes the correct values of an arbitrary number of bits beyond — is known. Here are the rules:

• If bit 54 is 0, round down
• If bit 54 is 1
  • If any bit beyond bit 54 is 1, round up
  • If all bits beyond bit 54 are 0 (meaning the number is halfway between two floating-point numbers)
    • If bit 53 is 0, round down
    • If bit 53 is 1, round up

The rules seem simple enough, but there’s a problem: the calculation needs to be precise enough so that the true value of bit 54 and the appropriate number of bits beyond is known. Decimal numbers that are close to halfway between two floating-point numbers require the most precision to get right; the three examples that follow are such numbers.

Example 1: 7.8459735791271921 x 10⁶⁵

7.8459735791271921 x 10⁶⁵ is shorthand for this number:

784597357912719210000000000000000000000000000000000000000000000000

In binary, it equals

11101110011010000000010001001110000010011000101001110011111111111111
11111111111111111111111111111111111111111111111111000111100101001110
11111110111110111000110111111101010000000000000000000000000000000000
000000000000000

Rounded correctly to 53 bits it equals

11101110011010000000010001001110000010011000101001110000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000000000000000000000000000000000
000000000000000

which in binary scientific notation equals

1.110111001101000000001000100111000001001100010100111 x 2²¹⁸

Why This Conversion is Difficult in Floating-Point

The full-precision binary value above was obtained using my decimal to binary converter. The converter uses arbitrary precision decimal arithmetic, so every bit in the result is correct. Since bit 54 (the first highlighted bit above) is 0, it is rounded down to get the correctly rounded 53-bit result.

When the conversion is done in floating-point, things get a bit murky. With insufficient precision, it’s impossible to tell whether bit 54 is really 0, since it’s a tiny fraction away from being 1. The conversion algorithm cannot make the distinction unless the calculation is known to be precise through the 119th bit (the second highlighted bit above).

Converting with Varying Floating-Point Precision

To demonstrate the effect of floating-point precision, I modified my quick and dirty conversion program to use GMP floating-point arithmetic. I increased the precision of the GMP conversions incrementally until I reached the intermediate result with sufficient precision; for each increment, I printed the result in binary and rounded correctly by hand.

I converted 7.8459735791271921 x 10⁶⁵ to double precision using the quick and dirty program (53-bit precision) and GMP with 64, 96, and 128 bits of precision (on my system, GMP changes precision internally at 32-bit granularity, plus a few extra bits); here are the results, expressed in binary scientific notation:

• Quick and Dirty (53-bit precision):

  1.110111001101000000001000100111000001001100010101 x 2²¹⁸

• GMP, with 64 and 96 bit precision:

  1.110111001101000000001000100111000001001100010100111 x 2²¹⁸

• GMP, with 128-bit precision:

  1.110111001101000000001000100111000001001100010100111 x 2²¹⁸

The quick and dirty result is incorrect, and is off by two units in the last place (ULPs). All three GMP values are correct, but the 64 and 96 bit values are correct by chance; before rounding, the intermediate result (not shown) has bit 53 equal to 0 and bit 54 equal to 1, with all bits beyond equal to 0. That makes it look like a halfway case, so it rounds down to the correct value — by sheer luck. The 128-bit GMP calculation is precise enough that the 119th bit is known to be 0, meaning bit 54 is 0 and that the result should be rounded down.

Example 2: 3.571 x 10²⁶⁶

3.571 x 10²⁶⁶ is an 886 bit binary number, so I’ll only show the relevant portion — its first 77 bits:

10110001001100100010011000110000111010100000110101001100000000000000
000000001…

Rounded correctly to 53 bits it equals

1.011000100110010001001100011000011101010000011010101 x 2⁸⁸⁵

Here are the quick and dirty and GMP conversions:

• Quick and Dirty (53-bit precision):

  1.011000100110010001001100011000011101010000011011001 x 2⁸⁸⁵

• GMP, with 64-bit precision:

  1.0110001001100100010011000110000111010100000110101001 x 2⁸⁸⁵

• GMP, with 96-bit precision:

  1.011000100110010001001100011000011101010000011010101 x 2⁸⁸⁵

The quick and dirty and 64-bit GMP results are both incorrect; the quick and dirty result is off by eight ULPs, and the 64-bit GMP result is off by one ULP. The 64-bit conversion goes wrong because GMP computes this intermediate value (it’s 66 bits, since GMP can carry extra bits beyond the expected precision):

101100010011001000100110001100001110101000001101010010111111111110

At 64-bit precision, bit 54 appears to be 0, and hence the result is rounded down.

As expected, the 96-bit GMP result is correct, since only 77 bits of precision are needed to decide the fate of bit 54.

Example 3: 3.08984926168550152811 x 10^-32

3.08984926168550152811 x 10^-32 is a decimal fraction which, when written in binary, is non-terminating; for our purposes, we only need to see its first 234 places (104 leading 0s, and significant bits 1-53 and 54-130):

0.000000000000000000000000000000000000000000000000000000000000000000
00000000000000000000000000000000000000101000000110111100100100001100
11101100110010100111010100000000000000000000000000000000000000000000
00000000000000000000000000000001…

Rounded correctly to 53 bits it equals

1.0100000011011110010010000110011101100110010100111011 x 2^-105

Here are the quick and dirty and GMP conversions:

• Quick and Dirty (53-bit precision):

  1.010000001101111001001000011001110110011001010011101 x 2^-105

• GMP, with 64 and 96 bit precision:

  1.010000001101111001001000011001110110011001010011101 x 2^-105

• GMP, with 128-bit precision:

  1.0100000011011110010010000110011101100110010100111011 x 2^-105

• GMP, with 160-bit precision:

  1.0100000011011110010010000110011101100110010100111011 x 2^-105

The quick and dirty, 64-bit GMP, and 96-bit GMP results are all incorrect; they are all the same, and are off by 1 ULP. (The Visual C++ compiler gets this conversion wrong too, getting the same wrong answer!)

The 128-bit GMP result is correct by accident. The intermediate value used to compute it (not shown) has a 0 at bit 53, a 1 at bit 54, and then a string of 74 0s followed by a 1 (yes, 129 bits). The full-precision result has 75 0s followed by a 1. The difference does not matter though, since both bit patterns indicate that the result should be rounded up.

The 160-bit GMP result is correct, which is expected since 130 bits of the full-precision binary result are needed to guarantee the right rounding decision is made.

How Much Precision is Enough?

In the examples above, I knew how much precision was needed for correct rounding because I already knew the full-precision binary equivalents of the decimal inputs — by external means. But how would a floating-point conversion routine know how much precision is enough? Naively, it would have to assume the worst case, computing with enough precision to match the input (for integer values, this would mean computing the full-precision binary equivalent). That would be overkill, since most of the time much less precision is needed to make a correct rounding decision.

David Gay’s strtod() function in dtoa.c — the de facto standard conversion routine — avoids this problem by avoiding high precision floating-point. It uses standard double-precision floating-point for most conversions, and arbitrary precision integer arithmetic for “tough” conversions; it has the smarts to know when to use one or the other. (I will discuss that code in a future article.)

Acknowledgement

Thanks to Vern Paxson for helping me understand the “very close to but less than 1/2 ULP case.”

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Decimal to Floating-Point Needs Arbitrary Precision

#include <string.h>

int main (void)
{
 double intPart = 0, fracPart = 0, conversion;
 unsigned int i;
 char decimal[] = "3.14159";

 i = 0; /* Left to right */
 while (decimal[i] != '.') {
    intPart = intPart*10 + (decimal[i] - '0');
    i++;
   }

 i = strlen(decimal)-1; /* Right to left */
 while (decimal[i] != '.') {
    fracPart = (fracPart + (decimal[i] - '0'))/10;
    i--;
   }

 conversion = intPart + fracPart;
}

The conversion is done using the elegant Horner’s method, summing each digit according to its decimal place value. So why do I call it “quick and dirty?” Because the binary floating-point value it produces is not necessarily the closest approximation to the input decimal value — the so-called correctly rounded result. (Remember that most real numbers cannot be represented exactly in floating-point.) Most of the time it will produce the correctly rounded result, but sometimes it won’t — the result will be off in its least significant bit(s). There’s just not enough precision in floating-point to guarantee the result is correct every time.

I will demonstrate this algorithm with different input values, some of which convert correctly, and some of which don’t. In the end, you’ll appreciate one reason why library functions like strtod() exist — to perform efficient, correctly rounded conversion.

Testing the Conversions

How do you test whether a quick and dirty conversion is correct? First you need to know what the correctly rounded value is, and then you need to compare it to the quick and dirty value.

Here are two ways to find the correctly rounded value:

• Assign the decimal constant to a double and let your compiler compute it for you.
• Enter the decimal constant in my arbitrary-precision decimal/binary converter and then round it by hand to 53 significant bits (the number of significant bits in a double).

I did both, because I know that sometimes even compilers do the conversions incorrectly (see my articles on Visual C++ and gcc/glibc). In any case, I verified that all the compiler conversions in the examples below are correct.

Next we need to compare the quick and dirty value to the correctly rounded value — but how? Comparing floating-point values for equality is not reliable. That means we should display the values and compare them manually. But we can’t depend on printf(), because in many implementations it won’t print all the significant digits, either before the decimal point or after the decimal point (this is true of Visual C++, the compiler I am using). This leaves two options:

• Print the raw IEEE representation of the double.
• Print the value of the double in binary, using my function fp2bin().

I chose to use fp2bin(), because seeing the binary representation is more enlightening.

The Test Code

I modified the code in the introduction of this article to allow for testing, as described; here it is

#include <string.h>
#include <stdio.h>
#include "fp2bin.h"

int main (void)
{
 char binString[FP2BIN_STRING_MAX];
 double intPart = 0, fracPart = 0, conversion;
 unsigned int i;
 char decimal[] = "3.14159"; /* Converted manually */
 double compiler_conversion = 3.14159; /* Converted by compiler */

 i = 0; /* Left to right */
 while (decimal[i] != '.')
   {
    intPart = intPart*10 + (decimal[i] - '0');
    i++;
   }

 i = strlen(decimal)-1; /* Right to left */
 while (decimal[i] != '.')
   {
    fracPart = (fracPart + (decimal[i] - '0'))/10;
    i--;
   }

 conversion = intPart + fracPart;

 /* Print the converted values, in binary */
 fp2bin(compiler_conversion,binString);
 printf("Correct = %s\n",binString);
 fp2bin(conversion,binString);
 printf("Q and D = %s\n",binString);
}

When run (I ran on an Intel Core Duo processor), the output shows that the compiler and quick and dirty conversions of 3.14159 are both correct:

Correct = 11.00100100001111110011111000000011011100001100110111
Q and D = 11.00100100001111110011111000000011011100001100110111

Their value in decimal, obtained from my binary/decimal converter, is
3.14158999999999988261834005243144929409027099609375.

Incorrect Conversions

Incorrect conversions are due to the limitations of floating-point arithmetic; specifically, these three things affect the conversion calculations:

• Large integers require more bits than can fit in a floating-point variable.

  Multiplying a binary integer by ten is exact, but once the integer exceeds the largest integer that a double can represent exactly, rounding errors can occur.

• Division by ten in binary arithmetic is inexact.

  Each time the fractional value is divided by ten, there is a chance of roundoff error.

• Addition of integers and fractions can incur rounding error.

  When an integer and fraction are added, some of the bits of the fraction can be lost.

The example values below — which I discovered by both random testing and trial and error — are all rounded incorrectly by my quick and dirty algorithm.

The Integer 18014398509481993

The number 18,014,398,509,481,993 (2⁵⁴ + 9) is not converted correctly by the quick and dirty algorithm; this is the output of the program (the lines I changed are char decimal[] = “18014398509481993.0″; and double compiler_conversion = 18014398509481993.0;):

Correct = 1000000000000000000000000000000000000000000000000001000.0
Q and D = 1000000000000000000000000000000000000000000000000001100.0

Their values, in decimal, are

• Correct = 18014398509481992 = 2⁵⁴ + 8
• Q and D = 18014398509481996 = 2⁵⁴ + 12

These values differ by 4, which is one unit in the last place (the last place for a double is the 53rd significant bit).

The Fraction 0.9199

The number 0.9199 is not converted correctly by the quick and dirty algorithm; this is the output of the program (the lines I changed are char decimal[] = “0.9199″; and double compiler_conversion = 0.9199;):

Correct = 0.11101011011111101001000011111111100101110010010001111
Q and D = 0.1110101101111110100100001111111110010111001001000111 

Their values, in decimal, are

• Correct = 0.91990000000000005098144129078718833625316619873046875

  = 8285722594436239/2⁵³

• Q and D = 0.9198999999999999399591388282715342938899993896484375

  = 4142861297218119/2⁵² = 8285722594436238/2⁵³

These values differ by 2^-53, which is one unit in the last place.

The Combined Integer and Fraction 1.89

The number 1.89 is not converted correctly by the quick and dirty algorithm; this is the output of the program (the lines I changed are char decimal[] = “1.89″; and double compiler_conversion = 1.89;):

Correct = 1.1110001111010111000010100011110101110000101000111101
Q and D = 1.111000111101011100001010001111010111000010100011111 

Their values, in decimal, are

• Correct = 1.8899999999999999023003738329862244427204132080078125

  = 1 + 4008203668359741/2⁵²

• Q and D = 1.890000000000000124344978758017532527446746826171875

  = 1 + 2004101834179871/2⁵¹ = 1 + 4008203668359742/2⁵²

These values differ by 2^-52, which is two units in the last place.

The interesting thing is that 0.89 by itself converts correctly; it’s only when it’s added to 1 that the final result becomes incorrect.

Discussion

In practice, you’ll depend on your compiler or interpreter or a library call for correctly rounded conversion. The de facto standard code used by many implementations is David Gay’s strtod() function in dtoa.c. Take a look at it to get a sense of how involved correctly rounded conversion is. (dtoa.c handles scientific notation, and also includes code for correct conversion in the other direction — floating-point to decimal).

For a little more detailed discussion of Horner’s method, see my PHP functions bin2dec_i() and bin2dec_f(). Those two functions implement the same quick and dirty algorithm, only using decimal arithmetic and powers of two (there’s nothing “dirty” about them though — they’re binary to decimal and arbitrary precision).

Conversion to Single-Precision

Conversion to floats works similarly, also producing correctly and incorrectly rounded results. There’s one thing I observed, which you might find surprising: a given decimal value may convert incorrectly to a double but correctly to a float!

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Quick and Dirty Decimal to Floating-Point Conversion

I started my search with my program that finds multiple-base palindromes. I generated palindromes in binary, and then checked them to see if they were also palindromes in hexadecimal and decimal. I looked for decimal/binary/hexadecimal palindromes up to 16 hex digits long, but did not find any.

To continue my search into bigger numbers, I wrote a program that uses arbitrary-precision integer arithmetic and a more efficient algorithm. Despite being able to search much further, I still have not found any.

In this article, I’ll analyze the size of the palindrome “search space”, explain my improved algorithm, and discuss the state of my search.

Search Space for Decimal/Binary/Hexadecimal Palindromes

The biggest obstacle in the search for palindromes is the growth in the amount of numbers that must be checked — the search space. As the number of digits increases, the search space grows exponentially. While this exponential growth can’t be eliminated, it can be delayed — by pruning the search space.

Consider the search space for 32 hexadecimal digit decimal/binary/hexadecimal palindromes. With no pruning, 15·16³¹ integers would have to be checked — 319,014,718,988,379,809,496,913,694,467,282,698,240 numbers, an impossible task. If only hexadecimal palindromes are considered, the search space is reduced to 15·16¹⁵ = 17,293,822,569,102,704,640 numbers — much smaller but still intractable. If only binary/hexadecimal palindromes are considered, then there are 2¹⁶ + 2³² = 4,295,032,832 numbers to be checked — quite manageable!

The binary/hexadecimal palindrome space can be partitioned by starting digit so that subsets can be searched independently. Consider the 32 hexadecimal digit binary/hexadecimal palindromes. There are 2¹⁶ = 65,536 that begin with the digits 1 and 3, and 2³² = 4,294,967,296 that begin with the digits 5, 7, 9, and F. The number of palindromes that start with 1 and 3 is markedly lower, allowing for a deeper search within them.

The following diagrams illustrate how the search space grows, showing the search spaces for the 16, 32, and 64 digit hexadecimal palindromes:

Number of 16-Digit Hexadecimal Palindromes, by Subset (not to scale)

Number of 32-Digit Hexadecimal Palindromes, by Subset (not to scale)

Number of 64-Digit Hexadecimal Palindromes, by Subset (not to scale)

Program to Find Decimal/Binary/Hexadecimal Palindromes

I wrote a C program to search for decimal/binary/hexadecimal palindromes — it is based on two things:

• Arbitrary-precision integers.

  This allows for palindrome numbers that are too big to fit in an unsigned long long integer variable; that is, numbers bigger than 64 bits, or 16 hex digits. I used GMP in Visual C++ on Windows.

• Direct generation of binary/hexadecimal palindromes.

  This reduces the search space dramatically, so that only known binary/hexadecimal palindromes are checked — for being palindromic in decimal. I used a simple algorithm derived from my analysis of binary/hexadecimal palindromes.

I will outline the algorithm instead of providing the source code.

The Algorithm

Binary/hexadecimal palindromes could be generated in two ways. One way would be incrementally, by generating the two-digit palindromes, then the three-digit palindromes, then the four-digit palindromes, etc. A set of m-digit palindromes would be built from a set of m-1 digit palindromes, by inserting a digit in the middle of each m-1 digit palindrome. This requires storage of palindromes, which will not scale — the number of palindromes grows exponentially.

Another way to generate binary/hexadecimal palindromes — the way I chose — is to construct them from permutations of hex digits. Each of the six subsets of binary/hexadecimal palindromes has its own rules for valid digits:

• Starts/ends with 1: each digit between can be 0 or 1
• Starts/ends with 3: each digit between can be 0 or 3
• Starts/ends with 5: each digit between can be 0, 2, 5, or 7
• Starts/ends with 7: each digit between can be 0, 2, 5, or 7
• Starts/ends with 9: each digit between can be 0, 6, 9, or F
• Starts/ends with F: each digit between can be 0, 6, 9, or F

You construct palindromes of n digits from permutations of floor(n/2) digits. For example, to generate a 32 digit palindrome that starts with 5, create a number that starts with 5 followed by a 15 digit permutation of the digits 0, 2, 5, and 7. Then, reverse that number and append it to itself to create a 32 digit palindrome. (The process for odd length palindromes is similar, except you insert an extra middle digit between the two halves.)

Permutations are easily generated, with a series of nested loops. You can nest loops statically, by coding a loop for each place of the permutation, or you can nest them dynamically, by calling a single loop recursively. The advantage of dynamic nesting is not having to add loops manually for increasing digit lengths; the advantage of static nesting is the ability to “checkpoint” the state of the search easily — an important feature for long running programs. I chose static nesting for that reason.

Results

I’ve searched the “1/3 space” through all 64 hex digit numbers, and the “5/7/9/F space” through all 32 hex digit numbers. I have not found a single multi-digit decimal/binary/hexadecimal palindrome.

I have suspended my search, since the running time of my program is starting to explode. For example, it took almost 60 hours on my Intel Core Duo processor to test all 2³² 64 hex digit binary/hexadecimal palindromes in the 1/3 space. Resuming the search would probably require a GIMPS-like distributive computing approach, which I have no plans to undertake.

Discussion

Are there any tricks to reduce the search space further, delaying the inevitable exponential explosion even more? Charlton’s binary/decimal palindrome search algorithm eliminates swaths of numbers as it runs, recognizing conditions under which copalindromic numbers can’t exist. Can similar ideas be applied to the search for decimal/binary/hexadecimal palindromes?

Another Way to Count Binary/Hexadecimal Palindromes

In my article “Counting Binary/Hexadecimal Palindromes” I derived formulas to count binary/hexadecimal palindromes using geometric series. An alternate way would be based on counting permutations. For example, the number of 32 digit palindromes starting with 5 is the number of 15 digit permutations of the digits 0, 2, 5, and 7; that is, 4¹⁵ = 2³⁰ = 1,073,741,824.

(You would still need geometric series to count the number of palindromes of n-digits or less.)

Are There Any?

Are there any decimal/binary/hexadecimal palindromes? Can it be proved either way?

I have no reason to believe they don’t exist. After all, there are at least two multi-digit decimal/binary/octal palindromes. And it’s always possible I made an error in my search, although Charlton’s results confirm my results through 32 hex digits plus.

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

In Search of Decimal/Binary/Hexadecimal Palindromes

Binary Code on Barbie's Laptop

It spells “Barbie” — repeatedly, in ASCII code.

B-a-r-b-i-e

The screen shows six full lines of code — three copies of these two lines:

01000010 01100001 01110010
01100010 01101001 01100101

Those happen to be six ASCII character codes, written in binary:

01000010 = B
01100001 = a
01110010 = r
01100010 = b
01101001 = i
01100101 = e

I’ve overlaid the characters on the screen so you can see the correspondence better:

Barbie's Laptop Overlaid with Corresponding ASCII Characters

(Now if instead those were numbers encoded in IEEE floating-point I’d have been really impressed .)

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

Barbie Goes Binary

In this article, I’ll do a similar analysis for double-precision floating-point variables, showing how similar anomalies arise when extended precision calculations are done. I modified my two example programs to use doubles instead of floats. Interestingly, the doubles version of program 2 does not exhibit the anomaly. I’ll explain.

The Environment in Which the Anomalies Occur

To see the described anomalies, floating-point calculations must be done in precision greater than double-precision. For Intel processors, this means that x87 FPU instructions must be used (as opposed to SSE instructions), and the processor must be calculating in extended-precision mode.

Your compiler determines which instructions to generate, although you may influence it with compiler options. On my system, Visual C++ generates x87 FPU instructions by default. As for the precision mode, it has a default — which is double-precision on my system — that can be changed by a call to the controlfp function. I call this function in the programs below.

Program 1

#include "stdio.h"
#include "float.h"
int main (void)
{
 double d1 = 0.3333333333333333, d2 = 3.0, d3;
 unsigned int cur_cw;

 _controlfp_s(&cur_cw,_PC_64,MCW_PC); 

 d3 = d1 * d2;
 if (d3 != d1 * d2)
   printf("Not equal\n");
}

Program 1 prints “Not equal”, for the same reason its float counterpart does: greater precision is being used under the covers.

Computing d3

Using my function fp2bin, I printed the binary values of the three floating point variables in the computation:

• d1 = 0.010101010101010101010101010101010101010101010101010101
• d2 = 11
• d3 = 1

The value of d3 is only an approximation to d1 * d2. Here’s the true value of d1 * d2, computed by hand:

  0.010101010101010101010101010101010101010101010101010101
  x                                                     11
  --------------------------------------------------------
     10101010101010101010101010101010101010101010101010101
    10101010101010101010101010101010101010101010101010101
  --------------------------------------------------------
  0.111111111111111111111111111111111111111111111111111111

The product is 54 bits, which can be represented exactly in extended precision, but only approximately in double-precision. Since bit 54 is a 1, the value is rounded up, to 53 bits. The result is d3 = 1.

Here’s the relevant assembler code generated by the compiler:

 d3 = d1 * d2;
0041131F  fld         qword ptr [d1]
00411322  fmul        qword ptr [d2]
00411325  fstp        qword ptr [d3]

Comparing d3 and d1 * d2

The assembler code shows what is happening:

 if (d3 != d1 * d2)
00411328  fld         qword ptr [d1]
0041132B  fmul        qword ptr [d2]
0041132E  fcomp       qword ptr [d3]

d1 * d2 is computed in extended precision and left on the stack. It is then compared to d3, which still has 53 bits of precision even though it’s “promoted” to extended precision. The two compared values will differ.

Program 2

#include "stdio.h"
#include "float.h"
int main (void)
{
 double d1 = 0.3333333333333333, d2 = 3.0, d3;
 long long i1, i2;
 unsigned int cur_cw;

 _controlfp_s(&cur_cw,_PC_64,MCW_PC); 

 d3 = d1 * d2;
 i1 = (int)d3;
 i2 = (int)(d1 * d2);
 if (i1 != i2)
   printf("Not equal\n");
}

Program 2 does not print anything, which is surprising knowing that its float counterpart does. According to our analysis of program 1 — program 2 uses the same values for d1 and d2 — d1 * d2 rounds to 1 in double precision, but is slightly less than 1 in extended precision. Therefore, we would expect i1 = 1 and i2 = 0. Instead, both are 1.

The explanation is simple: the extended precision calculation of d1 * d2 is “spilled” into a double before the conversion to integer, so it will be equal to d3:

Computing i2

 i2 = (int)(d1 * d2);
004114A5  fld         qword ptr [d1]
004114A8  fmul        qword ptr [d2]
004114AB  call        @ILT+245(__ftol2_sse) (4110FAh)
           ||
           ||
           \/
--- f:\dd\vctools\crt_bld\SELF_X86\crt\prebuild\tran\i386\ftol2.asm
00411B40  cmp         dword ptr [___sse2_available (417598h)],0
00411B47  je          _ftol2 (411B76h)
00411B49  push        ebp
00411B4A  mov         ebp,esp
00411B4C  sub         esp,8
00411B4F  and         esp,0FFFFFFF8h
00411B52  fstp        qword ptr [esp]
00411B55  cvttsd2si   eax,mmword ptr [esp]
00411B5A  leave
00411B5B  ret

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

When Doubles Don’t Behave Like Doubles

Program 1

#include "stdio.h"
int main (void)
{
 float f1 = 0.1f, f2 = 3.0f, f3;

 f3 = f1 * f2;
 if (f3 != f1 * f2)
   printf("Not equal\n");
}

Prints “Not equal”.

Program 2

#include "stdio.h"
int main (void)
{
 float f1 = 0.7f, f2 = 10.0f, f3;
 int i1, i2;

 f3 = f1 * f2;
 i1 = (int)f3;
 i2 = (int)(f1 * f2);
 if (i1 != i2)
   printf("Not equal\n");
}

Prints “Not equal”.

In each case, f3 and f1 * f2 differ. But why? I’ll explain what’s going on.

(This article was inspired by this question and several related questions on stackoverflow.com.)

Analyzing Program 1

Program 1 compares two floating point values for equality — something conventional wisdom says not to do. But why should that rule apply in this case? We’re just checking that a variable holds the value we just gave it. How could it not?

The root of the “problem” is this: the compiler generates instructions that do floating-point calculations in extended precision. Whereas floats are 4 bytes long and have 24 significant bits of precision, extended precision values — which are stored in registers on the floating point stack — are 10 bytes long and have 64 bits of precision.

Computing f3

Using my function fp2bin, I printed the binary values of the three floating point variables in the computation:

• f1 = 0.000110011001100110011001101
• f2 = 11
• f3 = 0.010011001100110011001101

The value of f3 is only an approximation to f1 * f2. To understand why, let’s calculate the true value of f1 * f2; that is, let’s calculate it by hand, using binary multiplication:

  0.000110011001100110011001101
  x                          11
  -----------------------------
       110011001100110011001101
      110011001100110011001101
-------------------------------
  0.010011001100110011001100111

So f1 * f2 = 0.010011001100110011001100111. It fits comfortably within extended precision. But to assign it to f3 — a float — it must be rounded. f1 * f2 has 26 significant bits, but a float holds only 24. Rounding it to the nearest 24 bit value makes it 0.010011001100110011001101.

The assembler code generated by the compiler confirms what we’re seeing:

 f3 = f1 * f2;
0041130B  fld         dword ptr [f1]
0041130E  fmul        dword ptr [f2]
00411311  fstp        dword ptr [f3]

f1 * f2 is computed in extended precision — enough bits to hold its true value — but that extra precision is lost when stored in f3.

Comparing f3 and f1 * f2

The answer is in the assembler code, so let’s get right to it:

 if (f3 != f1 * f2)
00411314  fld         dword ptr [f1]
00411317  fmul        dword ptr [f2]
0041131A  fld         dword ptr [f3]
0041131D  fucompp

Again, f1 * f2 is computed in extended precision, but this time its true value is retained — it is left on the stack. f3 is then loaded onto the stack (of course it still has only 24 bits of precision, even though it’s been “promoted” to extended precision). The two values on top of the stack are then compared. Clearly, they differ.

Analyzing Program 2

Program 2 “fails” for the same reason as program 1, except that the “error” is magnified by the conversion of f3 and f1 * f2 to integers: the integer part of f3 is 7, and the integer part of f1 * f2 is 6.

Computing i1

The binary values of the three floating point variables in the computation are:

• f1 = 0.101100110011001100110011
• f2 = 1010
• f3 = 111

The value of f3 is an integer. To understand why, let’s calculate the true value of f1 * f2:

  0.101100110011001100110011
  x                     1010
  --------------------------
                           0
   101100110011001100110011
                         0
 101100110011001100110011
----------------------------
110.11111111111111111111111

So f1 * f2 = 110.11111111111111111111111. To assign it to f3 it must be rounded. f1 * f2 has 26 significant bits, and rounding it to the nearest 24 bit value makes it 111, or 7 decimal. Clearly, this means i1 will be 7 as well.

Here’s the assembler code:

 f3 = f1 * f2;
0041130B  fld         dword ptr [f1]
0041130E  fmul        dword ptr [f2]
00411311  fstp        dword ptr [f3]
 i1 = (int)f3;
00411314  fld         dword ptr [f3]
00411317  call        @ILT+155(__ftol2_sse) (4110A0h)
           ||
           ||
           \/
--- f:\dd\vctools\crt_bld\SELF_X86\crt\prebuild\tran\i386\ftol2.asm
00411600  cmp         dword ptr [___sse2_available (416554h)],0
00411607  je          _ftol2 (411636h)
00411609  push        ebp
0041160A  mov         ebp,esp
0041160C  sub         esp,8
0041160F  and         esp,0FFFFFFF8h
00411612  fstp        qword ptr [esp]
00411615  cvttsd2si   eax,mmword ptr [esp]
0041161A  leave
0041161B  ret

Computing i2

Here’s the assembler code for computing i2:

 i2 = (int)(f1 * f2);
0041131F  fld         dword ptr [f1]
00411322  fmul        dword ptr [f2]
00411325  call        @ILT+155(__ftol2_sse) (4110A0h)
           ||
           ||
           \/
--- f:\dd\vctools\crt_bld\SELF_X86\crt\prebuild\tran\i386\ftol2.asm
00411600  cmp         dword ptr [___sse2_available (416554h)],0
00411607  je          _ftol2 (411636h)
00411609  push        ebp
0041160A  mov         ebp,esp
0041160C  sub         esp,8
0041160F  and         esp,0FFFFFFF8h
00411612  fstp        qword ptr [esp]
00411615  cvttsd2si   eax,mmword ptr [esp]
0041161A  leave
0041161B  ret

The interesting thing in the calculation of i2 is that f1 * f2 is stored in a double before being cast to an integer (see the highlighted fstp instruction). f1 * f2 in double precision is 110.11111111111111111111111; that is, its true 26 bit value. Casting this to an integer results in 110, or decimal 6. If the compiler had stored f1 * f2 in a float and used the cvttss2si instruction instead, the answer would have been 7 — just like in the first calculation.

Discussion

In program 2, it would appear that single-precision is more accurate than extended precision. After all, forcing the value into a float gives the expected answer. This is just a happy coincidence. Two losses of precision — the conversion of 0.7 to floating-point and the rounding up of the product — have effectively canceled each other out.

Behavior Depends on the Processor and How It Is Used

This anomaly may not occur on your machine. It depends on your processor, and in particular, the instructions used and the mode that it’s in. My example programs were compiled into Intel x87 FPU instructions, and ran with the x87 precision control field set to 53-bits. (I said above that my programs were using extended precision; technically they weren’t, but double precision is sufficient to make them “fail”. The true values of f1 * f2 are less than 53 bits.)

Intel processors can also do floating point using SSE (Streaming SIMD Extensions) instructions. I recompiled my programs using the Visual C++ compiler options /arch:SSE (single-precision floating-point) and /arch:SSE2 (double-precision floating-point). For the SSE option, there was no change; the compiler, at its discretion, still decided to generate x87 instructions. (See the comment below about the Mac’s use of SSE instructions making the anomaly disappear.)

Recompiling with /arch:SSE2 I got different assembler code, but the same output; here’s the assembler code for program 1:

f3 = f1 * f2;
00411313  cvtss2sd    xmm0,dword ptr [f1]
00411318  cvtss2sd    xmm1,dword ptr [f2]
0041131D  mulsd       xmm0,xmm1
00411321  cvtsd2ss    xmm0,xmm0
00411325  movss       dword ptr [f3],xmm0
if (f3 != f1 * f2)
0041132A  cvtss2sd    xmm0,dword ptr [f1]
0041132F  cvtss2sd    xmm1,dword ptr [f2]
00411334  mulsd       xmm0,xmm1
00411338  cvtss2sd    xmm1,dword ptr [f3]
0041133D  ucomisd     xmm1,xmm0

The SSE double-precision instructions are used. f3 — single precision — is compared to the double-precision intermediate result f1 * f2, resulting in a mismatch.

Please Try it Out

If you have access to a different compiler or processor, please try these programs out. Let me know what you find!

The Message

You can’t count on floats being handled as single precision values; they can be processed in double or extended precision, as dictated by your compiler and CPU.

(I have written a companion article called “When Doubles Don’t Behave Like Doubles”.)

By Rick Regan (Copyright © 2008-2010 Exploring Binary)

When Floats Don’t Behave Like Floats