1*2fa3546cSRichard Henderson /* 2*2fa3546cSRichard Henderson * fp-test-log2.c - test QEMU's softfloat log2 3*2fa3546cSRichard Henderson * 4*2fa3546cSRichard Henderson * Copyright (C) 2020, Linaro, Ltd. 5*2fa3546cSRichard Henderson * 6*2fa3546cSRichard Henderson * License: GNU GPL, version 2 or later. 7*2fa3546cSRichard Henderson * See the COPYING file in the top-level directory. 8*2fa3546cSRichard Henderson */ 9*2fa3546cSRichard Henderson #ifndef HW_POISON_H 10*2fa3546cSRichard Henderson #error Must define HW_POISON_H to work around TARGET_* poisoning 11*2fa3546cSRichard Henderson #endif 12*2fa3546cSRichard Henderson 13*2fa3546cSRichard Henderson #include "qemu/osdep.h" 14*2fa3546cSRichard Henderson #include "qemu/cutils.h" 15*2fa3546cSRichard Henderson #include <math.h> 16*2fa3546cSRichard Henderson #include "fpu/softfloat.h" 17*2fa3546cSRichard Henderson 18*2fa3546cSRichard Henderson typedef union { 19*2fa3546cSRichard Henderson double d; 20*2fa3546cSRichard Henderson float64 i; 21*2fa3546cSRichard Henderson } ufloat64; 22*2fa3546cSRichard Henderson 23*2fa3546cSRichard Henderson static int errors; 24*2fa3546cSRichard Henderson 25*2fa3546cSRichard Henderson static void compare(ufloat64 test, ufloat64 real, ufloat64 soft, bool exact) 26*2fa3546cSRichard Henderson { 27*2fa3546cSRichard Henderson int msb; 28*2fa3546cSRichard Henderson uint64_t ulp = UINT64_MAX; 29*2fa3546cSRichard Henderson 30*2fa3546cSRichard Henderson if (real.i == soft.i) { 31*2fa3546cSRichard Henderson return; 32*2fa3546cSRichard Henderson } 33*2fa3546cSRichard Henderson msb = 63 - __builtin_clzll(real.i ^ soft.i); 34*2fa3546cSRichard Henderson 35*2fa3546cSRichard Henderson if (msb < 52) { 36*2fa3546cSRichard Henderson if (real.i > soft.i) { 37*2fa3546cSRichard Henderson ulp = real.i - soft.i; 38*2fa3546cSRichard Henderson } else { 39*2fa3546cSRichard Henderson ulp = soft.i - real.i; 40*2fa3546cSRichard Henderson } 41*2fa3546cSRichard Henderson } 42*2fa3546cSRichard Henderson 43*2fa3546cSRichard Henderson /* glibc allows 3 ulp error in its libm-test-ulps; allow 4 here */ 44*2fa3546cSRichard Henderson if (!exact && ulp <= 4) { 45*2fa3546cSRichard Henderson return; 46*2fa3546cSRichard Henderson } 47*2fa3546cSRichard Henderson 48*2fa3546cSRichard Henderson printf("test: %016" PRIx64 " %+.13a\n" 49*2fa3546cSRichard Henderson " sf: %016" PRIx64 " %+.13a\n" 50*2fa3546cSRichard Henderson "libm: %016" PRIx64 " %+.13a\n", 51*2fa3546cSRichard Henderson test.i, test.d, soft.i, soft.d, real.i, real.d); 52*2fa3546cSRichard Henderson 53*2fa3546cSRichard Henderson if (msb == 63) { 54*2fa3546cSRichard Henderson printf("Error in sign!\n\n"); 55*2fa3546cSRichard Henderson } else if (msb >= 52) { 56*2fa3546cSRichard Henderson printf("Error in exponent: %d\n\n", 57*2fa3546cSRichard Henderson (int)(soft.i >> 52) - (int)(real.i >> 52)); 58*2fa3546cSRichard Henderson } else { 59*2fa3546cSRichard Henderson printf("Error in fraction: %" PRIu64 " ulp\n\n", ulp); 60*2fa3546cSRichard Henderson } 61*2fa3546cSRichard Henderson 62*2fa3546cSRichard Henderson if (++errors == 20) { 63*2fa3546cSRichard Henderson exit(1); 64*2fa3546cSRichard Henderson } 65*2fa3546cSRichard Henderson } 66*2fa3546cSRichard Henderson 67*2fa3546cSRichard Henderson int main(int ac, char **av) 68*2fa3546cSRichard Henderson { 69*2fa3546cSRichard Henderson ufloat64 test, real, soft; 70*2fa3546cSRichard Henderson float_status qsf = {0}; 71*2fa3546cSRichard Henderson int i; 72*2fa3546cSRichard Henderson 73*2fa3546cSRichard Henderson set_float_rounding_mode(float_round_nearest_even, &qsf); 74*2fa3546cSRichard Henderson 75*2fa3546cSRichard Henderson test.d = 0.0; 76*2fa3546cSRichard Henderson real.d = -__builtin_inf(); 77*2fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf); 78*2fa3546cSRichard Henderson compare(test, real, soft, true); 79*2fa3546cSRichard Henderson 80*2fa3546cSRichard Henderson test.d = 1.0; 81*2fa3546cSRichard Henderson real.d = 0.0; 82*2fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf); 83*2fa3546cSRichard Henderson compare(test, real, soft, true); 84*2fa3546cSRichard Henderson 85*2fa3546cSRichard Henderson test.d = 2.0; 86*2fa3546cSRichard Henderson real.d = 1.0; 87*2fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf); 88*2fa3546cSRichard Henderson compare(test, real, soft, true); 89*2fa3546cSRichard Henderson 90*2fa3546cSRichard Henderson test.d = 4.0; 91*2fa3546cSRichard Henderson real.d = 2.0; 92*2fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf); 93*2fa3546cSRichard Henderson compare(test, real, soft, true); 94*2fa3546cSRichard Henderson 95*2fa3546cSRichard Henderson test.d = 0x1p64; 96*2fa3546cSRichard Henderson real.d = 64.0; 97*2fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf); 98*2fa3546cSRichard Henderson compare(test, real, soft, true); 99*2fa3546cSRichard Henderson 100*2fa3546cSRichard Henderson test.d = __builtin_inf(); 101*2fa3546cSRichard Henderson real.d = __builtin_inf(); 102*2fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf); 103*2fa3546cSRichard Henderson compare(test, real, soft, true); 104*2fa3546cSRichard Henderson 105*2fa3546cSRichard Henderson for (i = 0; i < 10000; ++i) { 106*2fa3546cSRichard Henderson test.d = drand48() + 1.0; /* [1.0, 2.0) */ 107*2fa3546cSRichard Henderson real.d = log2(test.d); 108*2fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf); 109*2fa3546cSRichard Henderson compare(test, real, soft, false); 110*2fa3546cSRichard Henderson 111*2fa3546cSRichard Henderson test.d = drand48() * 100; /* [0.0, 100) */ 112*2fa3546cSRichard Henderson real.d = log2(test.d); 113*2fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf); 114*2fa3546cSRichard Henderson compare(test, real, soft, false); 115*2fa3546cSRichard Henderson } 116*2fa3546cSRichard Henderson 117*2fa3546cSRichard Henderson return 0; 118*2fa3546cSRichard Henderson } 119