Mercurial > hg > CbC > CbC_gcc

comparison libgcc/config/libbid/bid64_to_int64.c @ 0:a06113de4d67

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