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