gcc operational discrepancy

[Dale Williamson <Dale dot Williamson at trw dot com>]

/* File badbit1.c

June 5, 2000

On May 30, I reported this operational discrepancy from a program using
the gcc compiler:

   Using the gnu C compiler, it has been found that under certain cir-
   cumstances the bits of a double number returned from a function may
   not be reliably transferred to a variable across an equal sign, as 
   in the expression: (double)X = (double)func().

   The problem does not occur on other platforms and appears to be re-
   lated to the gcc compiler and x86 floating point hardware. 

   The program on this file shows a problem where bit 52 (counting from
   0 at leftmost) is changed from 0 to 1 under certain conditions: 

      gcc 2.7.2.3: when a double variable is equated to a function 
      that returns a double word value

      gcc 2.91.66: both when a double variable is equated to a func-
      tion that returns a double word value, and when a double vari-
      able is equated to another double variable

Responses from those taking an interest in this problem, allowing it to
be understood and circumvented, are greatly appreciated.  Thank you.

From all this it has been concluded that the original program submitted,
and the equivalent but simpler one given below, both contain the code 
Y=X that will exhibit an operational discrepancy that is unique to the 
gcc compiler on x86 machines.  

The operational discrepancy is one in which Y does not contain the bits
of X after the equate Y=X and both X and Y are words of matching type, 
double.

The upshot is that users of the x86 gcc compiler need to be aware that 
bit patterns carried in double words cannot simply be words.  All 64 
bits of words in double words must also match valid IEEE floating point 
numbers to avoid this operational discrepancy.

This note is being written to summarize my understanding of the problem,
to pose an ISO standards question that it raises, and to share lessons
learned with others who might also benefit.

Others benefiting from this are likely to be those who are moving pro-
grams onto x86 machines now that operating systems that work like unix 
are becoming available.  

To summarize again the discrepancy, results from a simple program (pro-
vided below) show cases A.2.2 and B.2 where bit 52 changes from 0 to 1 
in what is intended to be a simple memory-to-memory transfer of bits:

Case A.

  A.1 Original bits X:
 11111111 11111111 11111111 11111111 11111111 11111111 11110111 11111111

  A.2.1 (gcc 2.7.2.3) Y bits after Y=X, equate double Y to double X:
 11111111 11111111 11111111 11111111 11111111 11111111 11110111 11111111

  A.2.2 (gcc 2.91.66) Y bits after Y=X, equate double Y to double X:
 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111

Case B (gcc 2.7.2.3 and gcc 2.91.66).

  B.1 Original fX()=X bits:
 11111111 11111111 11111111 11111111 11111111 11111111 11110111 11111111

  B.2 Y bits after Y=fX(), equate double Y to double function fX():
 11111111 11111111 11111111 11111111 11111111 11111111 11111111 11111111

Case C (gcc 2.7.2.3 and gcc 2.91.66).

  C.1 Original fL()=X bits:
 11111111 11111111 11111111 11111111 11111111 11111111 11110111 11111111

  C.2 U.d bits after U.l=fL(), equate long long U.l to long long fL():
 11111111 11111111 11111111 11111111 11111111 11111111 11110111 11111111

From testing done and inputs received from others, it is believed that
the problem stems from certain bit patterns stored in double words--spe-
cifically ones that do not match valid IEEE floating point numbers--and
linking of the gcc compiler to the x87 numeric processor--an important 
part of the x86 machine--for all processes involving double words.

Here is some background:

Floating point numbers of 64-bits under the IEEE standard, familiar to
those working on a variety of machines, probably began jointly with the
efforts of Intel and development of the x87 numeric processor.

Bit pattern format--numbers of bits to define mantissa and exponent--of
machines using the IEEE standard is identical in x87 and other proces-
sors although the mantissa, exponent, and sign bits might be held in
differing regions of a 64-bit word.  In the x87 bit patterns shown above
and counting from 0 at the left, bits 48-51 and 57-63 are the 11 expo-
nent bits, bit 56 is the sign bit, and the remaining 52 are mantissa 
bits.

For those unfamiliar with the x87 processor, it is distinguished by a 
stack of eight 80-bit registers.  When a 64-bit IEEE floating point 
number is loaded into a register, the bits are converted to an 80-bit 
type called temporary real.  

When a temporary real is stored back into 64-bit memory, rounding is
done to fit the mantissa into 52 bits, and the exponent is biased to 
cover 11 bits in the range 2^1024 (INF) when all 11 bits are 1, and 
2^-1022 when just the lowest bit is set.

Processors using IEEE fp, not just x87, also recognize certain bit pat-
terns as not-a-number, NAN.  When these special patterns--which feature
all exponent bits of 1--are stored into 64-bit memory from an x87 reg-
ister, the mantissa bits are not rounded but simply chopped to fit the 
52-bit width, and the exponent is stored as-is (all 1s), to flag the 
value as NAN.

With this background in mind, here is the important x86 difference that
Alan Modra makes very clear:

=> and it's part of the x86 ABI to return a double from a function in a
=> floating point register (what a revolutionary concept!)
=> Thus you need to know something about the underlying hardware if
=> you want to examine bit patterns.  This is definitely not a gcc bug.

Since the compiler implementation always returns a double from an x87 
register, the following occurs for the expression (double)Y=(double)X:

   X is sent to the numeric processor and converted into an 
      80-bit temporary real

   Y receives X from the 80-bit register under one of the two 
      processes for returning 64-bits described above

Saying Y=X leads to more than a simple memory-to-memory transfer and 
Mr. Modra's advice concerning this operational discrepancy--it ain't
a bug--is right-on: you need to know something about the underlying 
hardware.

In fact, if the bit pattern of double word X happens to match a NAN--it
ain't a number--then sending it first to an 80-bit x87 register will re-
sult in the bit chopping described above when it is sent back to memory.

This is likely to cause Y to receive bits that are different from the 
original bits of X, and it is no surprise that the result of Y=X can be 
Y not-equal-to X as some of the bit patterns above demonstrate.

Whether other implementations of C compilers on x86 share this feature
for memory-to-memory transfers of double words is not known.

Nor is it known if allowing bits of the same type of word to change 
across a simple equate conforms with ISO standards.  The problem here 
raises this issue, and I would look to Martin v. Loewis, whose input 
regarding gcc and a different ISO topic is appreciated, for guidance on
this.  

A standards question for this problem might be:

   While one should expect an equate such as 
      (long long int)Y=(double)X 
   to yield Y with bits not the same as those of X, since X and Y 
   are of different types, under what conditions should one expect 
   an equate of identical types, such as 
      (double)Y=(double)X 
   to produce bits of Y that do not match bits of X?

For some the answer would be: never.

For x86 and gcc and doubles, at least one answer is: when X is NAN.  
Perhaps the standards say something like this about doubles too.

Regardless of standards, the fact is that the widely used gcc compiler 
works this way and probably for good reasons.  

Utilizing the x87 is a good thing: It is a jewel that has been buried
for a long time under an operating system evolving toward making the 
internals of these machines less accessible rather than more.  Now this
is changing, and programs from other IEEE machines can easily get at 
the x87.

The caution is: When working with general bit patterns--specifically 
double word types that actually contain words and not just floating 
point numbers (not a very revolutionary concept)--the lesson here is to
ensure that they are not sent to the numeric processor where they will 
undergo the transformation to 80-bits and may be mistaken for NANs.

As an example, this adds complication to an implementation that uses a 
64-bit virtual stack typed for double words and containing at various 
times elements that truly are words and not floating point numbers.

Why build such a stack of double words first place?  Because whenever 
numbers are upon it their patterns are in ready form for IEEE floating 
point processing--on any IEEE machine, not just x86.

For such a stack to encompass x86, equates of elements when they really
are not fp numbers would need to be done more carefully to avoid having
them, and their by-products, sent to the x87.

One way to make such a stack work properly on the x86 platform, an idea
suggested by Al Danial, is to masquerade double words as long long ints
when they do not hold floating point numbers.

In the program below, case C uses a union to accomplish this, as results
shown above confirm. 

Thanks to all for the prompt attention that allowed this problem to be 
quickly resolved.  The question it raises about ISO conformity remains.

Sincerely yours,

Dale R. Williamson

dale.williamson@trw.com

*/

#include <stdio.h>

static double X=0;
static long long int L=0; 

char *bitpat(unsigned char c)
/* Converting byte c into an 8 character bit pattern for display. */
{
   static unsigned char buf[9]="        \0";
   unsigned char c1[8]={'0','0','0','0','0','0','0','0'};
   unsigned char mask[8]={128,64,32,16,8,4,2,1};
   register int i=0;

   for(;i<8;i++) if(c&*(mask+i)) *(c1+i)='1';
   memcpy(buf,c1,8);
   return buf;
}

void bitprint(char* c, int n)
/* Displaying the bits of n bytes in c. */
{
   register int k=0;

   printf("\n\r");
   for(;k<n;k++) printf(" %s",bitpat(*(c+k)));
   printf("\n\r");
}

double fX()
/* Function returning a double number. */
{
   printf("\n\r  B.1 Original fX()=X bits:");
   bitprint((char *)&X,8);
   return X; 
} 

long long int fL() 
/* Function returning a long long int. */
{ 
   memcpy((char *)&L,(char *)&X,8);
   printf("\n\r  C.1 Original fL()=X bits:");
   bitprint((char *)&L,8);
   return L; 
}

void main()
/* Examining bit patterns of simple equates. */
{ 
   int I[2]={0xFFFFFFFF,0xFFF7FFFF}; /* a NAN bit pattern */

   union { /* 64-bit shared memory location: */
      double d;
      long long int l;
   } U;

   double Y;

   memcpy((char *)&X,(char *)&I,8);

/* Case A: X into Y */
   Y=X; 

   printf("\n\r  A.1 Original bits X:");
   bitprint((char *)&X,8);
   printf("  A.2 Y bits after Y=X, equate double Y to double X:");
   bitprint((char *)&Y,8);

/* Case B: fX() into Y */
   Y=fX(); 

   printf("  B.2 Y bits after Y=fX(),");
   printf(" equate double Y to double function fX():");
   bitprint((char *)&Y,8);

/* Case C: fL() into long long U.l shared by double U.d */
   U.l=fL(); 

   printf("  C.2 U.d bits after U.l=fL(),");
   printf(" equate long long U.l to long long fL():");
   bitprint((char *)&U.d,8);

   printf("\n\r");
}