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
compare(ufloat64 test,ufloat64 real,ufloat64 soft,bool exact)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
main(int ac,char ** av)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