2 * Copyright (c) 1998 Doug Rabson
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions
8 * 1. Redistributions of source code must retain the above copyright
9 * notice, this list of conditions and the following disclaimer.
10 * 2. Redistributions in binary form must reproduce the above copyright
11 * notice, this list of conditions and the following disclaimer in the
12 * documentation and/or other materials provided with the distribution.
14 * THIS SOFTWARE IS PROVIDED BY THE AUTHOR AND CONTRIBUTORS ``AS IS'' AND
15 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
16 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
17 * ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE
18 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
19 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
20 * OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
21 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
22 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
23 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
30 * An implementation of IEEE 754 floating point arithmetic supporting
31 * multiply, divide, addition, subtraction and conversion to and from
32 * integer. Probably not the fastest floating point code in the world
33 * but it should be pretty accurate.
35 * A special thanks to John Polstra for pointing out some problems
36 * with an earlier version of this code and for educating me as to the
37 * correct use of sticky bits.
40 #include <sys/types.h>
42 #include "../include/fpu.h"
43 #include "ieee_float.h"
45 #include <sys/param.h>
46 #include <sys/systm.h>
47 #include <sys/sysproto.h>
48 #include <sys/sysent.h>
50 #include <machine/fpu.h>
51 #include <alpha/alpha/ieee_float.h>
55 * The number of fraction bits in a T format float.
60 * The number of fraction bits in a S format float.
65 * Mask the fraction part of a float to contain only those bits which
66 * should be in single precision number.
68 #define S_FRACMASK ((1ULL << 52) - (1ULL << 29))
71 * The number of extra zero bits we shift into the fraction part
72 * to gain accuracy. Two guard bits and one sticky bit are required
73 * to ensure accurate rounding.
78 * Values for 1.0 and 2.0 fractions (including the extra FRAC_SHIFT
81 #define ONE (1ULL << (T_FRACBITS + FRAC_SHIFT))
82 #define TWO (ONE + ONE)
85 * The maximum and minimum values for S and T format exponents.
87 #define T_MAXEXP 0x3ff
88 #define T_MINEXP -0x3fe
90 #define S_MINEXP -0x7e
93 * Exponent values in registers are biased by adding this value.
95 #define BIAS_EXP 0x3ff
98 * Exponent value for INF and NaN.
100 #define NAN_EXP 0x7ff
103 * If this bit is set in the fraction part of a NaN, then the number
104 * is a quiet NaN, i.e. no traps are generated.
106 #define QNAN_BIT (1ULL << 51)
109 * Return true if the number is any kind of NaN.
112 isNaN(fp_register_t f)
114 return f.t.exponent == NAN_EXP && f.t.fraction != 0;
118 * Return true if the number is a quiet NaN.
121 isQNaN(fp_register_t f)
123 return f.t.exponent == NAN_EXP && (f.t.fraction & QNAN_BIT);
127 * Return true if the number is a signalling NaN.
130 isSNaN(fp_register_t f)
132 return isNaN(f) && !isQNaN(f);
136 * Return true if the number is +/- INF.
139 isINF(fp_register_t f)
141 return f.t.exponent == NAN_EXP && f.t.fraction == 0;
145 * Return true if the number is +/- 0.
148 isZERO(fp_register_t f)
150 return f.t.exponent == 0 && f.t.fraction == 0;
154 * Return true if the number is denormalised.
157 isDENORM(fp_register_t f)
159 return f.t.exponent == 0 && f.t.fraction != 0;
163 * Extract the exponent part of a float register. If the exponent is
164 * zero, the number may be denormalised (if the fraction is nonzero).
165 * If so, return the minimum exponent for the source datatype.
168 getexp(fp_register_t f, int src)
170 int minexp[] = { S_MINEXP, 0, T_MINEXP, 0 };
171 if (f.t.exponent == 0) {
177 return f.t.exponent - BIAS_EXP;
181 * Extract the fraction part of a float register, shift it up a bit
182 * to give extra accuracy and add in the implicit 1 bit. Must be
183 * careful to handle denormalised numbers and zero correctly.
185 static __inline u_int64_t
186 getfrac(fp_register_t f)
188 if (f.t.exponent == 0)
189 return f.t.fraction << FRAC_SHIFT;
191 return (f.t.fraction << FRAC_SHIFT) | ONE;
195 * Make a float (in register format) from a sign, exponent and
196 * fraction, normalising and rounding as necessary.
197 * Return the float and set *status if any traps are generated.
200 makefloat(int sign, int exp, u_int64_t frac,
202 u_int64_t control, u_int64_t *status)
205 int minexp = 0, maxexp = 0, alpha = 0;
206 u_int64_t epsilon = 0, max = 0;
217 * Fraction is >= 2.0.
218 * Shift the fraction down, preserving the 'sticky'
221 while (frac >= TWO) {
222 frac = (frac >> 1) | (frac & 1);
225 } else if (frac < ONE) {
227 * Fraction is < 1.0. Shift it up.
230 frac = (frac << 1) | (frac & 1);
240 epsilon = (1ULL << (T_FRACBITS - S_FRACBITS + FRAC_SHIFT));
248 epsilon = (1ULL << FRAC_SHIFT);
254 * Handle underflow before rounding so that denormalised
255 * numbers are rounded correctly.
259 if (control & IEEE_TRAP_ENABLE_UNF) {
264 while (exp < minexp) {
266 frac = (frac >> 1) | (frac & 1);
273 * Round the fraction according to the rounding mode.
275 if (frac & (epsilon - 1)) {
276 u_int64_t fraclo, frachi;
277 u_int64_t difflo, diffhi;
280 frachi = fraclo + epsilon;
285 case ROUND_MINUS_INF:
292 difflo = frac - fraclo;
293 diffhi = frachi - frac;
296 else if (diffhi < difflo)
300 if (fraclo & epsilon)
314 * Rounding up may take us to TWO if
315 * fraclo == (TWO - epsilon). Also If fraclo has been
316 * denormalised to (ONE - epsilon) then there is a
317 * possibility that we will round to ONE exactly.
320 frac = (frac >> 1) & ~(epsilon - 1);
322 } else if (exp == minexp - 1 && frac == ONE) {
323 /* Renormalise to ONE * 2^minexp */
331 * Check for overflow and round to the correct INF as needed.
334 *status |= FPCR_OVF | FPCR_INE;
335 if (control & IEEE_TRAP_ENABLE_OVF) {
343 case ROUND_MINUS_INF:
345 exp = maxexp + 1; /* INF */
353 exp = maxexp + 1; /* INF */
361 exp = maxexp + 1; /* INF */
370 if (exp > maxexp) /* NaN, INF */
371 f.t.exponent = NAN_EXP;
372 else if (exp < minexp) /* denorm, zero */
375 f.t.exponent = exp + BIAS_EXP;
376 f.t.fraction = (frac & ~ONE) >> FRAC_SHIFT;
381 * Return the canonical quiet NaN in register format.
388 f.t.exponent = NAN_EXP;
389 f.t.fraction = QNAN_BIT;
401 f.t.exponent = NAN_EXP;
420 ieee_add(fp_register_t fa, fp_register_t fb,
422 u_int64_t control, u_int64_t *status)
426 u_int64_t fraca, fracb, frac;
429 /* First handle NaNs */
430 if (isNaN(fa) || isNaN(fb)) {
431 fp_register_t result;
433 /* Instructions Descriptions (I) section 4.7.10.4 */
436 else if (isSNaN(fb)) {
438 result.t.fraction |= QNAN_BIT;
439 } else if (isQNaN(fa))
441 else if (isSNaN(fa)) {
443 result.t.fraction |= QNAN_BIT;
446 /* If either operand is a signalling NaN, trap. */
447 if (isSNaN(fa) || isSNaN(fb))
456 if (fa.t.sign != fb.t.sign) {
457 /* If adding -INF to +INF, generate a trap. */
468 * Unpack the registers.
470 expa = getexp(fa, src);
471 expb = getexp(fb, src);
478 sticky = (fraca & ((1ULL << shift) - 1)) != 0;
482 fraca = (fraca >> shift) | sticky;
483 } else if (shift > 0) {
485 sticky = (fracb & ((1ULL << shift) - 1)) != 0;
489 fracb = (fracb >> shift) | sticky;
492 if (fa.t.sign) fraca = -fraca;
493 if (fb.t.sign) fracb = -fracb;
494 frac = fraca + fracb;
502 if (fa.t.exponent == 0 && fa.t.fraction == 0
503 && fb.t.exponent == 0 && fb.t.fraction == 0)
504 sign = fa.t.sign && fb.t.sign;
506 return makefloat(sign, exp, frac, src, rnd, control, status);
510 ieee_sub(fp_register_t fa, fp_register_t fb,
512 u_int64_t control, u_int64_t *status)
514 fb.t.sign = !fb.t.sign;
515 return ieee_add(fa, fb, src, rnd, control, status);
523 #define SRL128(x, b) \
526 x.lo |= x.hi << (64 - b); \
530 #define SLL128(x, b) \
533 x.hi = x.lo << (b - 64); \
537 x.hi |= x.lo >> (64 - b); \
542 #define SUB128(a, b) \
544 int borrow = a.lo < b.lo; \
545 a.lo = a.lo - b.lo; \
546 a.hi = a.hi - b.hi - borrow; \
549 #define LESS128(a, b) (a.hi < b.hi || (a.hi == b.hi && a.lo < b.lo))
552 ieee_mul(fp_register_t fa, fp_register_t fb,
554 u_int64_t control, u_int64_t *status)
557 u_int64_t fraca, fracb, tmp;
561 /* First handle NaNs */
562 if (isNaN(fa) || isNaN(fb)) {
563 fp_register_t result;
565 /* Instructions Descriptions (I) section 4.7.10.4 */
568 else if (isSNaN(fb)) {
570 result.t.fraction |= QNAN_BIT;
571 } else if (isQNaN(fa))
573 else if (isSNaN(fa)) {
575 result.t.fraction |= QNAN_BIT;
578 /* If either operand is a signalling NaN, trap. */
579 if (isSNaN(fa) || isSNaN(fb))
585 /* Handle INF and 0 */
586 if ((isINF(fa) && isZERO(fb)) || (isZERO(fa) && isINF(fb))) {
591 /* If either is INF or zero, get the sign right */
592 if (isINF(fa) || isINF(fb))
593 return makeINF(fa.t.sign ^ fb.t.sign);
594 else if (isZERO(fa) || isZERO(fb))
595 return makeZERO(fa.t.sign ^ fb.t.sign);
598 * Unpack the registers.
600 expa = getexp(fa, src);
601 expb = getexp(fb, src);
604 sign = fa.t.sign ^ fb.t.sign;
606 #define LO32(x) ((x) & ((1ULL << 32) - 1))
607 #define HI32(x) ((x) >> 32)
610 * Calculate the 128bit result of multiplying fraca and fracb.
612 frac.lo = fraca * fracb;
615 * The alpha has a handy instruction to find the high word.
617 __asm__ __volatile__ ("umulh %1,%2,%0"
619 : "r"(fraca), "r"(fracb));
623 * Do the multiply longhand otherwise.
625 frac.hi = HI32(LO32(fraca) * HI32(fracb)
626 + HI32(fraca) * LO32(fracb)
627 + HI32(LO32(fraca) * LO32(fracb)))
628 + HI32(fraca) * HI32(fracb);
630 exp = expa + expb - (T_FRACBITS + FRAC_SHIFT);
632 while (frac.hi > 0) {
635 sticky = frac.lo & 1;
640 return makefloat(sign, exp, frac.lo, src, rnd, control, status);
644 divide_128(u_int128_t a, u_int128_t b)
651 * Make a couple of assumptions on the numbers passed in. The
652 * value in 'a' will have bits set in the upper 64 bits only
653 * and the number in 'b' will have zeros in the upper 64 bits.
654 * Also, 'b' will not be zero.
657 if (a.hi == 0 || b.hi != 0 || b.lo == 0)
662 * Find out how many bits of zeros are at the beginning of the divisor.
674 * Find out how much to shift the divisor so that its msb
675 * matches the msb of the dividend.
690 * Calculate the result in two parts to avoid keeping a 128bit
691 * value for the result bit.
694 bit = 1ULL << (i - 64);
696 if (!LESS128(a, b)) {
709 if (!LESS128(a, b)) {
723 ieee_div(fp_register_t fa, fp_register_t fb,
725 u_int64_t control, u_int64_t *status)
728 u_int128_t fraca, fracb, frac;
731 /* First handle NaNs, INFs and ZEROs */
732 if (isNaN(fa) || isNaN(fb)) {
733 fp_register_t result;
735 /* Instructions Descriptions (I) section 4.7.10.4 */
738 else if (isSNaN(fb)) {
740 result.t.fraction |= QNAN_BIT;
741 } else if (isQNaN(fa))
743 else if (isSNaN(fa)) {
745 result.t.fraction |= QNAN_BIT;
748 /* If either operand is a signalling NaN, trap. */
749 if (isSNaN(fa) || isSNaN(fb))
755 /* Handle INF and 0 */
756 if (isINF(fa) && isINF(fb)) {
759 } else if (isZERO(fb))
765 return makeINF(fa.t.sign ^ fb.t.sign);
768 return makeZERO(fa.t.sign ^ fb.t.sign);
771 * Unpack the registers.
773 expa = getexp(fa, src);
774 expb = getexp(fb, src);
775 fraca.hi = getfrac(fa);
777 fracb.lo = getfrac(fb);
779 sign = fa.t.sign ^ fb.t.sign;
781 frac = divide_128(fraca, fracb);
783 exp = expa - expb - (64 - T_FRACBITS - FRAC_SHIFT);
784 while (frac.hi > 0) {
787 sticky = frac.lo & 1;
793 return makefloat(sign, exp, frac.lo, src, rnd, control, status);
796 #define IEEE_TRUE 0x4000000000000000ULL
800 ieee_cmpun(fp_register_t fa, fp_register_t fb, u_int64_t *status)
802 fp_register_t result;
803 if (isNaN(fa) || isNaN(fb)) {
804 if (isSNaN(fa) || isSNaN(fb))
806 result.q = IEEE_TRUE;
808 result.q = IEEE_FALSE;
814 ieee_cmpeq(fp_register_t fa, fp_register_t fb, u_int64_t *status)
816 fp_register_t result;
817 if (isNaN(fa) || isNaN(fb)) {
818 if (isSNaN(fa) || isSNaN(fb))
820 result.q = IEEE_FALSE;
822 if (isZERO(fa) && isZERO(fb))
823 result.q = IEEE_TRUE;
824 else if (fa.q == fb.q)
825 result.q = IEEE_TRUE;
827 result.q = IEEE_FALSE;
834 ieee_cmplt(fp_register_t fa, fp_register_t fb, u_int64_t *status)
836 fp_register_t result;
837 if (isNaN(fa) || isNaN(fb)) {
838 if (isSNaN(fa) || isSNaN(fb))
840 result.q = IEEE_FALSE;
842 if (isZERO(fa) && isZERO(fb))
843 result.q = IEEE_FALSE;
844 else if (fa.t.sign) {
847 /* fb is positive, return true */
848 result.q = IEEE_TRUE;
849 else if (fa.t.exponent > fb.t.exponent)
850 /* fa has a larger exponent, return true */
851 result.q = IEEE_TRUE;
852 else if (fa.t.exponent == fb.t.exponent
853 && fa.t.fraction > fb.t.fraction)
854 /* compare fractions */
855 result.q = IEEE_TRUE;
857 result.q = IEEE_FALSE;
861 /* fb is negative, return false */
862 result.q = IEEE_FALSE;
863 else if (fb.t.exponent > fa.t.exponent)
864 /* fb has a larger exponent, return true */
865 result.q = IEEE_TRUE;
866 else if (fa.t.exponent == fb.t.exponent
867 && fa.t.fraction < fb.t.fraction)
868 /* compare fractions */
869 result.q = IEEE_TRUE;
871 result.q = IEEE_FALSE;
879 ieee_cmple(fp_register_t fa, fp_register_t fb, u_int64_t *status)
881 fp_register_t result;
882 if (isNaN(fa) || isNaN(fb)) {
883 if (isSNaN(fa) || isSNaN(fb))
885 result.q = IEEE_FALSE;
887 if (isZERO(fa) && isZERO(fb))
888 result.q = IEEE_TRUE;
889 else if (fa.t.sign) {
892 /* fb is positive, return true */
893 result.q = IEEE_TRUE;
894 else if (fa.t.exponent > fb.t.exponent)
895 /* fa has a larger exponent, return true */
896 result.q = IEEE_TRUE;
897 else if (fa.t.exponent == fb.t.exponent
898 && fa.t.fraction >= fb.t.fraction)
899 /* compare fractions */
900 result.q = IEEE_TRUE;
902 result.q = IEEE_FALSE;
906 /* fb is negative, return false */
907 result.q = IEEE_FALSE;
908 else if (fb.t.exponent > fa.t.exponent)
909 /* fb has a larger exponent, return true */
910 result.q = IEEE_TRUE;
911 else if (fa.t.exponent == fb.t.exponent
912 && fa.t.fraction <= fb.t.fraction)
913 /* compare fractions */
914 result.q = IEEE_TRUE;
916 result.q = IEEE_FALSE;
924 ieee_convert_S_T(fp_register_t f, int rnd,
925 u_int64_t control, u_int64_t *status)
928 * Handle exceptional values.
931 /* Instructions Descriptions (I) section 4.7.10.1 */
932 f.t.fraction |= QNAN_BIT;
935 if (isQNaN(f) || isINF(f))
939 * If the number is a denormalised float, renormalise.
942 return makefloat(f.t.sign,
945 T_FORMAT, rnd, control, status);
951 ieee_convert_T_S(fp_register_t f, int rnd,
952 u_int64_t control, u_int64_t *status)
955 * Handle exceptional values.
958 /* Instructions Descriptions (I) section 4.7.10.1 */
959 f.t.fraction |= QNAN_BIT;
960 f.t.fraction &= ~S_FRACMASK;
963 if (isQNaN(f) || isINF(f))
966 return makefloat(f.t.sign,
969 S_FORMAT, rnd, control, status);
973 ieee_convert_Q_S(fp_register_t f, int rnd,
974 u_int64_t control, u_int64_t *status)
976 u_int64_t frac = f.q;
986 * We shift up one bit to leave the sticky bit clear. This is
987 * possible unless frac == (1<<63), in which case the sticky
988 * bit is already clear.
990 exponent = T_FRACBITS + FRAC_SHIFT;
991 if (frac < (1ULL << 63)) {
996 return makefloat(sign, exponent, frac, S_FORMAT, rnd,
1001 ieee_convert_Q_T(fp_register_t f, int rnd,
1002 u_int64_t control, u_int64_t *status)
1004 u_int64_t frac = f.q;
1014 * We shift up one bit to leave the sticky bit clear. This is
1015 * possible unless frac == (1<<63), in which case the sticky
1016 * bit is already clear.
1018 exponent = T_FRACBITS + FRAC_SHIFT;
1019 if (frac < (1ULL << 63)) {
1024 return makefloat(sign, exponent, frac, T_FORMAT, rnd,
1029 ieee_convert_T_Q(fp_register_t f, int rnd,
1030 u_int64_t control, u_int64_t *status)
1036 * Handle exceptional values.
1039 /* Instructions Descriptions (I) section 4.7.10.1 */
1041 *status |= FPCR_INV;
1046 /* Instructions Descriptions (I) section 4.7.10.1 */
1047 *status |= FPCR_INV;
1052 exp = getexp(f, T_FORMAT) - (T_FRACBITS + FRAC_SHIFT);
1056 if (exp > 64 || frac >= (1 << (64 - exp)))
1057 *status |= FPCR_IOV | FPCR_INE;
1062 } else if (exp < 0) {
1064 u_int64_t fraclo, frachi;
1065 u_int64_t diffhi, difflo;
1073 *status |= FPCR_INE;
1076 } else if (exp == 64) {
1082 *status |= FPCR_INE;
1089 mask = (1 << exp) - 1;
1090 fraclo = frac >> exp;
1092 frachi = fraclo + 1;
1093 difflo = frac - (fraclo << exp);
1094 diffhi = (frachi << exp) - frac;
1095 *status |= FPCR_INE;
1106 case ROUND_MINUS_INF:
1117 if (difflo < diffhi)
1119 else if (diffhi > difflo)
1121 else if (fraclo & 1)
1132 case ROUND_PLUS_INF:
1142 if (frac > (1ULL << 63))
1143 *status |= FPCR_IOV | FPCR_INE;
1146 if (frac > (1ULL << 63) - 1)
1147 *status |= FPCR_IOV | FPCR_INE;
1155 ieee_convert_S_Q(fp_register_t f, int rnd,
1156 u_int64_t control, u_int64_t *status)
1158 f = ieee_convert_S_T(f, rnd, control, status);
1159 return ieee_convert_T_Q(f, rnd, control, status);
1180 a.r.t.fraction = ((long long)random() & (1ULL << 20) - 1) << 32
1182 exp = random() & 0x7ff;
1185 exp = 1; /* no denorms */
1186 else if (exp == 0x7ff)
1187 exp = 0x7fe; /* no NaNs and INFs */
1190 a.r.t.exponent = exp;
1191 a.r.t.sign = random() & 1;
1201 a.r.t.fraction = ((long)random() & (1ULL << 23) - 1) << 29;
1202 exp = random() & 0xff;
1205 exp = 1; /* no denorms */
1206 else if (exp == 0xff)
1207 exp = 0xfe; /* no NaNs and INFs */
1210 /* map exponent from S to T format */
1212 a.r.t.exponent = 0x7ff;
1213 else if (exp & 0x80)
1214 a.r.t.exponent = 0x400 + (exp & 0x7f);
1216 a.r.t.exponent = 0x380 + exp;
1219 a.r.t.sign = random() & 1;
1225 * Ignore epsilon errors
1228 equal_T(union value a, union value b)
1230 if (isZERO(a.r) && isZERO(b.r))
1232 if (a.r.t.sign != b.r.t.sign)
1234 if (a.r.t.exponent != b.r.t.exponent)
1237 return a.r.t.fraction == b.r.t.fraction;
1241 equal_S(union value a, union value b)
1243 int64_t epsilon = 1ULL << 29;
1245 if (isZERO(a.r) && isZERO(b.r))
1247 if (a.r.t.sign != b.r.t.sign)
1249 if (a.r.t.exponent != b.r.t.exponent)
1252 return ((a.r.t.fraction & ~(epsilon-1))
1253 == (b.r.t.fraction & ~(epsilon-1)));
1256 #define ITER 1000000
1261 union value a, b, c, x;
1262 u_int64_t status = 0;
1265 for (i = 0; i < ITER; i++) {
1266 a.d = random_double();
1267 b.d = random_double();
1269 c.r = ieee_add(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1271 /* ignore NaN and INF */
1272 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1275 if (!equal_T(c, x)) {
1276 printf("bad double add, %g + %g = %g (should be %g)\n",
1277 a.d, b.d, c.d, x.d);
1278 c.r = ieee_add(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1287 union value a, b, c, x, t;
1289 u_int64_t status = 0;
1292 for (i = 0; i < ITER; i++) {
1295 a.r.q = 0xb33acf292ca49700ULL;
1296 b.r.q = 0xcad3191058a693aeULL;
1299 a.d = random_float();
1300 b.d = random_float();
1302 c.r = ieee_add(a.r, b.r, S_FORMAT, ROUND_NORMAL,
1304 /* ignore NaN and INF */
1305 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1309 t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
1310 if (!equal_S(t, x)) {
1311 printf("bad single add, %g + %g = %g (should be %g)\n",
1312 a.d, b.d, t.d, x.d);
1313 c.r = ieee_add(a.r, b.r, S_FORMAT, ROUND_NORMAL,
1322 union value a, b, c, x;
1323 u_int64_t status = 0;
1326 for (i = 0; i < ITER; i++) {
1327 a.d = random_double();
1328 b.d = random_double();
1330 c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1332 /* ignore NaN and INF */
1333 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1336 if (!equal_T(c, x)) {
1337 printf("bad double mul, %g * %g = %g (should be %g)\n",
1338 a.d, b.d, c.d, x.d);
1339 c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1348 union value a, b, c, x, t;
1350 u_int64_t status = 0;
1353 for (i = 0; i < ITER; i++) {
1354 a.d = random_double();
1355 b.d = random_double();
1357 c.r = ieee_mul(a.r, b.r, S_FORMAT, ROUND_NORMAL,
1359 /* ignore NaN and INF */
1360 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1364 t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
1365 if (!equal_S(t, x)) {
1366 printf("bad single mul, %g * %g = %g (should be %g)\n",
1367 a.d, b.d, t.d, x.d);
1368 c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1377 union value a, b, c, x;
1378 u_int64_t status = 0;
1381 for (i = 0; i < ITER; i++) {
1382 a.d = random_double();
1383 b.d = random_double();
1385 c.r = ieee_div(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1387 /* ignore NaN and INF */
1388 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1391 if (!equal_T(c, x) && !isZERO(x.r)) {
1392 printf("bad double div, %g / %g = %g (should be %g)\n",
1393 a.d, b.d, c.d, x.d);
1394 c.r = ieee_div(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1403 union value a, b, c, x, t;
1405 u_int64_t status = 0;
1408 for (i = 0; i < ITER; i++) {
1409 a.d = random_double();
1410 b.d = random_double();
1412 c.r = ieee_div(a.r, b.r, S_FORMAT, ROUND_NORMAL,
1414 /* ignore NaN and INF */
1415 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1419 t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
1420 if (!equal_S(t, x)) {
1421 printf("bad single div, %g / %g = %g (should be %g)\n",
1422 a.d, b.d, t.d, x.d);
1423 c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1430 test_convert_int_to_double()
1432 union value a, c, x;
1433 u_int64_t status = 0;
1436 for (i = 0; i < ITER; i++) {
1437 a.r.q = (u_int64_t)random() << 32
1440 c.r = ieee_convert_Q_T(a.r, ROUND_NORMAL, 0, &status);
1441 /* ignore NaN and INF */
1442 if (isNaN(c.r) || isINF(c.r))
1444 x.d = (double) a.r.q;
1446 printf("bad convert double, (double)%qx = %g (should be %g)\n",
1448 c.r = ieee_convert_Q_T(a.r, ROUND_NORMAL, 0, &status);
1454 test_convert_int_to_single()
1456 union value a, c, x, t;
1458 u_int64_t status = 0;
1461 for (i = 0; i < ITER; i++) {
1462 a.r.q = (unsigned long long)random() << 32
1465 c.r = ieee_convert_Q_S(a.r, ROUND_NORMAL, 0, &status);
1466 /* ignore NaN and INF */
1467 if (isNaN(c.r) || isINF(c.r))
1471 t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
1473 printf("bad convert single, (double)%qx = %g (should be %g)\n",
1475 c.r = ieee_convert_Q_S(a.r, ROUND_NORMAL, 0, &status);
1481 test_convert_double_to_int()
1484 u_int64_t status = 0;
1487 for (i = 0; i < ITER; i++) {
1488 a.d = random_double();
1490 c.r = ieee_convert_T_Q(a.r, ROUND_NORMAL, 0, &status);
1491 if ((int)c.r.q != (int)a.d) {
1492 printf("bad convert double, (int)%g = %d (should be %d)\n",
1493 a.d, (int)c.r.q, (int)a.d);
1494 c.r = ieee_convert_T_Q(a.r, ROUND_NORMAL, 0, &status);
1500 main(int argc, char* argv[])
1510 test_convert_int_to_double();
1511 test_convert_int_to_single();
1513 /* x86 generates SIGFPE on overflows. */
1514 test_convert_double_to_int();