Why does isnormal () say the value is normal when it is not?
#include <stdlib.h>
#include <stdio.h>
#include <float.h>
#include <math.h>
void PrintBytes( const float value )
{
const char* const byte = ( const char* )&value ;
for( size_t i = 0 ; i < sizeof( value ) ; i++ )
{
printf( "%02hhx" , byte[i] );
}
}
int main(void)
{
float value = FLT_MIN;
while( 1 )
{
printf( "%e %d " , value , isnormal( value ) );
PrintBytes( value );
puts( "" );
if( !isnormal( value ) )
{
break;
}
value /= 2.0F;
}
return 0;
}
Output:
1.175494e-038 1 00008000
5.877472e-039 1 00004000
2.938736e-039 1 00002000
1.469368e-039 1 00001000
7.346840e-040 1 00000800
3.673420e-040 1 00000400
1.836710e-040 1 00000200
9.183550e-041 1 00000100
4.591775e-041 1 00800000
2.295887e-041 1 00400000
1.147944e-041 1 00200000
5.739719e-042 1 00100000
2.869859e-042 1 00080000
1.434930e-042 1 00040000
7.174648e-043 1 00020000
3.587324e-043 1 00010000
1.793662e-043 1 80000000
8.968310e-044 1 40000000
4.484155e-044 1 20000000
2.242078e-044 1 10000000
1.121039e-044 1 08000000
5.605194e-045 1 04000000
2.802597e-045 1 02000000
1.401298e-045 1 01000000
0.000000e+000 0 00000000
It is clear that the second value 5.877472e-039
is subnormal, as its index becomes 0 00004000
.
Ideone gives the correct result:
1.175494e-38 1 00008000
5.877472e-39 0 00004000
I am compiling my code using gcc (MinGW-w64) on Windows.
source to share
It seems to be a bug in MinGW-w64 with 32-bit target.
My conclusion:
- MinGW-w64 x86_64-4.9.2-win32-seh-rt_v3-rev1: Correct
- MinGW-w64 i686-4.9.2-win32-dwarf-rt_v4-rev2: Wrong
- Cygwin i686-pc-cygwin (gcc 4.9.2): Correct
Looking at the MinGW-W64 headers:
#define isnormal(x) (fpclassify(x) == FP_NORMAL)
Then we have:
#define fpclassify(x) (sizeof (x) == sizeof (float) ? __fpclassifyf (x) \
: sizeof (x) == sizeof (double) ? __fpclassify (x) \
: __fpclassifyl (x))
Function implementation __fpclassifyf
:
__CRT_INLINE int __cdecl __fpclassifyf (float x) {
#ifdef __x86_64__
__mingw_flt_type_t hlp;
hlp.x = x;
hlp.val &= 0x7fffffff;
if (hlp.val == 0)
return FP_ZERO;
if (hlp.val < 0x800000)
return FP_SUBNORMAL;
if (hlp.val >= 0x7f800000)
return (hlp.val > 0x7f800000 ? FP_NAN : FP_INFINITE);
return FP_NORMAL;
#else
unsigned short sw;
__asm__ __volatile__ ("fxam; fstsw %%ax;" : "=a" (sw): "t" (x));
return sw & (FP_NAN | FP_NORMAL | FP_ZERO );
#endif
}
Since the x64 version works, but the i686 version does not work, I think we should blame the above line fxam; fstsw
. I started a new new thread for this question; an answer has so far been given that fxam
checks a version that has been expanded to 80 bit precision.
The generated assembly in -O1
the broken version:
call ___main
flds LC1
fstps 24(%esp)
L7:
flds 24(%esp)
/APP
# 484 "F:/Prog/mingw-w64/i686-4.9.2-win32-dwarf-rt_v4-rev2/mingw32/i686-w64-mingw32/include/math.h" 1
fxam; fstsw %ax;
# 0 "" 2
/NO_APP
andw $17664, %ax
cmpw $1024, %ax
sete %al
movzbl %al, %eax
movl %eax, 12(%esp)
fstpl 4(%esp)
movl $LC2, (%esp)
call _printf
flds 24(%esp)
fstps (%esp)
call _PrintBytes
movl $LC3, (%esp)
call _puts
flds 24(%esp)
/APP
# 484 "F:/Prog/mingw-w64/i686-4.9.2-win32-dwarf-rt_v4-rev2/mingw32/i686-w64-mingw32/include/math.h" 1
fxam; fstsw %ax;
# 0 "" 2
/NO_APP
andw $17664, %ax
cmpw $1024, %ax
jne L6
fmuls LC4
fstps 24(%esp)
jmp L7
L6:
fstp %st(0)
movl $0, %eax
leave
.cfi_restore 5
.cfi_def_cfa 4, 4
ret
.cfi_endproc
LFE41:
.section .rdata,"dr"
.align 4
LC1:
.long 8388608
.align 4
LC4:
.long 1056964608
source to share