detecting underflow exception on Linux.

[alok]

I am getting inconsistent behvior on Linux and Solaris platfors while
computing doule ( 64 bit precision ) multiplications.
I have following two double numbers whose integer representation
is as following

I have a union

typedef union {
double double_val;
unsigned long long uint_val;
} double_reg;

double_reg a, b;
a.uint_val = 0x000ffffffffffff8ULL;
b.uint_val = 0xbff0000000000001ULL;

c = a.double_val * b.double_val;

Solaris raises UNDERFLOW Exception on this but Linux does not.
On Linux i am using following function from fenv.h to check the
exceptions on solaris.
fetestexcept (FE_ALL_EXCEPT)

The rounding mode for the above multiplication is round to Minus
Infinity which is set by
fesetround(FE_DOWNWARD)

Please suggest what am i doing wrong or i am missing any compile switch
such as to get Underflow exception on Linux as well. I am using gcc
version 3.3.3

THANKS IN ADVANCE.

Following is the Linux Program dmul_linux.c
I use
gcc dmul_linux.c -lm
command to compile the program
=======================================
#include <stdio.h>
#include <math.h>
#include <fpu_control.h>
#include <fenv.h>

// LINUX

// FE_INEXACT 0x20 inexact result
// FE_DIVBYZERO 0x4 division by zero
// FE_UNDERFLOW 0x10 result not representable due to underflow
// FE_OVERFLOW 0x8 result not representable due to overflow
// FE_INVALID 0x1 invalid operation

//FE TONEAREST IS 0
//FE UPWARD IS 800
//FE DOWNWARD IS 400
//FE TOWARDZERO IS c00

typedef union {
double double_val;
unsigned long long uint_val;
} double_reg;

int main()
{

#ifdef linux
/*
This puts the X86 FPU in 64-bit precision mode. The default
under Linux is to use 80-bit mode, which produces subtle
differences from FreeBSD and other systems, eg,
(int)(1000*atof("0.3")) is 300 in 64-bit mode, 299 in 80-bit
mode.
*/
fpu_control_t cw;
_FPU_GETCW(cw);
cw &= ~_FPU_EXTENDED;
cw |= _FPU_DOUBLE;
_FPU_SETCW(cw);
#endif

double_reg a,b,c;
int d;
a.uint_val = 0x000ffffffffffff8ULL;
b.uint_val = 0xbff0000000000001ULL;
fesetround(FE_DOWNWARD);
printf(" ROUNDING is %x \n", d);
feclearexcept(FE_ALL_EXCEPT);
d = fetestexcept (FE_ALL_EXCEPT);
printf (" D is %d \n", d);
c.double_val = (a.double_val) * (b.double_val);
d = fetestexcept (FE_ALL_EXCEPT);
printf (" D is %x \n", d);
printf (" DOUBLE RESULT = %f HEX RESULT = %llx \n", c.double_val,
c.uint_val);

};


Following is the Solaris dmul_sol.c I use
gcc dmul_sol.c
to compile the program
===================================================================

#include <ieeefp.h>

//#define FP_X_INV 0x10 /* invalid operation exception */
//#define FP_X_OFL 0x08 /* overflow exception */
//#define FP_X_UFL 0x04 /* underflow exception */
//#define FP_X_DZ 0x02 /* divide-by-zero exception */
//#define FP_X_IMP 0x01 /* imprecise (loss of precision) */

typedef union {
double double_val;
unsigned long long uint_val;
} double_reg;

int main()
{
double_reg a,b,c;
unsigned long long a_int, b_int;
int d;
a.uint_val = 0xffffffffffff8ULL;
b.uint_val = 0xbff0000000000001ULL;
a_int = 0xffffffffffff8ULL;
b_int = 0xbff0000000000001ULL;
fpsetsticky(0);
fpsetround(FP_RM);

c.double_val = (a.double_val) * (b.double_val);
//c.double_val = (double) (a_int) * (double) (b_int);
d = fpgetsticky();
printf (" EXCEPTIONS are %d \n", d);
printf (" DOUBLE RESULT = %lf HEX RESULT = %llx \n", c.double_val,
*(unsigned long long *)&(c.double_val));

};

 

[Chris Croughton]

I am getting inconsistent behvior on Linux and Solaris platfors while
computing doule ( 64 bit precision ) multiplications.
I have following two double numbers whose integer representation
is as following

I have a union

typedef union {
double double_val;
unsigned long long uint_val;
} double_reg;

double_reg a, b;
a.uint_val = 0x000ffffffffffff8ULL;
b.uint_val = 0xbff0000000000001ULL;

c = a.double_val * b.double_val;

You have just caused undefined behaviour. That means that anything at
all could happen, including nothing, underflow, overflow, the processor
going up in a puff of smoke, or the start of World War 3. Oh, or the
traditional demons flying out of your nose.

Solaris raises UNDERFLOW Exception on this but Linux does not.

That's fun. Although not as interesting as nasal demons.

Please suggest what am i doing wrong or i am missing any compile switch
such as to get Underflow exception on Linux as well. I am using gcc
version 3.3.3

You are writing a long long value into a union, and reading it as a
double. That is what you are doing wrong, doing that is Undefined
Behaviour. Messing about with the FPU in a non-portable way probably
isn't helping either (if it isn't in the mode the compiler expects
anything could happen).

Chris C

 

[Chris Torek]

I am getting inconsistent behvior on Linux and Solaris platfors while
computing doule ( 64 bit precision ) multiplications. ...

As someone else noted already, you cannot depend on this sort of
thing in the first place -- there are no guarantees in portable C
about the representation of floating-point numbers. (Even if you
know your implementation uses IEEE float, you cannot predict the
bit and byte order if you attempt to map them to integral types.)

You have several other problems:

#ifdef linux
/*
This puts the X86 FPU in 64-bit precision mode. The default
under Linux is to use 80-bit mode, which produces subtle
differences from FreeBSD and other systems, eg,
(int)(1000*atof("0.3")) is 300 in 64-bit mode, 299 in 80-bit
mode.
*/
fpu_control_t cw;
_FPU_GETCW(cw);
cw &= ~_FPU_EXTENDED;
cw |= _FPU_DOUBLE;
_FPU_SETCW(cw);
#endif

First, "Linux" does not necessarily mean "x86". (You will find
Linux running on PowerPC, SPARC, and other architectures these
days.)

Second, and perhaps most important, if you *are* on an x86, you
need to be aware that the FPU control word does not solve the
problem. It does limit the mantissa (fraction precision), but it
has no effect on the exponent (range), which always runs in full
80-bit extended mode. This means that some operations do not
overflow or underflow until "too late" (conversion from internal
FPU representation to in-memory representation). You need to obtain
and peruse the Intel CPU documentation, and then find out how your
compiler(s) use it, none of which has much to do with Standard C.
(The only real connection here is how common it is to find that C
compilers do not actually implement Standard C on the x86, because
doing so is too slow -- moving temporary results out of the FPU
and then back in is quite slow.)