/* == on floating-point numbers compares the values actually stored, and 0.1,
   0.2 and 0.3 are not stored exactly. In float the rounding happens to land
   0.1f + 0.2f on the same bits as 0.3f; in double it lands one step away.
   memcpy into a same-sized unsigned integer shows the bits without breaking
   the aliasing rules. */
#include <float.h>
#include <inttypes.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>
#include <string.h>

_Static_assert(sizeof(float) == sizeof(uint32_t), "float is 32 bits");
_Static_assert(sizeof(double) == sizeof(uint64_t), "double is 64 bits");

static uint32_t bits32(float f)
{
    uint32_t u;
    memcpy(&u, &f, sizeof u);
    return u;
}

static uint64_t bits64(double d)
{
    uint64_t u;
    memcpy(&u, &d, sizeof u);
    return u;
}

int main(void)
{
    float sum_f = 0.1f + 0.2f;
    printf("0.1f + 0.2f == 0.3f  -> %d\n", sum_f == 0.3f);
    printf("  0.1f + 0.2f bits   =  %08" PRIx32 "\n", bits32(sum_f));
    printf("  0.3f bits          =  %08" PRIx32 "\n", bits32(0.3f));

    double sum_d = 0.1 + 0.2;
    printf("0.1 + 0.2 == 0.3     -> %d\n", sum_d == 0.3);
    printf("  0.1 + 0.2 bits     =  %016" PRIx64 "\n", bits64(sum_d));
    printf("  0.3 bits           =  %016" PRIx64 "\n", bits64(0.3));
    printf("  0.1 + 0.2          =  %.17g\n", sum_d);

    /* An epsilon test. DBL_EPSILON is the gap between 1.0 and the next double,
       so this is a tolerance sized for numbers near 1. */
    printf("fabs(0.1 + 0.2 - 0.3) = %.17g\n", fabs(sum_d - 0.3));
    printf("DBL_EPSILON           = %.17g\n", DBL_EPSILON);
    printf("FLT_EPSILON           = %.9g\n", FLT_EPSILON);
    printf("fabs(0.1 + 0.2 - 0.3) <= DBL_EPSILON -> %d\n", fabs(sum_d - 0.3) <= DBL_EPSILON);

    /* NaN is the one value not equal to itself, so ask isnan instead. */
    double nan = NAN, inf = INFINITY;
    printf("NAN == NAN           -> %d\n", nan == nan);
    printf("NAN != NAN           -> %d\n", nan != nan);
    printf("isnan(NAN)           -> %d\n", isnan(nan) != 0);
    printf("isfinite(INFINITY)   -> %d\n", isfinite(inf) != 0);
    printf("isfinite(DBL_MAX)    -> %d\n", isfinite(DBL_MAX) != 0);
    printf("isnan(INFINITY - INFINITY) -> %d\n", isnan(inf - inf) != 0);
    return 0;
}