1/* Copyright (C) 2007-2020 Free Software Foundation, Inc.
2
3This file is part of GCC.
4
5GCC is free software; you can redistribute it and/or modify it under
6the terms of the GNU General Public License as published by the Free
7Software Foundation; either version 3, or (at your option) any later
8version.
9
10GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11WARRANTY; without even the implied warranty of MERCHANTABILITY or
12FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
13for more details.
14
15Under Section 7 of GPL version 3, you are granted additional
16permissions described in the GCC Runtime Library Exception, version
173.1, as published by the Free Software Foundation.
18
19You should have received a copy of the GNU General Public License and
20a copy of the GCC Runtime Library Exception along with this program;
21see 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
32void
33bid128_minnum (UINT128 * pres, UINT128 * px,
34	       UINT128 * py _EXC_FLAGS_PARAM) {
35  UINT128 x = *px;
36  UINT128 y = *py;
37#else
38UINT128
39bid128_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
296void
297bid128_minnum_mag (UINT128 * pres, UINT128 * px,
298		   UINT128 * py _EXC_FLAGS_PARAM) {
299  UINT128 x = *px;
300  UINT128 y = *py;
301#else
302UINT128
303bid128_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
568void
569bid128_maxnum (UINT128 * pres, UINT128 * px,
570	       UINT128 * py _EXC_FLAGS_PARAM) {
571  UINT128 x = *px;
572  UINT128 y = *py;
573#else
574UINT128
575bid128_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
831void
832bid128_maxnum_mag (UINT128 * pres, UINT128 * px,
833		   UINT128 * py _EXC_FLAGS_PARAM) {
834  UINT128 x = *px;
835  UINT128 y = *py;
836#else
837UINT128
838bid128_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