libgcc/config/libbid/bid64_minmax.c

   1 /* Copyright (C) 2007-2019 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 3, or (at your option) any later
   8 version.
   9
  10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  12 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  13 for more details.
  14
  15 Under Section 7 of GPL version 3, you are granted additional
  16 permissions described in the GCC Runtime Library Exception, version
  17 3.1, as published by the Free Software Foundation.
  18
  19 You should have received a copy of the GNU General Public License and
  20 a copy of the GCC Runtime Library Exception along with this program;
  21 see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
  22 <http://www.gnu.org/licenses/>.  */
  23
  24 #include "bid_internal.h"
  25
  26 /*****************************************************************************
  27  *  BID64 minimum function - returns greater of two numbers
  28  *****************************************************************************/
  29
  30 static const UINT64 mult_factor[16] = {
  31   1ull, 10ull, 100ull, 1000ull,
  32   10000ull, 100000ull, 1000000ull, 10000000ull,
  33   100000000ull, 1000000000ull, 10000000000ull, 100000000000ull,
  34   1000000000000ull, 10000000000000ull,
  35   100000000000000ull, 1000000000000000ull
  36 };
  37
  38 #if DECIMAL_CALL_BY_REFERENCE
  39 void
  40 bid64_minnum (UINT64 * pres, UINT64 * px, UINT64 * py _EXC_FLAGS_PARAM) {
  41   UINT64 x = *px;
  42   UINT64 y = *py;
  43 #else
  44 UINT64
  45 bid64_minnum (UINT64 x, UINT64 y _EXC_FLAGS_PARAM) {
  46 #endif
  47
  48   UINT64 res;
  49   int exp_x, exp_y;
  50   UINT64 sig_x, sig_y;
  51   UINT128 sig_n_prime;
  52   char x_is_zero = 0, y_is_zero = 0;
  53
  54   // check for non-canonical x
  55   if ((x & MASK_NAN) == MASK_NAN) {     // x is NaN
  56     x = x & 0xfe03ffffffffffffull;      // clear G6-G12
  57     if ((x & 0x0003ffffffffffffull) > 999999999999999ull) {
  58       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
  59     }
  60   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
  61     x = x & (MASK_SIGN | MASK_INF);
  62   } else {      // x is not special
  63     // check for non-canonical values - treated as zero
  64     if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
  65       // if the steering bits are 11, then the exponent is G[0:w+1]
  66       if (((x & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
  67           9999999999999999ull) {
  68         // non-canonical
  69         x = (x & MASK_SIGN) | ((x & MASK_BINARY_EXPONENT2) << 2);
  70       } // else canonical
  71     }   // else canonical
  72   }
  73
  74   // check for non-canonical y
  75   if ((y & MASK_NAN) == MASK_NAN) {     // y is NaN
  76     y = y & 0xfe03ffffffffffffull;      // clear G6-G12
  77     if ((y & 0x0003ffffffffffffull) > 999999999999999ull) {
  78       y = y & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
  79     }
  80   } else if ((y & MASK_INF) == MASK_INF) {      // check for Infinity
  81     y = y & (MASK_SIGN | MASK_INF);
  82   } else {      // y is not special
  83     // check for non-canonical values - treated as zero
  84     if ((y & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
  85       // if the steering bits are 11, then the exponent is G[0:w+1]
  86       if (((y & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
  87           9999999999999999ull) {
  88         // non-canonical
  89         y = (y & MASK_SIGN) | ((y & MASK_BINARY_EXPONENT2) << 2);
  90       } // else canonical
  91     }   // else canonical
  92   }
  93
  94   // NaN (CASE1)
  95   if ((x & MASK_NAN) == MASK_NAN) {     // x is NAN
  96     if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNaN
  97       // if x is SNAN, then return quiet (x)
  98       *pfpsf |= INVALID_EXCEPTION;      // set exception if SNaN
  99       x = x & 0xfdffffffffffffffull;    // quietize x
 100       res = x;
 101     } else {    // x is QNaN
 102       if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
 103         if ((y & MASK_SNAN) == MASK_SNAN) {     // y is SNAN
 104           *pfpsf |= INVALID_EXCEPTION;  // set invalid flag
 105         }
 106         res = x;
 107       } else {
 108         res = y;
 109       }
 110     }
 111     BID_RETURN (res);
 112   } else if ((y & MASK_NAN) == MASK_NAN) {      // y is NaN, but x is not
 113     if ((y & MASK_SNAN) == MASK_SNAN) {
 114       *pfpsf |= INVALID_EXCEPTION;      // set exception if SNaN
 115       y = y & 0xfdffffffffffffffull;    // quietize y
 116       res = y;
 117     } else {
 118       // will return x (which is not NaN)
 119       res = x;
 120     }
 121     BID_RETURN (res);
 122   }
 123   // SIMPLE (CASE2)
 124   // if all the bits are the same, these numbers are equal, return either number
 125   if (x == y) {
 126     res = x;
 127     BID_RETURN (res);
 128   }
 129   // INFINITY (CASE3)
 130   if ((x & MASK_INF) == MASK_INF) {
 131     // if x is neg infinity, there is no way it is greater than y, return x
 132     if (((x & MASK_SIGN) == MASK_SIGN)) {
 133       res = x;
 134       BID_RETURN (res);
 135     }
 136     // x is pos infinity, return y
 137     else {
 138       res = y;
 139       BID_RETURN (res);
 140     }
 141   } else if ((y & MASK_INF) == MASK_INF) {
 142     // x is finite, so if y is positive infinity, then x is less, return y
 143     //                 if y is negative infinity, then x is greater, return x
 144     res = ((y & MASK_SIGN) == MASK_SIGN) ? y : x;
 145     BID_RETURN (res);
 146   }
 147   // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
 148   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 149     exp_x = (x & MASK_BINARY_EXPONENT2) >> 51;
 150     sig_x = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 151   } else {
 152     exp_x = (x & MASK_BINARY_EXPONENT1) >> 53;
 153     sig_x = (x & MASK_BINARY_SIG1);
 154   }
 155
 156   // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
 157   if ((y & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 158     exp_y = (y & MASK_BINARY_EXPONENT2) >> 51;
 159     sig_y = (y & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 160   } else {
 161     exp_y = (y & MASK_BINARY_EXPONENT1) >> 53;
 162     sig_y = (y & MASK_BINARY_SIG1);
 163   }
 164
 165   // ZERO (CASE4)
 166   // some properties:
 167   //    (+ZERO == -ZERO) => therefore
 168   //        ignore the sign, and neither number is greater
 169   //    (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
 170   //        ignore the exponent field
 171   //    (Any non-canonical # is considered 0)
 172   if (sig_x == 0) {
 173     x_is_zero = 1;
 174   }
 175   if (sig_y == 0) {
 176     y_is_zero = 1;
 177   }
 178
 179   if (x_is_zero && y_is_zero) {
 180     // if both numbers are zero, neither is greater => return either
 181     res = y;
 182     BID_RETURN (res);
 183   } else if (x_is_zero) {
 184     // is x is zero, it is greater if Y is negative
 185     res = ((y & MASK_SIGN) == MASK_SIGN) ? y : x;
 186     BID_RETURN (res);
 187   } else if (y_is_zero) {
 188     // is y is zero, X is greater if it is positive
 189     res = ((x & MASK_SIGN) != MASK_SIGN) ? y : x;;
 190     BID_RETURN (res);
 191   }
 192   // OPPOSITE SIGN (CASE5)
 193   // now, if the sign bits differ, x is greater if y is negative
 194   if (((x ^ y) & MASK_SIGN) == MASK_SIGN) {
 195     res = ((y & MASK_SIGN) == MASK_SIGN) ? y : x;
 196     BID_RETURN (res);
 197   }
 198   // REDUNDANT REPRESENTATIONS (CASE6)
 199
 200   // if both components are either bigger or smaller,
 201   // it is clear what needs to be done
 202   if (sig_x > sig_y && exp_x >= exp_y) {
 203     res = ((x & MASK_SIGN) != MASK_SIGN) ? y : x;
 204     BID_RETURN (res);
 205   }
 206   if (sig_x < sig_y && exp_x <= exp_y) {
 207     res = ((x & MASK_SIGN) == MASK_SIGN) ? y : x;
 208     BID_RETURN (res);
 209   }
 210   // if exp_x is 15 greater than exp_y, no need for compensation
 211   if (exp_x - exp_y > 15) {
 212     res = ((x & MASK_SIGN) != MASK_SIGN) ? y : x;       // difference cannot be >10^15
 213     BID_RETURN (res);
 214   }
 215   // if exp_x is 15 less than exp_y, no need for compensation
 216   if (exp_y - exp_x > 15) {
 217     res = ((x & MASK_SIGN) == MASK_SIGN) ? y : x;
 218     BID_RETURN (res);
 219   }
 220   // if |exp_x - exp_y| < 15, it comes down to the compensated significand
 221   if (exp_x > exp_y) {  // to simplify the loop below,
 222
 223     // otherwise adjust the x significand upwards
 224     __mul_64x64_to_128MACH (sig_n_prime, sig_x,
 225                             mult_factor[exp_x - exp_y]);
 226     // if postitive, return whichever significand is larger
 227     // (converse if negative)
 228     if (sig_n_prime.w[1] == 0 && (sig_n_prime.w[0] == sig_y)) {
 229       res = y;
 230       BID_RETURN (res);
 231     }
 232
 233     res = (((sig_n_prime.w[1] > 0)
 234             || sig_n_prime.w[0] > sig_y) ^ ((x & MASK_SIGN) ==
 235                                             MASK_SIGN)) ? y : x;
 236     BID_RETURN (res);
 237   }
 238   // adjust the y significand upwards
 239   __mul_64x64_to_128MACH (sig_n_prime, sig_y,
 240                           mult_factor[exp_y - exp_x]);
 241
 242   // if postitive, return whichever significand is larger (converse if negative)
 243   if (sig_n_prime.w[1] == 0 && (sig_n_prime.w[0] == sig_x)) {
 244     res = y;
 245     BID_RETURN (res);
 246   }
 247   res = (((sig_n_prime.w[1] == 0)
 248           && (sig_x > sig_n_prime.w[0])) ^ ((x & MASK_SIGN) ==
 249                                             MASK_SIGN)) ? y : x;
 250   BID_RETURN (res);
 251 }
 252
 253 /*****************************************************************************
 254  *  BID64 minimum magnitude function - returns greater of two numbers
 255  *****************************************************************************/
 256
 257 #if DECIMAL_CALL_BY_REFERENCE
 258 void
 259 bid64_minnum_mag (UINT64 * pres, UINT64 * px,
 260                   UINT64 * py _EXC_FLAGS_PARAM) {
 261   UINT64 x = *px;
 262   UINT64 y = *py;
 263 #else
 264 UINT64
 265 bid64_minnum_mag (UINT64 x, UINT64 y _EXC_FLAGS_PARAM) {
 266 #endif
 267
 268   UINT64 res;
 269   int exp_x, exp_y;
 270   UINT64 sig_x, sig_y;
 271   UINT128 sig_n_prime;
 272
 273   // check for non-canonical x
 274   if ((x & MASK_NAN) == MASK_NAN) {     // x is NaN
 275     x = x & 0xfe03ffffffffffffull;      // clear G6-G12
 276     if ((x & 0x0003ffffffffffffull) > 999999999999999ull) {
 277       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 278     }
 279   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
 280     x = x & (MASK_SIGN | MASK_INF);
 281   } else {      // x is not special
 282     // check for non-canonical values - treated as zero
 283     if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 284       // if the steering bits are 11, then the exponent is G[0:w+1]
 285       if (((x & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
 286           9999999999999999ull) {
 287         // non-canonical
 288         x = (x & MASK_SIGN) | ((x & MASK_BINARY_EXPONENT2) << 2);
 289       } // else canonical
 290     }   // else canonical
 291   }
 292
 293   // check for non-canonical y
 294   if ((y & MASK_NAN) == MASK_NAN) {     // y is NaN
 295     y = y & 0xfe03ffffffffffffull;      // clear G6-G12
 296     if ((y & 0x0003ffffffffffffull) > 999999999999999ull) {
 297       y = y & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 298     }
 299   } else if ((y & MASK_INF) == MASK_INF) {      // check for Infinity
 300     y = y & (MASK_SIGN | MASK_INF);
 301   } else {      // y is not special
 302     // check for non-canonical values - treated as zero
 303     if ((y & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 304       // if the steering bits are 11, then the exponent is G[0:w+1]
 305       if (((y & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
 306           9999999999999999ull) {
 307         // non-canonical
 308         y = (y & MASK_SIGN) | ((y & MASK_BINARY_EXPONENT2) << 2);
 309       } // else canonical
 310     }   // else canonical
 311   }
 312
 313   // NaN (CASE1)
 314   if ((x & MASK_NAN) == MASK_NAN) {     // x is NAN
 315     if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNaN
 316       // if x is SNAN, then return quiet (x)
 317       *pfpsf |= INVALID_EXCEPTION;      // set exception if SNaN
 318       x = x & 0xfdffffffffffffffull;    // quietize x
 319       res = x;
 320     } else {    // x is QNaN
 321       if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
 322         if ((y & MASK_SNAN) == MASK_SNAN) {     // y is SNAN
 323           *pfpsf |= INVALID_EXCEPTION;  // set invalid flag
 324         }
 325         res = x;
 326       } else {
 327         res = y;
 328       }
 329     }
 330     BID_RETURN (res);
 331   } else if ((y & MASK_NAN) == MASK_NAN) {      // y is NaN, but x is not
 332     if ((y & MASK_SNAN) == MASK_SNAN) {
 333       *pfpsf |= INVALID_EXCEPTION;      // set exception if SNaN
 334       y = y & 0xfdffffffffffffffull;    // quietize y
 335       res = y;
 336     } else {
 337       // will return x (which is not NaN)
 338       res = x;
 339     }
 340     BID_RETURN (res);
 341   }
 342   // SIMPLE (CASE2)
 343   // if all the bits are the same, these numbers are equal, return either number
 344   if (x == y) {
 345     res = x;
 346     BID_RETURN (res);
 347   }
 348   // INFINITY (CASE3)
 349   if ((x & MASK_INF) == MASK_INF) {
 350     // x is infinity, its magnitude is greater than or equal to y
 351     // return x only if y is infinity and x is negative
 352     res = ((x & MASK_SIGN) == MASK_SIGN
 353            && (y & MASK_INF) == MASK_INF) ? x : y;
 354     BID_RETURN (res);
 355   } else if ((y & MASK_INF) == MASK_INF) {
 356     // y is infinity, then it must be greater in magnitude, return x
 357     res = x;
 358     BID_RETURN (res);
 359   }
 360   // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
 361   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 362     exp_x = (x & MASK_BINARY_EXPONENT2) >> 51;
 363     sig_x = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 364   } else {
 365     exp_x = (x & MASK_BINARY_EXPONENT1) >> 53;
 366     sig_x = (x & MASK_BINARY_SIG1);
 367   }
 368
 369   // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
 370   if ((y & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 371     exp_y = (y & MASK_BINARY_EXPONENT2) >> 51;
 372     sig_y = (y & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 373   } else {
 374     exp_y = (y & MASK_BINARY_EXPONENT1) >> 53;
 375     sig_y = (y & MASK_BINARY_SIG1);
 376   }
 377
 378   // ZERO (CASE4)
 379   // some properties:
 380   //    (+ZERO == -ZERO) => therefore
 381   //        ignore the sign, and neither number is greater
 382   //    (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
 383   //        ignore the exponent field
 384   //    (Any non-canonical # is considered 0)
 385   if (sig_x == 0) {
 386     res = x;    // x_is_zero, its magnitude must be smaller than y
 387     BID_RETURN (res);
 388   }
 389   if (sig_y == 0) {
 390     res = y;    // y_is_zero, its magnitude must be smaller than x
 391     BID_RETURN (res);
 392   }
 393   // REDUNDANT REPRESENTATIONS (CASE6)
 394   // if both components are either bigger or smaller,
 395   // it is clear what needs to be done
 396   if (sig_x > sig_y && exp_x >= exp_y) {
 397     res = y;
 398     BID_RETURN (res);
 399   }
 400   if (sig_x < sig_y && exp_x <= exp_y) {
 401     res = x;
 402     BID_RETURN (res);
 403   }
 404   // if exp_x is 15 greater than exp_y, no need for compensation
 405   if (exp_x - exp_y > 15) {
 406     res = y;    // difference cannot be greater than 10^15
 407     BID_RETURN (res);
 408   }
 409   // if exp_x is 15 less than exp_y, no need for compensation
 410   if (exp_y - exp_x > 15) {
 411     res = x;
 412     BID_RETURN (res);
 413   }
 414   // if |exp_x - exp_y| < 15, it comes down to the compensated significand
 415   if (exp_x > exp_y) {  // to simplify the loop below,
 416     // otherwise adjust the x significand upwards
 417     __mul_64x64_to_128MACH (sig_n_prime, sig_x,
 418                             mult_factor[exp_x - exp_y]);
 419     // now, sig_n_prime has: sig_x * 10^(exp_x-exp_y), this is
 420     // the compensated signif.
 421     if (sig_n_prime.w[1] == 0 && (sig_n_prime.w[0] == sig_y)) {
 422       // two numbers are equal, return minNum(x,y)
 423       res = ((y & MASK_SIGN) == MASK_SIGN) ? y : x;
 424       BID_RETURN (res);
 425     }
 426     // now, if compensated_x (sig_n_prime) is greater than y, return y,
 427     // otherwise return x
 428     res = ((sig_n_prime.w[1] != 0) || sig_n_prime.w[0] > sig_y) ? y : x;
 429     BID_RETURN (res);
 430   }
 431   // exp_y must be greater than exp_x, thus adjust the y significand upwards
 432   __mul_64x64_to_128MACH (sig_n_prime, sig_y,
 433                           mult_factor[exp_y - exp_x]);
 434
 435   if (sig_n_prime.w[1] == 0 && (sig_n_prime.w[0] == sig_x)) {
 436     res = ((y & MASK_SIGN) == MASK_SIGN) ? y : x;
 437     // two numbers are equal, return either
 438     BID_RETURN (res);
 439   }
 440
 441   res = ((sig_n_prime.w[1] == 0) && (sig_x > sig_n_prime.w[0])) ? y : x;
 442   BID_RETURN (res);
 443 }
 444
 445 /*****************************************************************************
 446  *  BID64 maximum function - returns greater of two numbers
 447  *****************************************************************************/
 448
 449 #if DECIMAL_CALL_BY_REFERENCE
 450 void
 451 bid64_maxnum (UINT64 * pres, UINT64 * px, UINT64 * py _EXC_FLAGS_PARAM) {
 452   UINT64 x = *px;
 453   UINT64 y = *py;
 454 #else
 455 UINT64
 456 bid64_maxnum (UINT64 x, UINT64 y _EXC_FLAGS_PARAM) {
 457 #endif
 458
 459   UINT64 res;
 460   int exp_x, exp_y;
 461   UINT64 sig_x, sig_y;
 462   UINT128 sig_n_prime;
 463   char x_is_zero = 0, y_is_zero = 0;
 464
 465   // check for non-canonical x
 466   if ((x & MASK_NAN) == MASK_NAN) {     // x is NaN
 467     x = x & 0xfe03ffffffffffffull;      // clear G6-G12
 468     if ((x & 0x0003ffffffffffffull) > 999999999999999ull) {
 469       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 470     }
 471   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
 472     x = x & (MASK_SIGN | MASK_INF);
 473   } else {      // x is not special
 474     // check for non-canonical values - treated as zero
 475     if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 476       // if the steering bits are 11, then the exponent is G[0:w+1]
 477       if (((x & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
 478           9999999999999999ull) {
 479         // non-canonical
 480         x = (x & MASK_SIGN) | ((x & MASK_BINARY_EXPONENT2) << 2);
 481       } // else canonical
 482     }   // else canonical
 483   }
 484
 485   // check for non-canonical y
 486   if ((y & MASK_NAN) == MASK_NAN) {     // y is NaN
 487     y = y & 0xfe03ffffffffffffull;      // clear G6-G12
 488     if ((y & 0x0003ffffffffffffull) > 999999999999999ull) {
 489       y = y & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 490     }
 491   } else if ((y & MASK_INF) == MASK_INF) {      // check for Infinity
 492     y = y & (MASK_SIGN | MASK_INF);
 493   } else {      // y is not special
 494     // check for non-canonical values - treated as zero
 495     if ((y & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 496       // if the steering bits are 11, then the exponent is G[0:w+1]
 497       if (((y & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
 498           9999999999999999ull) {
 499         // non-canonical
 500         y = (y & MASK_SIGN) | ((y & MASK_BINARY_EXPONENT2) << 2);
 501       } // else canonical
 502     }   // else canonical
 503   }
 504
 505   // NaN (CASE1)
 506   if ((x & MASK_NAN) == MASK_NAN) {     // x is NAN
 507     if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNaN
 508       // if x is SNAN, then return quiet (x)
 509       *pfpsf |= INVALID_EXCEPTION;      // set exception if SNaN
 510       x = x & 0xfdffffffffffffffull;    // quietize x
 511       res = x;
 512     } else {    // x is QNaN
 513       if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
 514         if ((y & MASK_SNAN) == MASK_SNAN) {     // y is SNAN
 515           *pfpsf |= INVALID_EXCEPTION;  // set invalid flag
 516         }
 517         res = x;
 518       } else {
 519         res = y;
 520       }
 521     }
 522     BID_RETURN (res);
 523   } else if ((y & MASK_NAN) == MASK_NAN) {      // y is NaN, but x is not
 524     if ((y & MASK_SNAN) == MASK_SNAN) {
 525       *pfpsf |= INVALID_EXCEPTION;      // set exception if SNaN
 526       y = y & 0xfdffffffffffffffull;    // quietize y
 527       res = y;
 528     } else {
 529       // will return x (which is not NaN)
 530       res = x;
 531     }
 532     BID_RETURN (res);
 533   }
 534   // SIMPLE (CASE2)
 535   // if all the bits are the same, these numbers are equal (not Greater).
 536   if (x == y) {
 537     res = x;
 538     BID_RETURN (res);
 539   }
 540   // INFINITY (CASE3)
 541   if ((x & MASK_INF) == MASK_INF) {
 542     // if x is neg infinity, there is no way it is greater than y, return y
 543     // x is pos infinity, it is greater, unless y is positive infinity =>
 544     // return y!=pos_infinity
 545     if (((x & MASK_SIGN) == MASK_SIGN)) {
 546       res = y;
 547       BID_RETURN (res);
 548     } else {
 549       res = (((y & MASK_INF) != MASK_INF)
 550              || ((y & MASK_SIGN) == MASK_SIGN)) ? x : y;
 551       BID_RETURN (res);
 552     }
 553   } else if ((y & MASK_INF) == MASK_INF) {
 554     // x is finite, so if y is positive infinity, then x is less, return y
 555     //                 if y is negative infinity, then x is greater, return x
 556     res = ((y & MASK_SIGN) == MASK_SIGN) ? x : y;
 557     BID_RETURN (res);
 558   }
 559   // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
 560   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 561     exp_x = (x & MASK_BINARY_EXPONENT2) >> 51;
 562     sig_x = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 563   } else {
 564     exp_x = (x & MASK_BINARY_EXPONENT1) >> 53;
 565     sig_x = (x & MASK_BINARY_SIG1);
 566   }
 567
 568   // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
 569   if ((y & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 570     exp_y = (y & MASK_BINARY_EXPONENT2) >> 51;
 571     sig_y = (y & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 572   } else {
 573     exp_y = (y & MASK_BINARY_EXPONENT1) >> 53;
 574     sig_y = (y & MASK_BINARY_SIG1);
 575   }
 576
 577   // ZERO (CASE4)
 578   // some properties:
 579   //    (+ZERO == -ZERO) => therefore
 580   //        ignore the sign, and neither number is greater
 581   //    (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
 582   //        ignore the exponent field
 583   //    (Any non-canonical # is considered 0)
 584   if (sig_x == 0) {
 585     x_is_zero = 1;
 586   }
 587   if (sig_y == 0) {
 588     y_is_zero = 1;
 589   }
 590
 591   if (x_is_zero && y_is_zero) {
 592     // if both numbers are zero, neither is greater => return NOTGREATERTHAN
 593     res = y;
 594     BID_RETURN (res);
 595   } else if (x_is_zero) {
 596     // is x is zero, it is greater if Y is negative
 597     res = ((y & MASK_SIGN) == MASK_SIGN) ? x : y;
 598     BID_RETURN (res);
 599   } else if (y_is_zero) {
 600     // is y is zero, X is greater if it is positive
 601     res = ((x & MASK_SIGN) != MASK_SIGN) ? x : y;;
 602     BID_RETURN (res);
 603   }
 604   // OPPOSITE SIGN (CASE5)
 605   // now, if the sign bits differ, x is greater if y is negative
 606   if (((x ^ y) & MASK_SIGN) == MASK_SIGN) {
 607     res = ((y & MASK_SIGN) == MASK_SIGN) ? x : y;
 608     BID_RETURN (res);
 609   }
 610   // REDUNDANT REPRESENTATIONS (CASE6)
 611
 612   // if both components are either bigger or smaller,
 613   //     it is clear what needs to be done
 614   if (sig_x > sig_y && exp_x >= exp_y) {
 615     res = ((x & MASK_SIGN) != MASK_SIGN) ? x : y;
 616     BID_RETURN (res);
 617   }
 618   if (sig_x < sig_y && exp_x <= exp_y) {
 619     res = ((x & MASK_SIGN) == MASK_SIGN) ? x : y;
 620     BID_RETURN (res);
 621   }
 622   // if exp_x is 15 greater than exp_y, no need for compensation
 623   if (exp_x - exp_y > 15) {
 624     res = ((x & MASK_SIGN) != MASK_SIGN) ? x : y;
 625     // difference cannot be > 10^15
 626     BID_RETURN (res);
 627   }
 628   // if exp_x is 15 less than exp_y, no need for compensation
 629   if (exp_y - exp_x > 15) {
 630     res = ((x & MASK_SIGN) == MASK_SIGN) ? x : y;
 631     BID_RETURN (res);
 632   }
 633   // if |exp_x - exp_y| < 15, it comes down to the compensated significand
 634   if (exp_x > exp_y) {  // to simplify the loop below,
 635     // otherwise adjust the x significand upwards
 636     __mul_64x64_to_128MACH (sig_n_prime, sig_x,
 637                             mult_factor[exp_x - exp_y]);
 638     // if postitive, return whichever significand is larger
 639     // (converse if negative)
 640     if (sig_n_prime.w[1] == 0 && (sig_n_prime.w[0] == sig_y)) {
 641       res = y;
 642       BID_RETURN (res);
 643     }
 644     res = (((sig_n_prime.w[1] > 0)
 645             || sig_n_prime.w[0] > sig_y) ^ ((x & MASK_SIGN) ==
 646                                             MASK_SIGN)) ? x : y;
 647     BID_RETURN (res);
 648   }
 649   // adjust the y significand upwards
 650   __mul_64x64_to_128MACH (sig_n_prime, sig_y,
 651                           mult_factor[exp_y - exp_x]);
 652
 653   // if postitive, return whichever significand is larger (converse if negative)
 654   if (sig_n_prime.w[1] == 0 && (sig_n_prime.w[0] == sig_x)) {
 655     res = y;
 656     BID_RETURN (res);
 657   }
 658   res = (((sig_n_prime.w[1] == 0)
 659           && (sig_x > sig_n_prime.w[0])) ^ ((x & MASK_SIGN) ==
 660                                             MASK_SIGN)) ? x : y;
 661   BID_RETURN (res);
 662 }
 663
 664 /*****************************************************************************
 665  *  BID64 maximum magnitude function - returns greater of two numbers
 666  *****************************************************************************/
 667
 668 #if DECIMAL_CALL_BY_REFERENCE
 669 void
 670 bid64_maxnum_mag (UINT64 * pres, UINT64 * px,
 671                   UINT64 * py _EXC_FLAGS_PARAM) {
 672   UINT64 x = *px;
 673   UINT64 y = *py;
 674 #else
 675 UINT64
 676 bid64_maxnum_mag (UINT64 x, UINT64 y _EXC_FLAGS_PARAM) {
 677 #endif
 678
 679   UINT64 res;
 680   int exp_x, exp_y;
 681   UINT64 sig_x, sig_y;
 682   UINT128 sig_n_prime;
 683
 684   // check for non-canonical x
 685   if ((x & MASK_NAN) == MASK_NAN) {     // x is NaN
 686     x = x & 0xfe03ffffffffffffull;      // clear G6-G12
 687     if ((x & 0x0003ffffffffffffull) > 999999999999999ull) {
 688       x = x & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 689     }
 690   } else if ((x & MASK_INF) == MASK_INF) {      // check for Infinity
 691     x = x & (MASK_SIGN | MASK_INF);
 692   } else {      // x is not special
 693     // check for non-canonical values - treated as zero
 694     if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 695       // if the steering bits are 11, then the exponent is G[0:w+1]
 696       if (((x & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
 697           9999999999999999ull) {
 698         // non-canonical
 699         x = (x & MASK_SIGN) | ((x & MASK_BINARY_EXPONENT2) << 2);
 700       } // else canonical
 701     }   // else canonical
 702   }
 703
 704   // check for non-canonical y
 705   if ((y & MASK_NAN) == MASK_NAN) {     // y is NaN
 706     y = y & 0xfe03ffffffffffffull;      // clear G6-G12
 707     if ((y & 0x0003ffffffffffffull) > 999999999999999ull) {
 708       y = y & 0xfe00000000000000ull;    // clear G6-G12 and the payload bits
 709     }
 710   } else if ((y & MASK_INF) == MASK_INF) {      // check for Infinity
 711     y = y & (MASK_SIGN | MASK_INF);
 712   } else {      // y is not special
 713     // check for non-canonical values - treated as zero
 714     if ((y & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 715       // if the steering bits are 11, then the exponent is G[0:w+1]
 716       if (((y & MASK_BINARY_SIG2) | MASK_BINARY_OR2) >
 717           9999999999999999ull) {
 718         // non-canonical
 719         y = (y & MASK_SIGN) | ((y & MASK_BINARY_EXPONENT2) << 2);
 720       } // else canonical
 721     }   // else canonical
 722   }
 723
 724   // NaN (CASE1)
 725   if ((x & MASK_NAN) == MASK_NAN) {     // x is NAN
 726     if ((x & MASK_SNAN) == MASK_SNAN) { // x is SNaN
 727       // if x is SNAN, then return quiet (x)
 728       *pfpsf |= INVALID_EXCEPTION;      // set exception if SNaN
 729       x = x & 0xfdffffffffffffffull;    // quietize x
 730       res = x;
 731     } else {    // x is QNaN
 732       if ((y & MASK_NAN) == MASK_NAN) { // y is NAN
 733         if ((y & MASK_SNAN) == MASK_SNAN) {     // y is SNAN
 734           *pfpsf |= INVALID_EXCEPTION;  // set invalid flag
 735         }
 736         res = x;
 737       } else {
 738         res = y;
 739       }
 740     }
 741     BID_RETURN (res);
 742   } else if ((y & MASK_NAN) == MASK_NAN) {      // y is NaN, but x is not
 743     if ((y & MASK_SNAN) == MASK_SNAN) {
 744       *pfpsf |= INVALID_EXCEPTION;      // set exception if SNaN
 745       y = y & 0xfdffffffffffffffull;    // quietize y
 746       res = y;
 747     } else {
 748       // will return x (which is not NaN)
 749       res = x;
 750     }
 751     BID_RETURN (res);
 752   }
 753   // SIMPLE (CASE2)
 754   // if all the bits are the same, these numbers are equal, return either number
 755   if (x == y) {
 756     res = x;
 757     BID_RETURN (res);
 758   }
 759   // INFINITY (CASE3)
 760   if ((x & MASK_INF) == MASK_INF) {
 761     // x is infinity, its magnitude is greater than or equal to y
 762     // return y as long as x isn't negative infinity
 763     res = ((x & MASK_SIGN) == MASK_SIGN
 764            && (y & MASK_INF) == MASK_INF) ? y : x;
 765     BID_RETURN (res);
 766   } else if ((y & MASK_INF) == MASK_INF) {
 767     // y is infinity, then it must be greater in magnitude
 768     res = y;
 769     BID_RETURN (res);
 770   }
 771   // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
 772   if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 773     exp_x = (x & MASK_BINARY_EXPONENT2) >> 51;
 774     sig_x = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 775   } else {
 776     exp_x = (x & MASK_BINARY_EXPONENT1) >> 53;
 777     sig_x = (x & MASK_BINARY_SIG1);
 778   }
 779
 780   // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
 781   if ((y & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
 782     exp_y = (y & MASK_BINARY_EXPONENT2) >> 51;
 783     sig_y = (y & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
 784   } else {
 785     exp_y = (y & MASK_BINARY_EXPONENT1) >> 53;
 786     sig_y = (y & MASK_BINARY_SIG1);
 787   }
 788
 789   // ZERO (CASE4)
 790   // some properties:
 791   //    (+ZERO == -ZERO) => therefore
 792   //        ignore the sign, and neither number is greater
 793   //    (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
 794   //        ignore the exponent field
 795   //    (Any non-canonical # is considered 0)
 796   if (sig_x == 0) {
 797     res = y;    // x_is_zero, its magnitude must be smaller than y
 798     BID_RETURN (res);
 799   }
 800   if (sig_y == 0) {
 801     res = x;    // y_is_zero, its magnitude must be smaller than x
 802     BID_RETURN (res);
 803   }
 804   // REDUNDANT REPRESENTATIONS (CASE6)
 805   // if both components are either bigger or smaller,
 806   // it is clear what needs to be done
 807   if (sig_x > sig_y && exp_x >= exp_y) {
 808     res = x;
 809     BID_RETURN (res);
 810   }
 811   if (sig_x < sig_y && exp_x <= exp_y) {
 812     res = y;
 813     BID_RETURN (res);
 814   }
 815   // if exp_x is 15 greater than exp_y, no need for compensation
 816   if (exp_x - exp_y > 15) {
 817     res = x;    // difference cannot be greater than 10^15
 818     BID_RETURN (res);
 819   }
 820   // if exp_x is 15 less than exp_y, no need for compensation
 821   if (exp_y - exp_x > 15) {
 822     res = y;
 823     BID_RETURN (res);
 824   }
 825   // if |exp_x - exp_y| < 15, it comes down to the compensated significand
 826   if (exp_x > exp_y) {  // to simplify the loop below,
 827     // otherwise adjust the x significand upwards
 828     __mul_64x64_to_128MACH (sig_n_prime, sig_x,
 829                             mult_factor[exp_x - exp_y]);
 830     // now, sig_n_prime has: sig_x * 10^(exp_x-exp_y),
 831     // this is the compensated signif.
 832     if (sig_n_prime.w[1] == 0 && (sig_n_prime.w[0] == sig_y)) {
 833       // two numbers are equal, return maxNum(x,y)
 834       res = ((y & MASK_SIGN) == MASK_SIGN) ? x : y;
 835       BID_RETURN (res);
 836     }
 837     // now, if compensated_x (sig_n_prime) is greater than y return y,
 838     // otherwise return x
 839     res = ((sig_n_prime.w[1] != 0) || sig_n_prime.w[0] > sig_y) ? x : y;
 840     BID_RETURN (res);
 841   }
 842   // exp_y must be greater than exp_x, thus adjust the y significand upwards
 843   __mul_64x64_to_128MACH (sig_n_prime, sig_y,
 844                           mult_factor[exp_y - exp_x]);
 845
 846   if (sig_n_prime.w[1] == 0 && (sig_n_prime.w[0] == sig_x)) {
 847     res = ((y & MASK_SIGN) == MASK_SIGN) ? x : y;
 848     // two numbers are equal, return either
 849     BID_RETURN (res);
 850   }
 851
 852   res = ((sig_n_prime.w[1] == 0) && (sig_x > sig_n_prime.w[0])) ? x : y;
 853   BID_RETURN (res);
 854 }