Home | History | Annotate | Line # | Download | only in libbid
      1 /* Copyright (C) 2007-2022 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 #define BID_128RES
     25 #include "bid_internal.h"
     26 
     27 /*****************************************************************************
     28  *  BID128 minimum number
     29  *****************************************************************************/
     30 
     31 #if DECIMAL_CALL_BY_REFERENCE
     32 void
     33 bid128_minnum (UINT128 * pres, UINT128 * px,
     34 	       UINT128 * py _EXC_FLAGS_PARAM) {
     35   UINT128 x = *px;
     36   UINT128 y = *py;
     37 #else
     38 UINT128
     39 bid128_minnum (UINT128 x, UINT128 y _EXC_FLAGS_PARAM) {
     40 #endif
     41 
     42   UINT128 res;
     43   int exp_x, exp_y;
     44   int diff;
     45   UINT128 sig_x, sig_y;
     46   UINT192 sig_n_prime192;
     47   UINT256 sig_n_prime256;
     48   char x_is_zero = 0, y_is_zero = 0;
     49 
     50   BID_SWAP128 (x);
     51   BID_SWAP128 (y);
     52 
     53   // check for non-canonical x
     54   if ((x.w[1] & MASK_NAN) == MASK_NAN) {	// x is NAN
     55     x.w[1] = x.w[1] & 0xfe003fffffffffffull;	// clear out G[6]-G[16]
     56     // check for non-canonical NaN payload
     57     if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
     58 	(((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
     59 	 (x.w[0] > 0x38c15b09ffffffffull))) {
     60       x.w[1] = x.w[1] & 0xffffc00000000000ull;
     61       x.w[0] = 0x0ull;
     62     }
     63   } else if ((x.w[1] & MASK_ANY_INF) == MASK_INF) {	// x = inf
     64     x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
     65     x.w[0] = 0x0ull;
     66   } else {	// x is not special
     67     // check for non-canonical values - treated as zero
     68     if ((x.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) {	// G0_G1=11
     69       // non-canonical
     70       x.w[1] = (x.w[1] & MASK_SIGN) | ((x.w[1] << 2) & MASK_EXP);
     71       x.w[0] = 0x0ull;
     72     } else {	// G0_G1 != 11
     73       if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
     74 	  ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
     75 	   && x.w[0] > 0x378d8e63ffffffffull)) {
     76 	// x is non-canonical if coefficient is larger than 10^34 -1
     77 	x.w[1] = (x.w[1] & MASK_SIGN) | (x.w[1] & MASK_EXP);
     78 	x.w[0] = 0x0ull;
     79       } else {	// canonical
     80 	;
     81       }
     82     }
     83   }
     84   // check for non-canonical y
     85   if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NAN
     86     y.w[1] = y.w[1] & 0xfe003fffffffffffull;	// clear out G[6]-G[16]
     87     // check for non-canonical NaN payload
     88     if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
     89 	(((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
     90 	 (y.w[0] > 0x38c15b09ffffffffull))) {
     91       y.w[1] = y.w[1] & 0xffffc00000000000ull;
     92       y.w[0] = 0x0ull;
     93     }
     94   } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) {	// y = inf
     95     y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
     96     y.w[0] = 0x0ull;
     97   } else {	// y is not special
     98     // check for non-canonical values - treated as zero
     99     if ((y.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) {	// G0_G1=11
    100       // non-canonical
    101       y.w[1] = (y.w[1] & MASK_SIGN) | ((y.w[1] << 2) & MASK_EXP);
    102       y.w[0] = 0x0ull;
    103     } else {	// G0_G1 != 11
    104       if ((y.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
    105 	  ((y.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
    106 	   && y.w[0] > 0x378d8e63ffffffffull)) {
    107 	// y is non-canonical if coefficient is larger than 10^34 -1
    108 	y.w[1] = (y.w[1] & MASK_SIGN) | (y.w[1] & MASK_EXP);
    109 	y.w[0] = 0x0ull;
    110       } else {	// canonical
    111 	;
    112       }
    113     }
    114   }
    115 
    116   // NaN (CASE1)
    117   if ((x.w[1] & MASK_NAN) == MASK_NAN) {	// x is NAN
    118     if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {	// x is SNaN
    119       // if x is SNAN, then return quiet (x)
    120       *pfpsf |= INVALID_EXCEPTION;	// set exception if SNaN
    121       x.w[1] = x.w[1] & 0xfdffffffffffffffull;	// quietize x
    122       res = x;
    123     } else {	// x is QNaN
    124       if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NAN
    125 	if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {	// y is SNAN
    126 	  *pfpsf |= INVALID_EXCEPTION;	// set invalid flag
    127 	}
    128 	res = x;
    129       } else {
    130 	res = y;
    131       }
    132     }
    133     BID_RETURN (res);
    134   } else if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NaN, but x is not
    135     if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {
    136       *pfpsf |= INVALID_EXCEPTION;	// set exception if SNaN
    137       y.w[1] = y.w[1] & 0xfdffffffffffffffull;	// quietize y
    138       res = y;
    139     } else {
    140       // will return x (which is not NaN)
    141       res = x;
    142     }
    143     BID_RETURN (res);
    144   }
    145   // SIMPLE (CASE2)
    146   // if all the bits are the same, these numbers are equal (not Greater).
    147   if (x.w[0] == y.w[0] && x.w[1] == y.w[1]) {
    148     res = x;
    149     BID_RETURN (res);
    150   }
    151   // INFINITY (CASE3)
    152   if ((x.w[1] & MASK_INF) == MASK_INF) {
    153     // if x is neg infinity, there is no way it is greater than y, return 0
    154     res = (((x.w[1] & MASK_SIGN) == MASK_SIGN)) ? x : y;
    155     BID_RETURN (res);
    156   } else if ((y.w[1] & MASK_INF) == MASK_INF) {
    157     // x is finite, so if y is positive infinity, then x is less, return 0
    158     //                 if y is negative infinity, then x is greater, return 1
    159     res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    160     BID_RETURN (res);
    161   }
    162   // CONVERT X
    163   sig_x.w[1] = x.w[1] & 0x0001ffffffffffffull;
    164   sig_x.w[0] = x.w[0];
    165   exp_x = (x.w[1] >> 49) & 0x000000000003fffull;
    166 
    167   // CONVERT Y
    168   exp_y = (y.w[1] >> 49) & 0x0000000000003fffull;
    169   sig_y.w[1] = y.w[1] & 0x0001ffffffffffffull;
    170   sig_y.w[0] = y.w[0];
    171 
    172   // ZERO (CASE4)
    173   // some properties:
    174   //    (+ZERO == -ZERO) => therefore ignore the sign
    175   //    (ZERO x 10^A == ZERO x 10^B) for any valid A, B => ignore the exponent
    176   //    field
    177   //    (Any non-canonical # is considered 0)
    178   if ((sig_x.w[1] == 0) && (sig_x.w[0] == 0)) {
    179     x_is_zero = 1;
    180   }
    181   if ((sig_y.w[1] == 0) && (sig_y.w[0] == 0)) {
    182     y_is_zero = 1;
    183   }
    184 
    185   if (x_is_zero && y_is_zero) {
    186     // if both numbers are zero, neither is greater => return either number
    187     res = x;
    188     BID_RETURN (res);
    189   } else if (x_is_zero) {
    190     // is x is zero, it is greater if Y is negative
    191     res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    192     BID_RETURN (res);
    193   } else if (y_is_zero) {
    194     // is y is zero, X is greater if it is positive
    195     res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? y : x;
    196     BID_RETURN (res);
    197   }
    198   // OPPOSITE SIGN (CASE5)
    199   // now, if the sign bits differ, x is greater if y is negative
    200   if (((x.w[1] ^ y.w[1]) & MASK_SIGN) == MASK_SIGN) {
    201     res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    202     BID_RETURN (res);
    203   }
    204   // REDUNDANT REPRESENTATIONS (CASE6)
    205   // if exponents are the same, then we have a simple comparison of
    206   //    the significands
    207   if (exp_y == exp_x) {
    208     res = (((sig_x.w[1] > sig_y.w[1])
    209 	    || (sig_x.w[1] == sig_y.w[1]
    210 		&& sig_x.w[0] >= sig_y.w[0])) ^ ((x.w[1] & MASK_SIGN) ==
    211 						 MASK_SIGN)) ? y : x;
    212     BID_RETURN (res);
    213   }
    214   // if both components are either bigger or smaller, it is clear what
    215   //    needs to be done
    216   if (sig_x.w[1] >= sig_y.w[1] && sig_x.w[0] >= sig_y.w[0]
    217       && exp_x > exp_y) {
    218     res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? y : x;
    219     BID_RETURN (res);
    220   }
    221   if (sig_x.w[1] <= sig_y.w[1] && sig_x.w[0] <= sig_y.w[0]
    222       && exp_x < exp_y) {
    223     res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    224     BID_RETURN (res);
    225   }
    226 
    227   diff = exp_x - exp_y;
    228 
    229   // if |exp_x - exp_y| < 33, it comes down to the compensated significand
    230   if (diff > 0) {	// to simplify the loop below,
    231     // if exp_x is 33 greater than exp_y, no need for compensation
    232     if (diff > 33) {
    233       // difference cannot be greater than 10^33
    234       res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? y : x;
    235       BID_RETURN (res);
    236     }
    237     if (diff > 19) {	//128 by 128 bit multiply -> 256 bits
    238       __mul_128x128_to_256 (sig_n_prime256, sig_x, ten2k128[diff - 20]);
    239       // if postitive, return whichever significand is larger
    240       // (converse if negative)
    241       res = ((((sig_n_prime256.w[3] > 0) || sig_n_prime256.w[2] > 0)
    242 	      || (sig_n_prime256.w[1] > sig_y.w[1])
    243 	      || (sig_n_prime256.w[1] == sig_y.w[1]
    244 		  && sig_n_prime256.w[0] >
    245 		  sig_y.w[0])) ^ ((y.w[1] & MASK_SIGN) ==
    246 				  MASK_SIGN)) ? y : x;
    247       BID_RETURN (res);
    248     }
    249     __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_x);
    250     // if postitive, return whichever significand is larger
    251     // (converse if negative)
    252     res =
    253       (((sig_n_prime192.w[2] > 0) || (sig_n_prime192.w[1] > sig_y.w[1])
    254 	|| (sig_n_prime192.w[1] == sig_y.w[1]
    255 	    && sig_n_prime192.w[0] >
    256 	    sig_y.w[0])) ^ ((y.w[1] & MASK_SIGN) == MASK_SIGN)) ? y : x;
    257     BID_RETURN (res);
    258   }
    259   diff = exp_y - exp_x;
    260   // if exp_x is 33 less than exp_y, no need for compensation
    261   if (diff > 33) {
    262     res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    263     BID_RETURN (res);
    264   }
    265   if (diff > 19) {	//128 by 128 bit multiply -> 256 bits
    266     // adjust the y significand upwards
    267     __mul_128x128_to_256 (sig_n_prime256, sig_y, ten2k128[diff - 20]);
    268     // if postitive, return whichever significand is larger
    269     // (converse if negative)
    270     res =
    271       ((sig_n_prime256.w[3] != 0 || sig_n_prime256.w[2] != 0
    272 	|| (sig_n_prime256.w[1] > sig_x.w[1]
    273 	    || (sig_n_prime256.w[1] == sig_x.w[1]
    274 		&& sig_n_prime256.w[0] >
    275 		sig_x.w[0]))) ^ ((x.w[1] & MASK_SIGN) ==
    276 				 MASK_SIGN)) ? x : y;
    277     BID_RETURN (res);
    278   }
    279   // adjust the y significand upwards
    280   __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_y);
    281   // if postitive, return whichever significand is larger (converse if negative)
    282   res =
    283     ((sig_n_prime192.w[2] != 0
    284       || (sig_n_prime192.w[1] > sig_x.w[1]
    285 	  || (sig_n_prime192.w[1] == sig_x.w[1]
    286 	      && sig_n_prime192.w[0] > sig_x.w[0])))
    287      ^ ((y.w[1] & MASK_SIGN) == MASK_SIGN)) ? x : y;
    288   BID_RETURN (res);
    289 }
    290 
    291 /*****************************************************************************
    292  *  BID128 minimum magnitude function - returns greater of two numbers
    293  *****************************************************************************/
    294 
    295 #if DECIMAL_CALL_BY_REFERENCE
    296 void
    297 bid128_minnum_mag (UINT128 * pres, UINT128 * px,
    298 		   UINT128 * py _EXC_FLAGS_PARAM) {
    299   UINT128 x = *px;
    300   UINT128 y = *py;
    301 #else
    302 UINT128
    303 bid128_minnum_mag (UINT128 x, UINT128 y _EXC_FLAGS_PARAM) {
    304 #endif
    305 
    306   UINT128 res;
    307   int exp_x, exp_y;
    308   int diff;
    309   UINT128 sig_x, sig_y;
    310   UINT192 sig_n_prime192;
    311   UINT256 sig_n_prime256;
    312 
    313   BID_SWAP128 (x);
    314   BID_SWAP128 (y);
    315 
    316   // check for non-canonical x
    317   if ((x.w[1] & MASK_NAN) == MASK_NAN) {	// x is NAN
    318     x.w[1] = x.w[1] & 0xfe003fffffffffffull;	// clear out G[6]-G[16]
    319     // check for non-canonical NaN payload
    320     if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
    321 	(((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
    322 	 (x.w[0] > 0x38c15b09ffffffffull))) {
    323       x.w[1] = x.w[1] & 0xffffc00000000000ull;
    324       x.w[0] = 0x0ull;
    325     }
    326   } else if ((x.w[1] & MASK_ANY_INF) == MASK_INF) {	// x = inf
    327     x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
    328     x.w[0] = 0x0ull;
    329   } else {	// x is not special
    330     // check for non-canonical values - treated as zero
    331     if ((x.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) {	// G0_G1=11
    332       // non-canonical
    333       x.w[1] = (x.w[1] & MASK_SIGN) | ((x.w[1] << 2) & MASK_EXP);
    334       x.w[0] = 0x0ull;
    335     } else {	// G0_G1 != 11
    336       if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
    337 	  ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
    338 	   && x.w[0] > 0x378d8e63ffffffffull)) {
    339 	// x is non-canonical if coefficient is larger than 10^34 -1
    340 	x.w[1] = (x.w[1] & MASK_SIGN) | (x.w[1] & MASK_EXP);
    341 	x.w[0] = 0x0ull;
    342       } else {	// canonical
    343 	;
    344       }
    345     }
    346   }
    347   // check for non-canonical y
    348   if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NAN
    349     y.w[1] = y.w[1] & 0xfe003fffffffffffull;	// clear out G[6]-G[16]
    350     // check for non-canonical NaN payload
    351     if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
    352 	(((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
    353 	 (y.w[0] > 0x38c15b09ffffffffull))) {
    354       y.w[1] = y.w[1] & 0xffffc00000000000ull;
    355       y.w[0] = 0x0ull;
    356     }
    357   } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) {	// y = inf
    358     y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
    359     y.w[0] = 0x0ull;
    360   } else {	// y is not special
    361     // check for non-canonical values - treated as zero
    362     if ((y.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) {	// G0_G1=11
    363       // non-canonical
    364       y.w[1] = (y.w[1] & MASK_SIGN) | ((y.w[1] << 2) & MASK_EXP);
    365       y.w[0] = 0x0ull;
    366     } else {	// G0_G1 != 11
    367       if ((y.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
    368 	  ((y.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
    369 	   && y.w[0] > 0x378d8e63ffffffffull)) {
    370 	// y is non-canonical if coefficient is larger than 10^34 -1
    371 	y.w[1] = (y.w[1] & MASK_SIGN) | (y.w[1] & MASK_EXP);
    372 	y.w[0] = 0x0ull;
    373       } else {	// canonical
    374 	;
    375       }
    376     }
    377   }
    378 
    379   // NaN (CASE1)
    380   if ((x.w[1] & MASK_NAN) == MASK_NAN) {	// x is NAN
    381     if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {	// x is SNaN
    382       // if x is SNAN, then return quiet (x)
    383       *pfpsf |= INVALID_EXCEPTION;	// set exception if SNaN
    384       x.w[1] = x.w[1] & 0xfdffffffffffffffull;	// quietize x
    385       res = x;
    386     } else {	// x is QNaN
    387       if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NAN
    388 	if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {	// y is SNAN
    389 	  *pfpsf |= INVALID_EXCEPTION;	// set invalid flag
    390 	}
    391 	res = x;
    392       } else {
    393 	res = y;
    394       }
    395     }
    396     BID_RETURN (res);
    397   } else if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NaN, but x is not
    398     if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {
    399       *pfpsf |= INVALID_EXCEPTION;	// set exception if SNaN
    400       y.w[1] = y.w[1] & 0xfdffffffffffffffull;	// quietize y
    401       res = y;
    402     } else {
    403       // will return x (which is not NaN)
    404       res = x;
    405     }
    406     BID_RETURN (res);
    407   }
    408   // SIMPLE (CASE2)
    409   // if all the bits are the same, these numbers are equal (not Greater).
    410   if (x.w[0] == y.w[0] && x.w[1] == y.w[1]) {
    411     res = y;
    412     BID_RETURN (res);
    413   }
    414   // INFINITY (CASE3)
    415   if ((x.w[1] & MASK_INF) == MASK_INF) {
    416     // if x infinity, it has maximum magnitude.
    417     // Check if magnitudes are equal.  If x is negative, return it.
    418     res = ((x.w[1] & MASK_SIGN) == MASK_SIGN
    419 	   && (y.w[1] & MASK_INF) == MASK_INF) ? x : y;
    420     BID_RETURN (res);
    421   } else if ((y.w[1] & MASK_INF) == MASK_INF) {
    422     // x is finite, so if y is infinity, then x is less in magnitude
    423     res = x;
    424     BID_RETURN (res);
    425   }
    426   // CONVERT X
    427   sig_x.w[1] = x.w[1] & 0x0001ffffffffffffull;
    428   sig_x.w[0] = x.w[0];
    429   exp_x = (x.w[1] >> 49) & 0x000000000003fffull;
    430 
    431   // CONVERT Y
    432   exp_y = (y.w[1] >> 49) & 0x0000000000003fffull;
    433   sig_y.w[1] = y.w[1] & 0x0001ffffffffffffull;
    434   sig_y.w[0] = y.w[0];
    435 
    436   // ZERO (CASE4)
    437   // some properties:
    438   //    (+ZERO == -ZERO) => therefore ignore the sign
    439   //    (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
    440   //        therefore ignore the exponent field
    441   //    (Any non-canonical # is considered 0)
    442   if ((sig_x.w[1] == 0) && (sig_x.w[0] == 0)) {
    443     res = x;
    444     BID_RETURN (res);
    445   }
    446   if ((sig_y.w[1] == 0) && (sig_y.w[0] == 0)) {
    447     res = y;
    448     BID_RETURN (res);
    449   }
    450   // REDUNDANT REPRESENTATIONS (CASE6)
    451   // check if exponents are the same and significands are the same
    452   if (exp_y == exp_x && sig_x.w[1] == sig_y.w[1]
    453       && sig_x.w[0] == sig_y.w[0]) {
    454     if (x.w[1] & 0x8000000000000000ull) {	// x is negative
    455       res = x;
    456       BID_RETURN (res);
    457     } else {
    458       res = y;
    459       BID_RETURN (res);
    460     }
    461   } else if (((sig_x.w[1] > sig_y.w[1] || (sig_x.w[1] == sig_y.w[1]
    462 					   && sig_x.w[0] > sig_y.w[0]))
    463 	      && exp_x == exp_y)
    464 	     || ((sig_x.w[1] > sig_y.w[1]
    465 		  || (sig_x.w[1] == sig_y.w[1]
    466 		      && sig_x.w[0] >= sig_y.w[0]))
    467 		 && exp_x > exp_y)) {
    468     // if both components are either bigger or smaller, it is clear what
    469     // needs to be done; also if the magnitudes are equal
    470     res = y;
    471     BID_RETURN (res);
    472   } else if (((sig_y.w[1] > sig_x.w[1] || (sig_y.w[1] == sig_x.w[1]
    473 					   && sig_y.w[0] > sig_x.w[0]))
    474 	      && exp_y == exp_x)
    475 	     || ((sig_y.w[1] > sig_x.w[1]
    476 		  || (sig_y.w[1] == sig_x.w[1]
    477 		      && sig_y.w[0] >= sig_x.w[0]))
    478 		 && exp_y > exp_x)) {
    479     res = x;
    480     BID_RETURN (res);
    481   } else {
    482     ;	// continue
    483   }
    484   diff = exp_x - exp_y;
    485   // if |exp_x - exp_y| < 33, it comes down to the compensated significand
    486   if (diff > 0) {	// to simplify the loop below,
    487     // if exp_x is 33 greater than exp_y, no need for compensation
    488     if (diff > 33) {
    489       res = y;	// difference cannot be greater than 10^33
    490       BID_RETURN (res);
    491     }
    492     if (diff > 19) {	//128 by 128 bit multiply -> 256 bits
    493       __mul_128x128_to_256 (sig_n_prime256, sig_x, ten2k128[diff - 20]);
    494       // if positive, return whichever significand is larger
    495       // (converse if negative)
    496       if (sig_n_prime256.w[3] == 0 && (sig_n_prime256.w[2] == 0)
    497 	  && sig_n_prime256.w[1] == sig_y.w[1]
    498 	  && (sig_n_prime256.w[0] == sig_y.w[0])) {
    499 	res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;	// if equal
    500 	BID_RETURN (res);
    501       }
    502       res = (((sig_n_prime256.w[3] > 0) || sig_n_prime256.w[2] > 0)
    503 	     || (sig_n_prime256.w[1] > sig_y.w[1])
    504 	     || (sig_n_prime256.w[1] == sig_y.w[1]
    505 		 && sig_n_prime256.w[0] > sig_y.w[0])) ? y : x;
    506       BID_RETURN (res);
    507     }
    508     __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_x);
    509     // if positive, return whichever significand is larger
    510     // (converse if negative)
    511     if ((sig_n_prime192.w[2] == 0) && sig_n_prime192.w[1] == sig_y.w[1]
    512 	&& (sig_n_prime192.w[0] == sig_y.w[0])) {
    513       // if = in magnitude, return +, (if possible)
    514       res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    515       BID_RETURN (res);
    516     }
    517     res = ((sig_n_prime192.w[2] > 0)
    518 	   || (sig_n_prime192.w[1] > sig_y.w[1])
    519 	   || (sig_n_prime192.w[1] == sig_y.w[1]
    520 	       && sig_n_prime192.w[0] > sig_y.w[0])) ? y : x;
    521     BID_RETURN (res);
    522   }
    523   diff = exp_y - exp_x;
    524   // if exp_x is 33 less than exp_y, no need for compensation
    525   if (diff > 33) {
    526     res = x;
    527     BID_RETURN (res);
    528   }
    529   if (diff > 19) {	//128 by 128 bit multiply -> 256 bits
    530     // adjust the y significand upwards
    531     __mul_128x128_to_256 (sig_n_prime256, sig_y, ten2k128[diff - 20]);
    532     // if positive, return whichever significand is larger
    533     // (converse if negative)
    534     if (sig_n_prime256.w[3] == 0 && (sig_n_prime256.w[2] == 0)
    535 	&& sig_n_prime256.w[1] == sig_x.w[1]
    536 	&& (sig_n_prime256.w[0] == sig_x.w[0])) {
    537       // if = in magnitude, return +, (if possible)
    538       res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    539       BID_RETURN (res);
    540     }
    541     res = (sig_n_prime256.w[3] == 0 && sig_n_prime256.w[2] == 0
    542 	   && (sig_n_prime256.w[1] < sig_x.w[1]
    543 	       || (sig_n_prime256.w[1] == sig_x.w[1]
    544 		   && sig_n_prime256.w[0] < sig_x.w[0]))) ? y : x;
    545     BID_RETURN (res);
    546   }
    547   // adjust the y significand upwards
    548   __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_y);
    549   // if positive, return whichever significand is larger (converse if negative)
    550   if ((sig_n_prime192.w[2] == 0) && sig_n_prime192.w[1] == sig_x.w[1]
    551       && (sig_n_prime192.w[0] == sig_x.w[0])) {
    552     // if = in magnitude, return +, if possible)
    553     res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    554     BID_RETURN (res);
    555   }
    556   res = (sig_n_prime192.w[2] == 0
    557 	 && (sig_n_prime192.w[1] < sig_x.w[1]
    558 	     || (sig_n_prime192.w[1] == sig_x.w[1]
    559 		 && sig_n_prime192.w[0] < sig_x.w[0]))) ? y : x;
    560   BID_RETURN (res);
    561 }
    562 
    563 /*****************************************************************************
    564  *  BID128 maximum function - returns greater of two numbers
    565  *****************************************************************************/
    566 
    567 #if DECIMAL_CALL_BY_REFERENCE
    568 void
    569 bid128_maxnum (UINT128 * pres, UINT128 * px,
    570 	       UINT128 * py _EXC_FLAGS_PARAM) {
    571   UINT128 x = *px;
    572   UINT128 y = *py;
    573 #else
    574 UINT128
    575 bid128_maxnum (UINT128 x, UINT128 y _EXC_FLAGS_PARAM) {
    576 #endif
    577 
    578   UINT128 res;
    579   int exp_x, exp_y;
    580   int diff;
    581   UINT128 sig_x, sig_y;
    582   UINT192 sig_n_prime192;
    583   UINT256 sig_n_prime256;
    584   char x_is_zero = 0, y_is_zero = 0;
    585 
    586   BID_SWAP128 (x);
    587   BID_SWAP128 (y);
    588 
    589   // check for non-canonical x
    590   if ((x.w[1] & MASK_NAN) == MASK_NAN) {	// x is NAN
    591     x.w[1] = x.w[1] & 0xfe003fffffffffffull;	// clear out G[6]-G[16]
    592     // check for non-canonical NaN payload
    593     if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
    594 	(((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
    595 	 (x.w[0] > 0x38c15b09ffffffffull))) {
    596       x.w[1] = x.w[1] & 0xffffc00000000000ull;
    597       x.w[0] = 0x0ull;
    598     }
    599   } else if ((x.w[1] & MASK_ANY_INF) == MASK_INF) {	// x = inf
    600     x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
    601     x.w[0] = 0x0ull;
    602   } else {	// x is not special
    603     // check for non-canonical values - treated as zero
    604     if ((x.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) {	// G0_G1=11
    605       // non-canonical
    606       x.w[1] = (x.w[1] & MASK_SIGN) | ((x.w[1] << 2) & MASK_EXP);
    607       x.w[0] = 0x0ull;
    608     } else {	// G0_G1 != 11
    609       if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
    610 	  ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
    611 	   && x.w[0] > 0x378d8e63ffffffffull)) {
    612 	// x is non-canonical if coefficient is larger than 10^34 -1
    613 	x.w[1] = (x.w[1] & MASK_SIGN) | (x.w[1] & MASK_EXP);
    614 	x.w[0] = 0x0ull;
    615       } else {	// canonical
    616 	;
    617       }
    618     }
    619   }
    620   // check for non-canonical y
    621   if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NAN
    622     y.w[1] = y.w[1] & 0xfe003fffffffffffull;	// clear out G[6]-G[16]
    623     // check for non-canonical NaN payload
    624     if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
    625 	(((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
    626 	 (y.w[0] > 0x38c15b09ffffffffull))) {
    627       y.w[1] = y.w[1] & 0xffffc00000000000ull;
    628       y.w[0] = 0x0ull;
    629     }
    630   } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) {	// y = inf
    631     y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
    632     y.w[0] = 0x0ull;
    633   } else {	// y is not special
    634     // check for non-canonical values - treated as zero
    635     if ((y.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) {	// G0_G1=11
    636       // non-canonical
    637       y.w[1] = (y.w[1] & MASK_SIGN) | ((y.w[1] << 2) & MASK_EXP);
    638       y.w[0] = 0x0ull;
    639     } else {	// G0_G1 != 11
    640       if ((y.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
    641 	  ((y.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
    642 	   && y.w[0] > 0x378d8e63ffffffffull)) {
    643 	// y is non-canonical if coefficient is larger than 10^34 -1
    644 	y.w[1] = (y.w[1] & MASK_SIGN) | (y.w[1] & MASK_EXP);
    645 	y.w[0] = 0x0ull;
    646       } else {	// canonical
    647 	;
    648       }
    649     }
    650   }
    651 
    652   // NaN (CASE1)
    653   if ((x.w[1] & MASK_NAN) == MASK_NAN) {	// x is NAN
    654     if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {	// x is SNaN
    655       // if x is SNAN, then return quiet (x)
    656       *pfpsf |= INVALID_EXCEPTION;	// set exception if SNaN
    657       x.w[1] = x.w[1] & 0xfdffffffffffffffull;	// quietize x
    658       res = x;
    659     } else {	// x is QNaN
    660       if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NAN
    661 	if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {	// y is SNAN
    662 	  *pfpsf |= INVALID_EXCEPTION;	// set invalid flag
    663 	}
    664 	res = x;
    665       } else {
    666 	res = y;
    667       }
    668     }
    669     BID_RETURN (res);
    670   } else if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NaN, but x is not
    671     if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {
    672       *pfpsf |= INVALID_EXCEPTION;	// set exception if SNaN
    673       y.w[1] = y.w[1] & 0xfdffffffffffffffull;	// quietize y
    674       res = y;
    675     } else {
    676       // will return x (which is not NaN)
    677       res = x;
    678     }
    679     BID_RETURN (res);
    680   }
    681   // SIMPLE (CASE2)
    682   // if all the bits are the same, these numbers are equal (not Greater).
    683   if (x.w[0] == y.w[0] && x.w[1] == y.w[1]) {
    684     res = x;
    685     BID_RETURN (res);
    686   }
    687   // INFINITY (CASE3)
    688   if ((x.w[1] & MASK_INF) == MASK_INF) {
    689     res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? y : x;
    690     BID_RETURN (res);
    691   } else if ((y.w[1] & MASK_INF) == MASK_INF) {
    692     // x is finite, so if y is positive infinity, then x is less, return 0
    693     //                 if y is negative infinity, then x is greater, return 1
    694     res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
    695     BID_RETURN (res);
    696   }
    697   // CONVERT X
    698   sig_x.w[1] = x.w[1] & 0x0001ffffffffffffull;
    699   sig_x.w[0] = x.w[0];
    700   exp_x = (x.w[1] >> 49) & 0x000000000003fffull;
    701 
    702   // CONVERT Y
    703   exp_y = (y.w[1] >> 49) & 0x0000000000003fffull;
    704   sig_y.w[1] = y.w[1] & 0x0001ffffffffffffull;
    705   sig_y.w[0] = y.w[0];
    706 
    707   // ZERO (CASE4)
    708   // some properties:
    709   //    (+ZERO == -ZERO) => therefore ignore the sign
    710   //    (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
    711   //        therefore ignore the exponent field
    712   //    (Any non-canonical # is considered 0)
    713   if ((sig_x.w[1] == 0) && (sig_x.w[0] == 0)) {
    714     x_is_zero = 1;
    715   }
    716   if ((sig_y.w[1] == 0) && (sig_y.w[0] == 0)) {
    717     y_is_zero = 1;
    718   }
    719 
    720   if (x_is_zero && y_is_zero) {
    721     // if both numbers are zero, neither is greater => return either number
    722     res = x;
    723     BID_RETURN (res);
    724   } else if (x_is_zero) {
    725     // is x is zero, it is greater if Y is negative
    726     res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
    727     BID_RETURN (res);
    728   } else if (y_is_zero) {
    729     // is y is zero, X is greater if it is positive
    730     res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? x : y;
    731     BID_RETURN (res);
    732   }
    733   // OPPOSITE SIGN (CASE5)
    734   // now, if the sign bits differ, x is greater if y is negative
    735   if (((x.w[1] ^ y.w[1]) & MASK_SIGN) == MASK_SIGN) {
    736     res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
    737     BID_RETURN (res);
    738   }
    739   // REDUNDANT REPRESENTATIONS (CASE6)
    740   // if exponents are the same, then we have a simple comparison of
    741   // the significands
    742   if (exp_y == exp_x) {
    743     res = (((sig_x.w[1] > sig_y.w[1]) || (sig_x.w[1] == sig_y.w[1] &&
    744 					  sig_x.w[0] >= sig_y.w[0])) ^
    745 	   ((x.w[1] & MASK_SIGN) == MASK_SIGN)) ? x : y;
    746     BID_RETURN (res);
    747   }
    748   // if both components are either bigger or smaller, it is clear what
    749   // needs to be done
    750   if ((sig_x.w[1] > sig_y.w[1]
    751        || (sig_x.w[1] == sig_y.w[1] && sig_x.w[0] > sig_y.w[0]))
    752       && exp_x >= exp_y) {
    753     res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? x : y;
    754     BID_RETURN (res);
    755   }
    756   if ((sig_x.w[1] < sig_y.w[1]
    757        || (sig_x.w[1] == sig_y.w[1] && sig_x.w[0] < sig_y.w[0]))
    758       && exp_x <= exp_y) {
    759     res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
    760     BID_RETURN (res);
    761   }
    762   diff = exp_x - exp_y;
    763   // if |exp_x - exp_y| < 33, it comes down to the compensated significand
    764   if (diff > 0) {	// to simplify the loop below,
    765     // if exp_x is 33 greater than exp_y, no need for compensation
    766     if (diff > 33) {
    767       // difference cannot be greater than 10^33
    768       res = ((x.w[1] & MASK_SIGN) != MASK_SIGN) ? x : y;
    769       BID_RETURN (res);
    770     }
    771     if (diff > 19) {	//128 by 128 bit multiply -> 256 bits
    772       __mul_128x128_to_256 (sig_n_prime256, sig_x, ten2k128[diff - 20]);
    773       // if postitive, return whichever significand is larger
    774       // (converse if negative)
    775       res = ((((sig_n_prime256.w[3] > 0) || sig_n_prime256.w[2] > 0)
    776 	      || (sig_n_prime256.w[1] > sig_y.w[1])
    777 	      || (sig_n_prime256.w[1] == sig_y.w[1]
    778 		  && sig_n_prime256.w[0] >
    779 		  sig_y.w[0])) ^ ((y.w[1] & MASK_SIGN) ==
    780 				  MASK_SIGN)) ? x : y;
    781       BID_RETURN (res);
    782     }
    783     __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_x);
    784     // if postitive, return whichever significand is larger
    785     // (converse if negative)
    786     res =
    787       (((sig_n_prime192.w[2] > 0) || (sig_n_prime192.w[1] > sig_y.w[1])
    788 	|| (sig_n_prime192.w[1] == sig_y.w[1]
    789 	    && sig_n_prime192.w[0] >
    790 	    sig_y.w[0])) ^ ((y.w[1] & MASK_SIGN) == MASK_SIGN)) ? x : y;
    791     BID_RETURN (res);
    792   }
    793   diff = exp_y - exp_x;
    794   // if exp_x is 33 less than exp_y, no need for compensation
    795   if (diff > 33) {
    796     res = ((x.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
    797     BID_RETURN (res);
    798   }
    799   if (diff > 19) {	//128 by 128 bit multiply -> 256 bits
    800     // adjust the y significand upwards
    801     __mul_128x128_to_256 (sig_n_prime256, sig_y, ten2k128[diff - 20]);
    802     // if postitive, return whichever significand is larger
    803     // (converse if negative)
    804     res =
    805       ((sig_n_prime256.w[3] != 0 || sig_n_prime256.w[2] != 0
    806 	|| (sig_n_prime256.w[1] > sig_x.w[1]
    807 	    || (sig_n_prime256.w[1] == sig_x.w[1]
    808 		&& sig_n_prime256.w[0] >
    809 		sig_x.w[0]))) ^ ((x.w[1] & MASK_SIGN) !=
    810 				 MASK_SIGN)) ? x : y;
    811     BID_RETURN (res);
    812   }
    813   // adjust the y significand upwards
    814   __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_y);
    815   // if postitive, return whichever significand is larger (converse if negative)
    816   res =
    817     ((sig_n_prime192.w[2] != 0
    818       || (sig_n_prime192.w[1] > sig_x.w[1]
    819 	  || (sig_n_prime192.w[1] == sig_x.w[1]
    820 	      && sig_n_prime192.w[0] >
    821 	      sig_x.w[0]))) ^ ((y.w[1] & MASK_SIGN) !=
    822 			       MASK_SIGN)) ? x : y;
    823   BID_RETURN (res);
    824 }
    825 
    826 /*****************************************************************************
    827  *  BID128 maximum magnitude function - returns greater of two numbers
    828  *****************************************************************************/
    829 
    830 #if DECIMAL_CALL_BY_REFERENCE
    831 void
    832 bid128_maxnum_mag (UINT128 * pres, UINT128 * px,
    833 		   UINT128 * py _EXC_FLAGS_PARAM) {
    834   UINT128 x = *px;
    835   UINT128 y = *py;
    836 #else
    837 UINT128
    838 bid128_maxnum_mag (UINT128 x, UINT128 y _EXC_FLAGS_PARAM) {
    839 #endif
    840 
    841   UINT128 res;
    842   int exp_x, exp_y;
    843   int diff;
    844   UINT128 sig_x, sig_y;
    845   UINT192 sig_n_prime192;
    846   UINT256 sig_n_prime256;
    847 
    848   BID_SWAP128 (x);
    849   BID_SWAP128 (y);
    850 
    851   // check for non-canonical x
    852   if ((x.w[1] & MASK_NAN) == MASK_NAN) {	// x is NAN
    853     x.w[1] = x.w[1] & 0xfe003fffffffffffull;	// clear out G[6]-G[16]
    854     // check for non-canonical NaN payload
    855     if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
    856 	(((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
    857 	 (x.w[0] > 0x38c15b09ffffffffull))) {
    858       x.w[1] = x.w[1] & 0xffffc00000000000ull;
    859       x.w[0] = 0x0ull;
    860     }
    861   } else if ((x.w[1] & MASK_ANY_INF) == MASK_INF) {	// x = inf
    862     x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
    863     x.w[0] = 0x0ull;
    864   } else {	// x is not special
    865     // check for non-canonical values - treated as zero
    866     if ((x.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) {	// G0_G1=11
    867       // non-canonical
    868       x.w[1] = (x.w[1] & MASK_SIGN) | ((x.w[1] << 2) & MASK_EXP);
    869       x.w[0] = 0x0ull;
    870     } else {	// G0_G1 != 11
    871       if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
    872 	  ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
    873 	   && x.w[0] > 0x378d8e63ffffffffull)) {
    874 	// x is non-canonical if coefficient is larger than 10^34 -1
    875 	x.w[1] = (x.w[1] & MASK_SIGN) | (x.w[1] & MASK_EXP);
    876 	x.w[0] = 0x0ull;
    877       } else {	// canonical
    878 	;
    879       }
    880     }
    881   }
    882   // check for non-canonical y
    883   if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NAN
    884     y.w[1] = y.w[1] & 0xfe003fffffffffffull;	// clear out G[6]-G[16]
    885     // check for non-canonical NaN payload
    886     if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
    887 	(((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
    888 	 (y.w[0] > 0x38c15b09ffffffffull))) {
    889       y.w[1] = y.w[1] & 0xffffc00000000000ull;
    890       y.w[0] = 0x0ull;
    891     }
    892   } else if ((y.w[1] & MASK_ANY_INF) == MASK_INF) {	// y = inf
    893     y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
    894     y.w[0] = 0x0ull;
    895   } else {	// y is not special
    896     // check for non-canonical values - treated as zero
    897     if ((y.w[1] & MASK_STEERING_BITS) == MASK_STEERING_BITS) {	// G0_G1=11
    898       // non-canonical
    899       y.w[1] = (y.w[1] & MASK_SIGN) | ((y.w[1] << 2) & MASK_EXP);
    900       y.w[0] = 0x0ull;
    901     } else {	// G0_G1 != 11
    902       if ((y.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
    903 	  ((y.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull &&
    904 	   y.w[0] > 0x378d8e63ffffffffull)) {
    905 	// y is non-canonical if coefficient is larger than 10^34 -1
    906 	y.w[1] = (y.w[1] & MASK_SIGN) | (y.w[1] & MASK_EXP);
    907 	y.w[0] = 0x0ull;
    908       } else {	// canonical
    909 	;
    910       }
    911     }
    912   }
    913 
    914   // NaN (CASE1)
    915   if ((x.w[1] & MASK_NAN) == MASK_NAN) {	// x is NAN
    916     if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {	// x is SNaN
    917       // if x is SNAN, then return quiet (x)
    918       *pfpsf |= INVALID_EXCEPTION;	// set exception if SNaN
    919       x.w[1] = x.w[1] & 0xfdffffffffffffffull;	// quietize x
    920       res = x;
    921     } else {	// x is QNaN
    922       if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NAN
    923 	if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {	// y is SNAN
    924 	  *pfpsf |= INVALID_EXCEPTION;	// set invalid flag
    925 	}
    926 	res = x;
    927       } else {
    928 	res = y;
    929       }
    930     }
    931     BID_RETURN (res);
    932   } else if ((y.w[1] & MASK_NAN) == MASK_NAN) {	// y is NaN, but x is not
    933     if ((y.w[1] & MASK_SNAN) == MASK_SNAN) {
    934       *pfpsf |= INVALID_EXCEPTION;	// set exception if SNaN
    935       y.w[1] = y.w[1] & 0xfdffffffffffffffull;	// quietize y
    936       res = y;
    937     } else {
    938       // will return x (which is not NaN)
    939       res = x;
    940     }
    941     BID_RETURN (res);
    942   }
    943   // SIMPLE (CASE2)
    944   // if all the bits are the same, these numbers are equal (not Greater).
    945   if (x.w[0] == y.w[0] && x.w[1] == y.w[1]) {
    946     res = y;
    947     BID_RETURN (res);
    948   }
    949   // INFINITY (CASE3)
    950   if ((x.w[1] & MASK_INF) == MASK_INF) {
    951     // if x infinity, it has maximum magnitude
    952     res = ((x.w[1] & MASK_SIGN) == MASK_SIGN
    953 	   && (y.w[1] & MASK_INF) == MASK_INF) ? y : x;
    954     BID_RETURN (res);
    955   } else if ((y.w[1] & MASK_INF) == MASK_INF) {
    956     // x is finite, so if y is positive infinity, then x is less, return 0
    957     //                 if y is negative infinity, then x is greater, return 1
    958     res = y;
    959     BID_RETURN (res);
    960   }
    961   // CONVERT X
    962   sig_x.w[1] = x.w[1] & 0x0001ffffffffffffull;
    963   sig_x.w[0] = x.w[0];
    964   exp_x = (x.w[1] >> 49) & 0x000000000003fffull;
    965 
    966   // CONVERT Y
    967   exp_y = (y.w[1] >> 49) & 0x0000000000003fffull;
    968   sig_y.w[1] = y.w[1] & 0x0001ffffffffffffull;
    969   sig_y.w[0] = y.w[0];
    970 
    971   // ZERO (CASE4)
    972   // some properties:
    973   //    (+ZERO == -ZERO) => therefore ignore the sign
    974   //    (ZERO x 10^A == ZERO x 10^B) for any valid A, B =>
    975   //         therefore ignore the exponent field
    976   //    (Any non-canonical # is considered 0)
    977   if ((sig_x.w[1] == 0) && (sig_x.w[0] == 0)) {
    978     res = y;
    979     BID_RETURN (res);
    980   }
    981   if ((sig_y.w[1] == 0) && (sig_y.w[0] == 0)) {
    982     res = x;
    983     BID_RETURN (res);
    984   }
    985   // REDUNDANT REPRESENTATIONS (CASE6)
    986   if (exp_y == exp_x && sig_x.w[1] == sig_y.w[1]
    987       && sig_x.w[0] == sig_y.w[0]) {
    988     // check if exponents are the same and significands are the same
    989     if (x.w[1] & 0x8000000000000000ull) {	// x is negative
    990       res = y;
    991       BID_RETURN (res);
    992     } else {
    993       res = x;
    994       BID_RETURN (res);
    995     }
    996   } else if (((sig_x.w[1] > sig_y.w[1] || (sig_x.w[1] == sig_y.w[1]
    997 					   && sig_x.w[0] > sig_y.w[0]))
    998 	      && exp_x == exp_y)
    999 	     || ((sig_x.w[1] > sig_y.w[1]
   1000 		  || (sig_x.w[1] == sig_y.w[1]
   1001 		      && sig_x.w[0] >= sig_y.w[0]))
   1002 		 && exp_x > exp_y)) {
   1003     // if both components are either bigger or smaller, it is clear what
   1004     // needs to be done; also if the magnitudes are equal
   1005     res = x;
   1006     BID_RETURN (res);
   1007   } else if (((sig_y.w[1] > sig_x.w[1] || (sig_y.w[1] == sig_x.w[1]
   1008 					   && sig_y.w[0] > sig_x.w[0]))
   1009 	      && exp_y == exp_x)
   1010 	     || ((sig_y.w[1] > sig_x.w[1]
   1011 		  || (sig_y.w[1] == sig_x.w[1]
   1012 		      && sig_y.w[0] >= sig_x.w[0]))
   1013 		 && exp_y > exp_x)) {
   1014     res = y;
   1015     BID_RETURN (res);
   1016   } else {
   1017     ;	// continue
   1018   }
   1019   diff = exp_x - exp_y;
   1020   // if |exp_x - exp_y| < 33, it comes down to the compensated significand
   1021   if (diff > 0) {	// to simplify the loop below,
   1022     // if exp_x is 33 greater than exp_y, no need for compensation
   1023     if (diff > 33) {
   1024       res = x;	// difference cannot be greater than 10^33
   1025       BID_RETURN (res);
   1026     }
   1027     if (diff > 19) {	//128 by 128 bit multiply -> 256 bits
   1028       __mul_128x128_to_256 (sig_n_prime256, sig_x, ten2k128[diff - 20]);
   1029       // if postitive, return whichever significand is larger
   1030       // (converse if negative)
   1031       if (sig_n_prime256.w[3] == 0 && (sig_n_prime256.w[2] == 0)
   1032 	  && sig_n_prime256.w[1] == sig_y.w[1]
   1033 	  && (sig_n_prime256.w[0] == sig_y.w[0])) {
   1034 	res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;	// if equal
   1035 	BID_RETURN (res);
   1036       }
   1037       res = (((sig_n_prime256.w[3] > 0) || sig_n_prime256.w[2] > 0)
   1038 	     || (sig_n_prime256.w[1] > sig_y.w[1])
   1039 	     || (sig_n_prime256.w[1] == sig_y.w[1]
   1040 		 && sig_n_prime256.w[0] > sig_y.w[0])) ? x : y;
   1041       BID_RETURN (res);
   1042     }
   1043     __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_x);
   1044     // if postitive, return whichever significand is larger (converse if negative)
   1045     if ((sig_n_prime192.w[2] == 0) && sig_n_prime192.w[1] == sig_y.w[1]
   1046 	&& (sig_n_prime192.w[0] == sig_y.w[0])) {
   1047       // if equal, return positive magnitude
   1048       res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
   1049       BID_RETURN (res);
   1050     }
   1051     res = ((sig_n_prime192.w[2] > 0)
   1052 	   || (sig_n_prime192.w[1] > sig_y.w[1])
   1053 	   || (sig_n_prime192.w[1] == sig_y.w[1]
   1054 	       && sig_n_prime192.w[0] > sig_y.w[0])) ? x : y;
   1055     BID_RETURN (res);
   1056   }
   1057   diff = exp_y - exp_x;
   1058   // if exp_x is 33 less than exp_y, no need for compensation
   1059   if (diff > 33) {
   1060     res = y;
   1061     BID_RETURN (res);
   1062   }
   1063   if (diff > 19) {	//128 by 128 bit multiply -> 256 bits
   1064     // adjust the y significand upwards
   1065     __mul_128x128_to_256 (sig_n_prime256, sig_y, ten2k128[diff - 20]);
   1066     // if postitive, return whichever significand is larger
   1067     // (converse if negative)
   1068     if (sig_n_prime256.w[3] == 0 && (sig_n_prime256.w[2] == 0)
   1069 	&& sig_n_prime256.w[1] == sig_x.w[1]
   1070 	&& (sig_n_prime256.w[0] == sig_x.w[0])) {
   1071       // if equal, return positive (if possible)
   1072       res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
   1073       BID_RETURN (res);
   1074     }
   1075     res = (sig_n_prime256.w[3] == 0 && sig_n_prime256.w[2] == 0
   1076 	   && (sig_n_prime256.w[1] < sig_x.w[1]
   1077 	       || (sig_n_prime256.w[1] == sig_x.w[1]
   1078 		   && sig_n_prime256.w[0] < sig_x.w[0]))) ? x : y;
   1079     BID_RETURN (res);
   1080   }
   1081   // adjust the y significand upwards
   1082   __mul_64x128_to_192 (sig_n_prime192, ten2k64[diff], sig_y);
   1083   // if postitive, return whichever significand is larger (converse if negative)
   1084   if ((sig_n_prime192.w[2] == 0) && sig_n_prime192.w[1] == sig_x.w[1]
   1085       && (sig_n_prime192.w[0] == sig_x.w[0])) {
   1086     // if equal, return positive (if possible)
   1087     res = ((y.w[1] & MASK_SIGN) == MASK_SIGN) ? x : y;
   1088     BID_RETURN (res);
   1089   }
   1090   res = (sig_n_prime192.w[2] == 0
   1091 	 && (sig_n_prime192.w[1] < sig_x.w[1]
   1092 	     || (sig_n_prime192.w[1] == sig_x.w[1]
   1093 		 && sig_n_prime192.w[0] < sig_x.w[0]))) ? x : y;
   1094   BID_RETURN (res);
   1095 }
   1096