float-commutativity.c

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

#if defined(__GNUC__)
#define NOINLINE __attribute__((noinline))
#else
#define NOINLINE
#endif

double u64_to_double(uint64_t u64)
{
    union {
        uint64_t u64;
        double dbl;
    } u;
    u.u64 = u64;
    return u.dbl;
}
uint64_t double_to_u64(double x)
{
    union {
        uint64_t u64;
        double dbl;
    } u;
    u.dbl = x;
    return u.u64;
}

NOINLINE
double add(double x, double y)
{
    return x + y;
}

NOINLINE
double mul(double x, double y)
{
    return x * y;
}

NOINLINE
void add2(double x, double y, double * restrict r1, double * restrict r2)
{
    *r1 = x + y;
    *r2 = y + x;
}

#if defined(__GNUC__) && defined(__x86_64__)
NOINLINE
double add_x87(double x, double y)
{
    asm volatile("fadd %1, %0\n"
                 : "+t"(x) : "f"(y));
    return x;
}

NOINLINE
double mul_x87(double x, double y)
{
    asm volatile("fmul %1, %0\n"
                 : "+t"(x) : "f"(y));
    return x;
}
#endif

int main(void)
{
    volatile double x = u64_to_double(UINT64_C(0x7ff800000000cafe));
    volatile double y = u64_to_double(UINT64_C(0xfff800000000beef));
    volatile double snan = u64_to_double(UINT64_C(0x7ff000000000deed));
    volatile double snan2 = u64_to_double(UINT64_C(0xfff000000000cafe));
    double r1, r2;
    add2(x, y, &r1, &r2);
    printf("          x = %" PRIx64 "\n", double_to_u64(x));
    printf("          y = %" PRIx64 "\n", double_to_u64(y));
    printf("      x + y = %" PRIx64 "\n", double_to_u64(add(x, y)));
    printf("      y + x = %" PRIx64 "\n", double_to_u64(add(y, x)));
    printf("(x + y, y + x) = (%" PRIx64 ", %" PRIx64 ")\n", double_to_u64(r1), double_to_u64(r2));
    printf("     x + -x = %" PRIx64 "\n", double_to_u64(add(x, -x)));
    printf("     -x + x = %" PRIx64 "\n", double_to_u64(add(-x, x)));
    printf("    -x + -x = %" PRIx64 "\n", double_to_u64(add(-x, -x)));
    printf(" -x + snan2 = %" PRIx64 "\n", double_to_u64(add(-x, snan2)));
    printf("  x + snan2 = %" PRIx64 "\n", double_to_u64(add(x, snan2)));
    printf("         -x = %" PRIx64 "\n", double_to_u64(- x));
    printf("   -1.0 * x = %" PRIx64 "\n", double_to_u64(mul(-1.0, x)));
    printf("  -1.0 * -x = %" PRIx64 "\n", double_to_u64(mul(-1.0, -x)));
    printf("   1.0 * -x = %" PRIx64 "\n", double_to_u64(mul(1.0, -x)));
    printf("      -snan = %" PRIx64 "\n", double_to_u64(- snan));
    printf("-1.0 * snan = %" PRIx64 "\n", double_to_u64(mul(-1.0, snan)));
    printf("      x * y = %" PRIx64 "\n", double_to_u64(mul(x, y)));
    printf("      y * x = %" PRIx64 "\n", double_to_u64(mul(y, x)));
    printf("     x * -x = %" PRIx64 "\n", double_to_u64(mul(x, -x)));
    printf("     -x * x = %" PRIx64 "\n", double_to_u64(mul(-x, x)));
    printf("    -x * -x = %" PRIx64 "\n", double_to_u64(mul(-x, -x)));
    printf(" -x * snan2 = %" PRIx64 "\n", double_to_u64(mul(-x, snan2)));
    printf("  x * snan2 = %" PRIx64 "\n", double_to_u64(mul(x, snan2)));
#if defined(__aarch64__)
    {
      unsigned int fpcr = __builtin_aarch64_get_fpcr();
      __builtin_aarch64_set_fpcr(fpcr | (1u << 25));
      fpcr = __builtin_aarch64_get_fpcr();
      if (fpcr & (1u << 25)) {
	puts("Default NaN mode set");
      } else {
	puts("Default NaN mode not set");
      }
    }
    printf("      x + y = %" PRIx64 "\n", double_to_u64(add(x, y)));
    printf("      y + x = %" PRIx64 "\n", double_to_u64(add(y, x)));
    printf("     x + -x = %" PRIx64 "\n", double_to_u64(add(x, -x)));
    printf("     -x + x = %" PRIx64 "\n", double_to_u64(add(-x, x)));
    printf("    -x + -x = %" PRIx64 "\n", double_to_u64(add(-x, -x)));
    printf(" -x + snan2 = %" PRIx64 "\n", double_to_u64(add(-x, snan2)));
    printf("  x + snan2 = %" PRIx64 "\n", double_to_u64(add(x, snan2)));
    printf("         -x = %" PRIx64 "\n", double_to_u64(- x));
    printf("   -1.0 * x = %" PRIx64 "\n", double_to_u64(mul(-1.0, x)));
    printf("  -1.0 * -x = %" PRIx64 "\n", double_to_u64(mul(-1.0, -x)));
    printf("   1.0 * -x = %" PRIx64 "\n", double_to_u64(mul(1.0, -x)));
    printf("      -snan = %" PRIx64 "\n", double_to_u64(- snan));
    printf("-1.0 * snan = %" PRIx64 "\n", double_to_u64(mul(-1.0, snan)));
    printf("      x * y = %" PRIx64 "\n", double_to_u64(mul(x, y)));
    printf("      y * x = %" PRIx64 "\n", double_to_u64(mul(y, x)));
    printf("     x * -x = %" PRIx64 "\n", double_to_u64(mul(x, -x)));
    printf("     -x * x = %" PRIx64 "\n", double_to_u64(mul(-x, x)));
    printf("    -x * -x = %" PRIx64 "\n", double_to_u64(mul(-x, -x)));
    printf(" -x * snan2 = %" PRIx64 "\n", double_to_u64(mul(-x, snan2)));
    printf("  x * snan2 = %" PRIx64 "\n", double_to_u64(mul(x, snan2)));
#elif defined(__GNUC__) && defined(__x86_64__)
    puts("x87 FPU:");
    printf("      x + y = %" PRIx64 "\n", double_to_u64(add_x87(x, y)));
    printf("      y + x = %" PRIx64 "\n", double_to_u64(add_x87(y, x)));
    printf("     x + -x = %" PRIx64 "\n", double_to_u64(add_x87(x, -x)));
    printf("     -x + x = %" PRIx64 "\n", double_to_u64(add_x87(-x, x)));
    printf("    -x + -x = %" PRIx64 "\n", double_to_u64(add_x87(-x, -x)));
    printf(" -x + snan2 = %" PRIx64 "\n", double_to_u64(add_x87(-x, snan2)));
    printf("  x + snan2 = %" PRIx64 "\n", double_to_u64(add_x87(x, snan2)));
    printf("   -1.0 * x = %" PRIx64 "\n", double_to_u64(mul_x87(-1.0, x)));
    printf("  -1.0 * -x = %" PRIx64 "\n", double_to_u64(mul_x87(-1.0, -x)));
    printf("   1.0 * -x = %" PRIx64 "\n", double_to_u64(mul_x87(1.0, -x)));
    printf("-1.0 * snan = %" PRIx64 "\n", double_to_u64(mul_x87(-1.0, snan)));
    printf("      x * y = %" PRIx64 "\n", double_to_u64(mul_x87(x, y)));
    printf("      y * x = %" PRIx64 "\n", double_to_u64(mul_x87(y, x)));
    printf("     x * -x = %" PRIx64 "\n", double_to_u64(mul_x87(x, -x)));
    printf("     -x * x = %" PRIx64 "\n", double_to_u64(mul_x87(-x, x)));
    printf("    -x * -x = %" PRIx64 "\n", double_to_u64(mul_x87(-x, -x)));
    printf(" -x * snan2 = %" PRIx64 "\n", double_to_u64(mul_x87(-x, snan2)));
    printf("  x * snan2 = %" PRIx64 "\n", double_to_u64(mul_x87(x, snan2)));
#endif
}