libgcc/config/libbid/bid128_next.c

   1 /* Copyright (C) 2007  Free Software Foundation, Inc.
   2
   3 This file is part of GCC.
   4
   5 GCC is free software; you can redistribute it and/or modify it under
   6 the terms of the GNU General Public License as published by the Free
   7 Software Foundation; either version 2, or (at your option) any later
   8 version.
   9
  10 In addition to the permissions in the GNU General Public License, the
  11 Free Software Foundation gives you unlimited permission to link the
  12 compiled version of this file into combinations with other programs,
  13 and to distribute those combinations without any restriction coming
  14 from the use of this file.  (The General Public License restrictions
  15 do apply in other respects; for example, they cover modification of
  16 the file, and distribution when not linked into a combine
  17 executable.)
  18
  19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  21 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  22 for more details.
  23
  24 You should have received a copy of the GNU General Public License
  25 along with GCC; see the file COPYING.  If not, write to the Free
  26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  27 02110-1301, USA.  */
  28
  29 #define BID_128RES
  30 #include "bid_internal.h"
  31
  32 /*****************************************************************************
  33  *  BID128 nextup
  34  ****************************************************************************/
  35
  36 #if DECIMAL_CALL_BY_REFERENCE
  37 void
  38 bid128_nextup (UINT128 * pres,
  39                UINT128 *
  40                px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  41   UINT128 x = *px;
  42 #else
  43 UINT128
  44 bid128_nextup (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  45                _EXC_INFO_PARAM) {
  46 #endif
  47
  48   UINT128 res;
  49   UINT64 x_sign;
  50   UINT64 x_exp;
  51   int exp;
  52   BID_UI64DOUBLE tmp1;
  53   int x_nr_bits;
  54   int q1, ind;
  55   UINT128 C1;                   // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
  56
  57   BID_SWAP128 (x);
  58   // unpack the argument
  59   x_sign = x.w[1] & MASK_SIGN;  // 0 for positive, MASK_SIGN for negative
  60   C1.w[1] = x.w[1] & MASK_COEFF;
  61   C1.w[0] = x.w[0];
  62
  63   // check for NaN or Infinity
  64   if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
  65     // x is special
  66     if ((x.w[1] & MASK_NAN) == MASK_NAN) {      // x is NAN
  67       // if x = NaN, then res = Q (x)
  68       // check first for non-canonical NaN payload
  69       if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
  70           (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
  71            && (x.w[0] > 0x38c15b09ffffffffull))) {
  72         x.w[1] = x.w[1] & 0xffffc00000000000ull;
  73         x.w[0] = 0x0ull;
  74       }
  75       if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {  // x is SNAN
  76         // set invalid flag
  77         *pfpsf |= INVALID_EXCEPTION;
  78         // return quiet (x)
  79         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out also G[6]-G[16]
  80         res.w[0] = x.w[0];
  81       } else {  // x is QNaN
  82         // return x
  83         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out G[6]-G[16]
  84         res.w[0] = x.w[0];
  85       }
  86     } else {    // x is not NaN, so it must be infinity
  87       if (!x_sign) {    // x is +inf
  88         res.w[1] = 0x7800000000000000ull;       // +inf
  89         res.w[0] = 0x0000000000000000ull;
  90       } else {  // x is -inf
  91         res.w[1] = 0xdfffed09bead87c0ull;       // -MAXFP = -999...99 * 10^emax
  92         res.w[0] = 0x378d8e63ffffffffull;
  93       }
  94     }
  95     BID_RETURN (res);
  96   }
  97   // check for non-canonical values (treated as zero)
  98   if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {      // G0_G1=11
  99     // non-canonical
 100     x_exp = (x.w[1] << 2) & MASK_EXP;   // biased and shifted left 49 bits
 101     C1.w[1] = 0;        // significand high
 102     C1.w[0] = 0;        // significand low
 103   } else {      // G0_G1 != 11
 104     x_exp = x.w[1] & MASK_EXP;  // biased and shifted left 49 bits
 105     if (C1.w[1] > 0x0001ed09bead87c0ull ||
 106         (C1.w[1] == 0x0001ed09bead87c0ull
 107          && C1.w[0] > 0x378d8e63ffffffffull)) {
 108       // x is non-canonical if coefficient is larger than 10^34 -1
 109       C1.w[1] = 0;
 110       C1.w[0] = 0;
 111     } else {    // canonical
 112       ;
 113     }
 114   }
 115
 116   if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
 117     // x is +/-0
 118     res.w[1] = 0x0000000000000000ull;   // +1 * 10^emin
 119     res.w[0] = 0x0000000000000001ull;
 120   } else {      // x is not special and is not zero
 121     if (x.w[1] == 0x5fffed09bead87c0ull
 122         && x.w[0] == 0x378d8e63ffffffffull) {
 123       // x = +MAXFP = 999...99 * 10^emax
 124       res.w[1] = 0x7800000000000000ull; // +inf
 125       res.w[0] = 0x0000000000000000ull;
 126     } else if (x.w[1] == 0x8000000000000000ull
 127                && x.w[0] == 0x0000000000000001ull) {
 128       // x = -MINFP = 1...99 * 10^emin
 129       res.w[1] = 0x8000000000000000ull; // -0
 130       res.w[0] = 0x0000000000000000ull;
 131     } else {    // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
 132       // can add/subtract 1 ulp to the significand
 133
 134       // Note: we could check here if x >= 10^34 to speed up the case q1 = 34
 135       // q1 = nr. of decimal digits in x
 136       // determine first the nr. of bits in x
 137       if (C1.w[1] == 0) {
 138         if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
 139           // split the 64-bit value in two 32-bit halves to avoid rnd errors
 140           if (C1.w[0] >= 0x0000000100000000ull) {       // x >= 2^32
 141             tmp1.d = (double) (C1.w[0] >> 32);  // exact conversion
 142             x_nr_bits =
 143               33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
 144                     0x3ff);
 145           } else {      // x < 2^32
 146             tmp1.d = (double) (C1.w[0]);        // exact conversion
 147             x_nr_bits =
 148               1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
 149                    0x3ff);
 150           }
 151         } else {        // if x < 2^53
 152           tmp1.d = (double) C1.w[0];    // exact conversion
 153           x_nr_bits =
 154             1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 155         }
 156       } else {  // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 157         tmp1.d = (double) C1.w[1];      // exact conversion
 158         x_nr_bits =
 159           65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 160       }
 161       q1 = nr_digits[x_nr_bits - 1].digits;
 162       if (q1 == 0) {
 163         q1 = nr_digits[x_nr_bits - 1].digits1;
 164         if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 165             || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
 166                 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 167           q1++;
 168       }
 169       // if q1 < P34 then pad the significand with zeros
 170       if (q1 < P34) {
 171         exp = (x_exp >> 49) - 6176;
 172         if (exp + 6176 > P34 - q1) {
 173           ind = P34 - q1;       // 1 <= ind <= P34 - 1
 174           // pad with P34 - q1 zeros, until exponent = emin
 175           // C1 = C1 * 10^ind
 176           if (q1 <= 19) {       // 64-bit C1
 177             if (ind <= 19) {    // 64-bit 10^ind and 64-bit C1
 178               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
 179             } else {    // 128-bit 10^ind and 64-bit C1
 180               __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
 181             }
 182           } else {      // C1 is (most likely) 128-bit
 183             if (ind <= 14) {    // 64-bit 10^ind and 128-bit C1 (most likely)
 184               __mul_128x64_to_128 (C1, ten2k64[ind], C1);
 185             } else if (ind <= 19) {     // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
 186               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
 187             } else {    // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
 188               __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
 189             }
 190           }
 191           x_exp = x_exp - ((UINT64) ind << 49);
 192         } else {        // pad with zeros until the exponent reaches emin
 193           ind = exp + 6176;
 194           // C1 = C1 * 10^ind
 195           if (ind <= 19) {      // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
 196             if (q1 <= 19) {     // 64-bit C1, 64-bit 10^ind
 197               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
 198             } else {    // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
 199               __mul_128x64_to_128 (C1, ten2k64[ind], C1);
 200             }
 201           } else {      // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
 202             // 64-bit C1, 128-bit 10^ind
 203             __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
 204           }
 205           x_exp = EXP_MIN;
 206         }
 207       }
 208       if (!x_sign) {    // x > 0
 209         // add 1 ulp (add 1 to the significand)
 210         C1.w[0]++;
 211         if (C1.w[0] == 0)
 212           C1.w[1]++;
 213         if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) {     // if  C1 = 10^34
 214           C1.w[1] = 0x0000314dc6448d93ull;      // C1 = 10^33
 215           C1.w[0] = 0x38c15b0a00000000ull;
 216           x_exp = x_exp + EXP_P1;
 217         }
 218       } else {  // x < 0
 219         // subtract 1 ulp (subtract 1 from the significand)
 220         C1.w[0]--;
 221         if (C1.w[0] == 0xffffffffffffffffull)
 222           C1.w[1]--;
 223         if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) {       // if  C1 = 10^33 - 1
 224           C1.w[1] = 0x0001ed09bead87c0ull;      // C1 = 10^34 - 1
 225           C1.w[0] = 0x378d8e63ffffffffull;
 226           x_exp = x_exp - EXP_P1;
 227         }
 228       }
 229       // assemble the result
 230       res.w[1] = x_sign | x_exp | C1.w[1];
 231       res.w[0] = C1.w[0];
 232     }   // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
 233   }     // end x is not special and is not zero
 234   BID_RETURN (res);
 235 }
 236
 237 /*****************************************************************************
 238  *  BID128 nextdown
 239  ****************************************************************************/
 240
 241 #if DECIMAL_CALL_BY_REFERENCE
 242 void
 243 bid128_nextdown (UINT128 * pres,
 244                  UINT128 *
 245                  px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
 246   UINT128 x = *px;
 247 #else
 248 UINT128
 249 bid128_nextdown (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 250                  _EXC_INFO_PARAM) {
 251 #endif
 252
 253   UINT128 res;
 254   UINT64 x_sign;
 255   UINT64 x_exp;
 256   int exp;
 257   BID_UI64DOUBLE tmp1;
 258   int x_nr_bits;
 259   int q1, ind;
 260   UINT128 C1;                   // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
 261
 262   BID_SWAP128 (x);
 263   // unpack the argument
 264   x_sign = x.w[1] & MASK_SIGN;  // 0 for positive, MASK_SIGN for negative
 265   C1.w[1] = x.w[1] & MASK_COEFF;
 266   C1.w[0] = x.w[0];
 267
 268   // check for NaN or Infinity
 269   if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
 270     // x is special
 271     if ((x.w[1] & MASK_NAN) == MASK_NAN) {      // x is NAN
 272       // if x = NaN, then res = Q (x)
 273       // check first for non-canonical NaN payload
 274       if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
 275           (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
 276            && (x.w[0] > 0x38c15b09ffffffffull))) {
 277         x.w[1] = x.w[1] & 0xffffc00000000000ull;
 278         x.w[0] = 0x0ull;
 279       }
 280       if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {  // x is SNAN
 281         // set invalid flag
 282         *pfpsf |= INVALID_EXCEPTION;
 283         // return quiet (x)
 284         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out also G[6]-G[16]
 285         res.w[0] = x.w[0];
 286       } else {  // x is QNaN
 287         // return x
 288         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out G[6]-G[16]
 289         res.w[0] = x.w[0];
 290       }
 291     } else {    // x is not NaN, so it must be infinity
 292       if (!x_sign) {    // x is +inf
 293         res.w[1] = 0x5fffed09bead87c0ull;       // +MAXFP = +999...99 * 10^emax
 294         res.w[0] = 0x378d8e63ffffffffull;
 295       } else {  // x is -inf
 296         res.w[1] = 0xf800000000000000ull;       // -inf
 297         res.w[0] = 0x0000000000000000ull;
 298       }
 299     }
 300     BID_RETURN (res);
 301   }
 302   // check for non-canonical values (treated as zero)
 303   if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {      // G0_G1=11
 304     // non-canonical
 305     x_exp = (x.w[1] << 2) & MASK_EXP;   // biased and shifted left 49 bits
 306     C1.w[1] = 0;        // significand high
 307     C1.w[0] = 0;        // significand low
 308   } else {      // G0_G1 != 11
 309     x_exp = x.w[1] & MASK_EXP;  // biased and shifted left 49 bits
 310     if (C1.w[1] > 0x0001ed09bead87c0ull ||
 311         (C1.w[1] == 0x0001ed09bead87c0ull
 312          && C1.w[0] > 0x378d8e63ffffffffull)) {
 313       // x is non-canonical if coefficient is larger than 10^34 -1
 314       C1.w[1] = 0;
 315       C1.w[0] = 0;
 316     } else {    // canonical
 317       ;
 318     }
 319   }
 320
 321   if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
 322     // x is +/-0
 323     res.w[1] = 0x8000000000000000ull;   // -1 * 10^emin
 324     res.w[0] = 0x0000000000000001ull;
 325   } else {      // x is not special and is not zero
 326     if (x.w[1] == 0xdfffed09bead87c0ull
 327         && x.w[0] == 0x378d8e63ffffffffull) {
 328       // x = -MAXFP = -999...99 * 10^emax
 329       res.w[1] = 0xf800000000000000ull; // -inf
 330       res.w[0] = 0x0000000000000000ull;
 331     } else if (x.w[1] == 0x0ull && x.w[0] == 0x0000000000000001ull) {   // +MINFP
 332       res.w[1] = 0x0000000000000000ull; // +0
 333       res.w[0] = 0x0000000000000000ull;
 334     } else {    // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
 335       // can add/subtract 1 ulp to the significand
 336
 337       // Note: we could check here if x >= 10^34 to speed up the case q1 = 34
 338       // q1 = nr. of decimal digits in x
 339       // determine first the nr. of bits in x
 340       if (C1.w[1] == 0) {
 341         if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
 342           // split the 64-bit value in two 32-bit halves to avoid rnd errors
 343           if (C1.w[0] >= 0x0000000100000000ull) {       // x >= 2^32
 344             tmp1.d = (double) (C1.w[0] >> 32);  // exact conversion
 345             x_nr_bits =
 346               33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
 347                     0x3ff);
 348           } else {      // x < 2^32
 349             tmp1.d = (double) (C1.w[0]);        // exact conversion
 350             x_nr_bits =
 351               1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
 352                    0x3ff);
 353           }
 354         } else {        // if x < 2^53
 355           tmp1.d = (double) C1.w[0];    // exact conversion
 356           x_nr_bits =
 357             1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 358         }
 359       } else {  // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
 360         tmp1.d = (double) C1.w[1];      // exact conversion
 361         x_nr_bits =
 362           65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
 363       }
 364       q1 = nr_digits[x_nr_bits - 1].digits;
 365       if (q1 == 0) {
 366         q1 = nr_digits[x_nr_bits - 1].digits1;
 367         if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
 368             || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
 369                 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
 370           q1++;
 371       }
 372       // if q1 < P then pad the significand with zeros
 373       if (q1 < P34) {
 374         exp = (x_exp >> 49) - 6176;
 375         if (exp + 6176 > P34 - q1) {
 376           ind = P34 - q1;       // 1 <= ind <= P34 - 1
 377           // pad with P34 - q1 zeros, until exponent = emin
 378           // C1 = C1 * 10^ind
 379           if (q1 <= 19) {       // 64-bit C1
 380             if (ind <= 19) {    // 64-bit 10^ind and 64-bit C1
 381               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
 382             } else {    // 128-bit 10^ind and 64-bit C1
 383               __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
 384             }
 385           } else {      // C1 is (most likely) 128-bit
 386             if (ind <= 14) {    // 64-bit 10^ind and 128-bit C1 (most likely)
 387               __mul_128x64_to_128 (C1, ten2k64[ind], C1);
 388             } else if (ind <= 19) {     // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
 389               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
 390             } else {    // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
 391               __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
 392             }
 393           }
 394           x_exp = x_exp - ((UINT64) ind << 49);
 395         } else {        // pad with zeros until the exponent reaches emin
 396           ind = exp + 6176;
 397           // C1 = C1 * 10^ind
 398           if (ind <= 19) {      // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
 399             if (q1 <= 19) {     // 64-bit C1, 64-bit 10^ind
 400               __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
 401             } else {    // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
 402               __mul_128x64_to_128 (C1, ten2k64[ind], C1);
 403             }
 404           } else {      // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
 405             // 64-bit C1, 128-bit 10^ind
 406             __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
 407           }
 408           x_exp = EXP_MIN;
 409         }
 410       }
 411       if (x_sign) {     // x < 0
 412         // add 1 ulp (add 1 to the significand)
 413         C1.w[0]++;
 414         if (C1.w[0] == 0)
 415           C1.w[1]++;
 416         if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) {     // if  C1 = 10^34
 417           C1.w[1] = 0x0000314dc6448d93ull;      // C1 = 10^33
 418           C1.w[0] = 0x38c15b0a00000000ull;
 419           x_exp = x_exp + EXP_P1;
 420         }
 421       } else {  // x > 0
 422         // subtract 1 ulp (subtract 1 from the significand)
 423         C1.w[0]--;
 424         if (C1.w[0] == 0xffffffffffffffffull)
 425           C1.w[1]--;
 426         if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) {       // if  C1 = 10^33 - 1
 427           C1.w[1] = 0x0001ed09bead87c0ull;      // C1 = 10^34 - 1
 428           C1.w[0] = 0x378d8e63ffffffffull;
 429           x_exp = x_exp - EXP_P1;
 430         }
 431       }
 432       // assemble the result
 433       res.w[1] = x_sign | x_exp | C1.w[1];
 434       res.w[0] = C1.w[0];
 435     }   // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
 436   }     // end x is not special and is not zero
 437   BID_RETURN (res);
 438 }
 439
 440 /*****************************************************************************
 441  *  BID128 nextafter
 442  ****************************************************************************/
 443
 444 #if DECIMAL_CALL_BY_REFERENCE
 445 void
 446 bid128_nextafter (UINT128 * pres, UINT128 * px,
 447                   UINT128 *
 448                   py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
 449 {
 450   UINT128 x = *px;
 451   UINT128 y = *py;
 452   UINT128 xnswp = *px;
 453   UINT128 ynswp = *py;
 454 #else
 455 UINT128
 456 bid128_nextafter (UINT128 x,
 457                   UINT128 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
 458                   _EXC_INFO_PARAM) {
 459   UINT128 xnswp = x;
 460   UINT128 ynswp = y;
 461 #endif
 462
 463   UINT128 res;
 464   UINT128 tmp1, tmp2, tmp3;
 465   FPSC tmp_fpsf = 0;            // dummy fpsf for calls to comparison functions
 466   int res1, res2;
 467   UINT64 x_exp;
 468
 469
 470   BID_SWAP128 (x);
 471   BID_SWAP128 (y);
 472   // check for NaNs
 473   if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL)
 474       || ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
 475     // x is special or y is special
 476     if ((x.w[1] & MASK_NAN) == MASK_NAN) {      // x is NAN
 477       // if x = NaN, then res = Q (x)
 478       // check first for non-canonical NaN payload
 479       if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
 480           (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
 481            && (x.w[0] > 0x38c15b09ffffffffull))) {
 482         x.w[1] = x.w[1] & 0xffffc00000000000ull;
 483         x.w[0] = 0x0ull;
 484       }
 485       if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {  // x is SNAN
 486         // set invalid flag
 487         *pfpsf |= INVALID_EXCEPTION;
 488         // return quiet (x)
 489         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out also G[6]-G[16]
 490         res.w[0] = x.w[0];
 491       } else {  // x is QNaN
 492         // return x
 493         res.w[1] = x.w[1] & 0xfc003fffffffffffull;      // clear out G[6]-G[16]
 494         res.w[0] = x.w[0];
 495         if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {        // y is SNAN
 496           // set invalid flag
 497           *pfpsf |= INVALID_EXCEPTION;
 498         }
 499       }
 500       BID_RETURN (res)
 501     } else if ((y.w[1] & MASK_NAN) == MASK_NAN) {       // y is NAN
 502       // if x = NaN, then res = Q (x)
 503       // check first for non-canonical NaN payload
 504       if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
 505           (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
 506            && (y.w[0] > 0x38c15b09ffffffffull))) {
 507         y.w[1] = y.w[1] & 0xffffc00000000000ull;
 508         y.w[0] = 0x0ull;
 509       }
 510       if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {  // y is SNAN
 511         // set invalid flag
 512         *pfpsf |= INVALID_EXCEPTION;
 513         // return quiet (x)
 514         res.w[1] = y.w[1] & 0xfc003fffffffffffull;      // clear out also G[6]-G[16]
 515         res.w[0] = y.w[0];
 516       } else {  // x is QNaN
 517         // return x
 518         res.w[1] = y.w[1] & 0xfc003fffffffffffull;      // clear out G[6]-G[16]
 519         res.w[0] = y.w[0];
 520       }
 521       BID_RETURN (res)
 522     } else {    // at least one is infinity
 523       if ((x.w[1] & MASK_ANY_INF) == MASK_INF) {        // x = inf
 524         x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
 525         x.w[0] = 0x0ull;
 526       }
 527       if ((y.w[1] & MASK_ANY_INF) == MASK_INF) {        // y = inf
 528         y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
 529         y.w[0] = 0x0ull;
 530       }
 531     }
 532   }
 533   // neither x nor y is NaN
 534
 535   // if not infinity, check for non-canonical values x (treated as zero)
 536   if ((x.w[1] & MASK_ANY_INF) != MASK_INF) {    // x != inf
 537     if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {    // G0_G1=11
 538       // non-canonical
 539       x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
 540       x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
 541       x.w[0] = 0x0ull;
 542     } else {    // G0_G1 != 11
 543       x_exp = x.w[1] & MASK_EXP;        // biased and shifted left 49 bits
 544       if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
 545           ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
 546            && x.w[0] > 0x378d8e63ffffffffull)) {
 547         // x is non-canonical if coefficient is larger than 10^34 -1
 548         x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
 549         x.w[0] = 0x0ull;
 550       } else {  // canonical
 551         ;
 552       }
 553     }
 554   }
 555   // no need to check for non-canonical y
 556
 557   // neither x nor y is NaN
 558   tmp_fpsf = *pfpsf;    // save fpsf
 559 #if DECIMAL_CALL_BY_REFERENCE
 560   bid128_quiet_equal (&res1, &xnswp,
 561                       &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
 562                       _EXC_INFO_ARG);
 563   bid128_quiet_greater (&res2, &xnswp,
 564                         &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
 565                         _EXC_INFO_ARG);
 566 #else
 567   res1 =
 568     bid128_quiet_equal (xnswp,
 569                         ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
 570                         _EXC_INFO_ARG);
 571   res2 =
 572     bid128_quiet_greater (xnswp,
 573                           ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
 574                           _EXC_INFO_ARG);
 575 #endif
 576   *pfpsf = tmp_fpsf;    // restore fpsf
 577
 578   if (res1) {   // x = y
 579     // return x with the sign of y
 580     res.w[1] =
 581       (x.w[1] & 0x7fffffffffffffffull) | (y.
 582                                           w[1] & 0x8000000000000000ull);
 583     res.w[0] = x.w[0];
 584   } else if (res2) {    // x > y
 585 #if DECIMAL_CALL_BY_REFERENCE
 586     bid128_nextdown (&res,
 587                      &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
 588                      _EXC_INFO_ARG);
 589 #else
 590     res =
 591       bid128_nextdown (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
 592                        _EXC_INFO_ARG);
 593 #endif
 594     BID_SWAP128 (res);
 595   } else {      // x < y
 596 #if DECIMAL_CALL_BY_REFERENCE
 597     bid128_nextup (&res,
 598                    &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 599 #else
 600     res =
 601       bid128_nextup (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
 602 #endif
 603     BID_SWAP128 (res);
 604   }
 605   // if the operand x is finite but the result is infinite, signal
 606   // overflow and inexact
 607   if (((x.w[1] & MASK_SPECIAL) != MASK_SPECIAL)
 608       && ((res.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
 609     // set the inexact flag
 610     *pfpsf |= INEXACT_EXCEPTION;
 611     // set the overflow flag
 612     *pfpsf |= OVERFLOW_EXCEPTION;
 613   }
 614   // if the result is in (-10^emin, 10^emin), and is different from the
 615   // operand x, signal underflow and inexact
 616   tmp1.w[HIGH_128W] = 0x0000314dc6448d93ull;
 617   tmp1.w[LOW_128W] = 0x38c15b0a00000000ull;     // +100...0[34] * 10^emin
 618   tmp2.w[HIGH_128W] = res.w[1] & 0x7fffffffffffffffull;
 619   tmp2.w[LOW_128W] = res.w[0];
 620   tmp3.w[HIGH_128W] = res.w[1];
 621   tmp3.w[LOW_128W] = res.w[0];
 622   tmp_fpsf = *pfpsf;    // save fpsf
 623 #if DECIMAL_CALL_BY_REFERENCE
 624   bid128_quiet_greater (&res1, &tmp1,
 625                         &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
 626                         _EXC_INFO_ARG);
 627   bid128_quiet_not_equal (&res2, &xnswp,
 628                           &tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
 629                           _EXC_INFO_ARG);
 630 #else
 631   res1 =
 632     bid128_quiet_greater (tmp1,
 633                           tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
 634                           _EXC_INFO_ARG);
 635   res2 =
 636     bid128_quiet_not_equal (xnswp,
 637                             tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
 638                             _EXC_INFO_ARG);
 639 #endif
 640   *pfpsf = tmp_fpsf;    // restore fpsf
 641   if (res1 && res2) {
 642     // set the inexact flag
 643     *pfpsf |= INEXACT_EXCEPTION;
 644     // set the underflow flag
 645     *pfpsf |= UNDERFLOW_EXCEPTION;
 646   }
 647   BID_RETURN (res);
 648 }