libmのsinの実装が x=2^n （-1000≤n≤1000）の形の入力に対して「正しく丸められた値」からどのぐらいずれているか検査するやつ

sin-test.c

// Written by @mod_poppo
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <mpfr.h>
#include <gmp.h>

static void to_normalized_hex(char *buf, mpfr_t a, int mant)
{
    // assume a is positive
    mpz_t z;
    mpz_init(z);
    mpfr_exp_t e = mpfr_get_z_2exp(z, a);
    size_t size_in_base_2 = mpz_sizeinbase(z, 2);
    // 2^(mant+k-1) <= 1XXXXXX(in hex) < 2^(mant+k), 0 <= k < 4
    // (mant+k-1) is a multiple of 4
    // k === 1-mant mod 4
    int k = (1 + 3 * mant) % 4;
    mpz_mul_2exp(z, z, mant + k - size_in_base_2);
    char buf2[(mant + 2) / 4 + 3];
    mpz_get_str(buf2, 16, z);
    // buf2: z * 2^(mant + k - size_in_base_2)
    // z * 2^(mant + k - size_in_base_2) * 2^(mant-1-k)
    sprintf(buf, "0x%c.%sp%+d", buf2[0] /* 0 */, buf2+1, (int)e + mant - 1);
    mpz_clear(z);
}

static void to_hex_str(char *buf, mpfr_t x)
{
    int prec = mpfr_get_prec(x);
    mpfr_t w;
    mpfr_init2(w, prec);
    int sign = mpfr_signbit(x);
    int t = mpfr_abs(w, x, MPFR_RNDN);
    if (t != 0) {
        fputs("Unexpected rounding in mpfr_abs\n", stderr);
        abort();
    }
    if (sign == 0) {
        to_normalized_hex(buf, w, prec);
    } else {
        buf[0] = '-';
        to_normalized_hex(buf + 1, w, prec);
    }
    mpfr_clear(w);
    // mpfr_printf("%s%Ra%s", s1, x, s2);
    // printf("%s%s%s", s1, buf, s2);
}

int main(void)
{
    for (int n = -1000; n <= 1000; n++) {
        double x = ldexp(2.0, n);
        double y = sin(x);
        double yc;
        int t;
        mpfr_t xm, ym, yp;
        mpfr_init2(xm, 53);
        mpfr_init2(ym, 53);
        mpfr_init2(yp, 113);
        int r = mpfr_set_d(xm, x, MPFR_RNDN);
        if (r != 0) {
            fprintf(stderr, "mpfr_set_d was inexact on %a!\n", x);
            exit(1);
        }
        t = mpfr_sin(ym, xm, MPFR_RNDN);
        mpfr_sin(yp, xm, MPFR_RNDN);
        yc = mpfr_get_d(ym, MPFR_RNDN);
        if (y != yc) {
            double yy = t > 0 ? nextafter(yc, -INFINITY) : t < 0 ? nextafter(yc, INFINITY) : yc;
            printf("MPFR: sin(%a) = %a\n", x, yc);
            char buf[1024];
            to_hex_str(buf, yp);
            printf("MPFR: sin(%a) = %s\n", x, buf);
            if (yy != y) {
                printf("libm: sin(%a) = %a <more than 1 ulp>\n", x, y);
            } else {
                printf("libm: sin(%a) = %a\n", x, y);
            }
        }
        mpfr_clear(xm);
        mpfr_clear(ym);
        mpfr_clear(yp);
    }
}