sys/alpha/alpha/ieee_float.c

   1 /*-
   2  * Copyright (c) 1998 Doug Rabson
   3  * All rights reserved.
   4  *
   5  * Redistribution and use in source and binary forms, with or without
   6  * modification, are permitted provided that the following conditions
   7  * are met:
   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.
  13  *
  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
  24  * SUCH DAMAGE.
  25  *
  26  * $FreeBSD$
  27  */
  28
  29 /*
  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.
  34  *
  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.
  38  */
  39
  40 #include <sys/types.h>
  41 #ifdef TEST
  42 #include "../include/fpu.h"
  43 #include "ieee_float.h"
  44 #else
  45 #include <sys/param.h>
  46 #include <sys/systm.h>
  47 #include <sys/sysproto.h>
  48 #include <sys/sysent.h>
  49 #include <sys/proc.h>
  50 #include <machine/fpu.h>
  51 #include <alpha/alpha/ieee_float.h>
  52 #endif
  53
  54 /*
  55  * The number of fraction bits in a T format float.
  56  */
  57 #define T_FRACBITS      52
  58
  59 /*
  60  * The number of fraction bits in a S format float.
  61  */
  62 #define S_FRACBITS      23
  63
  64 /*
  65  * Mask the fraction part of a float to contain only those bits which
  66  * should be in single precision number.
  67  */
  68 #define S_FRACMASK      ((1ULL << 52) - (1ULL << 29))
  69
  70 /*
  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.
  74  */
  75 #define FRAC_SHIFT      3
  76
  77 /*
  78  * Values for 1.0 and 2.0 fractions (including the extra FRAC_SHIFT
  79  * bits).
  80  */
  81 #define ONE             (1ULL << (T_FRACBITS + FRAC_SHIFT))
  82 #define TWO             (ONE + ONE)
  83
  84 /*
  85  * The maximum and minimum values for S and T format exponents.
  86  */
  87 #define T_MAXEXP        0x3ff
  88 #define T_MINEXP        -0x3fe
  89 #define S_MAXEXP        0x7f
  90 #define S_MINEXP        -0x7e
  91
  92 /*
  93  * Exponent values in registers are biased by adding this value.
  94  */
  95 #define BIAS_EXP        0x3ff
  96
  97 /*
  98  * Exponent value for INF and NaN.
  99  */
 100 #define NAN_EXP         0x7ff
 101
 102 /*
 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.
 105  */
 106 #define QNAN_BIT        (1ULL << 51)
 107
 108 /*
 109  * Return true if the number is any kind of NaN.
 110  */
 111 static __inline int
 112 isNaN(fp_register_t f)
 113 {
 114         return f.t.exponent == NAN_EXP && f.t.fraction != 0;
 115 }
 116
 117 /*
 118  * Return true if the number is a quiet NaN.
 119  */
 120 static __inline int
 121 isQNaN(fp_register_t f)
 122 {
 123         return f.t.exponent == NAN_EXP && (f.t.fraction & QNAN_BIT);
 124 }
 125
 126 /*
 127  * Return true if the number is a signalling NaN.
 128  */
 129 static __inline int
 130 isSNaN(fp_register_t f)
 131 {
 132         return isNaN(f) && !isQNaN(f);
 133 }
 134
 135 /*
 136  * Return true if the number is +/- INF.
 137  */
 138 static __inline int
 139 isINF(fp_register_t f)
 140 {
 141         return f.t.exponent == NAN_EXP && f.t.fraction == 0;
 142 }
 143
 144 /*
 145  * Return true if the number is +/- 0.
 146  */
 147 static __inline int
 148 isZERO(fp_register_t f)
 149 {
 150         return f.t.exponent == 0 && f.t.fraction == 0;
 151 }
 152
 153 /*
 154  * Return true if the number is denormalised.
 155  */
 156 static __inline int
 157 isDENORM(fp_register_t f)
 158 {
 159         return f.t.exponent == 0 && f.t.fraction != 0;
 160 }
 161
 162 /*
 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.
 166  */
 167 static __inline int
 168 getexp(fp_register_t f, int src)
 169 {
 170         int minexp[] = { S_MINEXP, 0, T_MINEXP, 0 };
 171         if (f.t.exponent == 0) {
 172                 if (f.t.fraction)
 173                         return minexp[src];
 174                 else
 175                         return 0;
 176         }
 177         return f.t.exponent - BIAS_EXP;
 178 }
 179
 180 /*
 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.
 184  */
 185 static __inline u_int64_t
 186 getfrac(fp_register_t f)
 187 {
 188         if (f.t.exponent == 0)
 189                 return f.t.fraction << FRAC_SHIFT;
 190         else
 191                 return (f.t.fraction << FRAC_SHIFT) | ONE;
 192 }
 193
 194 /*
 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.
 198  */
 199 static fp_register_t
 200 makefloat(int sign, int exp, u_int64_t frac,
 201           int src, int rnd,
 202           u_int64_t control, u_int64_t *status)
 203 {
 204         fp_register_t f;
 205         int minexp = 0, maxexp = 0, alpha = 0;
 206         u_int64_t epsilon = 0, max = 0;
 207
 208         if (frac == 0) {
 209                 f.t.sign = sign;
 210                 f.t.exponent = 0;
 211                 f.t.fraction = 0;
 212                 return f;
 213         }
 214
 215         if (frac >= TWO) {
 216                 /*
 217                  * Fraction is >= 2.0.
 218                  * Shift the fraction down, preserving the 'sticky'
 219                  * bit.
 220                  */
 221                 while (frac >= TWO) {
 222                         frac = (frac >> 1) | (frac & 1);
 223                         exp++;
 224                 }
 225         } else if (frac < ONE) {
 226                 /*
 227                  * Fraction is < 1.0. Shift it up.
 228                  */
 229                 while (frac < ONE) {
 230                         frac = (frac << 1) | (frac & 1);
 231                         exp--;
 232                 }
 233         }
 234
 235         switch (src) {
 236         case S_FORMAT:
 237                 minexp = S_MINEXP;
 238                 maxexp = S_MAXEXP;
 239                 alpha = 0xc0;
 240                 epsilon = (1ULL << (T_FRACBITS - S_FRACBITS + FRAC_SHIFT));
 241                 max = TWO - epsilon;
 242                 break;
 243
 244         case T_FORMAT:
 245                 minexp = T_MINEXP;
 246                 maxexp = T_MAXEXP;
 247                 alpha = 0x600;
 248                 epsilon = (1ULL << FRAC_SHIFT);
 249                 max = TWO - epsilon;
 250                 break;
 251         }
 252
 253         /*
 254          * Handle underflow before rounding so that denormalised
 255          * numbers are rounded correctly.
 256          */
 257         if (exp < minexp) {
 258                 *status |= FPCR_INE;
 259                 if (control & IEEE_TRAP_ENABLE_UNF) {
 260                         *status |= FPCR_UNF;
 261                         exp += alpha;
 262                 } else {
 263                         /* denormalise */
 264                         while (exp < minexp) {
 265                                 exp++;
 266                                 frac = (frac >> 1) | (frac & 1);
 267                         }
 268                         exp = minexp - 1;
 269                 }
 270         }
 271
 272         /*
 273          * Round the fraction according to the rounding mode.
 274          */
 275         if (frac & (epsilon - 1)) {
 276                 u_int64_t fraclo, frachi;
 277                 u_int64_t difflo, diffhi;
 278
 279                 fraclo = frac & max;
 280                 frachi = fraclo + epsilon;
 281                 switch (rnd) {
 282                 case ROUND_CHOP:
 283                         frac = fraclo;
 284                         break;
 285                 case ROUND_MINUS_INF:
 286                         if (f.t.sign)
 287                                 frac = frachi;
 288                         else
 289                                 frac = fraclo;
 290                         break;
 291                 case ROUND_NORMAL:
 292                         difflo = frac - fraclo;
 293                         diffhi = frachi - frac;
 294                         if (difflo < diffhi)
 295                                 frac = fraclo;
 296                         else if (diffhi < difflo)
 297                                 frac = frachi;
 298                         else
 299                                 /* round to even */
 300                                 if (fraclo & epsilon)
 301                                         frac = frachi;
 302                                 else
 303                                         frac = fraclo;
 304                         break;
 305                 case ROUND_PLUS_INF:
 306                         if (f.t.sign)
 307                                 frac = fraclo;
 308                         else
 309                                 frac = frachi;
 310                         break;
 311                 }
 312
 313                 /*
 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.
 318                  */
 319                 if (frac >= TWO) {
 320                         frac = (frac >> 1) & ~(epsilon - 1);
 321                         exp++;
 322                 } else if (exp == minexp - 1 && frac == ONE) {
 323                         /* Renormalise to ONE * 2^minexp */
 324                         exp = minexp;
 325                 }
 326
 327                 *status |= FPCR_INE;
 328         }
 329
 330         /*
 331          * Check for overflow and round to the correct INF as needed.
 332          */
 333         if (exp > maxexp) {
 334                 *status |= FPCR_OVF | FPCR_INE;
 335                 if (control & IEEE_TRAP_ENABLE_OVF) {
 336                         exp -= alpha;
 337                 } else {
 338                         switch (rnd) {
 339                         case ROUND_CHOP:
 340                                 exp = maxexp;
 341                                 frac = max;
 342                                 break;
 343                         case ROUND_MINUS_INF:
 344                                 if (sign) {
 345                                         exp = maxexp + 1; /* INF */
 346                                         frac = 0;
 347                                 } else {
 348                                         exp = maxexp;
 349                                         frac = max;
 350                                 }
 351                                 break;
 352                         case ROUND_NORMAL:
 353                                 exp = maxexp + 1; /* INF */
 354                                 frac = 0;
 355                                 break;
 356                         case ROUND_PLUS_INF:
 357                                 if (sign) {
 358                                         exp = maxexp;
 359                                         frac = max;
 360                                 } else {
 361                                         exp = maxexp + 1; /* INF */
 362                                         frac = 0;
 363                                 }
 364                                 break;
 365                         }
 366                 }
 367         }
 368
 369         f.t.sign = sign;
 370         if (exp > maxexp)       /* NaN, INF */
 371                 f.t.exponent = NAN_EXP;
 372         else if (exp < minexp)  /* denorm, zero */
 373                 f.t.exponent = 0;
 374         else
 375                 f.t.exponent = exp + BIAS_EXP;
 376         f.t.fraction = (frac & ~ONE) >> FRAC_SHIFT;
 377         return f;
 378 }
 379
 380 /*
 381  * Return the canonical quiet NaN in register format.
 382  */
 383 static fp_register_t
 384 makeQNaN(void)
 385 {
 386         fp_register_t f;
 387         f.t.sign = 0;
 388         f.t.exponent = NAN_EXP;
 389         f.t.fraction = QNAN_BIT;
 390         return f;
 391 }
 392
 393 /*
 394  * Return +/- INF.
 395  */
 396 static fp_register_t
 397 makeINF(int sign)
 398 {
 399         fp_register_t f;
 400         f.t.sign = sign;
 401         f.t.exponent = NAN_EXP;
 402         f.t.fraction = 0;
 403         return f;
 404 }
 405
 406 /*
 407  * Return +/- 0.
 408  */
 409 static fp_register_t
 410 makeZERO(int sign)
 411 {
 412         fp_register_t f;
 413         f.t.sign = sign;
 414         f.t.exponent = 0;
 415         f.t.fraction = 0;
 416         return f;
 417 }
 418
 419 fp_register_t
 420 ieee_add(fp_register_t fa, fp_register_t fb,
 421          int src, int rnd,
 422          u_int64_t control, u_int64_t *status)
 423 {
 424         int shift;
 425         int expa, expb, exp;
 426         u_int64_t fraca, fracb, frac;
 427         int sign, sticky;
 428
 429         /* First handle NaNs */
 430         if (isNaN(fa) || isNaN(fb)) {
 431                 fp_register_t result;
 432
 433                 /* Instructions Descriptions (I) section 4.7.10.4 */
 434                 if (isQNaN(fb))
 435                         result = fb;
 436                 else if (isSNaN(fb)) {
 437                         result = fb;
 438                         result.t.fraction |= QNAN_BIT;
 439                 } else if (isQNaN(fa))
 440                         result = fa;
 441                 else if (isSNaN(fa)) {
 442                         result = fa;
 443                         result.t.fraction |= QNAN_BIT;
 444                 }
 445
 446                 /* If either operand is a signalling NaN, trap. */
 447                 if (isSNaN(fa) || isSNaN(fb))
 448                         *status |= FPCR_INV;
 449
 450                 return result;
 451         }
 452
 453         /* Handle +/- INF */
 454         if (isINF(fa))
 455                 if (isINF(fb))
 456                         if (fa.t.sign != fb.t.sign) {
 457                                 /* If adding -INF to +INF, generate a trap. */
 458                                 *status |= FPCR_INV;
 459                                 return makeQNaN();
 460                         } else
 461                                 return fa;
 462                 else
 463                         return fa;
 464         else if (isINF(fb))
 465                 return fb;
 466
 467         /*
 468          * Unpack the registers.
 469          */
 470         expa = getexp(fa, src);
 471         expb = getexp(fb, src);
 472         fraca = getfrac(fa);
 473         fracb = getfrac(fb);
 474         shift = expa - expb;
 475         if (shift < 0) {
 476                 shift = -shift;
 477                 exp = expb;
 478                 sticky = (fraca & ((1ULL << shift) - 1)) != 0;
 479                 if (shift >= 64)
 480                         fraca = sticky;
 481                 else
 482                         fraca = (fraca >> shift) | sticky;
 483         } else if (shift > 0) {
 484                 exp = expa;
 485                 sticky = (fracb & ((1ULL << shift) - 1)) != 0;
 486                 if (shift >= 64)
 487                         fracb = sticky;
 488                 else
 489                         fracb = (fracb >> shift) | sticky;
 490         } else
 491                 exp = expa;
 492         if (fa.t.sign) fraca = -fraca;
 493         if (fb.t.sign) fracb = -fracb;
 494         frac = fraca + fracb;
 495         if (frac >> 63) {
 496                 sign = 1;
 497                 frac = -frac;
 498         } else
 499                 sign = 0;
 500
 501         /* -0 + -0 = -0 */
 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;
 505
 506         return makefloat(sign, exp, frac, src, rnd, control, status);
 507 }
 508
 509 fp_register_t
 510 ieee_sub(fp_register_t fa, fp_register_t fb,
 511          int src, int rnd,
 512          u_int64_t control, u_int64_t *status)
 513 {
 514         fb.t.sign = !fb.t.sign;
 515         return ieee_add(fa, fb, src, rnd, control, status);
 516 }
 517
 518 typedef struct {
 519         u_int64_t lo;
 520         u_int64_t hi;
 521 } u_int128_t;
 522
 523 #define SRL128(x, b)                            \
 524 do {                                            \
 525         x.lo >>= b;                             \
 526         x.lo |= x.hi << (64 - b);               \
 527         x.hi >>= b;                             \
 528 } while (0)
 529
 530 #define SLL128(x, b)                            \
 531 do {                                            \
 532         if (b >= 64) {                          \
 533                 x.hi = x.lo << (b - 64);        \
 534                 x.lo = 0;                       \
 535         } else {                                \
 536                 x.hi <<= b;                     \
 537                 x.hi |= x.lo >> (64 - b);       \
 538                 x.lo <<= b;                     \
 539         }                                       \
 540 } while (0)
 541
 542 #define SUB128(a, b)                            \
 543 do {                                            \
 544         int borrow = a.lo < b.lo;               \
 545         a.lo = a.lo - b.lo;                     \
 546         a.hi = a.hi - b.hi - borrow;            \
 547 } while (0)
 548
 549 #define LESS128(a, b) (a.hi < b.hi || (a.hi == b.hi && a.lo < b.lo))
 550
 551 fp_register_t
 552 ieee_mul(fp_register_t fa, fp_register_t fb,
 553          int src, int rnd,
 554          u_int64_t control, u_int64_t *status)
 555 {
 556         int expa, expb, exp;
 557         u_int64_t fraca, fracb, tmp;
 558         u_int128_t frac;
 559         int sign;
 560
 561         /* First handle NaNs */
 562         if (isNaN(fa) || isNaN(fb)) {
 563                 fp_register_t result;
 564
 565                 /* Instructions Descriptions (I) section 4.7.10.4 */
 566                 if (isQNaN(fb))
 567                         result = fb;
 568                 else if (isSNaN(fb)) {
 569                         result = fb;
 570                         result.t.fraction |= QNAN_BIT;
 571                 } else if (isQNaN(fa))
 572                         result = fa;
 573                 else if (isSNaN(fa)) {
 574                         result = fa;
 575                         result.t.fraction |= QNAN_BIT;
 576                 }
 577
 578                 /* If either operand is a signalling NaN, trap. */
 579                 if (isSNaN(fa) || isSNaN(fb))
 580                         *status |= FPCR_INV;
 581
 582                 return result;
 583         }
 584
 585         /* Handle INF and 0 */
 586         if ((isINF(fa) && isZERO(fb)) || (isZERO(fa) && isINF(fb))) {
 587                 /* INF * 0 = NaN */
 588                 *status |= FPCR_INV;
 589                 return makeQNaN();
 590         } else
 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);
 596
 597         /*
 598          * Unpack the registers.
 599          */
 600         expa = getexp(fa, src);
 601         expb = getexp(fb, src);
 602         fraca = getfrac(fa);
 603         fracb = getfrac(fb);
 604         sign = fa.t.sign ^ fb.t.sign;
 605
 606 #define LO32(x) ((x) & ((1ULL << 32) - 1))
 607 #define HI32(x) ((x) >> 32)
 608
 609         /*
 610          * Calculate the 128bit result of multiplying fraca and fracb.
 611          */
 612         frac.lo = fraca * fracb;
 613 #ifdef __alpha__
 614         /*
 615          * The alpha has a handy instruction to find the high word.
 616          */
 617         __asm__ __volatile__ ("umulh %1,%2,%0"
 618                               : "=r"(tmp)
 619                               : "r"(fraca), "r"(fracb));
 620         frac.hi = tmp;
 621 #else
 622         /*
 623          * Do the multiply longhand otherwise.
 624          */
 625         frac.hi = HI32(LO32(fraca) * HI32(fracb)
 626                        + HI32(fraca) * LO32(fracb)
 627                        + HI32(LO32(fraca) * LO32(fracb)))
 628                 + HI32(fraca) * HI32(fracb);
 629 #endif
 630         exp = expa + expb - (T_FRACBITS + FRAC_SHIFT);
 631
 632         while (frac.hi > 0) {
 633                 int sticky;
 634                 exp++;
 635                 sticky = frac.lo & 1;
 636                 SRL128(frac, 1);
 637                 frac.lo |= sticky;
 638         }
 639
 640         return makefloat(sign, exp, frac.lo, src, rnd, control, status);
 641 }
 642
 643 static u_int128_t
 644 divide_128(u_int128_t a, u_int128_t b)
 645 {
 646         u_int128_t result;
 647         u_int64_t bit;
 648         int i;
 649
 650         /*
 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.
 655          */
 656 #ifdef TEST
 657         if (a.hi == 0 || b.hi != 0 || b.lo == 0)
 658                 abort();
 659 #endif
 660
 661         /*
 662          * Find out how many bits of zeros are at the beginning of the divisor.
 663          */
 664         i = 64;
 665         bit = 1ULL << 63;
 666         while (i < 127) {
 667                 if (b.lo & bit)
 668                         break;
 669                 i++;
 670                 bit >>= 1;
 671         }
 672
 673         /*
 674          * Find out how much to shift the divisor so that its msb
 675          * matches the msb of the dividend.
 676          */
 677         bit = 1ULL << 63;
 678         while (i) {
 679                 if (a.hi & bit)
 680                         break;
 681                 i--;
 682                 bit >>= 1;
 683         }
 684
 685         result.lo = 0;
 686         result.hi = 0;
 687         SLL128(b, i);
 688
 689         /*
 690          * Calculate the result in two parts to avoid keeping a 128bit
 691          * value for the result bit.
 692          */
 693         if (i >= 64) {
 694                 bit = 1ULL << (i - 64);
 695                 while (bit) {
 696                         if (!LESS128(a, b)) {
 697                                 result.hi |= bit;
 698                                 SUB128(a, b);
 699                                 if (!a.lo && !a.hi)
 700                                         return result;
 701                         }
 702                         bit >>= 1;
 703                         SRL128(b, 1);
 704                 }
 705                 i = 63;
 706         }
 707         bit = 1ULL << i;
 708         while (bit) {
 709                 if (!LESS128(a, b)) {
 710                         result.lo |= bit;
 711                         SUB128(a, b);
 712                         if (!a.lo && !a.hi)
 713                                 return result;
 714                 }
 715                 bit >>= 1;
 716                 SRL128(b, 1);
 717         }
 718
 719         return result;
 720 }
 721
 722 fp_register_t
 723 ieee_div(fp_register_t fa, fp_register_t fb,
 724          int src, int rnd,
 725          u_int64_t control, u_int64_t *status)
 726 {
 727         int expa, expb, exp;
 728         u_int128_t fraca, fracb, frac;
 729         int sign;
 730
 731         /* First handle NaNs, INFs and ZEROs */
 732         if (isNaN(fa) || isNaN(fb)) {
 733                 fp_register_t result;
 734
 735                 /* Instructions Descriptions (I) section 4.7.10.4 */
 736                 if (isQNaN(fb))
 737                         result = fb;
 738                 else if (isSNaN(fb)) {
 739                         result = fb;
 740                         result.t.fraction |= QNAN_BIT;
 741                 } else if (isQNaN(fa))
 742                         result = fa;
 743                 else if (isSNaN(fa)) {
 744                         result = fa;
 745                         result.t.fraction |= QNAN_BIT;
 746                 }
 747
 748                 /* If either operand is a signalling NaN, trap. */
 749                 if (isSNaN(fa) || isSNaN(fb))
 750                         *status |= FPCR_INV;
 751
 752                 return result;
 753         }
 754
 755         /* Handle INF and 0 */
 756         if (isINF(fa) && isINF(fb)) {
 757                 *status |= FPCR_INV;
 758                 return makeQNaN();
 759         } else if (isZERO(fb))
 760                 if (isZERO(fa)) {
 761                         *status |= FPCR_INV;
 762                         return makeQNaN();
 763                 } else {
 764                         *status |= FPCR_DZE;
 765                         return makeINF(fa.t.sign ^ fb.t.sign);
 766                 }
 767         else if (isZERO(fa))
 768                 return makeZERO(fa.t.sign ^ fb.t.sign);
 769
 770         /*
 771          * Unpack the registers.
 772          */
 773         expa = getexp(fa, src);
 774         expb = getexp(fb, src);
 775         fraca.hi = getfrac(fa);
 776         fraca.lo = 0;
 777         fracb.lo = getfrac(fb);
 778         fracb.hi = 0;
 779         sign = fa.t.sign ^ fb.t.sign;
 780
 781         frac = divide_128(fraca, fracb);
 782
 783         exp = expa - expb - (64 - T_FRACBITS - FRAC_SHIFT);
 784         while (frac.hi > 0) {
 785                 int sticky;
 786                 exp++;
 787                 sticky = frac.lo & 1;
 788                 SRL128(frac, 1);
 789                 frac.lo |= sticky;
 790         }
 791         frac.lo |= 1;
 792
 793         return makefloat(sign, exp, frac.lo, src, rnd, control, status);
 794 }
 795
 796 #define IEEE_TRUE 0x4000000000000000ULL
 797 #define IEEE_FALSE 0
 798
 799 fp_register_t
 800 ieee_cmpun(fp_register_t fa, fp_register_t fb, u_int64_t *status)
 801 {
 802         fp_register_t result;
 803         if (isNaN(fa) || isNaN(fb)) {
 804                 if (isSNaN(fa) || isSNaN(fb))
 805                         *status |= FPCR_INV;
 806                 result.q = IEEE_TRUE;
 807         } else
 808                 result.q = IEEE_FALSE;
 809
 810         return result;
 811 }
 812
 813 fp_register_t
 814 ieee_cmpeq(fp_register_t fa, fp_register_t fb, u_int64_t *status)
 815 {
 816         fp_register_t result;
 817         if (isNaN(fa) || isNaN(fb)) {
 818                 if (isSNaN(fa) || isSNaN(fb))
 819                         *status |= FPCR_INV;
 820                 result.q = IEEE_FALSE;
 821         } else {
 822                 if (isZERO(fa) && isZERO(fb))
 823                         result.q = IEEE_TRUE;
 824                 else if (fa.q == fb.q)
 825                         result.q = IEEE_TRUE;
 826                 else
 827                         result.q = IEEE_FALSE;
 828         }
 829
 830         return result;
 831 }
 832
 833 fp_register_t
 834 ieee_cmplt(fp_register_t fa, fp_register_t fb, u_int64_t *status)
 835 {
 836         fp_register_t result;
 837         if (isNaN(fa) || isNaN(fb)) {
 838                 if (isSNaN(fa) || isSNaN(fb))
 839                         *status |= FPCR_INV;
 840                 result.q = IEEE_FALSE;
 841         } else {
 842                 if (isZERO(fa) && isZERO(fb))
 843                         result.q = IEEE_FALSE;
 844                 else if (fa.t.sign) {
 845                         /* fa is negative */
 846                         if (!fb.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;
 856                         else
 857                                 result.q = IEEE_FALSE;
 858                 } else {
 859                         /* fa is positive */
 860                         if (fb.t.sign)
 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;
 870                         else
 871                                 result.q = IEEE_FALSE;
 872                 }
 873         }
 874
 875         return result;
 876 }
 877
 878 fp_register_t
 879 ieee_cmple(fp_register_t fa, fp_register_t fb, u_int64_t *status)
 880 {
 881         fp_register_t result;
 882         if (isNaN(fa) || isNaN(fb)) {
 883                 if (isSNaN(fa) || isSNaN(fb))
 884                         *status |= FPCR_INV;
 885                 result.q = IEEE_FALSE;
 886         } else {
 887                 if (isZERO(fa) && isZERO(fb))
 888                         result.q = IEEE_TRUE;
 889                 else if (fa.t.sign) {
 890                         /* fa is negative */
 891                         if (!fb.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;
 901                         else
 902                                 result.q = IEEE_FALSE;
 903                 } else {
 904                         /* fa is positive */
 905                         if (fb.t.sign)
 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;
 915                         else
 916                                 result.q = IEEE_FALSE;
 917                 }
 918         }
 919
 920         return result;
 921 }
 922
 923 fp_register_t
 924 ieee_convert_S_T(fp_register_t f, int rnd,
 925                  u_int64_t control, u_int64_t *status)
 926 {
 927         /*
 928          * Handle exceptional values.
 929          */
 930         if (isNaN(f)) {
 931                 /* Instructions Descriptions (I) section 4.7.10.1 */
 932                 f.t.fraction |= QNAN_BIT;
 933                 *status |= FPCR_INV;
 934         }
 935         if (isQNaN(f) || isINF(f))
 936                 return f;
 937
 938         /*
 939          * If the number is a denormalised float, renormalise.
 940          */
 941         if (isDENORM(f))
 942                 return makefloat(f.t.sign,
 943                                  getexp(f, S_FORMAT),
 944                                  getfrac(f),
 945                                  T_FORMAT, rnd, control, status);
 946         else
 947                 return f;
 948 }
 949
 950 fp_register_t
 951 ieee_convert_T_S(fp_register_t f, int rnd,
 952                  u_int64_t control, u_int64_t *status)
 953 {
 954         /*
 955          * Handle exceptional values.
 956          */
 957         if (isNaN(f)) {
 958                 /* Instructions Descriptions (I) section 4.7.10.1 */
 959                 f.t.fraction |= QNAN_BIT;
 960                 f.t.fraction &= ~S_FRACMASK;
 961                 *status |= FPCR_INV;
 962         }
 963         if (isQNaN(f) || isINF(f))
 964                 return f;
 965
 966         return makefloat(f.t.sign,
 967                          getexp(f, T_FORMAT),
 968                          getfrac(f),
 969                          S_FORMAT, rnd, control, status);
 970 }
 971
 972 fp_register_t
 973 ieee_convert_Q_S(fp_register_t f, int rnd,
 974                  u_int64_t control, u_int64_t *status)
 975 {
 976         u_int64_t frac = f.q;
 977         int sign, exponent;
 978
 979         if (frac >> 63) {
 980                 sign = 1;
 981                 frac = -frac;
 982         } else
 983                 sign = 0;
 984
 985         /*
 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.
 989          */
 990         exponent = T_FRACBITS + FRAC_SHIFT;
 991         if (frac < (1ULL << 63)) {
 992                 frac <<= 1;
 993                 exponent--;
 994         }
 995
 996         return makefloat(sign, exponent, frac, S_FORMAT, rnd,
 997                          control, status);
 998 }
 999
1000 fp_register_t
1001 ieee_convert_Q_T(fp_register_t f, int rnd,
1002                  u_int64_t control, u_int64_t *status)
1003 {
1004         u_int64_t frac = f.q;
1005         int sign, exponent;
1006
1007         if (frac >> 63) {
1008                 sign = 1;
1009                 frac = -frac;
1010         } else
1011                 sign = 0;
1012
1013         /*
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.
1017          */
1018         exponent = T_FRACBITS + FRAC_SHIFT;
1019         if (frac < (1ULL << 63)) {
1020                 frac <<= 1;
1021                 exponent--;
1022         }
1023
1024         return makefloat(sign, exponent, frac, T_FORMAT, rnd,
1025                          control, status);
1026 }
1027
1028 fp_register_t
1029 ieee_convert_T_Q(fp_register_t f, int rnd,
1030                  u_int64_t control, u_int64_t *status)
1031 {
1032         u_int64_t frac;
1033         int exp;
1034
1035         /*
1036          * Handle exceptional values.
1037          */
1038         if (isNaN(f)) {
1039                 /* Instructions Descriptions (I) section 4.7.10.1 */
1040                 if (isSNaN(f))
1041                         *status |= FPCR_INV;
1042                 f.q = 0;
1043                 return f;
1044         }
1045         if (isINF(f)) {
1046                 /* Instructions Descriptions (I) section 4.7.10.1 */
1047                 *status |= FPCR_INV;
1048                 f.q = 0;
1049                 return f;
1050         }
1051
1052         exp = getexp(f, T_FORMAT) - (T_FRACBITS + FRAC_SHIFT);
1053         frac = getfrac(f);
1054
1055         if (exp > 0) {
1056                 if (exp > 64 || frac >= (1 << (64 - exp)))
1057                         *status |= FPCR_IOV | FPCR_INE;
1058                 if (exp < 64)
1059                         frac <<= exp;
1060                 else
1061                         frac = 0;
1062         } else if (exp < 0) {
1063                 u_int64_t mask;
1064                 u_int64_t fraclo, frachi;
1065                 u_int64_t diffhi, difflo;
1066                 exp = -exp;
1067                 if (exp > 64) {
1068                         fraclo = 0;
1069                         diffhi = 0;
1070                         difflo = 0;
1071                         if (frac) {
1072                                 frachi = 1;
1073                                 *status |= FPCR_INE;
1074                         } else
1075                                 frachi = 0;
1076                 } else if (exp == 64) {
1077                         fraclo = 0;
1078                         if (frac) {
1079                                 frachi = 1;
1080                                 difflo = frac;
1081                                 diffhi = -frac;
1082                                 *status |= FPCR_INE;
1083                         } else {
1084                                 frachi = 0;
1085                                 difflo = 0;
1086                                 diffhi = 0;
1087                         }
1088                 } else {
1089                         mask = (1 << exp) - 1;
1090                         fraclo = frac >> exp;
1091                         if (frac & mask) {
1092                                 frachi = fraclo + 1;
1093                                 difflo = frac - (fraclo << exp);
1094                                 diffhi = (frachi << exp) - frac;
1095                                 *status |= FPCR_INE;
1096                         } else {
1097                                 frachi = fraclo;
1098                                 difflo = 0;
1099                                 diffhi = 0;
1100                         }
1101                 }
1102                 switch (rnd) {
1103                 case ROUND_CHOP:
1104                         frac = fraclo;
1105                         break;
1106                 case ROUND_MINUS_INF:
1107                         if (f.t.sign)
1108                                 frac = frachi;
1109                         else
1110                                 frac = fraclo;
1111                         break;
1112                 case ROUND_NORMAL:
1113 #if 0
1114                         /*
1115                          * Round to nearest.
1116                          */
1117                         if (difflo < diffhi)
1118                                 frac = fraclo;
1119                         else if (diffhi > difflo)
1120                                 frac = frachi;
1121                         else if (fraclo & 1)
1122                                 frac = frachi;
1123                         else
1124                                 frac = fraclo;
1125 #else
1126                         /*
1127                          * Round to zero.
1128                          */
1129                         frac = fraclo;
1130 #endif
1131                         break;
1132                 case ROUND_PLUS_INF:
1133                         if (f.t.sign)
1134                                 frac = fraclo;
1135                         else
1136                                 frac = frachi;
1137                         break;
1138                 }
1139         }
1140
1141         if (f.t.sign) {
1142                 if (frac > (1ULL << 63))
1143                         *status |= FPCR_IOV | FPCR_INE;
1144                 frac = -frac;
1145         } else {
1146                 if (frac > (1ULL << 63) - 1)
1147                         *status |= FPCR_IOV | FPCR_INE;
1148         }
1149
1150         f.q = frac;
1151         return f;
1152 }
1153
1154 fp_register_t
1155 ieee_convert_S_Q(fp_register_t f, int rnd,
1156                  u_int64_t control, u_int64_t *status)
1157 {
1158         f = ieee_convert_S_T(f, rnd, control, status);
1159         return ieee_convert_T_Q(f, rnd, control, status);
1160 }
1161
1162 #ifndef _KERNEL
1163
1164 #include <stdio.h>
1165 #include <math.h>
1166 #include <stdlib.h>
1167
1168 union value {
1169         double d;
1170         fp_register_t r;
1171 };
1172
1173
1174 static double
1175 random_double()
1176 {
1177         union value a;
1178         int exp;
1179
1180         a.r.t.fraction = ((long long)random() & (1ULL << 20) - 1) << 32
1181                 | random();
1182         exp = random() & 0x7ff;
1183 #if 1
1184         if (exp == 0)
1185                 exp = 1;        /* no denorms */
1186         else if (exp == 0x7ff)
1187                 exp = 0x7fe;    /* no NaNs and INFs */
1188 #endif
1189
1190         a.r.t.exponent = exp;
1191         a.r.t.sign = random() & 1;
1192         return a.d;
1193 }
1194
1195 static float
1196 random_float()
1197 {
1198         union value a;
1199         int exp;
1200
1201         a.r.t.fraction = ((long)random() & (1ULL << 23) - 1) << 29;
1202         exp = random() & 0xff;
1203 #if 1
1204         if (exp == 0)
1205                 exp = 1;        /* no denorms */
1206         else if (exp == 0xff)
1207                 exp = 0xfe;     /* no NaNs and INFs */
1208 #endif
1209
1210         /* map exponent from S to T format */
1211         if (exp == 255)
1212                 a.r.t.exponent = 0x7ff;
1213         else if (exp & 0x80)
1214                 a.r.t.exponent = 0x400 + (exp & 0x7f);
1215         else if (exp)
1216                 a.r.t.exponent = 0x380 + exp;
1217         else
1218                 a.r.t.exponent = 0;
1219         a.r.t.sign = random() & 1;
1220
1221         return a.d;
1222 }
1223
1224 /*
1225  * Ignore epsilon errors
1226  */
1227 int
1228 equal_T(union value a, union value b)
1229 {
1230         if (isZERO(a.r) && isZERO(b.r))
1231                 return 1;
1232         if (a.r.t.sign != b.r.t.sign)
1233                 return 0;
1234         if (a.r.t.exponent != b.r.t.exponent)
1235                 return 0;
1236
1237         return a.r.t.fraction == b.r.t.fraction;
1238 }
1239
1240 int
1241 equal_S(union value a, union value b)
1242 {
1243         int64_t epsilon = 1ULL << 29;
1244
1245         if (isZERO(a.r) && isZERO(b.r))
1246                 return 1;
1247         if (a.r.t.sign != b.r.t.sign)
1248                 return 0;
1249         if (a.r.t.exponent != b.r.t.exponent)
1250                 return 0;
1251
1252         return ((a.r.t.fraction & ~(epsilon-1))
1253                 == (b.r.t.fraction & ~(epsilon-1)));
1254 }
1255
1256 #define ITER 1000000
1257
1258 static void
1259 test_double_add()
1260 {
1261         union value a, b, c, x;
1262         u_int64_t status = 0;
1263         int i;
1264
1265         for (i = 0; i < ITER; i++) {
1266                 a.d = random_double();
1267                 b.d = random_double();
1268                 status = 0;
1269                 c.r = ieee_add(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1270                                0, &status);
1271                 /* ignore NaN and INF */
1272                 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1273                         continue;
1274                 x.d = a.d + b.d;
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,
1279                                        0, &status);
1280                 }
1281         }
1282 }
1283
1284 static void
1285 test_single_add()
1286 {
1287         union value a, b, c, x, t;
1288         float xf;
1289         u_int64_t status = 0;
1290         int i;
1291
1292         for (i = 0; i < ITER; i++) {
1293 #if 0
1294                 if (i == 0) {
1295                         a.r.q = 0xb33acf292ca49700ULL;
1296                         b.r.q = 0xcad3191058a693aeULL;
1297                 }
1298 #endif
1299                 a.d = random_float();
1300                 b.d = random_float();
1301                 status = 0;
1302                 c.r = ieee_add(a.r, b.r, S_FORMAT, ROUND_NORMAL,
1303                                0, &status);
1304                 /* ignore NaN and INF */
1305                 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1306                         continue;
1307                 xf = a.d + b.d;
1308                 x.d = xf;
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,
1314                                        0, &status);
1315                 }
1316         }
1317 }
1318
1319 static void
1320 test_double_mul()
1321 {
1322         union value a, b, c, x;
1323         u_int64_t status = 0;
1324         int i;
1325
1326         for (i = 0; i < ITER; i++) {
1327                 a.d = random_double();
1328                 b.d = random_double();
1329                 status = 0;
1330                 c.r = ieee_mul(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1331                                0, &status);
1332                 /* ignore NaN and INF */
1333                 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1334                         continue;
1335                 x.d = a.d * b.d;
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,
1340                                        0, &status);
1341                 }
1342         }
1343 }
1344
1345 static void
1346 test_single_mul()
1347 {
1348         union value a, b, c, x, t;
1349         float xf;
1350         u_int64_t status = 0;
1351         int i;
1352
1353         for (i = 0; i < ITER; i++) {
1354                 a.d = random_double();
1355                 b.d = random_double();
1356                 status = 0;
1357                 c.r = ieee_mul(a.r, b.r, S_FORMAT, ROUND_NORMAL,
1358                                0, &status);
1359                 /* ignore NaN and INF */
1360                 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1361                         continue;
1362                 xf = a.d * b.d;
1363                 x.d = xf;
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,
1369                                        0, &status);
1370                 }
1371         }
1372 }
1373
1374 static void
1375 test_double_div()
1376 {
1377         union value a, b, c, x;
1378         u_int64_t status = 0;
1379         int i;
1380
1381         for (i = 0; i < ITER; i++) {
1382                 a.d = random_double();
1383                 b.d = random_double();
1384                 status = 0;
1385                 c.r = ieee_div(a.r, b.r, T_FORMAT, ROUND_NORMAL,
1386                                0, &status);
1387                 /* ignore NaN and INF */
1388                 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1389                         continue;
1390                 x.d = a.d / b.d;
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,
1395                                        0, &status);
1396                 }
1397         }
1398 }
1399
1400 static void
1401 test_single_div()
1402 {
1403         union value a, b, c, x, t;
1404         float xf;
1405         u_int64_t status = 0;
1406         int i;
1407
1408         for (i = 0; i < ITER; i++) {
1409                 a.d = random_double();
1410                 b.d = random_double();
1411                 status = 0;
1412                 c.r = ieee_div(a.r, b.r, S_FORMAT, ROUND_NORMAL,
1413                                0, &status);
1414                 /* ignore NaN and INF */
1415                 if (isNaN(c.r) || isINF(c.r) || isDENORM(c.r))
1416                         continue;
1417                 xf = a.d / b.d;
1418                 x.d = xf;
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,
1424                                        0, &status);
1425                 }
1426         }
1427 }
1428
1429 static void
1430 test_convert_int_to_double()
1431 {
1432         union value a, c, x;
1433         u_int64_t status = 0;
1434         int i;
1435
1436         for (i = 0; i < ITER; i++) {
1437                 a.r.q = (u_int64_t)random() << 32
1438                         | random();
1439                 status = 0;
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))
1443                         continue;
1444                 x.d = (double) a.r.q;
1445                 if (c.d != x.d) {
1446                         printf("bad convert double, (double)%qx = %g (should be %g)\n",
1447                                a.r.q, c.d, x.d);
1448                         c.r = ieee_convert_Q_T(a.r, ROUND_NORMAL, 0, &status);
1449                 }
1450         }
1451 }
1452
1453 static void
1454 test_convert_int_to_single()
1455 {
1456         union value a, c, x, t;
1457         float xf;
1458         u_int64_t status = 0;
1459         int i;
1460
1461         for (i = 0; i < ITER; i++) {
1462                 a.r.q = (unsigned long long)random() << 32
1463                         | random();
1464                 status = 0;
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))
1468                         continue;
1469                 xf = (float) a.r.q;
1470                 x.d = xf;
1471                 t.r = ieee_convert_S_T(c.r, ROUND_NORMAL, 0, &status);
1472                 if (t.d != x.d) {
1473                         printf("bad convert single, (double)%qx = %g (should be %g)\n",
1474                                a.r.q, c.d, x.d);
1475                         c.r = ieee_convert_Q_S(a.r, ROUND_NORMAL, 0, &status);
1476                 }
1477         }
1478 }
1479
1480 static void
1481 test_convert_double_to_int()
1482 {
1483         union value a, c;
1484         u_int64_t status = 0;
1485         int i;
1486
1487         for (i = 0; i < ITER; i++) {
1488                 a.d = random_double();
1489                 status = 0;
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);
1495                 }
1496         }
1497 }
1498
1499 int
1500 main(int argc, char* argv[])
1501 {
1502         srandom(0);
1503
1504         test_double_div();
1505         test_single_div();
1506         test_double_add();
1507         test_single_add();
1508         test_double_mul();
1509         test_single_mul();
1510         test_convert_int_to_double();
1511         test_convert_int_to_single();
1512 #if 0
1513         /* x86 generates SIGFPE on overflows. */
1514         test_convert_double_to_int();
1515 #endif
1516
1517         return 0;
1518 }
1519
1520 #endif