0
|
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 }
|