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