1 /* Copyright (C) 2007-2023 Free Software Foundation, Inc.
2
3 This file is part of GCC.
4
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 3, or (at your option) any later
8 version.
9
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13 for more details.
14
15 Under Section 7 of GPL version 3, you are granted additional
16 permissions described in the GCC Runtime Library Exception, version
17 3.1, as published by the Free Software Foundation.
18
19 You should have received a copy of the GNU General Public License and
20 a copy of the GCC Runtime Library Exception along with this program;
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22 <http://www.gnu.org/licenses/>. */
23
24 #define BID_128RES
25 #include "bid_internal.h"
26
27 /*****************************************************************************
28 * BID128 nextup
29 ****************************************************************************/
30
31 #if DECIMAL_CALL_BY_REFERENCE
32 void
33 bid128_nextup (UINT128 * pres,
34 UINT128 *
35 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
36 UINT128 x = *px;
37 #else
38 UINT128
39 bid128_nextup (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
40 _EXC_INFO_PARAM) {
41 #endif
42
43 UINT128 res;
44 UINT64 x_sign;
45 UINT64 x_exp;
46 int exp;
47 BID_UI64DOUBLE tmp1;
48 int x_nr_bits;
49 int q1, ind;
50 UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
51
52 BID_SWAP128 (x);
53 // unpack the argument
54 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
55 C1.w[1] = x.w[1] & MASK_COEFF;
56 C1.w[0] = x.w[0];
57
58 // check for NaN or Infinity
59 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
60 // x is special
61 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
62 // if x = NaN, then res = Q (x)
63 // check first for non-canonical NaN payload
64 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
65 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
66 && (x.w[0] > 0x38c15b09ffffffffull))) {
67 x.w[1] = x.w[1] & 0xffffc00000000000ull;
68 x.w[0] = 0x0ull;
69 }
70 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
71 // set invalid flag
72 *pfpsf |= INVALID_EXCEPTION;
73 // return quiet (x)
74 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
75 res.w[0] = x.w[0];
76 } else { // x is QNaN
77 // return x
78 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
79 res.w[0] = x.w[0];
80 }
81 } else { // x is not NaN, so it must be infinity
82 if (!x_sign) { // x is +inf
83 res.w[1] = 0x7800000000000000ull; // +inf
84 res.w[0] = 0x0000000000000000ull;
85 } else { // x is -inf
86 res.w[1] = 0xdfffed09bead87c0ull; // -MAXFP = -999...99 * 10^emax
87 res.w[0] = 0x378d8e63ffffffffull;
88 }
89 }
90 BID_RETURN (res);
91 }
92 // check for non-canonical values (treated as zero)
93 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
94 // non-canonical
95 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
96 C1.w[1] = 0; // significand high
97 C1.w[0] = 0; // significand low
98 } else { // G0_G1 != 11
99 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
100 if (C1.w[1] > 0x0001ed09bead87c0ull ||
101 (C1.w[1] == 0x0001ed09bead87c0ull
102 && C1.w[0] > 0x378d8e63ffffffffull)) {
103 // x is non-canonical if coefficient is larger than 10^34 -1
104 C1.w[1] = 0;
105 C1.w[0] = 0;
106 } else { // canonical
107 ;
108 }
109 }
110
111 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
112 // x is +/-0
113 res.w[1] = 0x0000000000000000ull; // +1 * 10^emin
114 res.w[0] = 0x0000000000000001ull;
115 } else { // x is not special and is not zero
116 if (x.w[1] == 0x5fffed09bead87c0ull
117 && x.w[0] == 0x378d8e63ffffffffull) {
118 // x = +MAXFP = 999...99 * 10^emax
119 res.w[1] = 0x7800000000000000ull; // +inf
120 res.w[0] = 0x0000000000000000ull;
121 } else if (x.w[1] == 0x8000000000000000ull
122 && x.w[0] == 0x0000000000000001ull) {
123 // x = -MINFP = 1...99 * 10^emin
124 res.w[1] = 0x8000000000000000ull; // -0
125 res.w[0] = 0x0000000000000000ull;
126 } else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
127 // can add/subtract 1 ulp to the significand
128
129 // Note: we could check here if x >= 10^34 to speed up the case q1 = 34
130 // q1 = nr. of decimal digits in x
131 // determine first the nr. of bits in x
132 if (C1.w[1] == 0) {
133 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
134 // split the 64-bit value in two 32-bit halves to avoid rnd errors
135 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
136 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
137 x_nr_bits =
138 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
139 0x3ff);
140 } else { // x < 2^32
141 tmp1.d = (double) (C1.w[0]); // exact conversion
142 x_nr_bits =
143 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
144 0x3ff);
145 }
146 } else { // if x < 2^53
147 tmp1.d = (double) C1.w[0]; // exact conversion
148 x_nr_bits =
149 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
150 }
151 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
152 tmp1.d = (double) C1.w[1]; // exact conversion
153 x_nr_bits =
154 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
155 }
156 q1 = nr_digits[x_nr_bits - 1].digits;
157 if (q1 == 0) {
158 q1 = nr_digits[x_nr_bits - 1].digits1;
159 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
160 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
161 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
162 q1++;
163 }
164 // if q1 < P34 then pad the significand with zeros
165 if (q1 < P34) {
166 exp = (x_exp >> 49) - 6176;
167 if (exp + 6176 > P34 - q1) {
168 ind = P34 - q1; // 1 <= ind <= P34 - 1
169 // pad with P34 - q1 zeros, until exponent = emin
170 // C1 = C1 * 10^ind
171 if (q1 <= 19) { // 64-bit C1
172 if (ind <= 19) { // 64-bit 10^ind and 64-bit C1
173 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
174 } else { // 128-bit 10^ind and 64-bit C1
175 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
176 }
177 } else { // C1 is (most likely) 128-bit
178 if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely)
179 __mul_128x64_to_128 (C1, ten2k64[ind], C1);
180 } else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
181 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
182 } else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
183 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
184 }
185 }
186 x_exp = x_exp - ((UINT64) ind << 49);
187 } else { // pad with zeros until the exponent reaches emin
188 ind = exp + 6176;
189 // C1 = C1 * 10^ind
190 if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
191 if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind
192 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
193 } else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
194 __mul_128x64_to_128 (C1, ten2k64[ind], C1);
195 }
196 } else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
197 // 64-bit C1, 128-bit 10^ind
198 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
199 }
200 x_exp = EXP_MIN;
201 }
202 }
203 if (!x_sign) { // x > 0
204 // add 1 ulp (add 1 to the significand)
205 C1.w[0]++;
206 if (C1.w[0] == 0)
207 C1.w[1]++;
208 if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34
209 C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33
210 C1.w[0] = 0x38c15b0a00000000ull;
211 x_exp = x_exp + EXP_P1;
212 }
213 } else { // x < 0
214 // subtract 1 ulp (subtract 1 from the significand)
215 C1.w[0]--;
216 if (C1.w[0] == 0xffffffffffffffffull)
217 C1.w[1]--;
218 if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1
219 C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1
220 C1.w[0] = 0x378d8e63ffffffffull;
221 x_exp = x_exp - EXP_P1;
222 }
223 }
224 // assemble the result
225 res.w[1] = x_sign | x_exp | C1.w[1];
226 res.w[0] = C1.w[0];
227 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
228 } // end x is not special and is not zero
229 BID_RETURN (res);
230 }
231
232 /*****************************************************************************
233 * BID128 nextdown
234 ****************************************************************************/
235
236 #if DECIMAL_CALL_BY_REFERENCE
237 void
238 bid128_nextdown (UINT128 * pres,
239 UINT128 *
240 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
241 UINT128 x = *px;
242 #else
243 UINT128
244 bid128_nextdown (UINT128 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
245 _EXC_INFO_PARAM) {
246 #endif
247
248 UINT128 res;
249 UINT64 x_sign;
250 UINT64 x_exp;
251 int exp;
252 BID_UI64DOUBLE tmp1;
253 int x_nr_bits;
254 int q1, ind;
255 UINT128 C1; // C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (UINT64)
256
257 BID_SWAP128 (x);
258 // unpack the argument
259 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
260 C1.w[1] = x.w[1] & MASK_COEFF;
261 C1.w[0] = x.w[0];
262
263 // check for NaN or Infinity
264 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
265 // x is special
266 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
267 // if x = NaN, then res = Q (x)
268 // check first for non-canonical NaN payload
269 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
270 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
271 && (x.w[0] > 0x38c15b09ffffffffull))) {
272 x.w[1] = x.w[1] & 0xffffc00000000000ull;
273 x.w[0] = 0x0ull;
274 }
275 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
276 // set invalid flag
277 *pfpsf |= INVALID_EXCEPTION;
278 // return quiet (x)
279 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
280 res.w[0] = x.w[0];
281 } else { // x is QNaN
282 // return x
283 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
284 res.w[0] = x.w[0];
285 }
286 } else { // x is not NaN, so it must be infinity
287 if (!x_sign) { // x is +inf
288 res.w[1] = 0x5fffed09bead87c0ull; // +MAXFP = +999...99 * 10^emax
289 res.w[0] = 0x378d8e63ffffffffull;
290 } else { // x is -inf
291 res.w[1] = 0xf800000000000000ull; // -inf
292 res.w[0] = 0x0000000000000000ull;
293 }
294 }
295 BID_RETURN (res);
296 }
297 // check for non-canonical values (treated as zero)
298 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
299 // non-canonical
300 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
301 C1.w[1] = 0; // significand high
302 C1.w[0] = 0; // significand low
303 } else { // G0_G1 != 11
304 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
305 if (C1.w[1] > 0x0001ed09bead87c0ull ||
306 (C1.w[1] == 0x0001ed09bead87c0ull
307 && C1.w[0] > 0x378d8e63ffffffffull)) {
308 // x is non-canonical if coefficient is larger than 10^34 -1
309 C1.w[1] = 0;
310 C1.w[0] = 0;
311 } else { // canonical
312 ;
313 }
314 }
315
316 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
317 // x is +/-0
318 res.w[1] = 0x8000000000000000ull; // -1 * 10^emin
319 res.w[0] = 0x0000000000000001ull;
320 } else { // x is not special and is not zero
321 if (x.w[1] == 0xdfffed09bead87c0ull
322 && x.w[0] == 0x378d8e63ffffffffull) {
323 // x = -MAXFP = -999...99 * 10^emax
324 res.w[1] = 0xf800000000000000ull; // -inf
325 res.w[0] = 0x0000000000000000ull;
326 } else if (x.w[1] == 0x0ull && x.w[0] == 0x0000000000000001ull) { // +MINFP
327 res.w[1] = 0x0000000000000000ull; // +0
328 res.w[0] = 0x0000000000000000ull;
329 } else { // -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
330 // can add/subtract 1 ulp to the significand
331
332 // Note: we could check here if x >= 10^34 to speed up the case q1 = 34
333 // q1 = nr. of decimal digits in x
334 // determine first the nr. of bits in x
335 if (C1.w[1] == 0) {
336 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
337 // split the 64-bit value in two 32-bit halves to avoid rnd errors
338 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
339 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
340 x_nr_bits =
341 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
342 0x3ff);
343 } else { // x < 2^32
344 tmp1.d = (double) (C1.w[0]); // exact conversion
345 x_nr_bits =
346 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) -
347 0x3ff);
348 }
349 } else { // if x < 2^53
350 tmp1.d = (double) C1.w[0]; // exact conversion
351 x_nr_bits =
352 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
353 }
354 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
355 tmp1.d = (double) C1.w[1]; // exact conversion
356 x_nr_bits =
357 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
358 }
359 q1 = nr_digits[x_nr_bits - 1].digits;
360 if (q1 == 0) {
361 q1 = nr_digits[x_nr_bits - 1].digits1;
362 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
363 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
364 && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
365 q1++;
366 }
367 // if q1 < P then pad the significand with zeros
368 if (q1 < P34) {
369 exp = (x_exp >> 49) - 6176;
370 if (exp + 6176 > P34 - q1) {
371 ind = P34 - q1; // 1 <= ind <= P34 - 1
372 // pad with P34 - q1 zeros, until exponent = emin
373 // C1 = C1 * 10^ind
374 if (q1 <= 19) { // 64-bit C1
375 if (ind <= 19) { // 64-bit 10^ind and 64-bit C1
376 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
377 } else { // 128-bit 10^ind and 64-bit C1
378 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
379 }
380 } else { // C1 is (most likely) 128-bit
381 if (ind <= 14) { // 64-bit 10^ind and 128-bit C1 (most likely)
382 __mul_128x64_to_128 (C1, ten2k64[ind], C1);
383 } else if (ind <= 19) { // 64-bit 10^ind and 64-bit C1 (q1 <= 19)
384 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
385 } else { // 128-bit 10^ind and 64-bit C1 (C1 must be 64-bit)
386 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
387 }
388 }
389 x_exp = x_exp - ((UINT64) ind << 49);
390 } else { // pad with zeros until the exponent reaches emin
391 ind = exp + 6176;
392 // C1 = C1 * 10^ind
393 if (ind <= 19) { // 1 <= P34 - q1 <= 19 <=> 15 <= q1 <= 33
394 if (q1 <= 19) { // 64-bit C1, 64-bit 10^ind
395 __mul_64x64_to_128MACH (C1, C1.w[0], ten2k64[ind]);
396 } else { // 20 <= q1 <= 33 => 128-bit C1, 64-bit 10^ind
397 __mul_128x64_to_128 (C1, ten2k64[ind], C1);
398 }
399 } else { // if 20 <= P34 - q1 <= 33 <=> 1 <= q1 <= 14 =>
400 // 64-bit C1, 128-bit 10^ind
401 __mul_128x64_to_128 (C1, C1.w[0], ten2k128[ind - 20]);
402 }
403 x_exp = EXP_MIN;
404 }
405 }
406 if (x_sign) { // x < 0
407 // add 1 ulp (add 1 to the significand)
408 C1.w[0]++;
409 if (C1.w[0] == 0)
410 C1.w[1]++;
411 if (C1.w[1] == 0x0001ed09bead87c0ull && C1.w[0] == 0x378d8e6400000000ull) { // if C1 = 10^34
412 C1.w[1] = 0x0000314dc6448d93ull; // C1 = 10^33
413 C1.w[0] = 0x38c15b0a00000000ull;
414 x_exp = x_exp + EXP_P1;
415 }
416 } else { // x > 0
417 // subtract 1 ulp (subtract 1 from the significand)
418 C1.w[0]--;
419 if (C1.w[0] == 0xffffffffffffffffull)
420 C1.w[1]--;
421 if (x_exp != 0 && C1.w[1] == 0x0000314dc6448d93ull && C1.w[0] == 0x38c15b09ffffffffull) { // if C1 = 10^33 - 1
422 C1.w[1] = 0x0001ed09bead87c0ull; // C1 = 10^34 - 1
423 C1.w[0] = 0x378d8e63ffffffffull;
424 x_exp = x_exp - EXP_P1;
425 }
426 }
427 // assemble the result
428 res.w[1] = x_sign | x_exp | C1.w[1];
429 res.w[0] = C1.w[0];
430 } // end -MAXFP <= x <= -MINFP - 1 ulp OR MINFP <= x <= MAXFP - 1 ulp
431 } // end x is not special and is not zero
432 BID_RETURN (res);
433 }
434
435 /*****************************************************************************
436 * BID128 nextafter
437 ****************************************************************************/
438
439 #if DECIMAL_CALL_BY_REFERENCE
440 void
441 bid128_nextafter (UINT128 * pres, UINT128 * px,
442 UINT128 *
443 py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
444 {
445 UINT128 x = *px;
446 UINT128 y = *py;
447 UINT128 xnswp = *px;
448 UINT128 ynswp = *py;
449 #else
450 UINT128
451 bid128_nextafter (UINT128 x,
452 UINT128 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
453 _EXC_INFO_PARAM) {
454 UINT128 xnswp = x;
455 UINT128 ynswp = y;
456 #endif
457
458 UINT128 res;
459 UINT128 tmp1, tmp2, tmp3;
460 FPSC tmp_fpsf = 0; // dummy fpsf for calls to comparison functions
461 int res1, res2;
462 UINT64 x_exp;
463
464
465 BID_SWAP128 (x);
466 BID_SWAP128 (y);
467 // check for NaNs
468 if (((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL)
469 || ((y.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
470 // x is special or y is special
471 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
472 // if x = NaN, then res = Q (x)
473 // check first for non-canonical NaN payload
474 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
475 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
476 && (x.w[0] > 0x38c15b09ffffffffull))) {
477 x.w[1] = x.w[1] & 0xffffc00000000000ull;
478 x.w[0] = 0x0ull;
479 }
480 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
481 // set invalid flag
482 *pfpsf |= INVALID_EXCEPTION;
483 // return quiet (x)
484 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
485 res.w[0] = x.w[0];
486 } else { // x is QNaN
487 // return x
488 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
489 res.w[0] = x.w[0];
490 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
491 // set invalid flag
492 *pfpsf |= INVALID_EXCEPTION;
493 }
494 }
495 BID_RETURN (res)
496 } else if ((y.w[1] & MASK_NAN) == MASK_NAN) { // y is NAN
497 // if x = NaN, then res = Q (x)
498 // check first for non-canonical NaN payload
499 if (((y.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
500 (((y.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull)
501 && (y.w[0] > 0x38c15b09ffffffffull))) {
502 y.w[1] = y.w[1] & 0xffffc00000000000ull;
503 y.w[0] = 0x0ull;
504 }
505 if ((y.w[1] & MASK_SNAN) == MASK_SNAN) { // y is SNAN
506 // set invalid flag
507 *pfpsf |= INVALID_EXCEPTION;
508 // return quiet (x)
509 res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
510 res.w[0] = y.w[0];
511 } else { // x is QNaN
512 // return x
513 res.w[1] = y.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
514 res.w[0] = y.w[0];
515 }
516 BID_RETURN (res)
517 } else { // at least one is infinity
518 if ((x.w[1] & MASK_ANY_INF) == MASK_INF) { // x = inf
519 x.w[1] = x.w[1] & (MASK_SIGN | MASK_INF);
520 x.w[0] = 0x0ull;
521 }
522 if ((y.w[1] & MASK_ANY_INF) == MASK_INF) { // y = inf
523 y.w[1] = y.w[1] & (MASK_SIGN | MASK_INF);
524 y.w[0] = 0x0ull;
525 }
526 }
527 }
528 // neither x nor y is NaN
529
530 // if not infinity, check for non-canonical values x (treated as zero)
531 if ((x.w[1] & MASK_ANY_INF) != MASK_INF) { // x != inf
532 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
533 // non-canonical
534 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
535 x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
536 x.w[0] = 0x0ull;
537 } else { // G0_G1 != 11
538 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
539 if ((x.w[1] & MASK_COEFF) > 0x0001ed09bead87c0ull ||
540 ((x.w[1] & MASK_COEFF) == 0x0001ed09bead87c0ull
541 && x.w[0] > 0x378d8e63ffffffffull)) {
542 // x is non-canonical if coefficient is larger than 10^34 -1
543 x.w[1] = (x.w[1] & MASK_SIGN) | x_exp;
544 x.w[0] = 0x0ull;
545 } else { // canonical
546 ;
547 }
548 }
549 }
550 // no need to check for non-canonical y
551
552 // neither x nor y is NaN
553 tmp_fpsf = *pfpsf; // save fpsf
554 #if DECIMAL_CALL_BY_REFERENCE
555 bid128_quiet_equal (&res1, &xnswp,
556 &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
557 _EXC_INFO_ARG);
558 bid128_quiet_greater (&res2, &xnswp,
559 &ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
560 _EXC_INFO_ARG);
561 #else
562 res1 =
563 bid128_quiet_equal (xnswp,
564 ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
565 _EXC_INFO_ARG);
566 res2 =
567 bid128_quiet_greater (xnswp,
568 ynswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
569 _EXC_INFO_ARG);
570 #endif
571 *pfpsf = tmp_fpsf; // restore fpsf
572
573 if (res1) { // x = y
574 // return x with the sign of y
575 res.w[1] =
576 (x.w[1] & 0x7fffffffffffffffull) | (y.
577 w[1] & 0x8000000000000000ull);
578 res.w[0] = x.w[0];
579 } else if (res2) { // x > y
580 #if DECIMAL_CALL_BY_REFERENCE
581 bid128_nextdown (&res,
582 &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
583 _EXC_INFO_ARG);
584 #else
585 res =
586 bid128_nextdown (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG
587 _EXC_INFO_ARG);
588 #endif
589 BID_SWAP128 (res);
590 } else { // x < y
591 #if DECIMAL_CALL_BY_REFERENCE
592 bid128_nextup (&res,
593 &xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
594 #else
595 res =
596 bid128_nextup (xnswp _EXC_FLAGS_ARG _EXC_MASKS_ARG _EXC_INFO_ARG);
597 #endif
598 BID_SWAP128 (res);
599 }
600 // if the operand x is finite but the result is infinite, signal
601 // overflow and inexact
602 if (((x.w[1] & MASK_SPECIAL) != MASK_SPECIAL)
603 && ((res.w[1] & MASK_SPECIAL) == MASK_SPECIAL)) {
604 // set the inexact flag
605 *pfpsf |= INEXACT_EXCEPTION;
606 // set the overflow flag
607 *pfpsf |= OVERFLOW_EXCEPTION;
608 }
609 // if the result is in (-10^emin, 10^emin), and is different from the
610 // operand x, signal underflow and inexact
611 tmp1.w[HIGH_128W] = 0x0000314dc6448d93ull;
612 tmp1.w[LOW_128W] = 0x38c15b0a00000000ull; // +100...0[34] * 10^emin
613 tmp2.w[HIGH_128W] = res.w[1] & 0x7fffffffffffffffull;
614 tmp2.w[LOW_128W] = res.w[0];
615 tmp3.w[HIGH_128W] = res.w[1];
616 tmp3.w[LOW_128W] = res.w[0];
617 tmp_fpsf = *pfpsf; // save fpsf
618 #if DECIMAL_CALL_BY_REFERENCE
619 bid128_quiet_greater (&res1, &tmp1,
620 &tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
621 _EXC_INFO_ARG);
622 bid128_quiet_not_equal (&res2, &xnswp,
623 &tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
624 _EXC_INFO_ARG);
625 #else
626 res1 =
627 bid128_quiet_greater (tmp1,
628 tmp2 _EXC_FLAGS_ARG _EXC_MASKS_ARG
629 _EXC_INFO_ARG);
630 res2 =
631 bid128_quiet_not_equal (xnswp,
632 tmp3 _EXC_FLAGS_ARG _EXC_MASKS_ARG
633 _EXC_INFO_ARG);
634 #endif
635 *pfpsf = tmp_fpsf; // restore fpsf
636 if (res1 && res2) {
637 // set the inexact flag
638 *pfpsf |= INEXACT_EXCEPTION;
639 // set the underflow flag
640 *pfpsf |= UNDERFLOW_EXCEPTION;
641 }
642 BID_RETURN (res);
643 }