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
26 #include "bid_internal.h"
27
28 /*****************************************************************************
29 * BID128_round_integral_exact
30 ****************************************************************************/
31
32 BID128_FUNCTION_ARG1 (bid128_round_integral_exact, x)
33
34 UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
35 };
36 UINT64 x_sign;
37 UINT64 x_exp;
38 int exp; // unbiased exponent
39 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
40 UINT64 tmp64;
41 BID_UI64DOUBLE tmp1;
42 unsigned int x_nr_bits;
43 int q, ind, shift;
44 UINT128 C1;
45 UINT256 fstar;
46 UINT256 P256;
47
48 // check for NaN or Infinity
49 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
50 // x is special
51 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
52 // if x = NaN, then res = Q (x)
53 // check first for non-canonical NaN payload
54 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
55 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
56 (x.w[0] > 0x38c15b09ffffffffull))) {
57 x.w[1] = x.w[1] & 0xffffc00000000000ull;
58 x.w[0] = 0x0ull;
59 }
60 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
61 // set invalid flag
62 *pfpsf |= INVALID_EXCEPTION;
63 // return quiet (x)
64 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
65 res.w[0] = x.w[0];
66 } else { // x is QNaN
67 // return x
68 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
69 res.w[0] = x.w[0];
70 }
71 BID_RETURN (res)
72 } else { // x is not a NaN, so it must be infinity
73 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
74 // return +inf
75 res.w[1] = 0x7800000000000000ull;
76 res.w[0] = 0x0000000000000000ull;
77 } else { // x is -inf
78 // return -inf
79 res.w[1] = 0xf800000000000000ull;
80 res.w[0] = 0x0000000000000000ull;
81 }
82 BID_RETURN (res);
83 }
84 }
85 // unpack x
86 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
87 C1.w[1] = x.w[1] & MASK_COEFF;
88 C1.w[0] = x.w[0];
89
90 // check for non-canonical values (treated as zero)
91 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
92 // non-canonical
93 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
94 C1.w[1] = 0; // significand high
95 C1.w[0] = 0; // significand low
96 } else { // G0_G1 != 11
97 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
98 if (C1.w[1] > 0x0001ed09bead87c0ull ||
99 (C1.w[1] == 0x0001ed09bead87c0ull
100 && C1.w[0] > 0x378d8e63ffffffffull)) {
101 // x is non-canonical if coefficient is larger than 10^34 -1
102 C1.w[1] = 0;
103 C1.w[0] = 0;
104 } else { // canonical
105 ;
106 }
107 }
108
109 // test for input equal to zero
110 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
111 // x is 0
112 // return 0 preserving the sign bit and the preferred exponent
113 // of MAX(Q(x), 0)
114 if (x_exp <= (0x1820ull << 49)) {
115 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
116 } else {
117 res.w[1] = x_sign | x_exp;
118 }
119 res.w[0] = 0x0000000000000000ull;
120 BID_RETURN (res);
121 }
122 // x is not special and is not zero
123
124 switch (rnd_mode) {
125 case ROUNDING_TO_NEAREST:
126 case ROUNDING_TIES_AWAY:
127 // if (exp <= -(p+1)) return 0.0
128 if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35
129 res.w[1] = x_sign | 0x3040000000000000ull;
130 res.w[0] = 0x0000000000000000ull;
131 *pfpsf |= INEXACT_EXCEPTION;
132 BID_RETURN (res);
133 }
134 break;
135 case ROUNDING_DOWN:
136 // if (exp <= -p) return -1.0 or +0.0
137 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffa000000000000ull == -34
138 if (x_sign) {
139 // if negative, return negative 1, because we know coefficient
140 // is non-zero (would have been caught above)
141 res.w[1] = 0xb040000000000000ull;
142 res.w[0] = 0x0000000000000001ull;
143 } else {
144 // if positive, return positive 0, because we know coefficient is
145 // non-zero (would have been caught above)
146 res.w[1] = 0x3040000000000000ull;
147 res.w[0] = 0x0000000000000000ull;
148 }
149 *pfpsf |= INEXACT_EXCEPTION;
150 BID_RETURN (res);
151 }
152 break;
153 case ROUNDING_UP:
154 // if (exp <= -p) return -0.0 or +1.0
155 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
156 if (x_sign) {
157 // if negative, return negative 0, because we know the coefficient
158 // is non-zero (would have been caught above)
159 res.w[1] = 0xb040000000000000ull;
160 res.w[0] = 0x0000000000000000ull;
161 } else {
162 // if positive, return positive 1, because we know coefficient is
163 // non-zero (would have been caught above)
164 res.w[1] = 0x3040000000000000ull;
165 res.w[0] = 0x0000000000000001ull;
166 }
167 *pfpsf |= INEXACT_EXCEPTION;
168 BID_RETURN (res);
169 }
170 break;
171 case ROUNDING_TO_ZERO:
172 // if (exp <= -p) return -0.0 or +0.0
173 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
174 res.w[1] = x_sign | 0x3040000000000000ull;
175 res.w[0] = 0x0000000000000000ull;
176 *pfpsf |= INEXACT_EXCEPTION;
177 BID_RETURN (res);
178 }
179 break;
180 }
181
182 // q = nr. of decimal digits in x
183 // determine first the nr. of bits in x
184 if (C1.w[1] == 0) {
185 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
186 // split the 64-bit value in two 32-bit halves to avoid rounding errors
187 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
188 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
189 x_nr_bits =
190 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
191 } else { // x < 2^32
192 tmp1.d = (double) (C1.w[0]); // exact conversion
193 x_nr_bits =
194 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
195 }
196 } else { // if x < 2^53
197 tmp1.d = (double) C1.w[0]; // exact conversion
198 x_nr_bits =
199 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200 }
201 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
202 tmp1.d = (double) C1.w[1]; // exact conversion
203 x_nr_bits =
204 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
205 }
206
207 q = nr_digits[x_nr_bits - 1].digits;
208 if (q == 0) {
209 q = nr_digits[x_nr_bits - 1].digits1;
210 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
211 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
212 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
213 q++;
214 }
215 exp = (x_exp >> 49) - 6176;
216 if (exp >= 0) { // -exp <= 0
217 // the argument is an integer already
218 res.w[1] = x.w[1];
219 res.w[0] = x.w[0];
220 BID_RETURN (res);
221 }
222 // exp < 0
223 switch (rnd_mode) {
224 case ROUNDING_TO_NEAREST:
225 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
226 // need to shift right -exp digits from the coefficient; exp will be 0
227 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
228 // chop off ind digits from the lower part of C1
229 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
230 tmp64 = C1.w[0];
231 if (ind <= 19) {
232 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
233 } else {
234 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
235 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
236 }
237 if (C1.w[0] < tmp64)
238 C1.w[1]++;
239 // calculate C* and f*
240 // C* is actually floor(C*) in this case
241 // C* and f* need shifting and masking, as shown by
242 // shiftright128[] and maskhigh128[]
243 // 1 <= x <= 34
244 // kx = 10^(-x) = ten2mk128[ind - 1]
245 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
246 // the approximation of 10^(-x) was rounded up to 118 bits
247 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
248 // determine the value of res and fstar
249
250 // determine inexactness of the rounding of C*
251 // if (0 < f* - 1/2 < 10^(-x)) then
252 // the result is exact
253 // else // if (f* - 1/2 > T*) then
254 // the result is inexact
255 // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[]
256
257 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
258 // redundant shift = shiftright128[ind - 1]; // shift = 0
259 res.w[1] = P256.w[3];
260 res.w[0] = P256.w[2];
261 // redundant fstar.w[3] = 0;
262 // redundant fstar.w[2] = 0;
263 fstar.w[1] = P256.w[1];
264 fstar.w[0] = P256.w[0];
265 // fraction f* < 10^(-x) <=> midpoint
266 // f* is in the right position to be compared with
267 // 10^(-x) from ten2mk128[]
268 // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even)
269 if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
270 ((fstar.w[1] < (ten2mk128[ind - 1].w[1]))
271 || ((fstar.w[1] == ten2mk128[ind - 1].w[1])
272 && (fstar.w[0] < ten2mk128[ind - 1].w[0])))) {
273 // subract 1 to make even
274 if (res.w[0]-- == 0) {
275 res.w[1]--;
276 }
277 }
278 if (fstar.w[1] > 0x8000000000000000ull ||
279 (fstar.w[1] == 0x8000000000000000ull
280 && fstar.w[0] > 0x0ull)) {
281 // f* > 1/2 and the result may be exact
282 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
283 if (tmp64 > ten2mk128[ind - 1].w[1] ||
284 (tmp64 == ten2mk128[ind - 1].w[1] &&
285 fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
286 // set the inexact flag
287 *pfpsf |= INEXACT_EXCEPTION;
288 } // else the result is exact
289 } else { // the result is inexact; f2* <= 1/2
290 // set the inexact flag
291 *pfpsf |= INEXACT_EXCEPTION;
292 }
293 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
294 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
295 res.w[1] = (P256.w[3] >> shift);
296 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
297 // redundant fstar.w[3] = 0;
298 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
299 fstar.w[1] = P256.w[1];
300 fstar.w[0] = P256.w[0];
301 // fraction f* < 10^(-x) <=> midpoint
302 // f* is in the right position to be compared with
303 // 10^(-x) from ten2mk128[]
304 if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
305 fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
306 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
307 fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
308 // subract 1 to make even
309 if (res.w[0]-- == 0) {
310 res.w[1]--;
311 }
312 }
313 if (fstar.w[2] > onehalf128[ind - 1] ||
314 (fstar.w[2] == onehalf128[ind - 1]
315 && (fstar.w[1] || fstar.w[0]))) {
316 // f2* > 1/2 and the result may be exact
317 // Calculate f2* - 1/2
318 tmp64 = fstar.w[2] - onehalf128[ind - 1];
319 if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
320 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
321 fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
322 // set the inexact flag
323 *pfpsf |= INEXACT_EXCEPTION;
324 } // else the result is exact
325 } else { // the result is inexact; f2* <= 1/2
326 // set the inexact flag
327 *pfpsf |= INEXACT_EXCEPTION;
328 }
329 } else { // 22 <= ind - 1 <= 33
330 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
331 res.w[1] = 0;
332 res.w[0] = P256.w[3] >> shift;
333 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
334 fstar.w[2] = P256.w[2];
335 fstar.w[1] = P256.w[1];
336 fstar.w[0] = P256.w[0];
337 // fraction f* < 10^(-x) <=> midpoint
338 // f* is in the right position to be compared with
339 // 10^(-x) from ten2mk128[]
340 if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
341 fstar.w[3] == 0 && fstar.w[2] == 0 &&
342 (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
343 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
344 fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
345 // subract 1 to make even
346 if (res.w[0]-- == 0) {
347 res.w[1]--;
348 }
349 }
350 if (fstar.w[3] > onehalf128[ind - 1] ||
351 (fstar.w[3] == onehalf128[ind - 1] &&
352 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
353 // f2* > 1/2 and the result may be exact
354 // Calculate f2* - 1/2
355 tmp64 = fstar.w[3] - onehalf128[ind - 1];
356 if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1]
357 || (fstar.w[1] == ten2mk128[ind - 1].w[1]
358 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
359 // set the inexact flag
360 *pfpsf |= INEXACT_EXCEPTION;
361 } // else the result is exact
362 } else { // the result is inexact; f2* <= 1/2
363 // set the inexact flag
364 *pfpsf |= INEXACT_EXCEPTION;
365 }
366 }
367 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
368 BID_RETURN (res);
369 } else { // if ((q + exp) < 0) <=> q < -exp
370 // the result is +0 or -0
371 res.w[1] = x_sign | 0x3040000000000000ull;
372 res.w[0] = 0x0000000000000000ull;
373 *pfpsf |= INEXACT_EXCEPTION;
374 BID_RETURN (res);
375 }
376 break;
377 case ROUNDING_TIES_AWAY:
378 if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
379 // need to shift right -exp digits from the coefficient; exp will be 0
380 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
381 // chop off ind digits from the lower part of C1
382 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
383 tmp64 = C1.w[0];
384 if (ind <= 19) {
385 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
386 } else {
387 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
388 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
389 }
390 if (C1.w[0] < tmp64)
391 C1.w[1]++;
392 // calculate C* and f*
393 // C* is actually floor(C*) in this case
394 // C* and f* need shifting and masking, as shown by
395 // shiftright128[] and maskhigh128[]
396 // 1 <= x <= 34
397 // kx = 10^(-x) = ten2mk128[ind - 1]
398 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
399 // the approximation of 10^(-x) was rounded up to 118 bits
400 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
401 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
402 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
403 // if (0 < f* < 10^(-x)) then the result is a midpoint
404 // if floor(C*) is even then C* = floor(C*) - logical right
405 // shift; C* has p decimal digits, correct by Prop. 1)
406 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
407 // shift; C* has p decimal digits, correct by Pr. 1)
408 // else
409 // C* = floor(C*) (logical right shift; C has p decimal digits,
410 // correct by Property 1)
411 // n = C* * 10^(e+x)
412
413 // determine also the inexactness of the rounding of C*
414 // if (0 < f* - 1/2 < 10^(-x)) then
415 // the result is exact
416 // else // if (f* - 1/2 > T*) then
417 // the result is inexact
418 // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[]
419 // shift right C* by Ex-128 = shiftright128[ind]
420 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
421 // redundant shift = shiftright128[ind - 1]; // shift = 0
422 res.w[1] = P256.w[3];
423 res.w[0] = P256.w[2];
424 // redundant fstar.w[3] = 0;
425 // redundant fstar.w[2] = 0;
426 fstar.w[1] = P256.w[1];
427 fstar.w[0] = P256.w[0];
428 if (fstar.w[1] > 0x8000000000000000ull ||
429 (fstar.w[1] == 0x8000000000000000ull
430 && fstar.w[0] > 0x0ull)) {
431 // f* > 1/2 and the result may be exact
432 tmp64 = fstar.w[1] - 0x8000000000000000ull; // f* - 1/2
433 if ((tmp64 > ten2mk128[ind - 1].w[1] ||
434 (tmp64 == ten2mk128[ind - 1].w[1] &&
435 fstar.w[0] >= ten2mk128[ind - 1].w[0]))) {
436 // set the inexact flag
437 *pfpsf |= INEXACT_EXCEPTION;
438 } // else the result is exact
439 } else { // the result is inexact; f2* <= 1/2
440 // set the inexact flag
441 *pfpsf |= INEXACT_EXCEPTION;
442 }
443 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
444 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
445 res.w[1] = (P256.w[3] >> shift);
446 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
447 // redundant fstar.w[3] = 0;
448 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
449 fstar.w[1] = P256.w[1];
450 fstar.w[0] = P256.w[0];
451 if (fstar.w[2] > onehalf128[ind - 1] ||
452 (fstar.w[2] == onehalf128[ind - 1]
453 && (fstar.w[1] || fstar.w[0]))) {
454 // f2* > 1/2 and the result may be exact
455 // Calculate f2* - 1/2
456 tmp64 = fstar.w[2] - onehalf128[ind - 1];
457 if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
458 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
459 fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
460 // set the inexact flag
461 *pfpsf |= INEXACT_EXCEPTION;
462 } // else the result is exact
463 } else { // the result is inexact; f2* <= 1/2
464 // set the inexact flag
465 *pfpsf |= INEXACT_EXCEPTION;
466 }
467 } else { // 22 <= ind - 1 <= 33
468 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
469 res.w[1] = 0;
470 res.w[0] = P256.w[3] >> shift;
471 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
472 fstar.w[2] = P256.w[2];
473 fstar.w[1] = P256.w[1];
474 fstar.w[0] = P256.w[0];
475 if (fstar.w[3] > onehalf128[ind - 1] ||
476 (fstar.w[3] == onehalf128[ind - 1] &&
477 (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
478 // f2* > 1/2 and the result may be exact
479 // Calculate f2* - 1/2
480 tmp64 = fstar.w[3] - onehalf128[ind - 1];
481 if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1]
482 || (fstar.w[1] == ten2mk128[ind - 1].w[1]
483 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
484 // set the inexact flag
485 *pfpsf |= INEXACT_EXCEPTION;
486 } // else the result is exact
487 } else { // the result is inexact; f2* <= 1/2
488 // set the inexact flag
489 *pfpsf |= INEXACT_EXCEPTION;
490 }
491 }
492 // if the result was a midpoint, it was already rounded away from zero
493 res.w[1] |= x_sign | 0x3040000000000000ull;
494 BID_RETURN (res);
495 } else { // if ((q + exp) < 0) <=> q < -exp
496 // the result is +0 or -0
497 res.w[1] = x_sign | 0x3040000000000000ull;
498 res.w[0] = 0x0000000000000000ull;
499 *pfpsf |= INEXACT_EXCEPTION;
500 BID_RETURN (res);
501 }
502 break;
503 case ROUNDING_DOWN:
504 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
505 // need to shift right -exp digits from the coefficient; exp will be 0
506 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
507 // (number of digits to be chopped off)
508 // chop off ind digits from the lower part of C1
509 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
510 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
511 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
512 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
513 // tmp64 = C1.w[0];
514 // if (ind <= 19) {
515 // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
516 // } else {
517 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
518 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
519 // }
520 // if (C1.w[0] < tmp64) C1.w[1]++;
521 // if carry-out from C1.w[0], increment C1.w[1]
522 // calculate C* and f*
523 // C* is actually floor(C*) in this case
524 // C* and f* need shifting and masking, as shown by
525 // shiftright128[] and maskhigh128[]
526 // 1 <= x <= 34
527 // kx = 10^(-x) = ten2mk128[ind - 1]
528 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
529 // the approximation of 10^(-x) was rounded up to 118 bits
530 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
531 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
532 res.w[1] = P256.w[3];
533 res.w[0] = P256.w[2];
534 // redundant fstar.w[3] = 0;
535 // redundant fstar.w[2] = 0;
536 // redundant fstar.w[1] = P256.w[1];
537 // redundant fstar.w[0] = P256.w[0];
538 // fraction f* > 10^(-x) <=> inexact
539 // f* is in the right position to be compared with
540 // 10^(-x) from ten2mk128[]
541 if ((P256.w[1] > ten2mk128[ind - 1].w[1])
542 || (P256.w[1] == ten2mk128[ind - 1].w[1]
543 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
544 *pfpsf |= INEXACT_EXCEPTION;
545 // if positive, the truncated value is already the correct result
546 if (x_sign) { // if negative
547 if (++res.w[0] == 0) {
548 res.w[1]++;
549 }
550 }
551 }
552 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
553 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
554 res.w[1] = (P256.w[3] >> shift);
555 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
556 // redundant fstar.w[3] = 0;
557 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
558 fstar.w[1] = P256.w[1];
559 fstar.w[0] = P256.w[0];
560 // fraction f* > 10^(-x) <=> inexact
561 // f* is in the right position to be compared with
562 // 10^(-x) from ten2mk128[]
563 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
564 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
565 fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
566 *pfpsf |= INEXACT_EXCEPTION;
567 // if positive, the truncated value is already the correct result
568 if (x_sign) { // if negative
569 if (++res.w[0] == 0) {
570 res.w[1]++;
571 }
572 }
573 }
574 } else { // 22 <= ind - 1 <= 33
575 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
576 res.w[1] = 0;
577 res.w[0] = P256.w[3] >> shift;
578 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
579 fstar.w[2] = P256.w[2];
580 fstar.w[1] = P256.w[1];
581 fstar.w[0] = P256.w[0];
582 // fraction f* > 10^(-x) <=> inexact
583 // f* is in the right position to be compared with
584 // 10^(-x) from ten2mk128[]
585 if (fstar.w[3] || fstar.w[2]
586 || fstar.w[1] > ten2mk128[ind - 1].w[1]
587 || (fstar.w[1] == ten2mk128[ind - 1].w[1]
588 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
589 *pfpsf |= INEXACT_EXCEPTION;
590 // if positive, the truncated value is already the correct result
591 if (x_sign) { // if negative
592 if (++res.w[0] == 0) {
593 res.w[1]++;
594 }
595 }
596 }
597 }
598 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
599 BID_RETURN (res);
600 } else { // if exp < 0 and q + exp <= 0
601 if (x_sign) { // negative rounds down to -1.0
602 res.w[1] = 0xb040000000000000ull;
603 res.w[0] = 0x0000000000000001ull;
604 } else { // positive rpunds down to +0.0
605 res.w[1] = 0x3040000000000000ull;
606 res.w[0] = 0x0000000000000000ull;
607 }
608 *pfpsf |= INEXACT_EXCEPTION;
609 BID_RETURN (res);
610 }
611 break;
612 case ROUNDING_UP:
613 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
614 // need to shift right -exp digits from the coefficient; exp will be 0
615 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
616 // (number of digits to be chopped off)
617 // chop off ind digits from the lower part of C1
618 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
619 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
620 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
621 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
622 // tmp64 = C1.w[0];
623 // if (ind <= 19) {
624 // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
625 // } else {
626 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
627 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
628 // }
629 // if (C1.w[0] < tmp64) C1.w[1]++;
630 // if carry-out from C1.w[0], increment C1.w[1]
631 // calculate C* and f*
632 // C* is actually floor(C*) in this case
633 // C* and f* need shifting and masking, as shown by
634 // shiftright128[] and maskhigh128[]
635 // 1 <= x <= 34
636 // kx = 10^(-x) = ten2mk128[ind - 1]
637 // C* = C1 * 10^(-x)
638 // the approximation of 10^(-x) was rounded up to 118 bits
639 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
640 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
641 res.w[1] = P256.w[3];
642 res.w[0] = P256.w[2];
643 // redundant fstar.w[3] = 0;
644 // redundant fstar.w[2] = 0;
645 // redundant fstar.w[1] = P256.w[1];
646 // redundant fstar.w[0] = P256.w[0];
647 // fraction f* > 10^(-x) <=> inexact
648 // f* is in the right position to be compared with
649 // 10^(-x) from ten2mk128[]
650 if ((P256.w[1] > ten2mk128[ind - 1].w[1])
651 || (P256.w[1] == ten2mk128[ind - 1].w[1]
652 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
653 *pfpsf |= INEXACT_EXCEPTION;
654 // if negative, the truncated value is already the correct result
655 if (!x_sign) { // if positive
656 if (++res.w[0] == 0) {
657 res.w[1]++;
658 }
659 }
660 }
661 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
662 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
663 res.w[1] = (P256.w[3] >> shift);
664 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
665 // redundant fstar.w[3] = 0;
666 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
667 fstar.w[1] = P256.w[1];
668 fstar.w[0] = P256.w[0];
669 // fraction f* > 10^(-x) <=> inexact
670 // f* is in the right position to be compared with
671 // 10^(-x) from ten2mk128[]
672 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
673 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
674 fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
675 *pfpsf |= INEXACT_EXCEPTION;
676 // if negative, the truncated value is already the correct result
677 if (!x_sign) { // if positive
678 if (++res.w[0] == 0) {
679 res.w[1]++;
680 }
681 }
682 }
683 } else { // 22 <= ind - 1 <= 33
684 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
685 res.w[1] = 0;
686 res.w[0] = P256.w[3] >> shift;
687 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
688 fstar.w[2] = P256.w[2];
689 fstar.w[1] = P256.w[1];
690 fstar.w[0] = P256.w[0];
691 // fraction f* > 10^(-x) <=> inexact
692 // f* is in the right position to be compared with
693 // 10^(-x) from ten2mk128[]
694 if (fstar.w[3] || fstar.w[2]
695 || fstar.w[1] > ten2mk128[ind - 1].w[1]
696 || (fstar.w[1] == ten2mk128[ind - 1].w[1]
697 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
698 *pfpsf |= INEXACT_EXCEPTION;
699 // if negative, the truncated value is already the correct result
700 if (!x_sign) { // if positive
701 if (++res.w[0] == 0) {
702 res.w[1]++;
703 }
704 }
705 }
706 }
707 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
708 BID_RETURN (res);
709 } else { // if exp < 0 and q + exp <= 0
710 if (x_sign) { // negative rounds up to -0.0
711 res.w[1] = 0xb040000000000000ull;
712 res.w[0] = 0x0000000000000000ull;
713 } else { // positive rpunds up to +1.0
714 res.w[1] = 0x3040000000000000ull;
715 res.w[0] = 0x0000000000000001ull;
716 }
717 *pfpsf |= INEXACT_EXCEPTION;
718 BID_RETURN (res);
719 }
720 break;
721 case ROUNDING_TO_ZERO:
722 if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
723 // need to shift right -exp digits from the coefficient; exp will be 0
724 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
725 // (number of digits to be chopped off)
726 // chop off ind digits from the lower part of C1
727 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
728 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
729 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
730 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
731 //tmp64 = C1.w[0];
732 // if (ind <= 19) {
733 // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
734 // } else {
735 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
736 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
737 // }
738 // if (C1.w[0] < tmp64) C1.w[1]++;
739 // if carry-out from C1.w[0], increment C1.w[1]
740 // calculate C* and f*
741 // C* is actually floor(C*) in this case
742 // C* and f* need shifting and masking, as shown by
743 // shiftright128[] and maskhigh128[]
744 // 1 <= x <= 34
745 // kx = 10^(-x) = ten2mk128[ind - 1]
746 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
747 // the approximation of 10^(-x) was rounded up to 118 bits
748 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
749 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
750 res.w[1] = P256.w[3];
751 res.w[0] = P256.w[2];
752 // redundant fstar.w[3] = 0;
753 // redundant fstar.w[2] = 0;
754 // redundant fstar.w[1] = P256.w[1];
755 // redundant fstar.w[0] = P256.w[0];
756 // fraction f* > 10^(-x) <=> inexact
757 // f* is in the right position to be compared with
758 // 10^(-x) from ten2mk128[]
759 if ((P256.w[1] > ten2mk128[ind - 1].w[1])
760 || (P256.w[1] == ten2mk128[ind - 1].w[1]
761 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
762 *pfpsf |= INEXACT_EXCEPTION;
763 }
764 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
765 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
766 res.w[1] = (P256.w[3] >> shift);
767 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
768 // redundant fstar.w[3] = 0;
769 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
770 fstar.w[1] = P256.w[1];
771 fstar.w[0] = P256.w[0];
772 // fraction f* > 10^(-x) <=> inexact
773 // f* is in the right position to be compared with
774 // 10^(-x) from ten2mk128[]
775 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
776 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
777 fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
778 *pfpsf |= INEXACT_EXCEPTION;
779 }
780 } else { // 22 <= ind - 1 <= 33
781 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
782 res.w[1] = 0;
783 res.w[0] = P256.w[3] >> shift;
784 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
785 fstar.w[2] = P256.w[2];
786 fstar.w[1] = P256.w[1];
787 fstar.w[0] = P256.w[0];
788 // fraction f* > 10^(-x) <=> inexact
789 // f* is in the right position to be compared with
790 // 10^(-x) from ten2mk128[]
791 if (fstar.w[3] || fstar.w[2]
792 || fstar.w[1] > ten2mk128[ind - 1].w[1]
793 || (fstar.w[1] == ten2mk128[ind - 1].w[1]
794 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
795 *pfpsf |= INEXACT_EXCEPTION;
796 }
797 }
798 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
799 BID_RETURN (res);
800 } else { // if exp < 0 and q + exp <= 0 the result is +0 or -0
801 res.w[1] = x_sign | 0x3040000000000000ull;
802 res.w[0] = 0x0000000000000000ull;
803 *pfpsf |= INEXACT_EXCEPTION;
804 BID_RETURN (res);
805 }
806 break;
807 }
808
809 BID_RETURN (res);
810 }
811
812 /*****************************************************************************
813 * BID128_round_integral_nearest_even
814 ****************************************************************************/
815
816 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_even, x)
817
818 UINT128 res;
819 UINT64 x_sign;
820 UINT64 x_exp;
821 int exp; // unbiased exponent
822 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
823 UINT64 tmp64;
824 BID_UI64DOUBLE tmp1;
825 unsigned int x_nr_bits;
826 int q, ind, shift;
827 UINT128 C1;
828 // UINT128 res is C* at first - represents up to 34 decimal digits ~ 113 bits
829 UINT256 fstar;
830 UINT256 P256;
831
832 // check for NaN or Infinity
833 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
834 // x is special
835 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
836 // if x = NaN, then res = Q (x)
837 // check first for non-canonical NaN payload
838 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
839 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
840 (x.w[0] > 0x38c15b09ffffffffull))) {
841 x.w[1] = x.w[1] & 0xffffc00000000000ull;
842 x.w[0] = 0x0ull;
843 }
844 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
845 // set invalid flag
846 *pfpsf |= INVALID_EXCEPTION;
847 // return quiet (x)
848 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
849 res.w[0] = x.w[0];
850 } else { // x is QNaN
851 // return x
852 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
853 res.w[0] = x.w[0];
854 }
855 BID_RETURN (res)
856 } else { // x is not a NaN, so it must be infinity
857 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
858 // return +inf
859 res.w[1] = 0x7800000000000000ull;
860 res.w[0] = 0x0000000000000000ull;
861 } else { // x is -inf
862 // return -inf
863 res.w[1] = 0xf800000000000000ull;
864 res.w[0] = 0x0000000000000000ull;
865 }
866 BID_RETURN (res);
867 }
868 }
869 // unpack x
870 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
871 C1.w[1] = x.w[1] & MASK_COEFF;
872 C1.w[0] = x.w[0];
873
874 // check for non-canonical values (treated as zero)
875 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
876 // non-canonical
877 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
878 C1.w[1] = 0; // significand high
879 C1.w[0] = 0; // significand low
880 } else { // G0_G1 != 11
881 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
882 if (C1.w[1] > 0x0001ed09bead87c0ull ||
883 (C1.w[1] == 0x0001ed09bead87c0ull
884 && C1.w[0] > 0x378d8e63ffffffffull)) {
885 // x is non-canonical if coefficient is larger than 10^34 -1
886 C1.w[1] = 0;
887 C1.w[0] = 0;
888 } else { // canonical
889 ;
890 }
891 }
892
893 // test for input equal to zero
894 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
895 // x is 0
896 // return 0 preserving the sign bit and the preferred exponent
897 // of MAX(Q(x), 0)
898 if (x_exp <= (0x1820ull << 49)) {
899 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
900 } else {
901 res.w[1] = x_sign | x_exp;
902 }
903 res.w[0] = 0x0000000000000000ull;
904 BID_RETURN (res);
905 }
906 // x is not special and is not zero
907
908 // if (exp <= -(p+1)) return 0
909 if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35
910 res.w[1] = x_sign | 0x3040000000000000ull;
911 res.w[0] = 0x0000000000000000ull;
912 BID_RETURN (res);
913 }
914 // q = nr. of decimal digits in x
915 // determine first the nr. of bits in x
916 if (C1.w[1] == 0) {
917 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
918 // split the 64-bit value in two 32-bit halves to avoid rounding errors
919 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
920 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
921 x_nr_bits =
922 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
923 } else { // x < 2^32
924 tmp1.d = (double) (C1.w[0]); // exact conversion
925 x_nr_bits =
926 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
927 }
928 } else { // if x < 2^53
929 tmp1.d = (double) C1.w[0]; // exact conversion
930 x_nr_bits =
931 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
932 }
933 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
934 tmp1.d = (double) C1.w[1]; // exact conversion
935 x_nr_bits =
936 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
937 }
938
939 q = nr_digits[x_nr_bits - 1].digits;
940 if (q == 0) {
941 q = nr_digits[x_nr_bits - 1].digits1;
942 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
943 || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
944 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
945 q++;
946 }
947 exp = (x_exp >> 49) - 6176;
948 if (exp >= 0) { // -exp <= 0
949 // the argument is an integer already
950 res.w[1] = x.w[1];
951 res.w[0] = x.w[0];
952 BID_RETURN (res);
953 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
954 // need to shift right -exp digits from the coefficient; the exp will be 0
955 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
956 // chop off ind digits from the lower part of C1
957 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
958 tmp64 = C1.w[0];
959 if (ind <= 19) {
960 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
961 } else {
962 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
963 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
964 }
965 if (C1.w[0] < tmp64)
966 C1.w[1]++;
967 // calculate C* and f*
968 // C* is actually floor(C*) in this case
969 // C* and f* need shifting and masking, as shown by
970 // shiftright128[] and maskhigh128[]
971 // 1 <= x <= 34
972 // kx = 10^(-x) = ten2mk128[ind - 1]
973 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
974 // the approximation of 10^(-x) was rounded up to 118 bits
975 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
976 // determine the value of res and fstar
977 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
978 // redundant shift = shiftright128[ind - 1]; // shift = 0
979 res.w[1] = P256.w[3];
980 res.w[0] = P256.w[2];
981 // redundant fstar.w[3] = 0;
982 // redundant fstar.w[2] = 0;
983 // redundant fstar.w[1] = P256.w[1];
984 // redundant fstar.w[0] = P256.w[0];
985 // fraction f* < 10^(-x) <=> midpoint
986 // f* is in the right position to be compared with
987 // 10^(-x) from ten2mk128[]
988 // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even)
989 if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
990 ((P256.w[1] < (ten2mk128[ind - 1].w[1]))
991 || ((P256.w[1] == ten2mk128[ind - 1].w[1])
992 && (P256.w[0] < ten2mk128[ind - 1].w[0])))) {
993 // subract 1 to make even
994 if (res.w[0]-- == 0) {
995 res.w[1]--;
996 }
997 }
998 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
999 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
1000 res.w[1] = (P256.w[3] >> shift);
1001 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1002 // redundant fstar.w[3] = 0;
1003 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1004 fstar.w[1] = P256.w[1];
1005 fstar.w[0] = P256.w[0];
1006 // fraction f* < 10^(-x) <=> midpoint
1007 // f* is in the right position to be compared with
1008 // 10^(-x) from ten2mk128[]
1009 if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
1010 fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
1011 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1012 fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
1013 // subract 1 to make even
1014 if (res.w[0]-- == 0) {
1015 res.w[1]--;
1016 }
1017 }
1018 } else { // 22 <= ind - 1 <= 33
1019 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
1020 res.w[1] = 0;
1021 res.w[0] = P256.w[3] >> shift;
1022 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1023 fstar.w[2] = P256.w[2];
1024 fstar.w[1] = P256.w[1];
1025 fstar.w[0] = P256.w[0];
1026 // fraction f* < 10^(-x) <=> midpoint
1027 // f* is in the right position to be compared with
1028 // 10^(-x) from ten2mk128[]
1029 if ((res.w[0] & 0x0000000000000001ull) && // is result odd?
1030 fstar.w[3] == 0 && fstar.w[2] == 0
1031 && (fstar.w[1] < ten2mk128[ind - 1].w[1]
1032 || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1033 && fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
1034 // subract 1 to make even
1035 if (res.w[0]-- == 0) {
1036 res.w[1]--;
1037 }
1038 }
1039 }
1040 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1041 BID_RETURN (res);
1042 } else { // if ((q + exp) < 0) <=> q < -exp
1043 // the result is +0 or -0
1044 res.w[1] = x_sign | 0x3040000000000000ull;
1045 res.w[0] = 0x0000000000000000ull;
1046 BID_RETURN (res);
1047 }
1048 }
1049
1050 /*****************************************************************************
1051 * BID128_round_integral_negative
1052 ****************************************************************************/
1053
1054 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_negative, x)
1055
1056 UINT128 res;
1057 UINT64 x_sign;
1058 UINT64 x_exp;
1059 int exp; // unbiased exponent
1060 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1061 // (all are UINT64)
1062 BID_UI64DOUBLE tmp1;
1063 unsigned int x_nr_bits;
1064 int q, ind, shift;
1065 UINT128 C1;
1066 // UINT128 res is C* at first - represents up to 34 decimal digits ~
1067 // 113 bits
1068 UINT256 fstar;
1069 UINT256 P256;
1070
1071 // check for NaN or Infinity
1072 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1073 // x is special
1074 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
1075 // if x = NaN, then res = Q (x)
1076 // check first for non-canonical NaN payload
1077 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1078 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1079 (x.w[0] > 0x38c15b09ffffffffull))) {
1080 x.w[1] = x.w[1] & 0xffffc00000000000ull;
1081 x.w[0] = 0x0ull;
1082 }
1083 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
1084 // set invalid flag
1085 *pfpsf |= INVALID_EXCEPTION;
1086 // return quiet (x)
1087 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
1088 res.w[0] = x.w[0];
1089 } else { // x is QNaN
1090 // return x
1091 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
1092 res.w[0] = x.w[0];
1093 }
1094 BID_RETURN (res)
1095 } else { // x is not a NaN, so it must be infinity
1096 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1097 // return +inf
1098 res.w[1] = 0x7800000000000000ull;
1099 res.w[0] = 0x0000000000000000ull;
1100 } else { // x is -inf
1101 // return -inf
1102 res.w[1] = 0xf800000000000000ull;
1103 res.w[0] = 0x0000000000000000ull;
1104 }
1105 BID_RETURN (res);
1106 }
1107 }
1108 // unpack x
1109 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1110 C1.w[1] = x.w[1] & MASK_COEFF;
1111 C1.w[0] = x.w[0];
1112
1113 // check for non-canonical values (treated as zero)
1114 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
1115 // non-canonical
1116 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
1117 C1.w[1] = 0; // significand high
1118 C1.w[0] = 0; // significand low
1119 } else { // G0_G1 != 11
1120 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
1121 if (C1.w[1] > 0x0001ed09bead87c0ull ||
1122 (C1.w[1] == 0x0001ed09bead87c0ull
1123 && C1.w[0] > 0x378d8e63ffffffffull)) {
1124 // x is non-canonical if coefficient is larger than 10^34 -1
1125 C1.w[1] = 0;
1126 C1.w[0] = 0;
1127 } else { // canonical
1128 ;
1129 }
1130 }
1131
1132 // test for input equal to zero
1133 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1134 // x is 0
1135 // return 0 preserving the sign bit and the preferred exponent
1136 // of MAX(Q(x), 0)
1137 if (x_exp <= (0x1820ull << 49)) {
1138 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1139 } else {
1140 res.w[1] = x_sign | x_exp;
1141 }
1142 res.w[0] = 0x0000000000000000ull;
1143 BID_RETURN (res);
1144 }
1145 // x is not special and is not zero
1146
1147 // if (exp <= -p) return -1.0 or +0.0
1148 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
1149 if (x_sign) {
1150 // if negative, return negative 1, because we know the coefficient
1151 // is non-zero (would have been caught above)
1152 res.w[1] = 0xb040000000000000ull;
1153 res.w[0] = 0x0000000000000001ull;
1154 } else {
1155 // if positive, return positive 0, because we know coefficient is
1156 // non-zero (would have been caught above)
1157 res.w[1] = 0x3040000000000000ull;
1158 res.w[0] = 0x0000000000000000ull;
1159 }
1160 BID_RETURN (res);
1161 }
1162 // q = nr. of decimal digits in x
1163 // determine first the nr. of bits in x
1164 if (C1.w[1] == 0) {
1165 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
1166 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1167 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
1168 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1169 x_nr_bits =
1170 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1171 } else { // x < 2^32
1172 tmp1.d = (double) (C1.w[0]); // exact conversion
1173 x_nr_bits =
1174 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1175 }
1176 } else { // if x < 2^53
1177 tmp1.d = (double) C1.w[0]; // exact conversion
1178 x_nr_bits =
1179 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1180 }
1181 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1182 tmp1.d = (double) C1.w[1]; // exact conversion
1183 x_nr_bits =
1184 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1185 }
1186
1187 q = nr_digits[x_nr_bits - 1].digits;
1188 if (q == 0) {
1189 q = nr_digits[x_nr_bits - 1].digits1;
1190 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1191 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1192 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1193 q++;
1194 }
1195 exp = (x_exp >> 49) - 6176;
1196 if (exp >= 0) { // -exp <= 0
1197 // the argument is an integer already
1198 res.w[1] = x.w[1];
1199 res.w[0] = x.w[0];
1200 BID_RETURN (res);
1201 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
1202 // need to shift right -exp digits from the coefficient; the exp will be 0
1203 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
1204 // (number of digits to be chopped off)
1205 // chop off ind digits from the lower part of C1
1206 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1207 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1208 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1209 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1210 //tmp64 = C1.w[0];
1211 // if (ind <= 19) {
1212 // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1213 // } else {
1214 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1215 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1216 // }
1217 // if (C1.w[0] < tmp64) C1.w[1]++;
1218 // if carry-out from C1.w[0], increment C1.w[1]
1219 // calculate C* and f*
1220 // C* is actually floor(C*) in this case
1221 // C* and f* need shifting and masking, as shown by
1222 // shiftright128[] and maskhigh128[]
1223 // 1 <= x <= 34
1224 // kx = 10^(-x) = ten2mk128[ind - 1]
1225 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1226 // the approximation of 10^(-x) was rounded up to 118 bits
1227 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1228 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
1229 res.w[1] = P256.w[3];
1230 res.w[0] = P256.w[2];
1231 // if positive, the truncated value is already the correct result
1232 if (x_sign) { // if negative
1233 // redundant fstar.w[3] = 0;
1234 // redundant fstar.w[2] = 0;
1235 // redundant fstar.w[1] = P256.w[1];
1236 // redundant fstar.w[0] = P256.w[0];
1237 // fraction f* > 10^(-x) <=> inexact
1238 // f* is in the right position to be compared with
1239 // 10^(-x) from ten2mk128[]
1240 if ((P256.w[1] > ten2mk128[ind - 1].w[1])
1241 || (P256.w[1] == ten2mk128[ind - 1].w[1]
1242 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
1243 if (++res.w[0] == 0) {
1244 res.w[1]++;
1245 }
1246 }
1247 }
1248 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1249 shift = shiftright128[ind - 1]; // 0 <= shift <= 102
1250 res.w[1] = (P256.w[3] >> shift);
1251 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1252 // if positive, the truncated value is already the correct result
1253 if (x_sign) { // if negative
1254 // redundant fstar.w[3] = 0;
1255 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1256 fstar.w[1] = P256.w[1];
1257 fstar.w[0] = P256.w[0];
1258 // fraction f* > 10^(-x) <=> inexact
1259 // f* is in the right position to be compared with
1260 // 10^(-x) from ten2mk128[]
1261 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
1262 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1263 fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1264 if (++res.w[0] == 0) {
1265 res.w[1]++;
1266 }
1267 }
1268 }
1269 } else { // 22 <= ind - 1 <= 33
1270 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
1271 res.w[1] = 0;
1272 res.w[0] = P256.w[3] >> shift;
1273 // if positive, the truncated value is already the correct result
1274 if (x_sign) { // if negative
1275 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1276 fstar.w[2] = P256.w[2];
1277 fstar.w[1] = P256.w[1];
1278 fstar.w[0] = P256.w[0];
1279 // fraction f* > 10^(-x) <=> inexact
1280 // f* is in the right position to be compared with
1281 // 10^(-x) from ten2mk128[]
1282 if (fstar.w[3] || fstar.w[2]
1283 || fstar.w[1] > ten2mk128[ind - 1].w[1]
1284 || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1285 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1286 if (++res.w[0] == 0) {
1287 res.w[1]++;
1288 }
1289 }
1290 }
1291 }
1292 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1293 BID_RETURN (res);
1294 } else { // if exp < 0 and q + exp <= 0
1295 if (x_sign) { // negative rounds down to -1.0
1296 res.w[1] = 0xb040000000000000ull;
1297 res.w[0] = 0x0000000000000001ull;
1298 } else { // positive rpunds down to +0.0
1299 res.w[1] = 0x3040000000000000ull;
1300 res.w[0] = 0x0000000000000000ull;
1301 }
1302 BID_RETURN (res);
1303 }
1304 }
1305
1306 /*****************************************************************************
1307 * BID128_round_integral_positive
1308 ****************************************************************************/
1309
1310 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_positive, x)
1311
1312 UINT128 res;
1313 UINT64 x_sign;
1314 UINT64 x_exp;
1315 int exp; // unbiased exponent
1316 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1317 // (all are UINT64)
1318 BID_UI64DOUBLE tmp1;
1319 unsigned int x_nr_bits;
1320 int q, ind, shift;
1321 UINT128 C1;
1322 // UINT128 res is C* at first - represents up to 34 decimal digits ~
1323 // 113 bits
1324 UINT256 fstar;
1325 UINT256 P256;
1326
1327 // check for NaN or Infinity
1328 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1329 // x is special
1330 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
1331 // if x = NaN, then res = Q (x)
1332 // check first for non-canonical NaN payload
1333 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1334 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1335 (x.w[0] > 0x38c15b09ffffffffull))) {
1336 x.w[1] = x.w[1] & 0xffffc00000000000ull;
1337 x.w[0] = 0x0ull;
1338 }
1339 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
1340 // set invalid flag
1341 *pfpsf |= INVALID_EXCEPTION;
1342 // return quiet (x)
1343 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
1344 res.w[0] = x.w[0];
1345 } else { // x is QNaN
1346 // return x
1347 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
1348 res.w[0] = x.w[0];
1349 }
1350 BID_RETURN (res)
1351 } else { // x is not a NaN, so it must be infinity
1352 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1353 // return +inf
1354 res.w[1] = 0x7800000000000000ull;
1355 res.w[0] = 0x0000000000000000ull;
1356 } else { // x is -inf
1357 // return -inf
1358 res.w[1] = 0xf800000000000000ull;
1359 res.w[0] = 0x0000000000000000ull;
1360 }
1361 BID_RETURN (res);
1362 }
1363 }
1364 // unpack x
1365 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1366 C1.w[1] = x.w[1] & MASK_COEFF;
1367 C1.w[0] = x.w[0];
1368
1369 // check for non-canonical values (treated as zero)
1370 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
1371 // non-canonical
1372 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
1373 C1.w[1] = 0; // significand high
1374 C1.w[0] = 0; // significand low
1375 } else { // G0_G1 != 11
1376 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
1377 if (C1.w[1] > 0x0001ed09bead87c0ull ||
1378 (C1.w[1] == 0x0001ed09bead87c0ull
1379 && C1.w[0] > 0x378d8e63ffffffffull)) {
1380 // x is non-canonical if coefficient is larger than 10^34 -1
1381 C1.w[1] = 0;
1382 C1.w[0] = 0;
1383 } else { // canonical
1384 ;
1385 }
1386 }
1387
1388 // test for input equal to zero
1389 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1390 // x is 0
1391 // return 0 preserving the sign bit and the preferred exponent
1392 // of MAX(Q(x), 0)
1393 if (x_exp <= (0x1820ull << 49)) {
1394 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1395 } else {
1396 res.w[1] = x_sign | x_exp;
1397 }
1398 res.w[0] = 0x0000000000000000ull;
1399 BID_RETURN (res);
1400 }
1401 // x is not special and is not zero
1402
1403 // if (exp <= -p) return -0.0 or +1.0
1404 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
1405 if (x_sign) {
1406 // if negative, return negative 0, because we know the coefficient
1407 // is non-zero (would have been caught above)
1408 res.w[1] = 0xb040000000000000ull;
1409 res.w[0] = 0x0000000000000000ull;
1410 } else {
1411 // if positive, return positive 1, because we know coefficient is
1412 // non-zero (would have been caught above)
1413 res.w[1] = 0x3040000000000000ull;
1414 res.w[0] = 0x0000000000000001ull;
1415 }
1416 BID_RETURN (res);
1417 }
1418 // q = nr. of decimal digits in x
1419 // determine first the nr. of bits in x
1420 if (C1.w[1] == 0) {
1421 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
1422 // split 64-bit value in two 32-bit halves to avoid rounding errors
1423 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
1424 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1425 x_nr_bits =
1426 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1427 } else { // x < 2^32
1428 tmp1.d = (double) (C1.w[0]); // exact conversion
1429 x_nr_bits =
1430 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1431 }
1432 } else { // if x < 2^53
1433 tmp1.d = (double) C1.w[0]; // exact conversion
1434 x_nr_bits =
1435 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1436 }
1437 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1438 tmp1.d = (double) C1.w[1]; // exact conversion
1439 x_nr_bits =
1440 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1441 }
1442
1443 q = nr_digits[x_nr_bits - 1].digits;
1444 if (q == 0) {
1445 q = nr_digits[x_nr_bits - 1].digits1;
1446 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1447 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1448 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1449 q++;
1450 }
1451 exp = (x_exp >> 49) - 6176;
1452 if (exp >= 0) { // -exp <= 0
1453 // the argument is an integer already
1454 res.w[1] = x.w[1];
1455 res.w[0] = x.w[0];
1456 BID_RETURN (res);
1457 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
1458 // need to shift right -exp digits from the coefficient; exp will be 0
1459 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
1460 // (number of digits to be chopped off)
1461 // chop off ind digits from the lower part of C1
1462 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1463 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1464 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1465 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1466 // tmp64 = C1.w[0];
1467 // if (ind <= 19) {
1468 // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1469 // } else {
1470 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1471 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1472 // }
1473 // if (C1.w[0] < tmp64) C1.w[1]++;
1474 // if carry-out from C1.w[0], increment C1.w[1]
1475 // calculate C* and f*
1476 // C* is actually floor(C*) in this case
1477 // C* and f* need shifting and masking, as shown by
1478 // shiftright128[] and maskhigh128[]
1479 // 1 <= x <= 34
1480 // kx = 10^(-x) = ten2mk128[ind - 1]
1481 // C* = C1 * 10^(-x)
1482 // the approximation of 10^(-x) was rounded up to 118 bits
1483 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1484 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
1485 res.w[1] = P256.w[3];
1486 res.w[0] = P256.w[2];
1487 // if negative, the truncated value is already the correct result
1488 if (!x_sign) { // if positive
1489 // redundant fstar.w[3] = 0;
1490 // redundant fstar.w[2] = 0;
1491 // redundant fstar.w[1] = P256.w[1];
1492 // redundant fstar.w[0] = P256.w[0];
1493 // fraction f* > 10^(-x) <=> inexact
1494 // f* is in the right position to be compared with
1495 // 10^(-x) from ten2mk128[]
1496 if ((P256.w[1] > ten2mk128[ind - 1].w[1])
1497 || (P256.w[1] == ten2mk128[ind - 1].w[1]
1498 && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
1499 if (++res.w[0] == 0) {
1500 res.w[1]++;
1501 }
1502 }
1503 }
1504 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1505 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
1506 res.w[1] = (P256.w[3] >> shift);
1507 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1508 // if negative, the truncated value is already the correct result
1509 if (!x_sign) { // if positive
1510 // redundant fstar.w[3] = 0;
1511 fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1512 fstar.w[1] = P256.w[1];
1513 fstar.w[0] = P256.w[0];
1514 // fraction f* > 10^(-x) <=> inexact
1515 // f* is in the right position to be compared with
1516 // 10^(-x) from ten2mk128[]
1517 if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
1518 (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1519 fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1520 if (++res.w[0] == 0) {
1521 res.w[1]++;
1522 }
1523 }
1524 }
1525 } else { // 22 <= ind - 1 <= 33
1526 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
1527 res.w[1] = 0;
1528 res.w[0] = P256.w[3] >> shift;
1529 // if negative, the truncated value is already the correct result
1530 if (!x_sign) { // if positive
1531 fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1532 fstar.w[2] = P256.w[2];
1533 fstar.w[1] = P256.w[1];
1534 fstar.w[0] = P256.w[0];
1535 // fraction f* > 10^(-x) <=> inexact
1536 // f* is in the right position to be compared with
1537 // 10^(-x) from ten2mk128[]
1538 if (fstar.w[3] || fstar.w[2]
1539 || fstar.w[1] > ten2mk128[ind - 1].w[1]
1540 || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1541 && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1542 if (++res.w[0] == 0) {
1543 res.w[1]++;
1544 }
1545 }
1546 }
1547 }
1548 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1549 BID_RETURN (res);
1550 } else { // if exp < 0 and q + exp <= 0
1551 if (x_sign) { // negative rounds up to -0.0
1552 res.w[1] = 0xb040000000000000ull;
1553 res.w[0] = 0x0000000000000000ull;
1554 } else { // positive rpunds up to +1.0
1555 res.w[1] = 0x3040000000000000ull;
1556 res.w[0] = 0x0000000000000001ull;
1557 }
1558 BID_RETURN (res);
1559 }
1560 }
1561
1562 /*****************************************************************************
1563 * BID128_round_integral_zero
1564 ****************************************************************************/
1565
1566 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_zero, x)
1567
1568 UINT128 res;
1569 UINT64 x_sign;
1570 UINT64 x_exp;
1571 int exp; // unbiased exponent
1572 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1573 // (all are UINT64)
1574 BID_UI64DOUBLE tmp1;
1575 unsigned int x_nr_bits;
1576 int q, ind, shift;
1577 UINT128 C1;
1578 // UINT128 res is C* at first - represents up to 34 decimal digits ~
1579 // 113 bits
1580 UINT256 P256;
1581
1582 // check for NaN or Infinity
1583 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1584 // x is special
1585 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
1586 // if x = NaN, then res = Q (x)
1587 // check first for non-canonical NaN payload
1588 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1589 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1590 (x.w[0] > 0x38c15b09ffffffffull))) {
1591 x.w[1] = x.w[1] & 0xffffc00000000000ull;
1592 x.w[0] = 0x0ull;
1593 }
1594 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
1595 // set invalid flag
1596 *pfpsf |= INVALID_EXCEPTION;
1597 // return quiet (x)
1598 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
1599 res.w[0] = x.w[0];
1600 } else { // x is QNaN
1601 // return x
1602 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
1603 res.w[0] = x.w[0];
1604 }
1605 BID_RETURN (res)
1606 } else { // x is not a NaN, so it must be infinity
1607 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1608 // return +inf
1609 res.w[1] = 0x7800000000000000ull;
1610 res.w[0] = 0x0000000000000000ull;
1611 } else { // x is -inf
1612 // return -inf
1613 res.w[1] = 0xf800000000000000ull;
1614 res.w[0] = 0x0000000000000000ull;
1615 }
1616 BID_RETURN (res);
1617 }
1618 }
1619 // unpack x
1620 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1621 C1.w[1] = x.w[1] & MASK_COEFF;
1622 C1.w[0] = x.w[0];
1623
1624 // check for non-canonical values (treated as zero)
1625 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
1626 // non-canonical
1627 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
1628 C1.w[1] = 0; // significand high
1629 C1.w[0] = 0; // significand low
1630 } else { // G0_G1 != 11
1631 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
1632 if (C1.w[1] > 0x0001ed09bead87c0ull ||
1633 (C1.w[1] == 0x0001ed09bead87c0ull
1634 && C1.w[0] > 0x378d8e63ffffffffull)) {
1635 // x is non-canonical if coefficient is larger than 10^34 -1
1636 C1.w[1] = 0;
1637 C1.w[0] = 0;
1638 } else { // canonical
1639 ;
1640 }
1641 }
1642
1643 // test for input equal to zero
1644 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1645 // x is 0
1646 // return 0 preserving the sign bit and the preferred exponent
1647 // of MAX(Q(x), 0)
1648 if (x_exp <= (0x1820ull << 49)) {
1649 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1650 } else {
1651 res.w[1] = x_sign | x_exp;
1652 }
1653 res.w[0] = 0x0000000000000000ull;
1654 BID_RETURN (res);
1655 }
1656 // x is not special and is not zero
1657
1658 // if (exp <= -p) return -0.0 or +0.0
1659 if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
1660 res.w[1] = x_sign | 0x3040000000000000ull;
1661 res.w[0] = 0x0000000000000000ull;
1662 BID_RETURN (res);
1663 }
1664 // q = nr. of decimal digits in x
1665 // determine first the nr. of bits in x
1666 if (C1.w[1] == 0) {
1667 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
1668 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1669 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
1670 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1671 x_nr_bits =
1672 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1673 } else { // x < 2^32
1674 tmp1.d = (double) (C1.w[0]); // exact conversion
1675 x_nr_bits =
1676 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1677 }
1678 } else { // if x < 2^53
1679 tmp1.d = (double) C1.w[0]; // exact conversion
1680 x_nr_bits =
1681 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1682 }
1683 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1684 tmp1.d = (double) C1.w[1]; // exact conversion
1685 x_nr_bits =
1686 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1687 }
1688
1689 q = nr_digits[x_nr_bits - 1].digits;
1690 if (q == 0) {
1691 q = nr_digits[x_nr_bits - 1].digits1;
1692 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1693 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1694 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1695 q++;
1696 }
1697 exp = (x_exp >> 49) - 6176;
1698 if (exp >= 0) { // -exp <= 0
1699 // the argument is an integer already
1700 res.w[1] = x.w[1];
1701 res.w[0] = x.w[0];
1702 BID_RETURN (res);
1703 } else if ((q + exp) > 0) { // exp < 0 and 1 <= -exp < q
1704 // need to shift right -exp digits from the coefficient; the exp will be 0
1705 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
1706 // (number of digits to be chopped off)
1707 // chop off ind digits from the lower part of C1
1708 // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1709 // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1710 // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1711 // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1712 //tmp64 = C1.w[0];
1713 // if (ind <= 19) {
1714 // C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1715 // } else {
1716 // C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1717 // C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1718 // }
1719 // if (C1.w[0] < tmp64) C1.w[1]++;
1720 // if carry-out from C1.w[0], increment C1.w[1]
1721 // calculate C* and f*
1722 // C* is actually floor(C*) in this case
1723 // C* and f* need shifting and masking, as shown by
1724 // shiftright128[] and maskhigh128[]
1725 // 1 <= x <= 34
1726 // kx = 10^(-x) = ten2mk128[ind - 1]
1727 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1728 // the approximation of 10^(-x) was rounded up to 118 bits
1729 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1730 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
1731 res.w[1] = P256.w[3];
1732 res.w[0] = P256.w[2];
1733 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1734 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
1735 res.w[1] = (P256.w[3] >> shift);
1736 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1737 } else { // 22 <= ind - 1 <= 33
1738 shift = shiftright128[ind - 1] - 64; // 2 <= shift <= 38
1739 res.w[1] = 0;
1740 res.w[0] = P256.w[3] >> shift;
1741 }
1742 res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1743 BID_RETURN (res);
1744 } else { // if exp < 0 and q + exp <= 0 the result is +0 or -0
1745 res.w[1] = x_sign | 0x3040000000000000ull;
1746 res.w[0] = 0x0000000000000000ull;
1747 BID_RETURN (res);
1748 }
1749 }
1750
1751 /*****************************************************************************
1752 * BID128_round_integral_nearest_away
1753 ****************************************************************************/
1754
1755 BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_away, x)
1756
1757 UINT128 res;
1758 UINT64 x_sign;
1759 UINT64 x_exp;
1760 int exp; // unbiased exponent
1761 // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1762 // (all are UINT64)
1763 UINT64 tmp64;
1764 BID_UI64DOUBLE tmp1;
1765 unsigned int x_nr_bits;
1766 int q, ind, shift;
1767 UINT128 C1;
1768 // UINT128 res is C* at first - represents up to 34 decimal digits ~
1769 // 113 bits
1770 // UINT256 fstar;
1771 UINT256 P256;
1772
1773 // check for NaN or Infinity
1774 if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1775 // x is special
1776 if ((x.w[1] & MASK_NAN) == MASK_NAN) { // x is NAN
1777 // if x = NaN, then res = Q (x)
1778 // check first for non-canonical NaN payload
1779 if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1780 (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1781 (x.w[0] > 0x38c15b09ffffffffull))) {
1782 x.w[1] = x.w[1] & 0xffffc00000000000ull;
1783 x.w[0] = 0x0ull;
1784 }
1785 if ((x.w[1] & MASK_SNAN) == MASK_SNAN) { // x is SNAN
1786 // set invalid flag
1787 *pfpsf |= INVALID_EXCEPTION;
1788 // return quiet (x)
1789 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out also G[6]-G[16]
1790 res.w[0] = x.w[0];
1791 } else { // x is QNaN
1792 // return x
1793 res.w[1] = x.w[1] & 0xfc003fffffffffffull; // clear out G[6]-G[16]
1794 res.w[0] = x.w[0];
1795 }
1796 BID_RETURN (res)
1797 } else { // x is not a NaN, so it must be infinity
1798 if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1799 // return +inf
1800 res.w[1] = 0x7800000000000000ull;
1801 res.w[0] = 0x0000000000000000ull;
1802 } else { // x is -inf
1803 // return -inf
1804 res.w[1] = 0xf800000000000000ull;
1805 res.w[0] = 0x0000000000000000ull;
1806 }
1807 BID_RETURN (res);
1808 }
1809 }
1810 // unpack x
1811 x_sign = x.w[1] & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1812 C1.w[1] = x.w[1] & MASK_COEFF;
1813 C1.w[0] = x.w[0];
1814
1815 // check for non-canonical values (treated as zero)
1816 if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) { // G0_G1=11
1817 // non-canonical
1818 x_exp = (x.w[1] << 2) & MASK_EXP; // biased and shifted left 49 bits
1819 C1.w[1] = 0; // significand high
1820 C1.w[0] = 0; // significand low
1821 } else { // G0_G1 != 11
1822 x_exp = x.w[1] & MASK_EXP; // biased and shifted left 49 bits
1823 if (C1.w[1] > 0x0001ed09bead87c0ull ||
1824 (C1.w[1] == 0x0001ed09bead87c0ull
1825 && C1.w[0] > 0x378d8e63ffffffffull)) {
1826 // x is non-canonical if coefficient is larger than 10^34 -1
1827 C1.w[1] = 0;
1828 C1.w[0] = 0;
1829 } else { // canonical
1830 ;
1831 }
1832 }
1833
1834 // test for input equal to zero
1835 if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1836 // x is 0
1837 // return 0 preserving the sign bit and the preferred exponent
1838 // of MAX(Q(x), 0)
1839 if (x_exp <= (0x1820ull << 49)) {
1840 res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1841 } else {
1842 res.w[1] = x_sign | x_exp;
1843 }
1844 res.w[0] = 0x0000000000000000ull;
1845 BID_RETURN (res);
1846 }
1847 // x is not special and is not zero
1848
1849 // if (exp <= -(p+1)) return 0.0
1850 if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35
1851 res.w[1] = x_sign | 0x3040000000000000ull;
1852 res.w[0] = 0x0000000000000000ull;
1853 BID_RETURN (res);
1854 }
1855 // q = nr. of decimal digits in x
1856 // determine first the nr. of bits in x
1857 if (C1.w[1] == 0) {
1858 if (C1.w[0] >= 0x0020000000000000ull) { // x >= 2^53
1859 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1860 if (C1.w[0] >= 0x0000000100000000ull) { // x >= 2^32
1861 tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1862 x_nr_bits =
1863 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1864 } else { // x < 2^32
1865 tmp1.d = (double) (C1.w[0]); // exact conversion
1866 x_nr_bits =
1867 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1868 }
1869 } else { // if x < 2^53
1870 tmp1.d = (double) C1.w[0]; // exact conversion
1871 x_nr_bits =
1872 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1873 }
1874 } else { // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1875 tmp1.d = (double) C1.w[1]; // exact conversion
1876 x_nr_bits =
1877 65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1878 }
1879
1880 q = nr_digits[x_nr_bits - 1].digits;
1881 if (q == 0) {
1882 q = nr_digits[x_nr_bits - 1].digits1;
1883 if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1884 (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1885 C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1886 q++;
1887 }
1888 exp = (x_exp >> 49) - 6176;
1889 if (exp >= 0) { // -exp <= 0
1890 // the argument is an integer already
1891 res.w[1] = x.w[1];
1892 res.w[0] = x.w[0];
1893 BID_RETURN (res);
1894 } else if ((q + exp) >= 0) { // exp < 0 and 1 <= -exp <= q
1895 // need to shift right -exp digits from the coefficient; the exp will be 0
1896 ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
1897 // chop off ind digits from the lower part of C1
1898 // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
1899 tmp64 = C1.w[0];
1900 if (ind <= 19) {
1901 C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1902 } else {
1903 C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1904 C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1905 }
1906 if (C1.w[0] < tmp64)
1907 C1.w[1]++;
1908 // calculate C* and f*
1909 // C* is actually floor(C*) in this case
1910 // C* and f* need shifting and masking, as shown by
1911 // shiftright128[] and maskhigh128[]
1912 // 1 <= x <= 34
1913 // kx = 10^(-x) = ten2mk128[ind - 1]
1914 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1915 // the approximation of 10^(-x) was rounded up to 118 bits
1916 __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1917 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1918 // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1919 // if (0 < f* < 10^(-x)) then the result is a midpoint
1920 // if floor(C*) is even then C* = floor(C*) - logical right
1921 // shift; C* has p decimal digits, correct by Prop. 1)
1922 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
1923 // shift; C* has p decimal digits, correct by Pr. 1)
1924 // else
1925 // C* = floor(C*) (logical right shift; C has p decimal digits,
1926 // correct by Property 1)
1927 // n = C* * 10^(e+x)
1928
1929 // shift right C* by Ex-128 = shiftright128[ind]
1930 if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
1931 res.w[1] = P256.w[3];
1932 res.w[0] = P256.w[2];
1933 } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1934 shift = shiftright128[ind - 1]; // 3 <= shift <= 63
1935 res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1936 res.w[1] = (P256.w[3] >> shift);
1937 } else { // 22 <= ind - 1 <= 33
1938 shift = shiftright128[ind - 1]; // 2 <= shift <= 38
1939 res.w[1] = 0;
1940 res.w[0] = (P256.w[3] >> (shift - 64)); // 2 <= shift - 64 <= 38
1941 }
1942 // if the result was a midpoint, it was already rounded away from zero
1943 res.w[1] |= x_sign | 0x3040000000000000ull;
1944 BID_RETURN (res);
1945 } else { // if ((q + exp) < 0) <=> q < -exp
1946 // the result is +0 or -0
1947 res.w[1] = x_sign | 0x3040000000000000ull;
1948 res.w[0] = 0x0000000000000000ull;
1949 BID_RETURN (res);
1950 }
1951 }