(root)/
gcc-13.2.0/
libgcc/
config/
libbid/
bid128_minmax.c
       1  /* Copyright (C) 2007-2023 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  }