?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Forgot ChangeLog of part1 pr25561.
[official-gcc.git] / libgcc / config / libbid / bid64_to_int32.c
blobb4e62371b6ae622ef4b875d124af5d7036083737
1 /* Copyright (C) 2007 Free Software Foundation, Inc.
2
3 This file is part of GCC.
4
5 GCC is free software; you can redistribute it and/or modify it under
6 the terms of the GNU General Public License as published by the Free
7 Software Foundation; either version 2, or (at your option) any later
8 version.
9
10 In addition to the permissions in the GNU General Public License, the
11 Free Software Foundation gives you unlimited permission to link the
12 compiled version of this file into combinations with other programs,
13 and to distribute those combinations without any restriction coming
14 from the use of this file. (The General Public License restrictions
15 do apply in other respects; for example, they cover modification of
16 the file, and distribution when not linked into a combine
17 executable.)
18
19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with GCC; see the file COPYING. If not, write to the Free
26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
27 02110-1301, USA. */
28
29 #include "bid_internal.h"
30
31 /*****************************************************************************
32 * BID64_to_int32_rnint
33 ****************************************************************************/
34
35 #if DECIMAL_CALL_BY_REFERENCE
36 void
37 bid64_to_int32_rnint (int *pres,
38 UINT64 *
39 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
40 _EXC_INFO_PARAM) {
41 UINT64 x = *px;
42 #else
43 int
44 bid64_to_int32_rnint (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
45 _EXC_INFO_PARAM) {
46 #endif
47 int res;
48 UINT64 x_sign;
49 UINT64 x_exp;
50 int exp; // unbiased exponent
51 // Note: C1 represents x_significand (UINT64)
52 UINT64 tmp64;
53 BID_UI64DOUBLE tmp1;
54 unsigned int x_nr_bits;
55 int q, ind, shift;
56 UINT64 C1;
57 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
58 UINT128 fstar;
59 UINT128 P128;
60
61 // check for NaN or Infinity
62 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
63 // set invalid flag
64 *pfpsf |= INVALID_EXCEPTION;
65 // return Integer Indefinite
66 res = 0x80000000;
67 BID_RETURN (res);
68 }
69 // unpack x
70 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
71 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
72 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
73 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
74 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
75 if (C1 > 9999999999999999ull) { // non-canonical
76 x_exp = 0;
77 C1 = 0;
78 }
79 } else {
80 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
81 C1 = x & MASK_BINARY_SIG1;
82 }
83
84 // check for zeros (possibly from non-canonical values)
85 if (C1 == 0x0ull) {
86 // x is 0
87 res = 0x00000000;
88 BID_RETURN (res);
89 }
90 // x is not special and is not zero
91
92 // q = nr. of decimal digits in x (1 <= q <= 54)
93 // determine first the nr. of bits in x
94 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
95 // split the 64-bit value in two 32-bit halves to avoid rounding errors
96 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
97 tmp1.d = (double) (C1 >> 32); // exact conversion
98 x_nr_bits =
99 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
100 } else { // x < 2^32
101 tmp1.d = (double) C1; // exact conversion
102 x_nr_bits =
103 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
104 }
105 } else { // if x < 2^53
106 tmp1.d = (double) C1; // exact conversion
107 x_nr_bits =
108 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
109 }
110 q = nr_digits[x_nr_bits - 1].digits;
111 if (q == 0) {
112 q = nr_digits[x_nr_bits - 1].digits1;
113 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
114 q++;
115 }
116 exp = x_exp - 398; // unbiased exponent
117
118 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
119 // set invalid flag
120 *pfpsf |= INVALID_EXCEPTION;
121 // return Integer Indefinite
122 res = 0x80000000;
123 BID_RETURN (res);
124 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
125 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
126 // so x rounded to an integer may or may not fit in a signed 32-bit int
127 // the cases that do not fit are identified here; the ones that fit
128 // fall through and will be handled with other cases further,
129 // under '1 <= q + exp <= 10'
130 if (x_sign) { // if n < 0 and q + exp = 10
131 // if n < -2^31 - 1/2 then n is too large
132 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31+1/2
133 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005, 1<=q<=16
134 // <=> C * 10^(11-q) > 0x500000005, 1<=q<=16
135 if (q <= 11) {
136 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
137 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
138 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
139 if (tmp64 > 0x500000005ull) {
140 // set invalid flag
141 *pfpsf |= INVALID_EXCEPTION;
142 // return Integer Indefinite
143 res = 0x80000000;
144 BID_RETURN (res);
145 }
146 // else cases that can be rounded to a 32-bit int fall through
147 // to '1 <= q + exp <= 10'
148 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
149 // C * 10^(11-q) > 0x500000005 <=>
150 // C > 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 5
151 // (scale 2^31+1/2 up)
152 // Note: 0x500000005*10^(q-11) has q-1 or q digits, where q <= 16
153 tmp64 = 0x500000005ull * ten2k64[q - 11];
154 if (C1 > tmp64) {
155 // set invalid flag
156 *pfpsf |= INVALID_EXCEPTION;
157 // return Integer Indefinite
158 res = 0x80000000;
159 BID_RETURN (res);
160 }
161 // else cases that can be rounded to a 32-bit int fall through
162 // to '1 <= q + exp <= 10'
163 }
164 } else { // if n > 0 and q + exp = 10
165 // if n >= 2^31 - 1/2 then n is too large
166 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
167 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=16
168 // <=> C * 10^(11-q) >= 0x4fffffffb, 1<=q<=16
169 if (q <= 11) {
170 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
171 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
172 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
173 if (tmp64 >= 0x4fffffffbull) {
174 // set invalid flag
175 *pfpsf |= INVALID_EXCEPTION;
176 // return Integer Indefinite
177 res = 0x80000000;
178 BID_RETURN (res);
179 }
180 // else cases that can be rounded to a 32-bit int fall through
181 // to '1 <= q + exp <= 10'
182 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
183 // C * 10^(11-q) >= 0x4fffffffb <=>
184 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
185 // (scale 2^31-1/2 up)
186 // Note: 0x4fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
187 tmp64 = 0x4fffffffbull * ten2k64[q - 11];
188 if (C1 >= tmp64) {
189 // set invalid flag
190 *pfpsf |= INVALID_EXCEPTION;
191 // return Integer Indefinite
192 res = 0x80000000;
193 BID_RETURN (res);
194 }
195 // else cases that can be rounded to a 32-bit int fall through
196 // to '1 <= q + exp <= 10'
197 }
198 }
199 }
200 // n is not too large to be converted to int32: -2^31 - 1/2 <= n < 2^31 - 1/2
201 // Note: some of the cases tested for above fall through to this point
202 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
203 // return 0
204 res = 0x00000000;
205 BID_RETURN (res);
206 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
207 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
208 // res = 0
209 // else
210 // res = +/-1
211 ind = q - 1;
212 if (C1 <= midpoint64[ind]) {
213 res = 0x00000000; // return 0
214 } else if (x_sign) { // n < 0
215 res = 0xffffffff; // return -1
216 } else { // n > 0
217 res = 0x00000001; // return +1
218 }
219 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
220 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
221 // to nearest to a 32-bit signed integer
222 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
223 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
224 // chop off ind digits from the lower part of C1
225 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
226 C1 = C1 + midpoint64[ind - 1];
227 // calculate C* and f*
228 // C* is actually floor(C*) in this case
229 // C* and f* need shifting and masking, as shown by
230 // shiftright128[] and maskhigh128[]
231 // 1 <= x <= 15
232 // kx = 10^(-x) = ten2mk64[ind - 1]
233 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
234 // the approximation of 10^(-x) was rounded up to 54 bits
235 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
236 Cstar = P128.w[1];
237 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
238 fstar.w[0] = P128.w[0];
239 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
240 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
241 // if (0 < f* < 10^(-x)) then the result is a midpoint
242 // if floor(C*) is even then C* = floor(C*) - logical right
243 // shift; C* has p decimal digits, correct by Prop. 1)
244 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
245 // shift; C* has p decimal digits, correct by Pr. 1)
246 // else
247 // C* = floor(C*) (logical right shift; C has p decimal digits,
248 // correct by Property 1)
249 // n = C* * 10^(e+x)
250
251 // shift right C* by Ex-64 = shiftright128[ind]
252 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
253 Cstar = Cstar >> shift;
254
255 // if the result was a midpoint it was rounded away from zero, so
256 // it will need a correction
257 // check for midpoints
258 if ((fstar.w[1] == 0) && fstar.w[0]
259 && (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) {
260 // ten2mk128trunc[ind -1].w[1] is identical to
261 // ten2mk128[ind -1].w[1]
262 // the result is a midpoint; round to nearest
263 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
264 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
265 Cstar--; // Cstar is now even
266 } // else MP in [ODD, EVEN]
267 }
268 if (x_sign)
269 res = -Cstar;
270 else
271 res = Cstar;
272 } else if (exp == 0) {
273 // 1 <= q <= 10
274 // res = +/-C (exact)
275 if (x_sign)
276 res = -C1;
277 else
278 res = C1;
279 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
280 // res = +/-C * 10^exp (exact)
281 if (x_sign)
282 res = -C1 * ten2k64[exp];
283 else
284 res = C1 * ten2k64[exp];
285 }
286 }
287 BID_RETURN (res);
288 }
289
290 /*****************************************************************************
291 * BID64_to_int32_xrnint
292 ****************************************************************************/
293
294 #if DECIMAL_CALL_BY_REFERENCE
295 void
296 bid64_to_int32_xrnint (int *pres,
297 UINT64 *
298 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
299 _EXC_INFO_PARAM) {
300 UINT64 x = *px;
301 #else
302 int
303 bid64_to_int32_xrnint (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
304 _EXC_INFO_PARAM) {
305 #endif
306 int res;
307 UINT64 x_sign;
308 UINT64 x_exp;
309 int exp; // unbiased exponent
310 // Note: C1 represents x_significand (UINT64)
311 UINT64 tmp64;
312 BID_UI64DOUBLE tmp1;
313 unsigned int x_nr_bits;
314 int q, ind, shift;
315 UINT64 C1;
316 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
317 UINT128 fstar;
318 UINT128 P128;
319
320 // check for NaN or Infinity
321 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
322 // set invalid flag
323 *pfpsf |= INVALID_EXCEPTION;
324 // return Integer Indefinite
325 res = 0x80000000;
326 BID_RETURN (res);
327 }
328 // unpack x
329 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
330 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
331 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
332 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
333 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
334 if (C1 > 9999999999999999ull) { // non-canonical
335 x_exp = 0;
336 C1 = 0;
337 }
338 } else {
339 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
340 C1 = x & MASK_BINARY_SIG1;
341 }
342
343 // check for zeros (possibly from non-canonical values)
344 if (C1 == 0x0ull) {
345 // x is 0
346 res = 0x00000000;
347 BID_RETURN (res);
348 }
349 // x is not special and is not zero
350
351 // q = nr. of decimal digits in x (1 <= q <= 54)
352 // determine first the nr. of bits in x
353 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
354 // split the 64-bit value in two 32-bit halves to avoid rounding errors
355 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
356 tmp1.d = (double) (C1 >> 32); // exact conversion
357 x_nr_bits =
358 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
359 } else { // x < 2^32
360 tmp1.d = (double) C1; // exact conversion
361 x_nr_bits =
362 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
363 }
364 } else { // if x < 2^53
365 tmp1.d = (double) C1; // exact conversion
366 x_nr_bits =
367 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
368 }
369 q = nr_digits[x_nr_bits - 1].digits;
370 if (q == 0) {
371 q = nr_digits[x_nr_bits - 1].digits1;
372 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
373 q++;
374 }
375 exp = x_exp - 398; // unbiased exponent
376
377 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
378 // set invalid flag
379 *pfpsf |= INVALID_EXCEPTION;
380 // return Integer Indefinite
381 res = 0x80000000;
382 BID_RETURN (res);
383 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
384 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
385 // so x rounded to an integer may or may not fit in a signed 32-bit int
386 // the cases that do not fit are identified here; the ones that fit
387 // fall through and will be handled with other cases further,
388 // under '1 <= q + exp <= 10'
389 if (x_sign) { // if n < 0 and q + exp = 10
390 // if n < -2^31 - 1/2 then n is too large
391 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31+1/2
392 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000005, 1<=q<=16
393 // <=> C * 10^(11-q) > 0x500000005, 1<=q<=16
394 if (q <= 11) {
395 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
396 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
397 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
398 if (tmp64 > 0x500000005ull) {
399 // set invalid flag
400 *pfpsf |= INVALID_EXCEPTION;
401 // return Integer Indefinite
402 res = 0x80000000;
403 BID_RETURN (res);
404 }
405 // else cases that can be rounded to a 32-bit int fall through
406 // to '1 <= q + exp <= 10'
407 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
408 // C * 10^(11-q) > 0x500000005 <=>
409 // C > 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 5
410 // (scale 2^31+1/2 up)
411 // Note: 0x500000005*10^(q-11) has q-1 or q digits, where q <= 16
412 tmp64 = 0x500000005ull * ten2k64[q - 11];
413 if (C1 > tmp64) {
414 // set invalid flag
415 *pfpsf |= INVALID_EXCEPTION;
416 // return Integer Indefinite
417 res = 0x80000000;
418 BID_RETURN (res);
419 }
420 // else cases that can be rounded to a 32-bit int fall through
421 // to '1 <= q + exp <= 10'
422 }
423 } else { // if n > 0 and q + exp = 10
424 // if n >= 2^31 - 1/2 then n is too large
425 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
426 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=16
427 // <=> C * 10^(11-q) >= 0x4fffffffb, 1<=q<=16
428 if (q <= 11) {
429 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
430 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
431 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
432 if (tmp64 >= 0x4fffffffbull) {
433 // set invalid flag
434 *pfpsf |= INVALID_EXCEPTION;
435 // return Integer Indefinite
436 res = 0x80000000;
437 BID_RETURN (res);
438 }
439 // else cases that can be rounded to a 32-bit int fall through
440 // to '1 <= q + exp <= 10'
441 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
442 // C * 10^(11-q) >= 0x4fffffffb <=>
443 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
444 // (scale 2^31-1/2 up)
445 // Note: 0x4fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
446 tmp64 = 0x4fffffffbull * ten2k64[q - 11];
447 if (C1 >= tmp64) {
448 // set invalid flag
449 *pfpsf |= INVALID_EXCEPTION;
450 // return Integer Indefinite
451 res = 0x80000000;
452 BID_RETURN (res);
453 }
454 // else cases that can be rounded to a 32-bit int fall through
455 // to '1 <= q + exp <= 10'
456 }
457 }
458 }
459 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
460 // Note: some of the cases tested for above fall through to this point
461 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
462 // set inexact flag
463 *pfpsf |= INEXACT_EXCEPTION;
464 // return 0
465 res = 0x00000000;
466 BID_RETURN (res);
467 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
468 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
469 // res = 0
470 // else
471 // res = +/-1
472 ind = q - 1;
473 if (C1 <= midpoint64[ind]) {
474 res = 0x00000000; // return 0
475 } else if (x_sign) { // n < 0
476 res = 0xffffffff; // return -1
477 } else { // n > 0
478 res = 0x00000001; // return +1
479 }
480 // set inexact flag
481 *pfpsf |= INEXACT_EXCEPTION;
482 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
483 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
484 // to nearest to a 32-bit signed integer
485 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
486 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
487 // chop off ind digits from the lower part of C1
488 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
489 C1 = C1 + midpoint64[ind - 1];
490 // calculate C* and f*
491 // C* is actually floor(C*) in this case
492 // C* and f* need shifting and masking, as shown by
493 // shiftright128[] and maskhigh128[]
494 // 1 <= x <= 15
495 // kx = 10^(-x) = ten2mk64[ind - 1]
496 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
497 // the approximation of 10^(-x) was rounded up to 54 bits
498 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
499 Cstar = P128.w[1];
500 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
501 fstar.w[0] = P128.w[0];
502 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
503 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
504 // if (0 < f* < 10^(-x)) then the result is a midpoint
505 // if floor(C*) is even then C* = floor(C*) - logical right
506 // shift; C* has p decimal digits, correct by Prop. 1)
507 // else if floor(C*) is odd C* = floor(C*)-1 (logical right
508 // shift; C* has p decimal digits, correct by Pr. 1)
509 // else
510 // C* = floor(C*) (logical right shift; C has p decimal digits,
511 // correct by Property 1)
512 // n = C* * 10^(e+x)
513
514 // shift right C* by Ex-64 = shiftright128[ind]
515 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
516 Cstar = Cstar >> shift;
517 // determine inexactness of the rounding of C*
518 // if (0 < f* - 1/2 < 10^(-x)) then
519 // the result is exact
520 // else // if (f* - 1/2 > T*) then
521 // the result is inexact
522 if (ind - 1 <= 2) {
523 if (fstar.w[0] > 0x8000000000000000ull) {
524 // f* > 1/2 and the result may be exact
525 tmp64 = fstar.w[0] - 0x8000000000000000ull; // f* - 1/2
526 if ((tmp64 > ten2mk128trunc[ind - 1].w[1])) {
527 // ten2mk128trunc[ind -1].w[1] is identical to
528 // ten2mk128[ind -1].w[1]
529 // set the inexact flag
530 *pfpsf |= INEXACT_EXCEPTION;
531 } // else the result is exact
532 } else { // the result is inexact; f2* <= 1/2
533 // set the inexact flag
534 *pfpsf |= INEXACT_EXCEPTION;
535 }
536 } else { // if 3 <= ind - 1 <= 14
537 if (fstar.w[1] > onehalf128[ind - 1] ||
538 (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) {
539 // f2* > 1/2 and the result may be exact
540 // Calculate f2* - 1/2
541 tmp64 = fstar.w[1] - onehalf128[ind - 1];
542 if (tmp64 || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
543 // ten2mk128trunc[ind -1].w[1] is identical to
544 // ten2mk128[ind -1].w[1]
545 // set the inexact flag
546 *pfpsf |= INEXACT_EXCEPTION;
547 } // else the result is exact
548 } else { // the result is inexact; f2* <= 1/2
549 // set the inexact flag
550 *pfpsf |= INEXACT_EXCEPTION;
551 }
552 }
553
554 // if the result was a midpoint it was rounded away from zero, so
555 // it will need a correction
556 // check for midpoints
557 if ((fstar.w[1] == 0) && fstar.w[0]
558 && (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) {
559 // ten2mk128trunc[ind -1].w[1] is identical to
560 // ten2mk128[ind -1].w[1]
561 // the result is a midpoint; round to nearest
562 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
563 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
564 Cstar--; // Cstar is now even
565 } // else MP in [ODD, EVEN]
566 }
567 if (x_sign)
568 res = -Cstar;
569 else
570 res = Cstar;
571 } else if (exp == 0) {
572 // 1 <= q <= 10
573 // res = +/-C (exact)
574 if (x_sign)
575 res = -C1;
576 else
577 res = C1;
578 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
579 // res = +/-C * 10^exp (exact)
580 if (x_sign)
581 res = -C1 * ten2k64[exp];
582 else
583 res = C1 * ten2k64[exp];
584 }
585 }
586 BID_RETURN (res);
587 }
588
589 /*****************************************************************************
590 * BID64_to_int32_floor
591 ****************************************************************************/
592
593 #if DECIMAL_CALL_BY_REFERENCE
594 void
595 bid64_to_int32_floor (int *pres,
596 UINT64 *
597 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
598 _EXC_INFO_PARAM) {
599 UINT64 x = *px;
600 #else
601 int
602 bid64_to_int32_floor (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
603 _EXC_INFO_PARAM) {
604 #endif
605 int res;
606 UINT64 x_sign;
607 UINT64 x_exp;
608 int exp; // unbiased exponent
609 // Note: C1 represents x_significand (UINT64)
610 UINT64 tmp64;
611 BID_UI64DOUBLE tmp1;
612 unsigned int x_nr_bits;
613 int q, ind, shift;
614 UINT64 C1;
615 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
616 UINT128 fstar;
617 UINT128 P128;
618
619 // check for NaN or Infinity
620 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
621 // set invalid flag
622 *pfpsf |= INVALID_EXCEPTION;
623 // return Integer Indefinite
624 res = 0x80000000;
625 BID_RETURN (res);
626 }
627 // unpack x
628 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
629 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
630 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
631 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
632 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
633 if (C1 > 9999999999999999ull) { // non-canonical
634 x_exp = 0;
635 C1 = 0;
636 }
637 } else {
638 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
639 C1 = x & MASK_BINARY_SIG1;
640 }
641
642 // check for zeros (possibly from non-canonical values)
643 if (C1 == 0x0ull) {
644 // x is 0
645 res = 0x00000000;
646 BID_RETURN (res);
647 }
648 // x is not special and is not zero
649
650 // q = nr. of decimal digits in x (1 <= q <= 54)
651 // determine first the nr. of bits in x
652 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
653 // split the 64-bit value in two 32-bit halves to avoid rounding errors
654 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
655 tmp1.d = (double) (C1 >> 32); // exact conversion
656 x_nr_bits =
657 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
658 } else { // x < 2^32
659 tmp1.d = (double) C1; // exact conversion
660 x_nr_bits =
661 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
662 }
663 } else { // if x < 2^53
664 tmp1.d = (double) C1; // exact conversion
665 x_nr_bits =
666 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
667 }
668 q = nr_digits[x_nr_bits - 1].digits;
669 if (q == 0) {
670 q = nr_digits[x_nr_bits - 1].digits1;
671 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
672 q++;
673 }
674 exp = x_exp - 398; // unbiased exponent
675
676 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
677 // set invalid flag
678 *pfpsf |= INVALID_EXCEPTION;
679 // return Integer Indefinite
680 res = 0x80000000;
681 BID_RETURN (res);
682 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
683 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
684 // so x rounded to an integer may or may not fit in a signed 32-bit int
685 // the cases that do not fit are identified here; the ones that fit
686 // fall through and will be handled with other cases further,
687 // under '1 <= q + exp <= 10'
688 if (x_sign) { // if n < 0 and q + exp = 10
689 // if n < -2^31 then n is too large
690 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31
691 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000, 1<=q<=16
692 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
693 if (q <= 11) {
694 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
695 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
696 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
697 if (tmp64 > 0x500000000ull) {
698 // set invalid flag
699 *pfpsf |= INVALID_EXCEPTION;
700 // return Integer Indefinite
701 res = 0x80000000;
702 BID_RETURN (res);
703 }
704 // else cases that can be rounded to a 32-bit int fall through
705 // to '1 <= q + exp <= 10'
706 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
707 // C * 10^(11-q) > 0x500000000 <=>
708 // C > 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
709 // (scale 2^31+1 up)
710 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
711 tmp64 = 0x500000000ull * ten2k64[q - 11];
712 if (C1 > tmp64) {
713 // set invalid flag
714 *pfpsf |= INVALID_EXCEPTION;
715 // return Integer Indefinite
716 res = 0x80000000;
717 BID_RETURN (res);
718 }
719 // else cases that can be rounded to a 32-bit int fall through
720 // to '1 <= q + exp <= 10'
721 }
722 } else { // if n > 0 and q + exp = 10
723 // if n >= 2^31 then n is too large
724 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
725 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=16
726 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
727 if (q <= 11) {
728 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
729 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
730 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
731 if (tmp64 >= 0x500000000ull) {
732 // set invalid flag
733 *pfpsf |= INVALID_EXCEPTION;
734 // return Integer Indefinite
735 res = 0x80000000;
736 BID_RETURN (res);
737 }
738 // else cases that can be rounded to a 32-bit int fall through
739 // to '1 <= q + exp <= 10'
740 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
741 // C * 10^(11-q) >= 0x500000000 <=>
742 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
743 // (scale 2^31-1 up)
744 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
745 tmp64 = 0x500000000ull * ten2k64[q - 11];
746 if (C1 >= tmp64) {
747 // set invalid flag
748 *pfpsf |= INVALID_EXCEPTION;
749 // return Integer Indefinite
750 res = 0x80000000;
751 BID_RETURN (res);
752 }
753 // else cases that can be rounded to a 32-bit int fall through
754 // to '1 <= q + exp <= 10'
755 }
756 }
757 }
758 // n is not too large to be converted to int32: -2^31 <= n < 2^31
759 // Note: some of the cases tested for above fall through to this point
760 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
761 // return -1 or 0
762 if (x_sign)
763 res = 0xffffffff;
764 else
765 res = 0x00000000;
766 BID_RETURN (res);
767 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
768 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
769 // to nearest to a 32-bit signed integer
770 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
771 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
772 // chop off ind digits from the lower part of C1
773 // C1 fits in 64 bits
774 // calculate C* and f*
775 // C* is actually floor(C*) in this case
776 // C* and f* need shifting and masking, as shown by
777 // shiftright128[] and maskhigh128[]
778 // 1 <= x <= 15
779 // kx = 10^(-x) = ten2mk64[ind - 1]
780 // C* = C1 * 10^(-x)
781 // the approximation of 10^(-x) was rounded up to 54 bits
782 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
783 Cstar = P128.w[1];
784 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
785 fstar.w[0] = P128.w[0];
786 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
787 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
788 // C* = floor(C*) (logical right shift; C has p decimal digits,
789 // correct by Property 1)
790 // n = C* * 10^(e+x)
791
792 // shift right C* by Ex-64 = shiftright128[ind]
793 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
794 Cstar = Cstar >> shift;
795 // determine inexactness of the rounding of C*
796 // if (0 < f* < 10^(-x)) then
797 // the result is exact
798 // else // if (f* > T*) then
799 // the result is inexact
800 if (ind - 1 <= 2) {
801 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
802 // ten2mk128trunc[ind -1].w[1] is identical to
803 // ten2mk128[ind -1].w[1]
804 if (x_sign) { // negative and inexact
805 Cstar++;
806 }
807 } // else the result is exact
808 } else { // if 3 <= ind - 1 <= 14
809 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
810 // ten2mk128trunc[ind -1].w[1] is identical to
811 // ten2mk128[ind -1].w[1]
812 if (x_sign) { // negative and inexact
813 Cstar++;
814 }
815 } // else the result is exact
816 }
817
818 if (x_sign)
819 res = -Cstar;
820 else
821 res = Cstar;
822 } else if (exp == 0) {
823 // 1 <= q <= 10
824 // res = +/-C (exact)
825 if (x_sign)
826 res = -C1;
827 else
828 res = C1;
829 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
830 // res = +/-C * 10^exp (exact)
831 if (x_sign)
832 res = -C1 * ten2k64[exp];
833 else
834 res = C1 * ten2k64[exp];
835 }
836 }
837 BID_RETURN (res);
838 }
839
840 /*****************************************************************************
841 * BID64_to_int32_xfloor
842 ****************************************************************************/
843
844 #if DECIMAL_CALL_BY_REFERENCE
845 void
846 bid64_to_int32_xfloor (int *pres,
847 UINT64 *
848 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
849 _EXC_INFO_PARAM) {
850 UINT64 x = *px;
851 #else
852 int
853 bid64_to_int32_xfloor (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
854 _EXC_INFO_PARAM) {
855 #endif
856 int res;
857 UINT64 x_sign;
858 UINT64 x_exp;
859 int exp; // unbiased exponent
860 // Note: C1 represents x_significand (UINT64)
861 UINT64 tmp64;
862 BID_UI64DOUBLE tmp1;
863 unsigned int x_nr_bits;
864 int q, ind, shift;
865 UINT64 C1;
866 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
867 UINT128 fstar;
868 UINT128 P128;
869
870 // check for NaN or Infinity
871 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
872 // set invalid flag
873 *pfpsf |= INVALID_EXCEPTION;
874 // return Integer Indefinite
875 res = 0x80000000;
876 BID_RETURN (res);
877 }
878 // unpack x
879 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
880 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
881 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
882 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
883 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
884 if (C1 > 9999999999999999ull) { // non-canonical
885 x_exp = 0;
886 C1 = 0;
887 }
888 } else {
889 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
890 C1 = x & MASK_BINARY_SIG1;
891 }
892
893 // check for zeros (possibly from non-canonical values)
894 if (C1 == 0x0ull) {
895 // x is 0
896 res = 0x00000000;
897 BID_RETURN (res);
898 }
899 // x is not special and is not zero
900
901 // q = nr. of decimal digits in x (1 <= q <= 54)
902 // determine first the nr. of bits in x
903 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
904 // split the 64-bit value in two 32-bit halves to avoid rounding errors
905 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
906 tmp1.d = (double) (C1 >> 32); // exact conversion
907 x_nr_bits =
908 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
909 } else { // x < 2^32
910 tmp1.d = (double) C1; // exact conversion
911 x_nr_bits =
912 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
913 }
914 } else { // if x < 2^53
915 tmp1.d = (double) C1; // exact conversion
916 x_nr_bits =
917 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
918 }
919 q = nr_digits[x_nr_bits - 1].digits;
920 if (q == 0) {
921 q = nr_digits[x_nr_bits - 1].digits1;
922 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
923 q++;
924 }
925 exp = x_exp - 398; // unbiased exponent
926
927 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
928 // set invalid flag
929 *pfpsf |= INVALID_EXCEPTION;
930 // return Integer Indefinite
931 res = 0x80000000;
932 BID_RETURN (res);
933 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
934 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
935 // so x rounded to an integer may or may not fit in a signed 32-bit int
936 // the cases that do not fit are identified here; the ones that fit
937 // fall through and will be handled with other cases further,
938 // under '1 <= q + exp <= 10'
939 if (x_sign) { // if n < 0 and q + exp = 10
940 // if n < -2^31 then n is too large
941 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31
942 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x500000000, 1<=q<=16
943 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
944 if (q <= 11) {
945 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
946 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
947 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
948 if (tmp64 > 0x500000000ull) {
949 // set invalid flag
950 *pfpsf |= INVALID_EXCEPTION;
951 // return Integer Indefinite
952 res = 0x80000000;
953 BID_RETURN (res);
954 }
955 // else cases that can be rounded to a 32-bit int fall through
956 // to '1 <= q + exp <= 10'
957 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
958 // C * 10^(11-q) > 0x500000000 <=>
959 // C > 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
960 // (scale 2^31+1 up)
961 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
962 tmp64 = 0x500000000ull * ten2k64[q - 11];
963 if (C1 > tmp64) {
964 // set invalid flag
965 *pfpsf |= INVALID_EXCEPTION;
966 // return Integer Indefinite
967 res = 0x80000000;
968 BID_RETURN (res);
969 }
970 // else cases that can be rounded to a 32-bit int fall through
971 // to '1 <= q + exp <= 10'
972 }
973 } else { // if n > 0 and q + exp = 10
974 // if n >= 2^31 then n is too large
975 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
976 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=16
977 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
978 if (q <= 11) {
979 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
980 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
981 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
982 if (tmp64 >= 0x500000000ull) {
983 // set invalid flag
984 *pfpsf |= INVALID_EXCEPTION;
985 // return Integer Indefinite
986 res = 0x80000000;
987 BID_RETURN (res);
988 }
989 // else cases that can be rounded to a 32-bit int fall through
990 // to '1 <= q + exp <= 10'
991 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
992 // C * 10^(11-q) >= 0x500000000 <=>
993 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
994 // (scale 2^31-1 up)
995 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
996 tmp64 = 0x500000000ull * ten2k64[q - 11];
997 if (C1 >= tmp64) {
998 // set invalid flag
999 *pfpsf |= INVALID_EXCEPTION;
1000 // return Integer Indefinite
1001 res = 0x80000000;
1002 BID_RETURN (res);
1003 }
1004 // else cases that can be rounded to a 32-bit int fall through
1005 // to '1 <= q + exp <= 10'
1006 }
1007 }
1008 }
1009 // n is not too large to be converted to int32: -2^31 <= n < 2^31
1010 // Note: some of the cases tested for above fall through to this point
1011 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1012 // set inexact flag
1013 *pfpsf |= INEXACT_EXCEPTION;
1014 // return -1 or 0
1015 if (x_sign)
1016 res = 0xffffffff;
1017 else
1018 res = 0x00000000;
1019 BID_RETURN (res);
1020 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1021 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
1022 // to nearest to a 32-bit signed integer
1023 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1024 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
1025 // chop off ind digits from the lower part of C1
1026 // C1 fits in 64 bits
1027 // calculate C* and f*
1028 // C* is actually floor(C*) in this case
1029 // C* and f* need shifting and masking, as shown by
1030 // shiftright128[] and maskhigh128[]
1031 // 1 <= x <= 15
1032 // kx = 10^(-x) = ten2mk64[ind - 1]
1033 // C* = C1 * 10^(-x)
1034 // the approximation of 10^(-x) was rounded up to 54 bits
1035 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
1036 Cstar = P128.w[1];
1037 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
1038 fstar.w[0] = P128.w[0];
1039 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1040 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1041 // C* = floor(C*) (logical right shift; C has p decimal digits,
1042 // correct by Property 1)
1043 // n = C* * 10^(e+x)
1044
1045 // shift right C* by Ex-64 = shiftright128[ind]
1046 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
1047 Cstar = Cstar >> shift;
1048 // determine inexactness of the rounding of C*
1049 // if (0 < f* < 10^(-x)) then
1050 // the result is exact
1051 // else // if (f* > T*) then
1052 // the result is inexact
1053 if (ind - 1 <= 2) {
1054 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1055 // ten2mk128trunc[ind -1].w[1] is identical to
1056 // ten2mk128[ind -1].w[1]
1057 if (x_sign) { // negative and inexact
1058 Cstar++;
1059 }
1060 // set the inexact flag
1061 *pfpsf |= INEXACT_EXCEPTION;
1062 } // else the result is exact
1063 } else { // if 3 <= ind - 1 <= 14
1064 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1065 // ten2mk128trunc[ind -1].w[1] is identical to
1066 // ten2mk128[ind -1].w[1]
1067 if (x_sign) { // negative and inexact
1068 Cstar++;
1069 }
1070 // set the inexact flag
1071 *pfpsf |= INEXACT_EXCEPTION;
1072 } // else the result is exact
1073 }
1074
1075 if (x_sign)
1076 res = -Cstar;
1077 else
1078 res = Cstar;
1079 } else if (exp == 0) {
1080 // 1 <= q <= 10
1081 // res = +/-C (exact)
1082 if (x_sign)
1083 res = -C1;
1084 else
1085 res = C1;
1086 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1087 // res = +/-C * 10^exp (exact)
1088 if (x_sign)
1089 res = -C1 * ten2k64[exp];
1090 else
1091 res = C1 * ten2k64[exp];
1092 }
1093 }
1094 BID_RETURN (res);
1095 }
1096
1097 /*****************************************************************************
1098 * BID64_to_int32_ceil
1099 ****************************************************************************/
1100
1101 #if DECIMAL_CALL_BY_REFERENCE
1102 void
1103 bid64_to_int32_ceil (int *pres,
1104 UINT64 *
1105 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1106 _EXC_INFO_PARAM) {
1107 UINT64 x = *px;
1108 #else
1109 int
1110 bid64_to_int32_ceil (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1111 _EXC_INFO_PARAM) {
1112 #endif
1113 int res;
1114 UINT64 x_sign;
1115 UINT64 x_exp;
1116 int exp; // unbiased exponent
1117 // Note: C1 represents x_significand (UINT64)
1118 UINT64 tmp64;
1119 BID_UI64DOUBLE tmp1;
1120 unsigned int x_nr_bits;
1121 int q, ind, shift;
1122 UINT64 C1;
1123 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
1124 UINT128 fstar;
1125 UINT128 P128;
1126
1127 // check for NaN or Infinity
1128 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
1129 // set invalid flag
1130 *pfpsf |= INVALID_EXCEPTION;
1131 // return Integer Indefinite
1132 res = 0x80000000;
1133 BID_RETURN (res);
1134 }
1135 // unpack x
1136 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1137 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1138 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1139 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
1140 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1141 if (C1 > 9999999999999999ull) { // non-canonical
1142 x_exp = 0;
1143 C1 = 0;
1144 }
1145 } else {
1146 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
1147 C1 = x & MASK_BINARY_SIG1;
1148 }
1149
1150 // check for zeros (possibly from non-canonical values)
1151 if (C1 == 0x0ull) {
1152 // x is 0
1153 res = 0x00000000;
1154 BID_RETURN (res);
1155 }
1156 // x is not special and is not zero
1157
1158 // q = nr. of decimal digits in x (1 <= q <= 54)
1159 // determine first the nr. of bits in x
1160 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1161 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1162 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1163 tmp1.d = (double) (C1 >> 32); // exact conversion
1164 x_nr_bits =
1165 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1166 } else { // x < 2^32
1167 tmp1.d = (double) C1; // exact conversion
1168 x_nr_bits =
1169 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1170 }
1171 } else { // if x < 2^53
1172 tmp1.d = (double) C1; // exact conversion
1173 x_nr_bits =
1174 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1175 }
1176 q = nr_digits[x_nr_bits - 1].digits;
1177 if (q == 0) {
1178 q = nr_digits[x_nr_bits - 1].digits1;
1179 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1180 q++;
1181 }
1182 exp = x_exp - 398; // unbiased exponent
1183
1184 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1185 // set invalid flag
1186 *pfpsf |= INVALID_EXCEPTION;
1187 // return Integer Indefinite
1188 res = 0x80000000;
1189 BID_RETURN (res);
1190 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1191 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1192 // so x rounded to an integer may or may not fit in a signed 32-bit int
1193 // the cases that do not fit are identified here; the ones that fit
1194 // fall through and will be handled with other cases further,
1195 // under '1 <= q + exp <= 10'
1196 if (x_sign) { // if n < 0 and q + exp = 10
1197 // if n <= -2^31 - 1 then n is too large
1198 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1199 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x50000000a, 1<=q<=16
1200 // <=> C * 10^(11-q) >= 0x50000000a, 1<=q<=16
1201 if (q <= 11) {
1202 // Note: C * 10^(11-q) has 10 or 11 digits; 0x50000000a has 11 digits
1203 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1204 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1205 if (tmp64 >= 0x50000000aull) {
1206 // set invalid flag
1207 *pfpsf |= INVALID_EXCEPTION;
1208 // return Integer Indefinite
1209 res = 0x80000000;
1210 BID_RETURN (res);
1211 }
1212 // else cases that can be rounded to a 32-bit int fall through
1213 // to '1 <= q + exp <= 10'
1214 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1215 // C * 10^(11-q) >= 0x50000000a <=>
1216 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 5
1217 // (scale 2^31+1 up)
1218 // Note: 0x50000000a*10^(q-11) has q-1 or q digits, where q <= 16
1219 tmp64 = 0x50000000aull * ten2k64[q - 11];
1220 if (C1 >= tmp64) {
1221 // set invalid flag
1222 *pfpsf |= INVALID_EXCEPTION;
1223 // return Integer Indefinite
1224 res = 0x80000000;
1225 BID_RETURN (res);
1226 }
1227 // else cases that can be rounded to a 32-bit int fall through
1228 // to '1 <= q + exp <= 10'
1229 }
1230 } else { // if n > 0 and q + exp = 10
1231 // if n > 2^31 - 1 then n is too large
1232 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 - 1
1233 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6, 1<=q<=16
1234 // <=> C * 10^(11-q) > 0x4fffffff6, 1<=q<=16
1235 if (q <= 11) {
1236 // Note: C * 10^(11-q) has 10 or 11 digits; 0x4fffffff6 has 11 digits
1237 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1238 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1239 if (tmp64 > 0x4fffffff6ull) {
1240 // set invalid flag
1241 *pfpsf |= INVALID_EXCEPTION;
1242 // return Integer Indefinite
1243 res = 0x80000000;
1244 BID_RETURN (res);
1245 }
1246 // else cases that can be rounded to a 32-bit int fall through
1247 // to '1 <= q + exp <= 10'
1248 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1249 // C * 10^(11-q) > 0x4fffffff6 <=>
1250 // C > 0x4fffffff6 * 10^(q-11) where 1 <= q - 11 <= 5
1251 // (scale 2^31-1 up)
1252 // Note: 0x4fffffff6*10^(q-11) has q-1 or q digits, where q <= 16
1253 tmp64 = 0x4fffffff6ull * ten2k64[q - 11];
1254 if (C1 > tmp64) {
1255 // set invalid flag
1256 *pfpsf |= INVALID_EXCEPTION;
1257 // return Integer Indefinite
1258 res = 0x80000000;
1259 BID_RETURN (res);
1260 }
1261 // else cases that can be rounded to a 32-bit int fall through
1262 // to '1 <= q + exp <= 10'
1263 }
1264 }
1265 }
1266 // n is not too large to be converted to int32: -2^31 - 1 < n <= 2^31 - 1
1267 // Note: some of the cases tested for above fall through to this point
1268 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1269 // return 0 or 1
1270 if (x_sign)
1271 res = 0x00000000;
1272 else
1273 res = 0x00000001;
1274 BID_RETURN (res);
1275 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1276 // -2^31-1 < x <= -1 or 1 <= x <= 2^31-1 so x can be rounded
1277 // to nearest to a 32-bit signed integer
1278 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1279 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
1280 // chop off ind digits from the lower part of C1
1281 // C1 fits in 64 bits
1282 // calculate C* and f*
1283 // C* is actually floor(C*) in this case
1284 // C* and f* need shifting and masking, as shown by
1285 // shiftright128[] and maskhigh128[]
1286 // 1 <= x <= 15
1287 // kx = 10^(-x) = ten2mk64[ind - 1]
1288 // C* = C1 * 10^(-x)
1289 // the approximation of 10^(-x) was rounded up to 54 bits
1290 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
1291 Cstar = P128.w[1];
1292 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
1293 fstar.w[0] = P128.w[0];
1294 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1295 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1296 // C* = floor(C*) (logical right shift; C has p decimal digits,
1297 // correct by Property 1)
1298 // n = C* * 10^(e+x)
1299
1300 // shift right C* by Ex-64 = shiftright128[ind]
1301 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
1302 Cstar = Cstar >> shift;
1303 // determine inexactness of the rounding of C*
1304 // if (0 < f* < 10^(-x)) then
1305 // the result is exact
1306 // else // if (f* > T*) then
1307 // the result is inexact
1308 if (ind - 1 <= 2) {
1309 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1310 // ten2mk128trunc[ind -1].w[1] is identical to
1311 // ten2mk128[ind -1].w[1]
1312 if (!x_sign) { // positive and inexact
1313 Cstar++;
1314 }
1315 } // else the result is exact
1316 } else { // if 3 <= ind - 1 <= 14
1317 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1318 // ten2mk128trunc[ind -1].w[1] is identical to
1319 // ten2mk128[ind -1].w[1]
1320 if (!x_sign) { // positive and inexact
1321 Cstar++;
1322 }
1323 } // else the result is exact
1324 }
1325
1326 if (x_sign)
1327 res = -Cstar;
1328 else
1329 res = Cstar;
1330 } else if (exp == 0) {
1331 // 1 <= q <= 10
1332 // res = +/-C (exact)
1333 if (x_sign)
1334 res = -C1;
1335 else
1336 res = C1;
1337 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1338 // res = +/-C * 10^exp (exact)
1339 if (x_sign)
1340 res = -C1 * ten2k64[exp];
1341 else
1342 res = C1 * ten2k64[exp];
1343 }
1344 }
1345 BID_RETURN (res);
1346 }
1347
1348 /*****************************************************************************
1349 * BID64_to_int32_xceil
1350 ****************************************************************************/
1351
1352 #if DECIMAL_CALL_BY_REFERENCE
1353 void
1354 bid64_to_int32_xceil (int *pres,
1355 UINT64 *
1356 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1357 _EXC_INFO_PARAM) {
1358 UINT64 x = *px;
1359 #else
1360 int
1361 bid64_to_int32_xceil (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1362 _EXC_INFO_PARAM) {
1363 #endif
1364 int res;
1365 UINT64 x_sign;
1366 UINT64 x_exp;
1367 int exp; // unbiased exponent
1368 // Note: C1 represents x_significand (UINT64)
1369 UINT64 tmp64;
1370 BID_UI64DOUBLE tmp1;
1371 unsigned int x_nr_bits;
1372 int q, ind, shift;
1373 UINT64 C1;
1374 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
1375 UINT128 fstar;
1376 UINT128 P128;
1377
1378 // check for NaN or Infinity
1379 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
1380 // set invalid flag
1381 *pfpsf |= INVALID_EXCEPTION;
1382 // return Integer Indefinite
1383 res = 0x80000000;
1384 BID_RETURN (res);
1385 }
1386 // unpack x
1387 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1388 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1389 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1390 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
1391 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1392 if (C1 > 9999999999999999ull) { // non-canonical
1393 x_exp = 0;
1394 C1 = 0;
1395 }
1396 } else {
1397 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
1398 C1 = x & MASK_BINARY_SIG1;
1399 }
1400
1401 // check for zeros (possibly from non-canonical values)
1402 if (C1 == 0x0ull) {
1403 // x is 0
1404 res = 0x00000000;
1405 BID_RETURN (res);
1406 }
1407 // x is not special and is not zero
1408
1409 // q = nr. of decimal digits in x (1 <= q <= 54)
1410 // determine first the nr. of bits in x
1411 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1412 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1413 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1414 tmp1.d = (double) (C1 >> 32); // exact conversion
1415 x_nr_bits =
1416 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1417 } else { // x < 2^32
1418 tmp1.d = (double) C1; // exact conversion
1419 x_nr_bits =
1420 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1421 }
1422 } else { // if x < 2^53
1423 tmp1.d = (double) C1; // exact conversion
1424 x_nr_bits =
1425 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1426 }
1427 q = nr_digits[x_nr_bits - 1].digits;
1428 if (q == 0) {
1429 q = nr_digits[x_nr_bits - 1].digits1;
1430 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1431 q++;
1432 }
1433 exp = x_exp - 398; // unbiased exponent
1434
1435 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1436 // set invalid flag
1437 *pfpsf |= INVALID_EXCEPTION;
1438 // return Integer Indefinite
1439 res = 0x80000000;
1440 BID_RETURN (res);
1441 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1442 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1443 // so x rounded to an integer may or may not fit in a signed 32-bit int
1444 // the cases that do not fit are identified here; the ones that fit
1445 // fall through and will be handled with other cases further,
1446 // under '1 <= q + exp <= 10'
1447 if (x_sign) { // if n < 0 and q + exp = 10
1448 // if n <= -2^31 - 1 then n is too large
1449 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1450 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x50000000a, 1<=q<=16
1451 // <=> C * 10^(11-q) >= 0x50000000a, 1<=q<=16
1452 if (q <= 11) {
1453 // Note: C * 10^(11-q) has 10 or 11 digits; 0x50000000a has 11 digits
1454 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1455 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1456 if (tmp64 >= 0x50000000aull) {
1457 // set invalid flag
1458 *pfpsf |= INVALID_EXCEPTION;
1459 // return Integer Indefinite
1460 res = 0x80000000;
1461 BID_RETURN (res);
1462 }
1463 // else cases that can be rounded to a 32-bit int fall through
1464 // to '1 <= q + exp <= 10'
1465 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1466 // C * 10^(11-q) >= 0x50000000a <=>
1467 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 5
1468 // (scale 2^31+1 up)
1469 // Note: 0x50000000a*10^(q-11) has q-1 or q digits, where q <= 16
1470 tmp64 = 0x50000000aull * ten2k64[q - 11];
1471 if (C1 >= tmp64) {
1472 // set invalid flag
1473 *pfpsf |= INVALID_EXCEPTION;
1474 // return Integer Indefinite
1475 res = 0x80000000;
1476 BID_RETURN (res);
1477 }
1478 // else cases that can be rounded to a 32-bit int fall through
1479 // to '1 <= q + exp <= 10'
1480 }
1481 } else { // if n > 0 and q + exp = 10
1482 // if n > 2^31 - 1 then n is too large
1483 // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^31 - 1
1484 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x4fffffff6, 1<=q<=16
1485 // <=> C * 10^(11-q) > 0x4fffffff6, 1<=q<=16
1486 if (q <= 11) {
1487 // Note: C * 10^(11-q) has 10 or 11 digits; 0x4fffffff6 has 11 digits
1488 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1489 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1490 if (tmp64 > 0x4fffffff6ull) {
1491 // set invalid flag
1492 *pfpsf |= INVALID_EXCEPTION;
1493 // return Integer Indefinite
1494 res = 0x80000000;
1495 BID_RETURN (res);
1496 }
1497 // else cases that can be rounded to a 32-bit int fall through
1498 // to '1 <= q + exp <= 10'
1499 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1500 // C * 10^(11-q) > 0x4fffffff6 <=>
1501 // C > 0x4fffffff6 * 10^(q-11) where 1 <= q - 11 <= 5
1502 // (scale 2^31-1 up)
1503 // Note: 0x4fffffff6*10^(q-11) has q-1 or q digits, where q <= 16
1504 tmp64 = 0x4fffffff6ull * ten2k64[q - 11];
1505 if (C1 > tmp64) {
1506 // set invalid flag
1507 *pfpsf |= INVALID_EXCEPTION;
1508 // return Integer Indefinite
1509 res = 0x80000000;
1510 BID_RETURN (res);
1511 }
1512 // else cases that can be rounded to a 32-bit int fall through
1513 // to '1 <= q + exp <= 10'
1514 }
1515 }
1516 }
1517 // n is not too large to be converted to int32: -2^31 - 1 < n <= 2^31 - 1
1518 // Note: some of the cases tested for above fall through to this point
1519 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1520 // set inexact flag
1521 *pfpsf |= INEXACT_EXCEPTION;
1522 // return 0 or 1
1523 if (x_sign)
1524 res = 0x00000000;
1525 else
1526 res = 0x00000001;
1527 BID_RETURN (res);
1528 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1529 // -2^31-1 < x <= -1 or 1 <= x <= 2^31-1 so x can be rounded
1530 // to nearest to a 32-bit signed integer
1531 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1532 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
1533 // chop off ind digits from the lower part of C1
1534 // C1 fits in 64 bits
1535 // calculate C* and f*
1536 // C* is actually floor(C*) in this case
1537 // C* and f* need shifting and masking, as shown by
1538 // shiftright128[] and maskhigh128[]
1539 // 1 <= x <= 15
1540 // kx = 10^(-x) = ten2mk64[ind - 1]
1541 // C* = C1 * 10^(-x)
1542 // the approximation of 10^(-x) was rounded up to 54 bits
1543 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
1544 Cstar = P128.w[1];
1545 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
1546 fstar.w[0] = P128.w[0];
1547 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1548 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1549 // C* = floor(C*) (logical right shift; C has p decimal digits,
1550 // correct by Property 1)
1551 // n = C* * 10^(e+x)
1552
1553 // shift right C* by Ex-64 = shiftright128[ind]
1554 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
1555 Cstar = Cstar >> shift;
1556 // determine inexactness of the rounding of C*
1557 // if (0 < f* < 10^(-x)) then
1558 // the result is exact
1559 // else // if (f* > T*) then
1560 // the result is inexact
1561 if (ind - 1 <= 2) {
1562 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1563 // ten2mk128trunc[ind -1].w[1] is identical to
1564 // ten2mk128[ind -1].w[1]
1565 if (!x_sign) { // positive and inexact
1566 Cstar++;
1567 }
1568 // set the inexact flag
1569 *pfpsf |= INEXACT_EXCEPTION;
1570 } // else the result is exact
1571 } else { // if 3 <= ind - 1 <= 14
1572 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
1573 // ten2mk128trunc[ind -1].w[1] is identical to
1574 // ten2mk128[ind -1].w[1]
1575 if (!x_sign) { // positive and inexact
1576 Cstar++;
1577 }
1578 // set the inexact flag
1579 *pfpsf |= INEXACT_EXCEPTION;
1580 } // else the result is exact
1581 }
1582
1583 if (x_sign)
1584 res = -Cstar;
1585 else
1586 res = Cstar;
1587 } else if (exp == 0) {
1588 // 1 <= q <= 10
1589 // res = +/-C (exact)
1590 if (x_sign)
1591 res = -C1;
1592 else
1593 res = C1;
1594 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1595 // res = +/-C * 10^exp (exact)
1596 if (x_sign)
1597 res = -C1 * ten2k64[exp];
1598 else
1599 res = C1 * ten2k64[exp];
1600 }
1601 }
1602 BID_RETURN (res);
1603 }
1604
1605 /*****************************************************************************
1606 * BID64_to_int32_int
1607 ****************************************************************************/
1608
1609 #if DECIMAL_CALL_BY_REFERENCE
1610 void
1611 bid64_to_int32_int (int *pres,
1612 UINT64 *
1613 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1614 _EXC_INFO_PARAM) {
1615 UINT64 x = *px;
1616 #else
1617 int
1618 bid64_to_int32_int (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1619 _EXC_INFO_PARAM) {
1620 #endif
1621 int res;
1622 UINT64 x_sign;
1623 UINT64 x_exp;
1624 int exp; // unbiased exponent
1625 // Note: C1 represents x_significand (UINT64)
1626 UINT64 tmp64;
1627 BID_UI64DOUBLE tmp1;
1628 unsigned int x_nr_bits;
1629 int q, ind, shift;
1630 UINT64 C1;
1631 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
1632 UINT128 P128;
1633
1634 // check for NaN or Infinity
1635 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
1636 // set invalid flag
1637 *pfpsf |= INVALID_EXCEPTION;
1638 // return Integer Indefinite
1639 res = 0x80000000;
1640 BID_RETURN (res);
1641 }
1642 // unpack x
1643 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1644 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1645 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1646 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
1647 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1648 if (C1 > 9999999999999999ull) { // non-canonical
1649 x_exp = 0;
1650 C1 = 0;
1651 }
1652 } else {
1653 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
1654 C1 = x & MASK_BINARY_SIG1;
1655 }
1656
1657 // check for zeros (possibly from non-canonical values)
1658 if (C1 == 0x0ull) {
1659 // x is 0
1660 res = 0x00000000;
1661 BID_RETURN (res);
1662 }
1663 // x is not special and is not zero
1664
1665 // q = nr. of decimal digits in x (1 <= q <= 54)
1666 // determine first the nr. of bits in x
1667 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1668 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1669 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1670 tmp1.d = (double) (C1 >> 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; // 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; // exact conversion
1680 x_nr_bits =
1681 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1682 }
1683 q = nr_digits[x_nr_bits - 1].digits;
1684 if (q == 0) {
1685 q = nr_digits[x_nr_bits - 1].digits1;
1686 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1687 q++;
1688 }
1689 exp = x_exp - 398; // unbiased exponent
1690
1691 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1692 // set invalid flag
1693 *pfpsf |= INVALID_EXCEPTION;
1694 // return Integer Indefinite
1695 res = 0x80000000;
1696 BID_RETURN (res);
1697 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1698 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1699 // so x rounded to an integer may or may not fit in a signed 32-bit int
1700 // the cases that do not fit are identified here; the ones that fit
1701 // fall through and will be handled with other cases further,
1702 // under '1 <= q + exp <= 10'
1703 if (x_sign) { // if n < 0 and q + exp = 10
1704 // if n <= -2^31 - 1 then n is too large
1705 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1706 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x50000000a, 1<=q<=16
1707 // <=> C * 10^(11-q) >= 0x50000000a, 1<=q<=16
1708 if (q <= 11) {
1709 // Note: C * 10^(11-q) has 10 or 11 digits; 0x50000000a has 11 digits
1710 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1711 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1712 if (tmp64 >= 0x50000000aull) {
1713 // set invalid flag
1714 *pfpsf |= INVALID_EXCEPTION;
1715 // return Integer Indefinite
1716 res = 0x80000000;
1717 BID_RETURN (res);
1718 }
1719 // else cases that can be rounded to a 32-bit int fall through
1720 // to '1 <= q + exp <= 10'
1721 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1722 // C * 10^(11-q) >= 0x50000000a <=>
1723 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 5
1724 // (scale 2^31+1 up)
1725 // Note: 0x50000000a*10^(q-11) has q-1 or q digits, where q <= 16
1726 tmp64 = 0x50000000aull * ten2k64[q - 11];
1727 if (C1 >= tmp64) {
1728 // set invalid flag
1729 *pfpsf |= INVALID_EXCEPTION;
1730 // return Integer Indefinite
1731 res = 0x80000000;
1732 BID_RETURN (res);
1733 }
1734 // else cases that can be rounded to a 32-bit int fall through
1735 // to '1 <= q + exp <= 10'
1736 }
1737 } else { // if n > 0 and q + exp = 10
1738 // if n >= 2^31 then n is too large
1739 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
1740 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=16
1741 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
1742 if (q <= 11) {
1743 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
1744 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1745 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1746 if (tmp64 >= 0x500000000ull) {
1747 // set invalid flag
1748 *pfpsf |= INVALID_EXCEPTION;
1749 // return Integer Indefinite
1750 res = 0x80000000;
1751 BID_RETURN (res);
1752 }
1753 // else cases that can be rounded to a 32-bit int fall through
1754 // to '1 <= q + exp <= 10'
1755 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1756 // C * 10^(11-q) >= 0x500000000 <=>
1757 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
1758 // (scale 2^31-1 up)
1759 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
1760 tmp64 = 0x500000000ull * ten2k64[q - 11];
1761 if (C1 >= tmp64) {
1762 // set invalid flag
1763 *pfpsf |= INVALID_EXCEPTION;
1764 // return Integer Indefinite
1765 res = 0x80000000;
1766 BID_RETURN (res);
1767 }
1768 // else cases that can be rounded to a 32-bit int fall through
1769 // to '1 <= q + exp <= 10'
1770 }
1771 }
1772 }
1773 // n is not too large to be converted to int32: -2^31 - 1 < n < 2^31
1774 // Note: some of the cases tested for above fall through to this point
1775 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1776 // return 0
1777 res = 0x00000000;
1778 BID_RETURN (res);
1779 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
1780 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
1781 // to nearest to a 32-bit signed integer
1782 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
1783 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
1784 // chop off ind digits from the lower part of C1
1785 // C1 fits in 64 bits
1786 // calculate C* and f*
1787 // C* is actually floor(C*) in this case
1788 // C* and f* need shifting and masking, as shown by
1789 // shiftright128[] and maskhigh128[]
1790 // 1 <= x <= 15
1791 // kx = 10^(-x) = ten2mk64[ind - 1]
1792 // C* = C1 * 10^(-x)
1793 // the approximation of 10^(-x) was rounded up to 54 bits
1794 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
1795 Cstar = P128.w[1];
1796 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
1797 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
1798 // C* = floor(C*) (logical right shift; C has p decimal digits,
1799 // correct by Property 1)
1800 // n = C* * 10^(e+x)
1801
1802 // shift right C* by Ex-64 = shiftright128[ind]
1803 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
1804 Cstar = Cstar >> shift;
1805 if (x_sign)
1806 res = -Cstar;
1807 else
1808 res = Cstar;
1809 } else if (exp == 0) {
1810 // 1 <= q <= 10
1811 // res = +/-C (exact)
1812 if (x_sign)
1813 res = -C1;
1814 else
1815 res = C1;
1816 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1817 // res = +/-C * 10^exp (exact)
1818 if (x_sign)
1819 res = -C1 * ten2k64[exp];
1820 else
1821 res = C1 * ten2k64[exp];
1822 }
1823 }
1824 BID_RETURN (res);
1825 }
1826
1827 /*****************************************************************************
1828 * BID64_to_int32_xint
1829 ****************************************************************************/
1830
1831 #if DECIMAL_CALL_BY_REFERENCE
1832 void
1833 bid64_to_int32_xint (int *pres,
1834 UINT64 *
1835 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1836 _EXC_INFO_PARAM) {
1837 UINT64 x = *px;
1838 #else
1839 int
1840 bid64_to_int32_xint (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1841 _EXC_INFO_PARAM) {
1842 #endif
1843 int res;
1844 UINT64 x_sign;
1845 UINT64 x_exp;
1846 int exp; // unbiased exponent
1847 // Note: C1 represents x_significand (UINT64)
1848 UINT64 tmp64;
1849 BID_UI64DOUBLE tmp1;
1850 unsigned int x_nr_bits;
1851 int q, ind, shift;
1852 UINT64 C1;
1853 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
1854 UINT128 fstar;
1855 UINT128 P128;
1856
1857 // check for NaN or Infinity
1858 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
1859 // set invalid flag
1860 *pfpsf |= INVALID_EXCEPTION;
1861 // return Integer Indefinite
1862 res = 0x80000000;
1863 BID_RETURN (res);
1864 }
1865 // unpack x
1866 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
1867 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
1868 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
1869 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
1870 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
1871 if (C1 > 9999999999999999ull) { // non-canonical
1872 x_exp = 0;
1873 C1 = 0;
1874 }
1875 } else {
1876 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
1877 C1 = x & MASK_BINARY_SIG1;
1878 }
1879
1880 // check for zeros (possibly from non-canonical values)
1881 if (C1 == 0x0ull) {
1882 // x is 0
1883 res = 0x00000000;
1884 BID_RETURN (res);
1885 }
1886 // x is not special and is not zero
1887
1888 // q = nr. of decimal digits in x (1 <= q <= 54)
1889 // determine first the nr. of bits in x
1890 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
1891 // split the 64-bit value in two 32-bit halves to avoid rounding errors
1892 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1893 tmp1.d = (double) (C1 >> 32); // exact conversion
1894 x_nr_bits =
1895 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1896 } else { // x < 2^32
1897 tmp1.d = (double) C1; // exact conversion
1898 x_nr_bits =
1899 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1900 }
1901 } else { // if x < 2^53
1902 tmp1.d = (double) C1; // exact conversion
1903 x_nr_bits =
1904 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1905 }
1906 q = nr_digits[x_nr_bits - 1].digits;
1907 if (q == 0) {
1908 q = nr_digits[x_nr_bits - 1].digits1;
1909 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
1910 q++;
1911 }
1912 exp = x_exp - 398; // unbiased exponent
1913
1914 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1915 // set invalid flag
1916 *pfpsf |= INVALID_EXCEPTION;
1917 // return Integer Indefinite
1918 res = 0x80000000;
1919 BID_RETURN (res);
1920 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1921 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1922 // so x rounded to an integer may or may not fit in a signed 32-bit int
1923 // the cases that do not fit are identified here; the ones that fit
1924 // fall through and will be handled with other cases further,
1925 // under '1 <= q + exp <= 10'
1926 if (x_sign) { // if n < 0 and q + exp = 10
1927 // if n <= -2^31 - 1 then n is too large
1928 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1
1929 // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x50000000a, 1<=q<=16
1930 // <=> C * 10^(11-q) >= 0x50000000a, 1<=q<=16
1931 if (q <= 11) {
1932 // Note: C * 10^(11-q) has 10 or 11 digits; 0x50000000a has 11 digits
1933 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1934 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1935 if (tmp64 >= 0x50000000aull) {
1936 // set invalid flag
1937 *pfpsf |= INVALID_EXCEPTION;
1938 // return Integer Indefinite
1939 res = 0x80000000;
1940 BID_RETURN (res);
1941 }
1942 // else cases that can be rounded to a 32-bit int fall through
1943 // to '1 <= q + exp <= 10'
1944 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1945 // C * 10^(11-q) >= 0x50000000a <=>
1946 // C >= 0x50000000a * 10^(q-11) where 1 <= q - 11 <= 5
1947 // (scale 2^31+1 up)
1948 // Note: 0x50000000a*10^(q-11) has q-1 or q digits, where q <= 16
1949 tmp64 = 0x50000000aull * ten2k64[q - 11];
1950 if (C1 >= tmp64) {
1951 // set invalid flag
1952 *pfpsf |= INVALID_EXCEPTION;
1953 // return Integer Indefinite
1954 res = 0x80000000;
1955 BID_RETURN (res);
1956 }
1957 // else cases that can be rounded to a 32-bit int fall through
1958 // to '1 <= q + exp <= 10'
1959 }
1960 } else { // if n > 0 and q + exp = 10
1961 // if n >= 2^31 then n is too large
1962 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31
1963 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000000, 1<=q<=16
1964 // <=> C * 10^(11-q) >= 0x500000000, 1<=q<=16
1965 if (q <= 11) {
1966 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000000 has 11 digits
1967 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
1968 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1969 if (tmp64 >= 0x500000000ull) {
1970 // set invalid flag
1971 *pfpsf |= INVALID_EXCEPTION;
1972 // return Integer Indefinite
1973 res = 0x80000000;
1974 BID_RETURN (res);
1975 }
1976 // else cases that can be rounded to a 32-bit int fall through
1977 // to '1 <= q + exp <= 10'
1978 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
1979 // C * 10^(11-q) >= 0x500000000 <=>
1980 // C >= 0x500000000 * 10^(q-11) where 1 <= q - 11 <= 5
1981 // (scale 2^31-1 up)
1982 // Note: 0x500000000*10^(q-11) has q-1 or q digits, where q <= 16
1983 tmp64 = 0x500000000ull * ten2k64[q - 11];
1984 if (C1 >= tmp64) {
1985 // set invalid flag
1986 *pfpsf |= INVALID_EXCEPTION;
1987 // return Integer Indefinite
1988 res = 0x80000000;
1989 BID_RETURN (res);
1990 }
1991 // else cases that can be rounded to a 32-bit int fall through
1992 // to '1 <= q + exp <= 10'
1993 }
1994 }
1995 }
1996 // n is not too large to be converted to int32: -2^31 - 1 < n < 2^31
1997 // Note: some of the cases tested for above fall through to this point
1998 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
1999 // set inexact flag
2000 *pfpsf |= INEXACT_EXCEPTION;
2001 // return 0
2002 res = 0x00000000;
2003 BID_RETURN (res);
2004 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
2005 // -2^31-1 < x <= -1 or 1 <= x < 2^31 so x can be rounded
2006 // to nearest to a 32-bit signed integer
2007 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
2008 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
2009 // chop off ind digits from the lower part of C1
2010 // C1 fits in 64 bits
2011 // calculate C* and f*
2012 // C* is actually floor(C*) in this case
2013 // C* and f* need shifting and masking, as shown by
2014 // shiftright128[] and maskhigh128[]
2015 // 1 <= x <= 15
2016 // kx = 10^(-x) = ten2mk64[ind - 1]
2017 // C* = C1 * 10^(-x)
2018 // the approximation of 10^(-x) was rounded up to 54 bits
2019 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
2020 Cstar = P128.w[1];
2021 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
2022 fstar.w[0] = P128.w[0];
2023 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2024 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2025 // C* = floor(C*) (logical right shift; C has p decimal digits,
2026 // correct by Property 1)
2027 // n = C* * 10^(e+x)
2028
2029 // shift right C* by Ex-64 = shiftright128[ind]
2030 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
2031 Cstar = Cstar >> shift;
2032 // determine inexactness of the rounding of C*
2033 // if (0 < f* < 10^(-x)) then
2034 // the result is exact
2035 // else // if (f* > T*) then
2036 // the result is inexact
2037 if (ind - 1 <= 2) {
2038 if (fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
2039 // ten2mk128trunc[ind -1].w[1] is identical to
2040 // ten2mk128[ind -1].w[1]
2041 // set the inexact flag
2042 *pfpsf |= INEXACT_EXCEPTION;
2043 } // else the result is exact
2044 } else { // if 3 <= ind - 1 <= 14
2045 if (fstar.w[1] || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
2046 // ten2mk128trunc[ind -1].w[1] is identical to
2047 // ten2mk128[ind -1].w[1]
2048 // set the inexact flag
2049 *pfpsf |= INEXACT_EXCEPTION;
2050 } // else the result is exact
2051 }
2052
2053 if (x_sign)
2054 res = -Cstar;
2055 else
2056 res = Cstar;
2057 } else if (exp == 0) {
2058 // 1 <= q <= 10
2059 // res = +/-C (exact)
2060 if (x_sign)
2061 res = -C1;
2062 else
2063 res = C1;
2064 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2065 // res = +/-C * 10^exp (exact)
2066 if (x_sign)
2067 res = -C1 * ten2k64[exp];
2068 else
2069 res = C1 * ten2k64[exp];
2070 }
2071 }
2072 BID_RETURN (res);
2073 }
2074
2075 /*****************************************************************************
2076 * BID64_to_int32_rninta
2077 ****************************************************************************/
2078
2079 #if DECIMAL_CALL_BY_REFERENCE
2080 void
2081 bid64_to_int32_rninta (int *pres,
2082 UINT64 *
2083 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2084 _EXC_INFO_PARAM) {
2085 UINT64 x = *px;
2086 #else
2087 int
2088 bid64_to_int32_rninta (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2089 _EXC_INFO_PARAM) {
2090 #endif
2091 int res;
2092 UINT64 x_sign;
2093 UINT64 x_exp;
2094 int exp; // unbiased exponent
2095 // Note: C1 represents x_significand (UINT64)
2096 UINT64 tmp64;
2097 BID_UI64DOUBLE tmp1;
2098 unsigned int x_nr_bits;
2099 int q, ind, shift;
2100 UINT64 C1;
2101 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
2102 UINT128 P128;
2103
2104 // check for NaN or Infinity
2105 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
2106 // set invalid flag
2107 *pfpsf |= INVALID_EXCEPTION;
2108 // return Integer Indefinite
2109 res = 0x80000000;
2110 BID_RETURN (res);
2111 }
2112 // unpack x
2113 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2114 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
2115 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
2116 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
2117 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
2118 if (C1 > 9999999999999999ull) { // non-canonical
2119 x_exp = 0;
2120 C1 = 0;
2121 }
2122 } else {
2123 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
2124 C1 = x & MASK_BINARY_SIG1;
2125 }
2126
2127 // check for zeros (possibly from non-canonical values)
2128 if (C1 == 0x0ull) {
2129 // x is 0
2130 res = 0x00000000;
2131 BID_RETURN (res);
2132 }
2133 // x is not special and is not zero
2134
2135 // q = nr. of decimal digits in x (1 <= q <= 54)
2136 // determine first the nr. of bits in x
2137 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
2138 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2139 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
2140 tmp1.d = (double) (C1 >> 32); // exact conversion
2141 x_nr_bits =
2142 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2143 } else { // x < 2^32
2144 tmp1.d = (double) C1; // exact conversion
2145 x_nr_bits =
2146 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2147 }
2148 } else { // if x < 2^53
2149 tmp1.d = (double) C1; // exact conversion
2150 x_nr_bits =
2151 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2152 }
2153 q = nr_digits[x_nr_bits - 1].digits;
2154 if (q == 0) {
2155 q = nr_digits[x_nr_bits - 1].digits1;
2156 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
2157 q++;
2158 }
2159 exp = x_exp - 398; // unbiased exponent
2160
2161 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2162 // set invalid flag
2163 *pfpsf |= INVALID_EXCEPTION;
2164 // return Integer Indefinite
2165 res = 0x80000000;
2166 BID_RETURN (res);
2167 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2168 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2169 // so x rounded to an integer may or may not fit in a signed 32-bit int
2170 // the cases that do not fit are identified here; the ones that fit
2171 // fall through and will be handled with other cases further,
2172 // under '1 <= q + exp <= 10'
2173 if (x_sign) { // if n < 0 and q + exp = 10
2174 // if n <= -2^31 - 1/2 then n is too large
2175 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1/2
2176 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005, 1<=q<=16
2177 // <=> C * 10^(11-q) >= 0x500000005, 1<=q<=16
2178 if (q <= 11) {
2179 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
2180 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
2181 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2182 if (tmp64 >= 0x500000005ull) {
2183 // set invalid flag
2184 *pfpsf |= INVALID_EXCEPTION;
2185 // return Integer Indefinite
2186 res = 0x80000000;
2187 BID_RETURN (res);
2188 }
2189 // else cases that can be rounded to a 32-bit int fall through
2190 // to '1 <= q + exp <= 10'
2191 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2192 // C * 10^(11-q) >= 0x500000005 <=>
2193 // C >= 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 5
2194 // (scale 2^31+1/2 up)
2195 // Note: 0x500000005*10^(q-11) has q-1 or q digits, where q <= 16
2196 tmp64 = 0x500000005ull * ten2k64[q - 11];
2197 if (C1 >= tmp64) {
2198 // set invalid flag
2199 *pfpsf |= INVALID_EXCEPTION;
2200 // return Integer Indefinite
2201 res = 0x80000000;
2202 BID_RETURN (res);
2203 }
2204 // else cases that can be rounded to a 32-bit int fall through
2205 // to '1 <= q + exp <= 10'
2206 }
2207 } else { // if n > 0 and q + exp = 10
2208 // if n >= 2^31 - 1/2 then n is too large
2209 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
2210 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=16
2211 // <=> C * 10^(11-q) >= 0x4fffffffb, 1<=q<=16
2212 if (q <= 11) {
2213 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
2214 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
2215 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2216 if (tmp64 >= 0x4fffffffbull) {
2217 // set invalid flag
2218 *pfpsf |= INVALID_EXCEPTION;
2219 // return Integer Indefinite
2220 res = 0x80000000;
2221 BID_RETURN (res);
2222 }
2223 // else cases that can be rounded to a 32-bit int fall through
2224 // to '1 <= q + exp <= 10'
2225 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2226 // C * 10^(11-q) >= 0x4fffffffb <=>
2227 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
2228 // (scale 2^31-1/2 up)
2229 // Note: 0x4fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
2230 tmp64 = 0x4fffffffbull * ten2k64[q - 11];
2231 if (C1 >= tmp64) {
2232 // set invalid flag
2233 *pfpsf |= INVALID_EXCEPTION;
2234 // return Integer Indefinite
2235 res = 0x80000000;
2236 BID_RETURN (res);
2237 }
2238 // else cases that can be rounded to a 32-bit int fall through
2239 // to '1 <= q + exp <= 10'
2240 }
2241 }
2242 }
2243 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
2244 // Note: some of the cases tested for above fall through to this point
2245 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2246 // return 0
2247 res = 0x00000000;
2248 BID_RETURN (res);
2249 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2250 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
2251 // res = 0
2252 // else
2253 // res = +/-1
2254 ind = q - 1;
2255 if (C1 < midpoint64[ind]) {
2256 res = 0x00000000; // return 0
2257 } else if (x_sign) { // n < 0
2258 res = 0xffffffff; // return -1
2259 } else { // n > 0
2260 res = 0x00000001; // return +1
2261 }
2262 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
2263 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
2264 // to nearest away to a 32-bit signed integer
2265 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
2266 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
2267 // chop off ind digits from the lower part of C1
2268 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
2269 C1 = C1 + midpoint64[ind - 1];
2270 // calculate C* and f*
2271 // C* is actually floor(C*) in this case
2272 // C* and f* need shifting and masking, as shown by
2273 // shiftright128[] and maskhigh128[]
2274 // 1 <= x <= 15
2275 // kx = 10^(-x) = ten2mk64[ind - 1]
2276 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2277 // the approximation of 10^(-x) was rounded up to 54 bits
2278 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
2279 Cstar = P128.w[1];
2280 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2281 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2282 // C* = floor(C*)-1 (logical right shift; C* has p decimal digits,
2283 // correct by Pr. 1)
2284 // n = C* * 10^(e+x)
2285
2286 // shift right C* by Ex-64 = shiftright128[ind]
2287 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
2288 Cstar = Cstar >> shift;
2289
2290 // if the result was a midpoint it was rounded away from zero
2291 if (x_sign)
2292 res = -Cstar;
2293 else
2294 res = Cstar;
2295 } else if (exp == 0) {
2296 // 1 <= q <= 10
2297 // res = +/-C (exact)
2298 if (x_sign)
2299 res = -C1;
2300 else
2301 res = C1;
2302 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2303 // res = +/-C * 10^exp (exact)
2304 if (x_sign)
2305 res = -C1 * ten2k64[exp];
2306 else
2307 res = C1 * ten2k64[exp];
2308 }
2309 }
2310 BID_RETURN (res);
2311 }
2312
2313 /*****************************************************************************
2314 * BID64_to_int32_xrninta
2315 ****************************************************************************/
2316
2317 #if DECIMAL_CALL_BY_REFERENCE
2318 void
2319 bid64_to_int32_xrninta (int *pres,
2320 UINT64 *
2321 px _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2322 _EXC_INFO_PARAM) {
2323 UINT64 x = *px;
2324 #else
2325 int
2326 bid64_to_int32_xrninta (UINT64 x _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2327 _EXC_INFO_PARAM) {
2328 #endif
2329 int res;
2330 UINT64 x_sign;
2331 UINT64 x_exp;
2332 int exp; // unbiased exponent
2333 // Note: C1 represents x_significand (UINT64)
2334 UINT64 tmp64;
2335 BID_UI64DOUBLE tmp1;
2336 unsigned int x_nr_bits;
2337 int q, ind, shift;
2338 UINT64 C1;
2339 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
2340 UINT128 fstar;
2341 UINT128 P128;
2342
2343 // check for NaN or Infinity
2344 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
2345 // set invalid flag
2346 *pfpsf |= INVALID_EXCEPTION;
2347 // return Integer Indefinite
2348 res = 0x80000000;
2349 BID_RETURN (res);
2350 }
2351 // unpack x
2352 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
2353 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
2354 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
2355 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
2356 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
2357 if (C1 > 9999999999999999ull) { // non-canonical
2358 x_exp = 0;
2359 C1 = 0;
2360 }
2361 } else {
2362 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
2363 C1 = x & MASK_BINARY_SIG1;
2364 }
2365
2366 // check for zeros (possibly from non-canonical values)
2367 if (C1 == 0x0ull) {
2368 // x is 0
2369 res = 0x00000000;
2370 BID_RETURN (res);
2371 }
2372 // x is not special and is not zero
2373
2374 // q = nr. of decimal digits in x (1 <= q <= 54)
2375 // determine first the nr. of bits in x
2376 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
2377 // split the 64-bit value in two 32-bit halves to avoid rounding errors
2378 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
2379 tmp1.d = (double) (C1 >> 32); // exact conversion
2380 x_nr_bits =
2381 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2382 } else { // x < 2^32
2383 tmp1.d = (double) C1; // exact conversion
2384 x_nr_bits =
2385 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2386 }
2387 } else { // if x < 2^53
2388 tmp1.d = (double) C1; // exact conversion
2389 x_nr_bits =
2390 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2391 }
2392 q = nr_digits[x_nr_bits - 1].digits;
2393 if (q == 0) {
2394 q = nr_digits[x_nr_bits - 1].digits1;
2395 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
2396 q++;
2397 }
2398 exp = x_exp - 398; // unbiased exponent
2399
2400 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2401 // set invalid flag
2402 *pfpsf |= INVALID_EXCEPTION;
2403 // return Integer Indefinite
2404 res = 0x80000000;
2405 BID_RETURN (res);
2406 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2407 // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2408 // so x rounded to an integer may or may not fit in a signed 32-bit int
2409 // the cases that do not fit are identified here; the ones that fit
2410 // fall through and will be handled with other cases further,
2411 // under '1 <= q + exp <= 10'
2412 if (x_sign) { // if n < 0 and q + exp = 10
2413 // if n <= -2^31 - 1/2 then n is too large
2414 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31+1/2
2415 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x500000005, 1<=q<=16
2416 // <=> C * 10^(11-q) >= 0x500000005, 1<=q<=16
2417 if (q <= 11) {
2418 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
2419 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
2420 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2421 if (tmp64 >= 0x500000005ull) {
2422 // set invalid flag
2423 *pfpsf |= INVALID_EXCEPTION;
2424 // return Integer Indefinite
2425 res = 0x80000000;
2426 BID_RETURN (res);
2427 }
2428 // else cases that can be rounded to a 32-bit int fall through
2429 // to '1 <= q + exp <= 10'
2430 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2431 // C * 10^(11-q) >= 0x500000005 <=>
2432 // C >= 0x500000005 * 10^(q-11) where 1 <= q - 11 <= 5
2433 // (scale 2^31+1/2 up)
2434 // Note: 0x500000005*10^(q-11) has q-1 or q digits, where q <= 16
2435 tmp64 = 0x500000005ull * ten2k64[q - 11];
2436 if (C1 >= tmp64) {
2437 // set invalid flag
2438 *pfpsf |= INVALID_EXCEPTION;
2439 // return Integer Indefinite
2440 res = 0x80000000;
2441 BID_RETURN (res);
2442 }
2443 // else cases that can be rounded to a 32-bit int fall through
2444 // to '1 <= q + exp <= 10'
2445 }
2446 } else { // if n > 0 and q + exp = 10
2447 // if n >= 2^31 - 1/2 then n is too large
2448 // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^31-1/2
2449 // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x4fffffffb, 1<=q<=16
2450 // <=> C * 10^(11-q) >= 0x4fffffffb, 1<=q<=16
2451 if (q <= 11) {
2452 // Note: C * 10^(11-q) has 10 or 11 digits; 0x500000005 has 11 digits
2453 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
2454 // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2455 if (tmp64 >= 0x4fffffffbull) {
2456 // set invalid flag
2457 *pfpsf |= INVALID_EXCEPTION;
2458 // return Integer Indefinite
2459 res = 0x80000000;
2460 BID_RETURN (res);
2461 }
2462 // else cases that can be rounded to a 32-bit int fall through
2463 // to '1 <= q + exp <= 10'
2464 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
2465 // C * 10^(11-q) >= 0x4fffffffb <=>
2466 // C >= 0x4fffffffb * 10^(q-11) where 1 <= q - 11 <= 5
2467 // (scale 2^31-1/2 up)
2468 // Note: 0x4fffffffb*10^(q-11) has q-1 or q digits, where q <= 16
2469 tmp64 = 0x4fffffffbull * ten2k64[q - 11];
2470 if (C1 >= tmp64) {
2471 // set invalid flag
2472 *pfpsf |= INVALID_EXCEPTION;
2473 // return Integer Indefinite
2474 res = 0x80000000;
2475 BID_RETURN (res);
2476 }
2477 // else cases that can be rounded to a 32-bit int fall through
2478 // to '1 <= q + exp <= 10'
2479 }
2480 }
2481 }
2482 // n is not too large to be converted to int32: -2^31 - 1/2 < n < 2^31 - 1/2
2483 // Note: some of the cases tested for above fall through to this point
2484 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
2485 // set inexact flag
2486 *pfpsf |= INEXACT_EXCEPTION;
2487 // return 0
2488 res = 0x00000000;
2489 BID_RETURN (res);
2490 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
2491 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
2492 // res = 0
2493 // else
2494 // res = +/-1
2495 ind = q - 1;
2496 if (C1 < midpoint64[ind]) {
2497 res = 0x00000000; // return 0
2498 } else if (x_sign) { // n < 0
2499 res = 0xffffffff; // return -1
2500 } else { // n > 0
2501 res = 0x00000001; // return +1
2502 }
2503 // set inexact flag
2504 *pfpsf |= INEXACT_EXCEPTION;
2505 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
2506 // -2^31-1/2 <= x <= -1 or 1 <= x < 2^31-1/2 so x can be rounded
2507 // to nearest away to a 32-bit signed integer
2508 if (exp < 0) { // 2 <= q <= 16, -15 <= exp <= -1, 1 <= q + exp <= 10
2509 ind = -exp; // 1 <= ind <= 15; ind is a synonym for 'x'
2510 // chop off ind digits from the lower part of C1
2511 // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 64 bits
2512 C1 = C1 + midpoint64[ind - 1];
2513 // calculate C* and f*
2514 // C* is actually floor(C*) in this case
2515 // C* and f* need shifting and masking, as shown by
2516 // shiftright128[] and maskhigh128[]
2517 // 1 <= x <= 15
2518 // kx = 10^(-x) = ten2mk64[ind - 1]
2519 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2520 // the approximation of 10^(-x) was rounded up to 54 bits
2521 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
2522 Cstar = P128.w[1];
2523 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
2524 fstar.w[0] = P128.w[0];
2525 // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind].w[0], e.g.
2526 // if x=1, T*=ten2mk128trunc[0].w[0]=0x1999999999999999
2527 // C* = floor(C*)-1 (logical right shift; C* has p decimal digits,
2528 // correct by Pr. 1)
2529 // n = C* * 10^(e+x)
2530
2531 // shift right C* by Ex-64 = shiftright128[ind]
2532 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
2533 Cstar = Cstar >> shift;
2534 // determine inexactness of the rounding of C*
2535 // if (0 < f* - 1/2 < 10^(-x)) then
2536 // the result is exact
2537 // else // if (f* - 1/2 > T*) then
2538 // the result is inexact
2539 if (ind - 1 <= 2) {
2540 if (fstar.w[0] > 0x8000000000000000ull) {
2541 // f* > 1/2 and the result may be exact
2542 tmp64 = fstar.w[0] - 0x8000000000000000ull; // f* - 1/2
2543 if ((tmp64 > ten2mk128trunc[ind - 1].w[1])) {
2544 // ten2mk128trunc[ind -1].w[1] is identical to
2545 // ten2mk128[ind -1].w[1]
2546 // set the inexact flag
2547 *pfpsf |= INEXACT_EXCEPTION;
2548 } // else the result is exact
2549 } else { // the result is inexact; f2* <= 1/2
2550 // set the inexact flag
2551 *pfpsf |= INEXACT_EXCEPTION;
2552 }
2553 } else { // if 3 <= ind - 1 <= 14
2554 if (fstar.w[1] > onehalf128[ind - 1] ||
2555 (fstar.w[1] == onehalf128[ind - 1] && fstar.w[0])) {
2556 // f2* > 1/2 and the result may be exact
2557 // Calculate f2* - 1/2
2558 tmp64 = fstar.w[1] - onehalf128[ind - 1];
2559 if (tmp64 || fstar.w[0] > ten2mk128trunc[ind - 1].w[1]) {
2560 // ten2mk128trunc[ind -1].w[1] is identical to
2561 // ten2mk128[ind -1].w[1]
2562 // set the inexact flag
2563 *pfpsf |= INEXACT_EXCEPTION;
2564 } // else the result is exact
2565 } else { // the result is inexact; f2* <= 1/2
2566 // set the inexact flag
2567 *pfpsf |= INEXACT_EXCEPTION;
2568 }
2569 }
2570
2571 // if the result was a midpoint it was rounded away from zero
2572 if (x_sign)
2573 res = -Cstar;
2574 else
2575 res = Cstar;
2576 } else if (exp == 0) {
2577 // 1 <= q <= 10
2578 // res = +/-C (exact)
2579 if (x_sign)
2580 res = -C1;
2581 else
2582 res = C1;
2583 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2584 // res = +/-C * 10^exp (exact)
2585 if (x_sign)
2586 res = -C1 * ten2k64[exp];
2587 else
2588 res = C1 * ten2k64[exp];
2589 }
2590 }
2591 BID_RETURN (res);
2592 }
Mirror of the offical gcc git repository hosted at gcc.gnu.org