While debugging some code, I came across x87 behavior that I don't understand. Consider the following program:
#include <iostream>
using namespace std;
double test1(double x,double y)
{
__asm {
fld x
fsin
fld y
fsin
fsubp st(1),st(0)
}
}
int main(int argc, char* argv[])
{
cout << test1( 1, 1)<<endl;
return 0;
}
Everything is fine here, the output is 0. Now we slightly modify the test function so that the intermediate result is stored in memory:
double test2(double x,double y)
{
double tmp;
__asm {
fld x
fsin
fld y
fsin
fstp tmp
fld tmp
fsubp st(1),st(0)
}
}
Now the output is 1.78893e-18. This result can be explained by rounding the intermediate value stored in memory. But I don’t understand why the result is non-zero even if the precision of x87 floating-point calculations is set to double or single and regardless of rounding mode:
double test3(double x,double y)
{
// Double precision,
// round toward zero
unsigned short cw = 0xe7f;
double tmp;
__asm {
fldcw cw
fld x
fsin
fld y
fsin
fstp tmp
fld tmp
fsubp st(1),st(0)
}
}
The fsubp instruction should consider the precision mode and ignore unnecessary bits, right? Why, then, is the result non-zero in any rounding mode? After all, in at least one of the rounding modes, rounded bits must match?
