test-x87-double-rounding.c

#include <stdint.h>
#include <float.h>
#include <stdlib.h>
#include <stdio.h>
#include <fenv.h>
#pragma STDC FENV_ACCESS ON

#if defined(__x86_64__) && defined(__GNUC__)
#define HAS_X87
static void set_x87_prec_24(void)
{
    uint16_t cword;
    asm volatile("fstcw %0" : "=m"(cword));
    uint16_t newcword = (cword & ~(3u << 8)) | (0u << 8); // precision: single
    asm volatile("fldcw %0" : : "m"(newcword));
}
static void set_x87_prec_53(void)
{
    uint16_t cword;
    asm volatile("fstcw %0" : "=m"(cword));
    uint16_t newcword = (cword & ~(3u << 8)) | (2u << 8); // precision: double
    asm volatile("fldcw %0" : : "m"(newcword));
}
static void set_x87_prec_64(void)
{
    uint16_t cword;
    asm volatile("fstcw %0" : "=m"(cword));
    uint16_t newcword = (cword & ~(3u << 8)) | (3u << 8); // precision: double extended
    asm volatile("fldcw %0" : : "m"(newcword));
}

static double add_x87(double x, double y)
{
    asm volatile("fadd %1, %0\n" : "+t"(x) : "f"(y));
    return x;
}
static double mul_x87(double x, double y)
{
    asm volatile("fmul %1, %0\n" : "+t"(x) : "f"(y));
    return x;
}
#elif (defined(_M_IX86) || defined(_M_AMD64)) && defined(_MSC_VER)
#define HAS_X87
static void set_x87_prec_24(void)
{
    _control87(_PC_24, _MCW_PC);
}
static void set_x87_prec_53(void)
{
    _control87(_PC_53, _MCW_PC);
}
static void set_x87_prec_64(void)
{
    _control87(_PC_64, _MCW_PC);
}
#if _M_IX86_FP == 2
#pragma message("Warning: mul_x87 will use SSE2")
#endif
static double mul_x87(double x, double y)
{
    return x * y;
}
#endif

int main(void)
{
    struct {
        double x, y;
        const char *exact_value_s;
    } cases[] = {
                 {0x1.fffe0effffffep-51, 0x1.0000000000001p-1000, "0x1.fffe0effffffffffe0effffffep-1051"}, // underflow
                 {0x1.fffe0effffffep0, 0x1.0000000000001p0, "0x1.fffe0effffffffffe0effffffep0"}, // without underflow
    };
    for (int i = 0; i < sizeof(cases)/sizeof(cases[0]); ++i) {
        printf("Case #%d:\n", i+1);
        volatile double x = cases[i].x;
        volatile double y = cases[i].y;
        double z = x * y;
        const char *exact_value_s = cases[i].exact_value_s;
        printf("x = %a\n", x);
        printf("y = %a\n", y);
        printf("x * y = %a\n", z);
        printf("Correctly-rounded value: %a\n", strtod(exact_value_s, NULL));
#if defined(HAS_X87)
        fenv_t fenv;
        fegetenv(&fenv);
        set_x87_prec_64();
        printf("[x87, 64 bits] x * y = %a\n", mul_x87(x, y));
        set_x87_prec_53();
        printf("[x87, 53 bits] x * y = %a\n", mul_x87(x, y));
        set_x87_prec_24();
        printf("[x87, 24 bits] x * y = %a\n", mul_x87(x, y));
        fesetenv(&fenv);
#endif
        puts("---");
    }
}