Subversion Repositories Kolibri OS

Rev

Blame | Last modification | View Log | RSS feed

  1. /* Copyright (C) 2007-2015 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.