libgcc/config/libbid/bid128_minmax.c

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