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