12fa3546cSRichard Henderson /*
22fa3546cSRichard Henderson * fp-test-log2.c - test QEMU's softfloat log2
32fa3546cSRichard Henderson *
42fa3546cSRichard Henderson * Copyright (C) 2020, Linaro, Ltd.
52fa3546cSRichard Henderson *
62fa3546cSRichard Henderson * License: GNU GPL, version 2 or later.
72fa3546cSRichard Henderson * See the COPYING file in the top-level directory.
82fa3546cSRichard Henderson */
92fa3546cSRichard Henderson #ifndef HW_POISON_H
102fa3546cSRichard Henderson #error Must define HW_POISON_H to work around TARGET_* poisoning
112fa3546cSRichard Henderson #endif
122fa3546cSRichard Henderson
132fa3546cSRichard Henderson #include "qemu/osdep.h"
142fa3546cSRichard Henderson #include "qemu/cutils.h"
152fa3546cSRichard Henderson #include <math.h>
162fa3546cSRichard Henderson #include "fpu/softfloat.h"
172fa3546cSRichard Henderson
182fa3546cSRichard Henderson typedef union {
192fa3546cSRichard Henderson double d;
202fa3546cSRichard Henderson float64 i;
212fa3546cSRichard Henderson } ufloat64;
222fa3546cSRichard Henderson
232fa3546cSRichard Henderson static int errors;
242fa3546cSRichard Henderson
compare(ufloat64 test,ufloat64 real,ufloat64 soft,bool exact)252fa3546cSRichard Henderson static void compare(ufloat64 test, ufloat64 real, ufloat64 soft, bool exact)
262fa3546cSRichard Henderson {
272fa3546cSRichard Henderson int msb;
282fa3546cSRichard Henderson uint64_t ulp = UINT64_MAX;
292fa3546cSRichard Henderson
302fa3546cSRichard Henderson if (real.i == soft.i) {
312fa3546cSRichard Henderson return;
322fa3546cSRichard Henderson }
332fa3546cSRichard Henderson msb = 63 - __builtin_clzll(real.i ^ soft.i);
342fa3546cSRichard Henderson
352fa3546cSRichard Henderson if (msb < 52) {
362fa3546cSRichard Henderson if (real.i > soft.i) {
372fa3546cSRichard Henderson ulp = real.i - soft.i;
382fa3546cSRichard Henderson } else {
392fa3546cSRichard Henderson ulp = soft.i - real.i;
402fa3546cSRichard Henderson }
412fa3546cSRichard Henderson }
422fa3546cSRichard Henderson
432fa3546cSRichard Henderson /* glibc allows 3 ulp error in its libm-test-ulps; allow 4 here */
442fa3546cSRichard Henderson if (!exact && ulp <= 4) {
452fa3546cSRichard Henderson return;
462fa3546cSRichard Henderson }
472fa3546cSRichard Henderson
482fa3546cSRichard Henderson printf("test: %016" PRIx64 " %+.13a\n"
492fa3546cSRichard Henderson " sf: %016" PRIx64 " %+.13a\n"
502fa3546cSRichard Henderson "libm: %016" PRIx64 " %+.13a\n",
512fa3546cSRichard Henderson test.i, test.d, soft.i, soft.d, real.i, real.d);
522fa3546cSRichard Henderson
532fa3546cSRichard Henderson if (msb == 63) {
542fa3546cSRichard Henderson printf("Error in sign!\n\n");
552fa3546cSRichard Henderson } else if (msb >= 52) {
562fa3546cSRichard Henderson printf("Error in exponent: %d\n\n",
572fa3546cSRichard Henderson (int)(soft.i >> 52) - (int)(real.i >> 52));
582fa3546cSRichard Henderson } else {
592fa3546cSRichard Henderson printf("Error in fraction: %" PRIu64 " ulp\n\n", ulp);
602fa3546cSRichard Henderson }
612fa3546cSRichard Henderson
622fa3546cSRichard Henderson if (++errors == 20) {
632fa3546cSRichard Henderson exit(1);
642fa3546cSRichard Henderson }
652fa3546cSRichard Henderson }
662fa3546cSRichard Henderson
main(int ac,char ** av)672fa3546cSRichard Henderson int main(int ac, char **av)
682fa3546cSRichard Henderson {
692fa3546cSRichard Henderson ufloat64 test, real, soft;
702fa3546cSRichard Henderson float_status qsf = {0};
712fa3546cSRichard Henderson int i;
722fa3546cSRichard Henderson
73*d22c9949SPeter Maydell set_float_2nan_prop_rule(float_2nan_prop_s_ab, &qsf);
742fa3546cSRichard Henderson set_float_rounding_mode(float_round_nearest_even, &qsf);
752fa3546cSRichard Henderson
762fa3546cSRichard Henderson test.d = 0.0;
772fa3546cSRichard Henderson real.d = -__builtin_inf();
782fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf);
792fa3546cSRichard Henderson compare(test, real, soft, true);
802fa3546cSRichard Henderson
812fa3546cSRichard Henderson test.d = 1.0;
822fa3546cSRichard Henderson real.d = 0.0;
832fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf);
842fa3546cSRichard Henderson compare(test, real, soft, true);
852fa3546cSRichard Henderson
862fa3546cSRichard Henderson test.d = 2.0;
872fa3546cSRichard Henderson real.d = 1.0;
882fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf);
892fa3546cSRichard Henderson compare(test, real, soft, true);
902fa3546cSRichard Henderson
912fa3546cSRichard Henderson test.d = 4.0;
922fa3546cSRichard Henderson real.d = 2.0;
932fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf);
942fa3546cSRichard Henderson compare(test, real, soft, true);
952fa3546cSRichard Henderson
962fa3546cSRichard Henderson test.d = 0x1p64;
972fa3546cSRichard Henderson real.d = 64.0;
982fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf);
992fa3546cSRichard Henderson compare(test, real, soft, true);
1002fa3546cSRichard Henderson
1012fa3546cSRichard Henderson test.d = __builtin_inf();
1022fa3546cSRichard Henderson real.d = __builtin_inf();
1032fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf);
1042fa3546cSRichard Henderson compare(test, real, soft, true);
1052fa3546cSRichard Henderson
1062fa3546cSRichard Henderson for (i = 0; i < 10000; ++i) {
1072fa3546cSRichard Henderson test.d = drand48() + 1.0; /* [1.0, 2.0) */
1082fa3546cSRichard Henderson real.d = log2(test.d);
1092fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf);
1102fa3546cSRichard Henderson compare(test, real, soft, false);
1112fa3546cSRichard Henderson
1122fa3546cSRichard Henderson test.d = drand48() * 100; /* [0.0, 100) */
1132fa3546cSRichard Henderson real.d = log2(test.d);
1142fa3546cSRichard Henderson soft.i = float64_log2(test.i, &qsf);
1152fa3546cSRichard Henderson compare(test, real, soft, false);
1162fa3546cSRichard Henderson }
1172fa3546cSRichard Henderson
1182fa3546cSRichard Henderson return 0;
1192fa3546cSRichard Henderson }
120