#602 __mingw_printf invalid output for nan value

[Duncan Ogilvie]

This code:

#include <math.h>
#include <stdio.h>
#include <string.h>

int main()
{
    long double ld = nanl("nan");
    __mingw_printf("nan (%d): %Lf\n", isnan(ld), ld);
    unsigned char data[10]={0x5b, 0x01, 0x04, 0x5e, 0x85, 0x00, 0x00, 0x00, 0xd8, 0x59};
    memcpy(&ld, data, 10);
    __mingw_printf("nan (%d): %Lf\n", isnan(ld), ld);
}

Will output the following on Windows 8.1 x64 with g++ (x86_64-win32-seh-rev2, Built by MinGW-W64 project) 6.3.0:

nan (1): nan
nan (1): 51121152321801233769528628351853187448384447522053840083788534546827429704643290381340679219510392362325477918202365924663653936157016779794625974445014197143562076732008008738128391263363163574965197975470359672288313174823442406254048743002805452006157937878331142146798072696855831972926249999727974322449691877753828930695850548966676286079757955984128191916481949201177743973903996572963765928618534081793280129105149856531260982551553469556429687623577078046545852198566211939155677028446451165609819075837450164249776102845298460469771918289548752998338223854607227934434765179856884186807766927334957254748394233630664914775000455646748666717708093934285269156086702020354011540591392665589375067807455867173656994975148894173206984974769181739207306862381317984266053473189056331126277978170416496979674173995459863769696678139373469071902631219029009882566470150137655735124177599243342172252211543365526650540919950000331153682021554718999864489566506339904693091924474536686798677359626831249371422914336837875642156815196149795381551612204397230483764946192534114535346234814185380840663925616737129762828956809594806005021218926647104971253221663691633803548451583289405163859070274525946904733996109568233540983324711005688985472233093988319080327085238864951989022529775137959073989229879976005224619428151507207204127577684739182214549342293405410680249028464282743001152759784716025069416795757606758511947646024601325845023201020928965906189624416406030848379505691504516615002554935359779127043369953207671477905398294201081775159338488386294932914644247061475060610844251702566871004568191800019003391839796982595646233069513712001461158765658211174522442116616289845700046396282353720595469749930911408197188628124217336603684988877899113364285268336103879082322140847309239491610769009092472379988705831151518195226501624576728449337637413597709101054377371751575601666492702420321446551893642494512728924667549191793875921344186506961309343308327172881506107392.000000

It works fine on linux (ideone uses gcc 6.3, but I also tested on lbunutu with gcc 6.2):

https://ideone.com/CqqtDe

I'll try to find the relevant source file but I have no clue where to start...

[Duncan Ogilvie]

Okay, the solution is to either fix __fpclassifyl (x64 uses a bad software implementation) or just use __builtin_fpclassify:

#include <math.h>
#include <stdio.h>
#include <string.h>

#  define fufu(x) __builtin_fpclassify (FP_NAN, FP_INFINITE,        \
     FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x)

const char* show_classification(long double x) {
    switch(fufu(x)) {
        case FP_INFINITE:  return "Inf";
        case FP_NAN:       return "NaN";
        case FP_NORMAL:    return "normal";
        case FP_SUBNORMAL: return "subnormal";
        case FP_ZERO:      return "zero";
        default:           return "unknown";
    }
}

int main()
{
    long double ld = nanl("nan");
    __mingw_printf("nan (%d, %s)\n", __isnanl(ld), show_classification(ld));
    unsigned char data[10]={0x5b, 0x01, 0x04, 0x5e, 0x85, 0x00, 0x00, 0x00, 0xd8, 0x59};
    memcpy(&ld, data, 10);
    __mingw_printf("nan (%d, %s)\n", __isnanl(ld), show_classification(ld));
}

Which prints:

nan (-1, NaN)
nan (0, NaN)

[Duncan Ogilvie]

I have confirmed this fix, by replacing the implementations of __fpclassifyl with __builtin_fpclassify everything works correctly, but I don't know how to properly fix this (there are various direct uses of __fpclassifyl in mingw_pformat.c so my solution would be to replace those with fpclassify and then replace the macro for fpclassify with a call to __builtin_fpclassify).

[LIU Hao]

The number in question is an invalid operand, because its exponent is neither all-zeroes nor all-ones and its bit 63 is zero.

[Duncan Ogilvie]

Hm okay, so is nan an incorrect classification or is that what should be looked at for the software implementation?

[Duncan Ogilvie]

Okay I found and the bug in the software implementation (it also appears to exist in glibc by the way). I still believe __builtin_fpclassify should be called in the fpclassify macro and the __fpclassifyX family of functions should never be used directly.

This is the code for math.h (I'm pretty sure the other software implementations are also broken, for long double but also the other floating point types):

#ifndef __CRT__NO_INLINE
  __CRT_INLINE int __cdecl __fpclassifyl (long double x) {
#if defined(__x86_64__) || defined(_AMD64_)
    __mingw_fp_types_t hlp;
    unsigned int e, h;
    hlp.ld = &x;
    e = hlp.ldt->lh.sign_exponent & 0x7fff;
    h = hlp.ldt->lh.high;
    if (!e)
      {
        if (!(hlp.ldt->lh.low | h))
          return FP_ZERO;
        else if (!(h & 0x80000000))
          return FP_SUBNORMAL;
      }
    else if (e == 0x7fff)
      return (((hlp.ldt->lh.high & 0x7fffffff) | hlp.ldt->lh.low) == 0 ?
              FP_INFINITE : FP_NAN);
    return ((h & 0x80000000) == 0) ? FP_NAN : FP_NORMAL;
#elif defined(__arm__) || defined(_ARM_)
    return __fpclassify(x);
#elif defined(__i386__) || defined(_X86_)
    unsigned short sw;
    __asm__ __volatile__ ("fxam; fstsw %%ax;" : "=a" (sw): "t" (x));
    return sw & (FP_NAN | FP_NORMAL | FP_ZERO );
#endif
  } 

[LIU Hao]

If we weren't emulating it, a hardware solution is available here (the code is purely in Intel syntax and requires a little modification if you are using AT&T syntax).

Nevertheless, fxam could still say the operand is unsupported by setting C3, C2 and C0 to all zeroes, which is impossible on 32-bit, single precision or 64-bit, double precision floating point numbers, both of which don't have the extra bit in their mantissas.

Also have a look at the C99 standard (ISO/IEC WG14 N1256):

7.12 Mathematics <math.h>
6 The number classification macros
FP_INFINITE
FP_NAN
FP_NORMAL
FP_SUBNORMAL
FP_ZERO
represent the mutually exclusive kinds of floating-point values. They expand to integer
constant expressions with distinct values. Additional implementation-defined floatingpoint
classifications, with macro definitions beginning with FP_ and an uppercase letter,
may also be specified by the implementation.

The set of standard number classification macros does not contain a value for invalid operands, I am looking forward to an implementation-defined one for it, as the 80-bit, extended precision long double has already been implementation-defined.

[Duncan Ogilvie]

The builtin classifier for gcc is very clear: that value is nan. The 32 bit
inline assembly is also very clear: nan. These macros are based on the
hardware implementation so following the hardware is definitely correct I
think. Another argument for using what the hardware tells you is that if
you don't you get the kind of garbage output I got because the value is
invalid but you try to do something with it anyway.

[Duncan Ogilvie]

To clarify: the hardware implementation I'm talking about is also in that
function for 32 bit.

[LIU Hao]

No that value is not a NaN.

A NaN can be copied without generating an exception, no matter whether it is a signaling NaN (SNaN) or a quiet NaN (QNaN). It is just that if you don't do arithmetic operations on it, you don't get an exception. If you do arithmetic operations on a QNaN, you get a QNaN. If you do arithmetic operations on a SNaN you get an #IA exception (the exception can be masked, but it is generated unconditionally).

This invalid number, on the other hand, can't be copied. you can FST it, but x87 will not actually store it but generate an #IA instead. This can be observed by unmasking the LSB of the x87 control word using FSTCW before FST.

The hardware isn't lying either. FXAM tells you that the number is unsupported. FXAM doesn't say the number is a NaN, in which case bytewise and'ing the x87 status word with 0x4500 would yield 0x0100 instead of 0x0000.

[LIU Hao]

FST does not generate an #IA exception for a SNaN in 80-bit format, but it does generate an exception for a SNaN if the memory location is in 64-bit or 32-bit format.

[Duncan Ogilvie]

Hm so then there is a bug in GCC... Or because it's implementation defined they chose to classify this as NAN :-D

http://x86.renejeschke.de/html/file_module_x86_id_125.html

I'm currently rather unclear about how this table relates to the values of FP_NAN, etc. but I can add FP_UNSUPPORTED. Could you check if my current implementation of that wikipedia table is correct?

int myclassify(long double ld)
{
    __mingw_ldbl_type_t fuck = *(__mingw_ldbl_type_t*)&ld;
    //https://en.wikipedia.org/wiki/Extended_precision#x86_extended_precision_format
    auto bit63 = (fuck.lh.high & 0x80000000) != 0;
    switch (fuck.lh.sign_exponent & 0x7fff) //Exponent
    {
        case 0: //Exponent: All Zeros
        {
            if (!bit63) //Bit 63: Zero
            {
                if (fuck.lh.low == 0 && (fuck.lh.high & 0x7fffffff) == 0) //Bits 62-0: Zero
                {
                    return FP_ZERO; //Zero. The sign bit gives the sign of the zero.
                }
                else //Bits 62-0: Non-zero
                {
                    return FP_SUBNORMAL; //Denormal. The value is (−1)^s * m * 2^(−16382)
                }
            }
            else //Bit 63: One
            {
                return FP_SUBNORMAL; //Pseudo Denormal. The 80387 and later properly interpret this value but will not generate it. The value is (−1)^s * m * 2^(−16382)
            }
        }

        case 0x7fff: //Exponent: All Ones
        {
            auto bit62 = (fuck.lh.high & 0x40000000) != 0;
            auto bit61to0_zero = (fuck.lh.high & 0x3fffffff) == 0 && fuck.lh.low == 0;
            if (!bit63 && !bit62) //Bits 63,62: 00
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    return FP_UNSUPPORTED; //Pseudo-Infinity. The sign bit gives the sign of the infinity. The 8087 and 80287 treat this as Infinity. The 80387 and later treat this as an invalid operand.
                }
                else //Bits 61-0: Non-zero
                {
                    return FP_UNSUPPORTED; //Pseudo Not a Number. The sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number. The 80387 and later treat this as an invalid operand.
                }
            }
            else if (!bit63 && bit62) //Bits 63,62: 01
            {
                //Bits 61-0: Anything
                return FP_UNSUPPORTED; //Pseudo Not a Number. The sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number. The 80387 and later treat this as an invalid operand.
            }
            else if (bit63 && !bit62) //Bits 63,62: 10
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    return FP_INFINITE; //Infinity. The sign bit gives the sign of the infinity. The 8087 and 80287 treat this as a Signaling Not a Number. The 8087 and 80287 coprocessors used the pseudo-infinity representation for infinities.
                }
                else //Bits 61-0: Non-zero
                {
                    return FP_NAN; //Signalling Not a Number, the sign bit is meaningless.
                }
            }
            else if (bit63 && bit62) //Bits 63,62: 11
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    return FP_NAN; //Floating-point Indefinite, the result of invalid calculations such as square root of a negative number, logarithm of a negative number, 0/0, infinity / infinity, infinity times 0, and others when the processor has been configured to not generate exceptions for invalid operands. The sign bit is meaningless. This is a special case of a Quiet Not a Number.
                }
                else //Bits 61-0: Non-zero
                {
                    return FP_NAN; //Quiet Not a Number, the sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number.
                }
            }
        }

        default: //Exponent: All other values
        {
            if (!bit63) //Bit 63: Zero
            {
                return FP_NAN;
            }
            else //Bit 63: One
            {
                return FP_NORMAL; //Normalized value. The value is (−1)^s * m * 2^(e−16383)
            }
        }
    }
}

[LIU Hao]

Keep an eye on https://gcc.gnu.org/ml/gcc-help/2017-04/msg00101.html.

BTW, why the fuck do you want to fuck.lh? :<

[LIU Hao]

default: //Exponent: All other values
        {
            if (!bit63) //Bit 63: Zero
            {
                return FP_NAN;
            }

^^^ This is the number in question and it should be FP_UNSUPPORTED in my opinion.

[Duncan Ogilvie]

Sorry about the swearing I forgot to clean that up lol. Also I agree, just missed that one. Does this look correct? I'm not particularly sure about the bit manipulation stuff...

int myclassify(long double ld)
{
    __mingw_ldbl_type_t lds = *(__mingw_ldbl_type_t*)&ld;
    //https://en.wikipedia.org/wiki/Extended_precision#x86_extended_precision_format
    auto bit63 = (lds.lh.high & 0x80000000) != 0;
    switch (lds.lh.sign_exponent & 0x7fff) //Exponent
    {
        case 0: //Exponent: All Zeros
        {
            if (!bit63) //Bit 63: Zero
            {
                if (lds.lh.low == 0 && (lds.lh.high & 0x7fffffff) == 0) //Bits 62-0: Zero
                {
                    return FP_ZERO; //Zero. The sign bit gives the sign of the zero.
                }
                else //Bits 62-0: Non-zero
                {
                    return FP_SUBNORMAL; //Denormal. The value is (−1)^s * m * 2^(−16382)
                }
            }
            else //Bit 63: One
            {
                return FP_SUBNORMAL; //Pseudo Denormal. The 80387 and later properly interpret this value but will not generate it. The value is (−1)^s * m * 2^(−16382)
            }
        }

        case 0x7fff: //Exponent: All Ones
        {
            auto bit62 = (lds.lh.high & 0x40000000) != 0;
            auto bit61to0_zero = (lds.lh.high & 0x3fffffff) == 0 && lds.lh.low == 0;
            if (!bit63 && !bit62) //Bits 63,62: 00
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    return FP_UNSUPPORTED; //Pseudo-Infinity. The sign bit gives the sign of the infinity. The 8087 and 80287 treat this as Infinity. The 80387 and later treat this as an invalid operand.
                }
                else //Bits 61-0: Non-zero
                {
                    return FP_UNSUPPORTED; //Pseudo Not a Number. The sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number. The 80387 and later treat this as an invalid operand.
                }
            }
            else if (!bit63 && bit62) //Bits 63,62: 01
            {
                //Bits 61-0: Anything
                return FP_UNSUPPORTED; //Pseudo Not a Number. The sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number. The 80387 and later treat this as an invalid operand.
            }
            else if (bit63 && !bit62) //Bits 63,62: 10
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    return FP_INFINITE; //Infinity. The sign bit gives the sign of the infinity. The 8087 and 80287 treat this as a Signaling Not a Number. The 8087 and 80287 coprocessors used the pseudo-infinity representation for infinities.
                }
                else //Bits 61-0: Non-zero
                {
                    return FP_NAN; //Signalling Not a Number, the sign bit is meaningless.
                }
            }
            else if (bit63 && bit62) //Bits 63,62: 11
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    return FP_NAN; //Floating-point Indefinite, the result of invalid calculations such as square root of a negative number, logarithm of a negative number, 0/0, infinity / infinity, infinity times 0, and others when the processor has been configured to not generate exceptions for invalid operands. The sign bit is meaningless. This is a special case of a Quiet Not a Number.
                }
                else //Bits 61-0: Non-zero
                {
                    return FP_NAN; //Quiet Not a Number, the sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number.
                }
            }
        }

        default: //Exponent: All other values
        {
            if (!bit63) //Bit 63: Zero
            {
                return FP_UNSUPPORTED; //Unnormal. Only generated on the 8087 and 80287. The 80387 and later treat this as an invalid operand. The value is (−1)^s * m * 2^(e−16383)
            }
            else //Bit 63: One
            {
                return FP_NORMAL; //Normalized value. The value is (−1)^s * m * 2^(e−16383)
            }
        }
    }
}

[LIU Hao]

Just get some numbers and test it. :>

Why don't you just declare the mantissa as a uint64_t ? You could just right-shift it by 62 and switch the result, so if bit 63 and 62 are both zeroes you get 0, and if both are ones you get 3, and so on.

[Duncan Ogilvie]

This code:

#include <math.h>
int main() {
  long double ld = 0;
  unsigned char data[10] = {
    0x5b, 0x01, 0x04, 0x5e, 0x85, 0x00, 0x00, 0x00, 0xd8, 0x59
  };
  __builtin_memcpy(&ld, data, 10);
  __builtin_printf("sizeof(long double) = %d\n", sizeof(long double));
  __builtin_printf("FP_NAN: 0x%X\nFP_NORMAL: 0x%X\n", FP_NAN, FP_NORMAL);
  __builtin_printf("fpclassify: 0x%X\n", fpclassify(ld), FP_NAN, FP_NORMAL);
}

Shows (on GCC 6.3):

sizeof(long double) = 16
FP_NAN: 0x0
FP_NORMAL: 0x4
fpclassify: 0x0

And on Apple LLVM version 8.1.0 (clang-802.0.38):

sizeof(long double) = 16
FP_NAN: 0x1
FP_NORMAL: 0x4
fpclassify: 0x4

So clang has a similar issue apparently...

[Duncan Ogilvie]

I did some more testing where I take all possible sign_exponent and high values that influence the classification according to wikipedia and compare my implementation (with #define FP_UNSUPPORTED FP_NAN) with __builtin_fpclassify and it appears that my implementation is correct (or at least matches the builtin classification): https://ideone.com/RB1sZm

I will do some more testing also on the 32 bit version of GCC and mingw-w64 later to see if the 32 bit builtin/software implementation is the same. And yes this code is dumb but it works and tests all the cases so I'm fine with it.

[Duncan Ogilvie]

Okay that code is broken I will try again...

[Duncan Ogilvie]

Here is my final code, the function fxclassify is a candidate to replace ___fxclassifyl on x64 (it is correct in the sense that the currently unimplementation FP_UNSUPPORTED flag is returned as FP_NAN).

The function hwclassify is inline assembly and this works, I added support for FP_UNSUPPORTED (and if you define unsupported as nan it works).

The macro btclassify calls __builtin_fpclassify and should not be used because it is incorrect (admittedly only in certain edge cases where you have a denormal zero value).

The function myclassify is an implementation of the Wikipedia article and has been tested against all possible representations of 10-byte floating point values that can affect the status of the number. This code is implemented to be as readable as possible and to reflect the article, not to be efficient.

The testing mentioned above is implemented by iterating over all possible exponent values (sign is set to 0 since it doesn't affect the classification) and then over bits 63,62,61 which allow you to represent all cases. The bits 60-0 have to be set to 0 for this to work! (bit 61/62 fulfil the conditions zero/nonzero of the whole range). The myclassify function also sets bits 0-11 (one bit for every case) of the visited global variable which allow you to test if all cases at least happened once.

Because I don't have any experience with standard libraries I will refrain from a submitting a patch, but to fix the original issue you just have to fake fxclassify and adapt it in all __fpclassifyl implementations. For my use case I call myclassify before actually printing the number to support the FP_UNSUPPORTED so my work is done.

#include <iostream>
#include <math.h>
#include <string.h>

#define FP_UNSUPPORTED 0

/* 7.12.3.1 */
/*
   Return values for fpclassify.
   These are based on Intel x87 fpu condition codes
   in the high byte of status word and differ from
   the return values for MS IEEE 754 extension _fpclass()
*/
/*
#define FP_NAN      0x0100
#define FP_NORMAL   0x0400
#define FP_INFINITE (FP_NAN | FP_NORMAL)
#define FP_ZERO     0x4000
#define FP_SUBNORMAL    (FP_NORMAL | FP_ZERO)
// 0x0200 is signbit mask
*/

typedef union my__mingw_ldbl_type_t
{
    long double x;
    struct
    {
        unsigned int low, high;
        int sign_exponent : 16;
        int res1 : 16;
        int res0 : 32;
    } lh;
} my__mingw_ldbl_type_t;

typedef union my__mingw_fp_types_t
{
    long double *ld;
    my__mingw_ldbl_type_t *ldt;
} my__mingw_fp_types_t;

int fxclassify(long double x)
{
    my__mingw_fp_types_t hlp;
    unsigned int e, h;
    hlp.ld = &x;
    e = hlp.ldt->lh.sign_exponent & 0x7fff;
    h = hlp.ldt->lh.high;
    if (!e)
        return (!(hlp.ldt->lh.low | h)) ? FP_ZERO : FP_SUBNORMAL;
    else if (e == 0x7fff)
        return (((h & 0x7fffffff) | hlp.ldt->lh.low) == 0 ?
                ((h & 0x80000000) == 0) ? FP_NAN : FP_INFINITE : FP_NAN);
    return ((h & 0x80000000) == 0) ? FP_NAN : FP_NORMAL;
}

unsigned int visited = 0;

//Based on: https://en.wikipedia.org/wiki/Extended_precision#x86_extended_precision_format
int myclassify(long double ld)
{
    auto putlbl = [](const char* l)
    {
//        puts(l);
    };
    my__mingw_ldbl_type_t lds = *(my__mingw_ldbl_type_t*)&ld;
    bool bit63 = (lds.lh.high & 0x80000000) != 0;
    switch (lds.lh.sign_exponent & 0x7fff) //Exponent
    {
        case 0: //Exponent: All Zeros
        {
            if (!bit63) //Bit 63: Zero
            {
                if (lds.lh.low == 0 && (lds.lh.high & 0x7fffffff) == 0) //Bits 62-0: Zero
                {
                    visited |= 1 << 0;
                    putlbl("FP_ZERO (Bits 62-0: Zero)");
                    return FP_ZERO; //Zero. The sign bit gives the sign of the zero.
                }
                else //Bits 62-0: Non-zero
                {
                    visited |= 1 << 1;
                    putlbl("FP_SUBNORMAL (Bits 62-0: Non-zero)");
                    return FP_SUBNORMAL; //Denormal. The value is (−1)^s * m * 2^(−16382)
                }
            }
            else //Bit 63: One
            {
                visited |= 1 << 2;
                putlbl("FP_SUBNORMAL (Bit 63: One)");
                return FP_SUBNORMAL; //Pseudo Denormal. The 80387 and later properly interpret this value but will not generate it. The value is (−1)^s * m * 2^(−16382)
            }
        }

        case 0x7fff: //Exponent: All Ones
        {
            bool bit62 = (lds.lh.high & 0x40000000) != 0;
            bool bit61to0_zero = (lds.lh.high & 0x3fffffff) == 0 && lds.lh.low == 0;
            if (!bit63 && !bit62) //Bits 63,62: 00
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    visited |= 1 << 3;
                    putlbl("FP_UNSUPPORTED (Bits 61-0: Zero)");
                    return FP_UNSUPPORTED; //Pseudo-Infinity. The sign bit gives the sign of the infinity. The 8087 and 80287 treat this as Infinity. The 80387 and later treat this as an invalid operand.
                }
                else //Bits 61-0: Non-zero
                {
                    visited |= 1 << 4;
                    putlbl("FP_UNSUPPORTED (Bits 61-0: Non-zero)");
                    return FP_UNSUPPORTED; //Pseudo Not a Number. The sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number. The 80387 and later treat this as an invalid operand.
                }
            }
            else if (!bit63/* && bit62*/) //Bits 63,62: 01
            {
                //Bits 61-0: Anything
                visited |= 1 << 5;
                putlbl("FP_UNSUPPORTED (Bits 61-0: Anything)");
                return FP_UNSUPPORTED; //Pseudo Not a Number. The sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number. The 80387 and later treat this as an invalid operand.
            }
            else if (/*bit63 && */!bit62) //Bits 63,62: 10
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    visited |= 1 << 6;
                    putlbl("FP_INFINITE (Bits 61-0: Zero)");
                    return FP_INFINITE; //Infinity. The sign bit gives the sign of the infinity. The 8087 and 80287 treat this as a Signaling Not a Number. The 8087 and 80287 coprocessors used the pseudo-infinity representation for infinities.
                }
                else //Bits 61-0: Non-zero
                {
                    visited |= 1 << 7;
                    putlbl("FP_NAN (Bits 61-0: Non-zero)");
                    return FP_NAN; //Signalling Not a Number, the sign bit is meaningless.
                }
            }
            else /*if (bit63 && bit62)*/ //Bits 63,62: 11
            {
                if (bit61to0_zero) //Bits 61-0: Zero
                {
                    visited |= 1 << 8;
                    putlbl("FP_NAN (Bits 61-0: Zero)");
                    return FP_NAN; //Floating-point Indefinite, the result of invalid calculations such as square root of a negative number, logarithm of a negative number, 0/0, infinity / infinity, infinity times 0, and others when the processor has been configured to not generate exceptions for invalid operands. The sign bit is meaningless. This is a special case of a Quiet Not a Number.
                }
                else //Bits 61-0: Non-zero
                {
                    visited |= 1 << 9;
                    putlbl("FP_NAN (Bits 61-0: Non-zero)");
                    return FP_NAN; //Quiet Not a Number, the sign bit is meaningless. The 8087 and 80287 treat this as a Signaling Not a Number.
                }
            }
        }

        default: //Exponent: All other values
        {
            if (!bit63) //Bit 63: Zero
            {
                visited |= 1 << 10;
                putlbl("FP_UNSUPPORTED (Bit 63: Zero)");
                return FP_UNSUPPORTED; //Unnormal. Only generated on the 8087 and 80287. The 80387 and later treat this as an invalid operand. The value is (−1)^s * m * 2^(e−16383)
            }
            else //Bit 63: One
            {
                visited |= 1 << 11;
                putlbl("FP_NORMAL (Bit 63: One)");
                return FP_NORMAL; //Normalized value. The value is (−1)^s * m * 2^(e−16383)
            }
        }
    }
}

int hwclassify(long double x)
{
    unsigned short sw;
    __asm__ __volatile__ ("fxam; fstsw %%ax;" : "=a" (sw): "t" (x));
    auto status = sw & (FP_NAN | FP_NORMAL | FP_ZERO);
    return status == 0 ? FP_UNSUPPORTED : status;
}

#define btclassify(x) __builtin_fpclassify(FP_NAN, FP_INFINITE, FP_NORMAL, FP_SUBNORMAL, FP_ZERO, x)

void printhex(const void* p, size_t s)
{
    printf("{ ");
    for(size_t i = 0; i < s; i++)
    {
        if(i)
            printf(", ");
        printf("0x%02X", ((unsigned char*)p)[i]);
    }
    puts(" }");
}

int main()
{
    long double ld;
    unsigned char data[10] = {0x5b, 0x01, 0x04, 0x5e, 0x85, 0x00, 0x00, 0x00, 0xd8, 0x59};
    memcpy(&ld, data, 10);
    ld=INFINITY;
    my__mingw_ldbl_type_t lds = *(my__mingw_ldbl_type_t*)&ld;

    unsigned int resultMap[6];
    memset(resultMap, 0, sizeof(resultMap));
    char unknown[32]="???";
    const char* fpnames[] = {"FP_INFINITE", "FP_NAN", "FP_NORMAL", "FP_SUBNORMAL", "FP_ZERO", unknown};

    auto fpindex = [&unknown](int t)
    {
        switch(t)
        {
            case FP_INFINITE: return 0;
            case FP_NAN: return 1;
            case FP_NORMAL: return 2;
            case FP_SUBNORMAL: return 3;
            case FP_ZERO: return 4;
            default: sprintf(unknown, "FP_UNSUPPORTED (0x%X)", t); return 5;
        }
    };
    auto fpname = [&fpindex, &fpnames](int t) -> const char*
    {
        return fpnames[fpindex(t)];
    };

    //printf("lh.low: 0x%08X\n", lds.lh.low);
    //printf("lh.high: 0x%08X\n", lds.lh.high);
    //printf("lh.sign_exponent: 0x%04X\n", lds.lh.sign_exponent);
    //printf("lh.res0: 0x%04X\n", lds.lh.res0);
    //printf("lh.res1: 0x%08X\n", lds.lh.res1);
    //printf("ld: %s == %s\n", fpname(btclassify(ld)), fpname(myclassify(ld)));

    puts("bt: __builtin_fpclassify\nfx: __fpclassifyl (my fixed version)\nmy: my own implementation\nhw: inline asm\n");

    memset(&lds, 0, sizeof(lds));

    for(unsigned int exponent = 0; exponent < 0x8000; exponent++)
    {
        for(unsigned int high = 0; high < 8; high++)
        {
            int builtin, my;
            lds.lh.sign_exponent = exponent;
            lds.lh.high = high << 29;
            //builtin = btclassify(lds.x);
            //my = fxclassify(lds.x);
            //my = btclassify(lds.x);
            builtin = hwclassify(lds.x);
            my = myclassify(lds.x);
            if(builtin != my)
            {
                printf("NAY (exponent: 0x%X, high: 0x%X, bt: %s, fx: %s, my: %s, hw: %s)\n",
                       (unsigned short)lds.lh.sign_exponent,
                       lds.lh.high,
                       fpname(btclassify(lds.x)),
                       fpname(fxclassify(lds.x)),
                       fpname(myclassify(lds.x)),
                       fpname(hwclassify(lds.x)));
                printf("    HEX: "); printhex(&lds.x, 10);
                //return 1;
            }
            else
                resultMap[fpindex(builtin)]++;
        }
    }

    for(int i = 0; i < sizeof(resultMap) / sizeof(*resultMap); i++)
        printf("%s: %u\n", fpnames[i], resultMap[i]);

    printf("YAY (visited: 0x%X == 0xFFF, %s)\n", visited, visited == 0xFFF ? "YAY" : "NAY");

    return 0;
}