]>
Commit | Line | Data |
---|---|---|
7adcbafe | 1 | /* Copyright (C) 2007-2022 Free Software Foundation, Inc. |
200359e8 L |
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 | |
748086b7 | 7 | Software Foundation; either version 3, or (at your option) any later |
200359e8 L |
8 | version. |
9 | ||
200359e8 L |
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 | ||
748086b7 JJ |
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/>. */ | |
200359e8 L |
23 | |
24 | #include "bid_internal.h" | |
25 | ||
26 | /***************************************************************************** | |
27 | * BID64_to_uint64_rnint | |
28 | ****************************************************************************/ | |
29 | ||
30 | #if DECIMAL_CALL_BY_REFERENCE | |
31 | void | |
b2a00c89 | 32 | bid64_to_uint64_rnint (UINT64 * pres, UINT64 * px |
200359e8 L |
33 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
34 | _EXC_INFO_PARAM) { | |
35 | UINT64 x = *px; | |
36 | #else | |
37 | UINT64 | |
b2a00c89 | 38 | bid64_to_uint64_rnint (UINT64 x |
200359e8 L |
39 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
40 | _EXC_INFO_PARAM) { | |
41 | #endif | |
42 | UINT64 res; | |
43 | UINT64 x_sign; | |
44 | UINT64 x_exp; | |
b2a00c89 | 45 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 52 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 65 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 68 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 69 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 70 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
71 | x_exp = 0; |
72 | C1 = 0; | |
73 | } | |
74 | } else { | |
b2a00c89 | 75 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | |
b2a00c89 | 89 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 90 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
91 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
92 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
93 | x_nr_bits = |
94 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
95 | } else { // x < 2^32 |
96 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
97 | x_nr_bits = |
98 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
99 | } | |
b2a00c89 L |
100 | } else { // if x < 2^53 |
101 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
102 | x_nr_bits = |
103 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
104 | } | |
b2a00c89 | 105 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 106 | if (q == 0) { |
b2a00c89 L |
107 | q = nr_digits[x_nr_bits - 1].digits1; |
108 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
109 | q++; |
110 | } | |
b2a00c89 | 111 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 112 | |
b2a00c89 | 113 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
114 | // set invalid flag |
115 | *pfpsf |= INVALID_EXCEPTION; | |
116 | // return Integer Indefinite | |
117 | res = 0x8000000000000000ull; | |
118 | BID_RETURN (res); | |
b2a00c89 | 119 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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' | |
b2a00c89 | 125 | if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1/2 |
200359e8 L |
126 | // => set invalid flag |
127 | *pfpsf |= INVALID_EXCEPTION; | |
128 | // return Integer Indefinite | |
129 | res = 0x8000000000000000ull; | |
130 | BID_RETURN (res); | |
b2a00c89 | 131 | } else { // if n > 0 and q + exp = 20 |
200359e8 L |
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 | |
b2a00c89 | 139 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 150 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
151 | // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb |
152 | // has 21 digits | |
b2a00c89 | 153 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 169 | if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) |
200359e8 L |
170 | // return 0 |
171 | res = 0x0000000000000000ull; | |
172 | BID_RETURN (res); | |
b2a00c89 | 173 | } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) |
200359e8 L |
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 | |
b2a00c89 L |
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 | |
200359e8 L |
186 | res = 0x8000000000000000ull; |
187 | *pfpsf |= INVALID_EXCEPTION; | |
188 | BID_RETURN (res); | |
189 | } | |
b2a00c89 | 190 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
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 | |
b2a00c89 | 193 | if (x_sign) { // x <= -1 |
200359e8 L |
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 | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 206 | C1 = C1 + midpoint64[ind - 1]; |
200359e8 L |
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 | |
b2a00c89 | 210 | // shiftright128[] and maskhigh128[] |
200359e8 | 211 | // 1 <= x <= 15 |
b2a00c89 | 212 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
213 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) |
214 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 215 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 216 | Cstar = P128.w[1]; |
b2a00c89 | 217 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 | 218 | fstar.w[0] = P128.w[0]; |
b2a00c89 L |
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 | |
200359e8 L |
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 | ||
b2a00c89 L |
231 | // shift right C* by Ex-64 = shiftright128[ind] |
232 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 L |
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] && | |
b2a00c89 L |
239 | (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) { |
240 | // ten2mk128trunc[ind -1].w[1] is identical to | |
241 | // ten2mk128[ind -1].w[1] | |
200359e8 | 242 | // the result is a midpoint; round to nearest |
b2a00c89 | 243 | if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD] |
200359e8 | 244 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 |
b2a00c89 L |
245 | Cstar--; // Cstar is now even |
246 | } // else MP in [ODD, EVEN] | |
200359e8 | 247 | } |
b2a00c89 | 248 | res = Cstar; // the result is positive |
200359e8 L |
249 | } else if (exp == 0) { |
250 | // 1 <= q <= 10 | |
251 | // res = +C (exact) | |
b2a00c89 L |
252 | res = C1; // the result is positive |
253 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 254 | // res = +C * 10^exp (exact) |
b2a00c89 | 255 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
256 | } |
257 | } | |
258 | BID_RETURN (res); | |
259 | } | |
260 | ||
261 | /***************************************************************************** | |
262 | * BID64_to_uint64_xrnint | |
263 | ****************************************************************************/ | |
264 | ||
265 | #if DECIMAL_CALL_BY_REFERENCE | |
266 | void | |
b2a00c89 | 267 | bid64_to_uint64_xrnint (UINT64 * pres, UINT64 * px |
200359e8 L |
268 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
269 | _EXC_INFO_PARAM) { | |
270 | UINT64 x = *px; | |
271 | #else | |
272 | UINT64 | |
b2a00c89 | 273 | bid64_to_uint64_xrnint (UINT64 x |
200359e8 L |
274 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
275 | _EXC_INFO_PARAM) { | |
276 | #endif | |
277 | UINT64 res; | |
278 | UINT64 x_sign; | |
279 | UINT64 x_exp; | |
b2a00c89 | 280 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 288 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 301 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 304 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 305 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 306 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
307 | x_exp = 0; |
308 | C1 = 0; | |
309 | } | |
310 | } else { | |
b2a00c89 | 311 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | |
b2a00c89 | 325 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 326 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
327 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
328 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
329 | x_nr_bits = |
330 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
331 | } else { // x < 2^32 |
332 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
333 | x_nr_bits = |
334 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
335 | } | |
b2a00c89 L |
336 | } else { // if x < 2^53 |
337 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
338 | x_nr_bits = |
339 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
340 | } | |
b2a00c89 | 341 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 342 | if (q == 0) { |
b2a00c89 L |
343 | q = nr_digits[x_nr_bits - 1].digits1; |
344 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
345 | q++; |
346 | } | |
b2a00c89 | 347 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 348 | |
b2a00c89 | 349 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
350 | // set invalid flag |
351 | *pfpsf |= INVALID_EXCEPTION; | |
352 | // return Integer Indefinite | |
353 | res = 0x8000000000000000ull; | |
354 | BID_RETURN (res); | |
b2a00c89 | 355 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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' | |
b2a00c89 | 361 | if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1/2 |
200359e8 L |
362 | // => set invalid flag |
363 | *pfpsf |= INVALID_EXCEPTION; | |
364 | // return Integer Indefinite | |
365 | res = 0x8000000000000000ull; | |
366 | BID_RETURN (res); | |
b2a00c89 | 367 | } else { // if n > 0 and q + exp = 20 |
200359e8 L |
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 | |
b2a00c89 | 375 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 386 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
387 | // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb |
388 | // has 21 digits | |
b2a00c89 | 389 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 405 | if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) |
200359e8 L |
406 | // set inexact flag |
407 | *pfpsf |= INEXACT_EXCEPTION; | |
408 | // return 0 | |
409 | res = 0x0000000000000000ull; | |
410 | BID_RETURN (res); | |
b2a00c89 | 411 | } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) |
200359e8 L |
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 | |
b2a00c89 L |
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 | |
200359e8 L |
424 | res = 0x8000000000000000ull; |
425 | *pfpsf |= INVALID_EXCEPTION; | |
426 | BID_RETURN (res); | |
427 | } | |
428 | // set inexact flag | |
429 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 | 430 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
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 | |
b2a00c89 | 433 | if (x_sign) { // x <= -1 |
200359e8 L |
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 | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 446 | C1 = C1 + midpoint64[ind - 1]; |
200359e8 L |
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 | |
b2a00c89 | 450 | // shiftright128[] and maskhigh128[] |
200359e8 | 451 | // 1 <= x <= 15 |
b2a00c89 | 452 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
453 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) |
454 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 455 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 456 | Cstar = P128.w[1]; |
b2a00c89 | 457 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 | 458 | fstar.w[0] = P128.w[0]; |
b2a00c89 L |
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 | |
200359e8 L |
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 | ||
b2a00c89 L |
471 | // shift right C* by Ex-64 = shiftright128[ind] |
472 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 L |
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 | |
b2a00c89 | 479 | if (ind - 1 <= 2) { // fstar.w[1] is 0 |
200359e8 L |
480 | if (fstar.w[0] > 0x8000000000000000ull) { |
481 | // f* > 1/2 and the result may be exact | |
b2a00c89 L |
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] | |
200359e8 L |
486 | // set the inexact flag |
487 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
488 | } // else the result is exact |
489 | } else { // the result is inexact; f2* <= 1/2 | |
200359e8 L |
490 | // set the inexact flag |
491 | *pfpsf |= INEXACT_EXCEPTION; | |
492 | } | |
b2a00c89 L |
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])) { | |
200359e8 L |
496 | // f2* > 1/2 and the result may be exact |
497 | // Calculate f2* - 1/2 | |
b2a00c89 L |
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] | |
200359e8 L |
502 | // set the inexact flag |
503 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
504 | } // else the result is exact |
505 | } else { // the result is inexact; f2* <= 1/2 | |
200359e8 L |
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] && | |
b2a00c89 L |
515 | (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) { |
516 | // ten2mk128trunc[ind -1].w[1] is identical to | |
517 | // ten2mk128[ind -1].w[1] | |
200359e8 | 518 | // the result is a midpoint; round to nearest |
b2a00c89 | 519 | if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD] |
200359e8 | 520 | // if floor(C*) is odd C = floor(C*) - 1; the result >= 1 |
b2a00c89 L |
521 | Cstar--; // Cstar is now even |
522 | } // else MP in [ODD, EVEN] | |
200359e8 | 523 | } |
b2a00c89 | 524 | res = Cstar; // the result is positive |
200359e8 L |
525 | } else if (exp == 0) { |
526 | // 1 <= q <= 10 | |
527 | // res = +C (exact) | |
b2a00c89 L |
528 | res = C1; // the result is positive |
529 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 530 | // res = +C * 10^exp (exact) |
b2a00c89 | 531 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
532 | } |
533 | } | |
534 | BID_RETURN (res); | |
535 | } | |
536 | ||
537 | /***************************************************************************** | |
538 | * BID64_to_uint64_floor | |
539 | ****************************************************************************/ | |
540 | ||
541 | #if DECIMAL_CALL_BY_REFERENCE | |
542 | void | |
b2a00c89 | 543 | bid64_to_uint64_floor (UINT64 * pres, UINT64 * px |
200359e8 L |
544 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
545 | _EXC_INFO_PARAM) { | |
546 | UINT64 x = *px; | |
547 | #else | |
548 | UINT64 | |
b2a00c89 | 549 | bid64_to_uint64_floor (UINT64 x |
200359e8 L |
550 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
551 | _EXC_INFO_PARAM) { | |
552 | #endif | |
553 | UINT64 res; | |
554 | UINT64 x_sign; | |
555 | UINT64 x_exp; | |
b2a00c89 | 556 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 563 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 575 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 578 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 579 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 580 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
581 | x_exp = 0; |
582 | C1 = 0; | |
583 | } | |
584 | } else { | |
b2a00c89 | 585 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | ||
b2a00c89 | 597 | if (x_sign) { // if n < 0 the conversion is invalid |
200359e8 L |
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 | |
b2a00c89 | 606 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 607 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
608 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
609 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
610 | x_nr_bits = |
611 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
612 | } else { // x < 2^32 |
613 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
614 | x_nr_bits = |
615 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
616 | } | |
b2a00c89 L |
617 | } else { // if x < 2^53 |
618 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
619 | x_nr_bits = |
620 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
621 | } | |
b2a00c89 | 622 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 623 | if (q == 0) { |
b2a00c89 L |
624 | q = nr_digits[x_nr_bits - 1].digits1; |
625 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
626 | q++; |
627 | } | |
b2a00c89 | 628 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 629 | |
b2a00c89 | 630 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
631 | // set invalid flag |
632 | *pfpsf |= INVALID_EXCEPTION; | |
633 | // return Integer Indefinite | |
634 | res = 0x8000000000000000ull; | |
635 | BID_RETURN (res); | |
b2a00c89 | 636 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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 | |
b2a00c89 | 650 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 661 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
662 | // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000 |
663 | // has 21 digits | |
b2a00c89 | 664 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 679 | if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1) |
200359e8 L |
680 | // return 0 |
681 | res = 0x0000000000000000ull; | |
682 | BID_RETURN (res); | |
b2a00c89 | 683 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
684 | // 1 <= x < 2^64 so x can be rounded |
685 | // to nearest to a 64-bit unsigned signed integer | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 693 | // shiftright128[] and maskhigh128[] |
200359e8 | 694 | // 1 <= x <= 15 |
b2a00c89 | 695 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
696 | // C* = C1 * 10^(-x) |
697 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 698 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 699 | Cstar = P128.w[1]; |
b2a00c89 L |
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 | |
200359e8 L |
702 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
703 | // correct by Property 1) | |
704 | // n = C* * 10^(e+x) | |
705 | ||
b2a00c89 L |
706 | // shift right C* by Ex-64 = shiftright128[ind] |
707 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 | 708 | Cstar = Cstar >> shift; |
b2a00c89 | 709 | res = Cstar; // the result is positive |
200359e8 L |
710 | } else if (exp == 0) { |
711 | // 1 <= q <= 10 | |
712 | // res = +C (exact) | |
b2a00c89 L |
713 | res = C1; // the result is positive |
714 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 715 | // res = +C * 10^exp (exact) |
b2a00c89 | 716 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
717 | } |
718 | } | |
719 | BID_RETURN (res); | |
720 | } | |
721 | ||
722 | /***************************************************************************** | |
723 | * BID64_to_uint64_xfloor | |
724 | ****************************************************************************/ | |
725 | ||
726 | #if DECIMAL_CALL_BY_REFERENCE | |
727 | void | |
b2a00c89 | 728 | bid64_to_uint64_xfloor (UINT64 * pres, UINT64 * px |
200359e8 L |
729 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
730 | _EXC_INFO_PARAM) { | |
731 | UINT64 x = *px; | |
732 | #else | |
733 | UINT64 | |
b2a00c89 | 734 | bid64_to_uint64_xfloor (UINT64 x |
200359e8 L |
735 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
736 | _EXC_INFO_PARAM) { | |
737 | #endif | |
738 | UINT64 res; | |
739 | UINT64 x_sign; | |
740 | UINT64 x_exp; | |
b2a00c89 | 741 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 748 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 761 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 764 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 765 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 766 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
767 | x_exp = 0; |
768 | C1 = 0; | |
769 | } | |
770 | } else { | |
b2a00c89 | 771 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | ||
b2a00c89 | 783 | if (x_sign) { // if n < 0 the conversion is invalid |
200359e8 L |
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 | |
b2a00c89 | 792 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 793 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
794 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
795 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
796 | x_nr_bits = |
797 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
798 | } else { // x < 2^32 |
799 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
800 | x_nr_bits = |
801 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
802 | } | |
b2a00c89 L |
803 | } else { // if x < 2^53 |
804 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
805 | x_nr_bits = |
806 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
807 | } | |
b2a00c89 | 808 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 809 | if (q == 0) { |
b2a00c89 L |
810 | q = nr_digits[x_nr_bits - 1].digits1; |
811 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
812 | q++; |
813 | } | |
b2a00c89 | 814 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 815 | |
b2a00c89 | 816 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
817 | // set invalid flag |
818 | *pfpsf |= INVALID_EXCEPTION; | |
819 | // return Integer Indefinite | |
820 | res = 0x8000000000000000ull; | |
821 | BID_RETURN (res); | |
b2a00c89 | 822 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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 | |
b2a00c89 | 836 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 847 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
848 | // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000 |
849 | // has 21 digits | |
b2a00c89 | 850 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 865 | if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1) |
200359e8 L |
866 | // set inexact flag |
867 | *pfpsf |= INEXACT_EXCEPTION; | |
868 | // return 0 | |
869 | res = 0x0000000000000000ull; | |
870 | BID_RETURN (res); | |
b2a00c89 | 871 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
872 | // 1 <= x < 2^64 so x can be rounded |
873 | // to nearest to a 64-bit unsigned signed integer | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 881 | // shiftright128[] and maskhigh128[] |
200359e8 | 882 | // 1 <= x <= 15 |
b2a00c89 | 883 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
884 | // C* = C1 * 10^(-x) |
885 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 886 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 887 | Cstar = P128.w[1]; |
b2a00c89 | 888 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 | 889 | fstar.w[0] = P128.w[0]; |
b2a00c89 L |
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 | |
200359e8 L |
892 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
893 | // correct by Property 1) | |
894 | // n = C* * 10^(e+x) | |
895 | ||
b2a00c89 L |
896 | // shift right C* by Ex-64 = shiftright128[ind] |
897 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 L |
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 | |
b2a00c89 L |
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] | |
200359e8 L |
908 | // set the inexact flag |
909 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
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] | |
200359e8 L |
915 | // set the inexact flag |
916 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 | 917 | } // else the result is exact |
200359e8 L |
918 | } |
919 | ||
b2a00c89 | 920 | res = Cstar; // the result is positive |
200359e8 L |
921 | } else if (exp == 0) { |
922 | // 1 <= q <= 10 | |
923 | // res = +C (exact) | |
b2a00c89 L |
924 | res = C1; // the result is positive |
925 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 926 | // res = +C * 10^exp (exact) |
b2a00c89 | 927 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
928 | } |
929 | } | |
930 | BID_RETURN (res); | |
931 | } | |
932 | ||
933 | /***************************************************************************** | |
934 | * BID64_to_uint64_ceil | |
935 | ****************************************************************************/ | |
936 | ||
937 | #if DECIMAL_CALL_BY_REFERENCE | |
938 | void | |
b2a00c89 | 939 | bid64_to_uint64_ceil (UINT64 * pres, UINT64 * px |
200359e8 L |
940 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
941 | _EXC_INFO_PARAM) { | |
942 | UINT64 x = *px; | |
943 | #else | |
944 | UINT64 | |
b2a00c89 | 945 | bid64_to_uint64_ceil (UINT64 x |
200359e8 L |
946 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
947 | _EXC_INFO_PARAM) { | |
948 | #endif | |
949 | UINT64 res; | |
950 | UINT64 x_sign; | |
951 | UINT64 x_exp; | |
b2a00c89 | 952 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 959 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 972 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 975 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 976 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 977 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
978 | x_exp = 0; |
979 | C1 = 0; | |
980 | } | |
981 | } else { | |
b2a00c89 | 982 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | |
b2a00c89 | 996 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 997 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
998 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
999 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
1000 | x_nr_bits = |
1001 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
1002 | } else { // x < 2^32 |
1003 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1004 | x_nr_bits = |
1005 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1006 | } | |
b2a00c89 L |
1007 | } else { // if x < 2^53 |
1008 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1009 | x_nr_bits = |
1010 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1011 | } | |
b2a00c89 | 1012 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 1013 | if (q == 0) { |
b2a00c89 L |
1014 | q = nr_digits[x_nr_bits - 1].digits1; |
1015 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
1016 | q++; |
1017 | } | |
b2a00c89 | 1018 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 1019 | |
b2a00c89 | 1020 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
1021 | // set invalid flag |
1022 | *pfpsf |= INVALID_EXCEPTION; | |
1023 | // return Integer Indefinite | |
1024 | res = 0x8000000000000000ull; | |
1025 | BID_RETURN (res); | |
b2a00c89 | 1026 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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' | |
b2a00c89 | 1032 | if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1 |
200359e8 L |
1033 | // => set invalid flag |
1034 | *pfpsf |= INVALID_EXCEPTION; | |
1035 | // return Integer Indefinite | |
1036 | res = 0x8000000000000000ull; | |
1037 | BID_RETURN (res); | |
b2a00c89 | 1038 | } else { // if n > 0 and q + exp = 20 |
200359e8 L |
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 | |
b2a00c89 | 1046 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 1057 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
1058 | // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffff6 |
1059 | // has 21 digits | |
b2a00c89 | 1060 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 1076 | if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) |
200359e8 L |
1077 | // return 0 or 1 |
1078 | if (x_sign) | |
1079 | res = 0x0000000000000000ull; | |
1080 | else | |
1081 | res = 0x0000000000000001ull; | |
1082 | BID_RETURN (res); | |
b2a00c89 | 1083 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
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 | |
b2a00c89 | 1086 | if (x_sign) { // x <= -1 |
200359e8 L |
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 | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 1102 | // shiftright128[] and maskhigh128[] |
200359e8 | 1103 | // 1 <= x <= 15 |
b2a00c89 | 1104 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
1105 | // C* = C1 * 10^(-x) |
1106 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 1107 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 1108 | Cstar = P128.w[1]; |
b2a00c89 | 1109 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 | 1110 | fstar.w[0] = P128.w[0]; |
b2a00c89 L |
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 | |
200359e8 L |
1113 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
1114 | // correct by Property 1) | |
1115 | // n = C* * 10^(e+x) | |
1116 | ||
b2a00c89 L |
1117 | // shift right C* by Ex-64 = shiftright128[ind] |
1118 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 L |
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 | |
b2a00c89 L |
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] | |
200359e8 | 1129 | Cstar++; |
b2a00c89 L |
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] | |
200359e8 | 1135 | Cstar++; |
b2a00c89 | 1136 | } // else the result is exact |
200359e8 L |
1137 | } |
1138 | ||
b2a00c89 | 1139 | res = Cstar; // the result is positive |
200359e8 L |
1140 | } else if (exp == 0) { |
1141 | // 1 <= q <= 10 | |
1142 | // res = +C (exact) | |
b2a00c89 L |
1143 | res = C1; // the result is positive |
1144 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 1145 | // res = +C * 10^exp (exact) |
b2a00c89 | 1146 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
1147 | } |
1148 | } | |
1149 | BID_RETURN (res); | |
1150 | } | |
1151 | ||
1152 | /***************************************************************************** | |
1153 | * BID64_to_uint64_xceil | |
1154 | ****************************************************************************/ | |
1155 | ||
1156 | #if DECIMAL_CALL_BY_REFERENCE | |
1157 | void | |
b2a00c89 | 1158 | bid64_to_uint64_xceil (UINT64 * pres, UINT64 * px |
200359e8 L |
1159 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
1160 | _EXC_INFO_PARAM) { | |
1161 | UINT64 x = *px; | |
1162 | #else | |
1163 | UINT64 | |
b2a00c89 | 1164 | bid64_to_uint64_xceil (UINT64 x |
200359e8 L |
1165 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
1166 | _EXC_INFO_PARAM) { | |
1167 | #endif | |
1168 | UINT64 res; | |
1169 | UINT64 x_sign; | |
1170 | UINT64 x_exp; | |
b2a00c89 | 1171 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 1178 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 1191 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 1194 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 1195 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 1196 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
1197 | x_exp = 0; |
1198 | C1 = 0; | |
1199 | } | |
1200 | } else { | |
b2a00c89 | 1201 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | |
b2a00c89 | 1215 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 1216 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
1217 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
1218 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
1219 | x_nr_bits = |
1220 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
1221 | } else { // x < 2^32 |
1222 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1223 | x_nr_bits = |
1224 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1225 | } | |
b2a00c89 L |
1226 | } else { // if x < 2^53 |
1227 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1228 | x_nr_bits = |
1229 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1230 | } | |
b2a00c89 | 1231 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 1232 | if (q == 0) { |
b2a00c89 L |
1233 | q = nr_digits[x_nr_bits - 1].digits1; |
1234 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
1235 | q++; |
1236 | } | |
b2a00c89 | 1237 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 1238 | |
b2a00c89 | 1239 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
1240 | // set invalid flag |
1241 | *pfpsf |= INVALID_EXCEPTION; | |
1242 | // return Integer Indefinite | |
1243 | res = 0x8000000000000000ull; | |
1244 | BID_RETURN (res); | |
b2a00c89 | 1245 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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' | |
b2a00c89 | 1251 | if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1 |
200359e8 L |
1252 | // => set invalid flag |
1253 | *pfpsf |= INVALID_EXCEPTION; | |
1254 | // return Integer Indefinite | |
1255 | res = 0x8000000000000000ull; | |
1256 | BID_RETURN (res); | |
b2a00c89 | 1257 | } else { // if n > 0 and q + exp = 20 |
200359e8 L |
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 | |
b2a00c89 | 1265 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 1276 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
1277 | // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffff6 |
1278 | // has 21 digits | |
b2a00c89 | 1279 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 1295 | if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) |
200359e8 L |
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); | |
b2a00c89 | 1304 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
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 | |
b2a00c89 | 1307 | if (x_sign) { // x <= -1 |
200359e8 L |
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 | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 1323 | // shiftright128[] and maskhigh128[] |
200359e8 | 1324 | // 1 <= x <= 15 |
b2a00c89 | 1325 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
1326 | // C* = C1 * 10^(-x) |
1327 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 1328 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 1329 | Cstar = P128.w[1]; |
b2a00c89 | 1330 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 | 1331 | fstar.w[0] = P128.w[0]; |
b2a00c89 L |
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 | |
200359e8 L |
1334 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
1335 | // correct by Property 1) | |
1336 | // n = C* * 10^(e+x) | |
1337 | ||
b2a00c89 L |
1338 | // shift right C* by Ex-64 = shiftright128[ind] |
1339 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 L |
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 | |
b2a00c89 L |
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] | |
200359e8 L |
1350 | Cstar++; |
1351 | // set the inexact flag | |
1352 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
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] | |
200359e8 L |
1358 | Cstar++; |
1359 | // set the inexact flag | |
1360 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 | 1361 | } // else the result is exact |
200359e8 L |
1362 | } |
1363 | ||
b2a00c89 | 1364 | res = Cstar; // the result is positive |
200359e8 L |
1365 | } else if (exp == 0) { |
1366 | // 1 <= q <= 10 | |
1367 | // res = +C (exact) | |
b2a00c89 L |
1368 | res = C1; // the result is positive |
1369 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 1370 | // res = +C * 10^exp (exact) |
b2a00c89 | 1371 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
1372 | } |
1373 | } | |
1374 | BID_RETURN (res); | |
1375 | } | |
1376 | ||
1377 | /***************************************************************************** | |
1378 | * BID64_to_uint64_int | |
1379 | ****************************************************************************/ | |
1380 | ||
1381 | #if DECIMAL_CALL_BY_REFERENCE | |
1382 | void | |
b2a00c89 | 1383 | bid64_to_uint64_int (UINT64 * pres, UINT64 * px |
200359e8 L |
1384 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) |
1385 | { | |
1386 | UINT64 x = *px; | |
1387 | #else | |
1388 | UINT64 | |
b2a00c89 | 1389 | bid64_to_uint64_int (UINT64 x |
200359e8 L |
1390 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) |
1391 | { | |
1392 | #endif | |
1393 | UINT64 res; | |
1394 | UINT64 x_sign; | |
1395 | UINT64 x_exp; | |
b2a00c89 | 1396 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 1403 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 1415 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 1418 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 1419 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 1420 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
1421 | x_exp = 0; |
1422 | C1 = 0; | |
1423 | } | |
1424 | } else { | |
b2a00c89 | 1425 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | |
b2a00c89 | 1439 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 1440 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
1441 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
1442 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
1443 | x_nr_bits = |
1444 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
1445 | } else { // x < 2^32 |
1446 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1447 | x_nr_bits = |
1448 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1449 | } | |
b2a00c89 L |
1450 | } else { // if x < 2^53 |
1451 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1452 | x_nr_bits = |
1453 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1454 | } | |
b2a00c89 | 1455 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 1456 | if (q == 0) { |
b2a00c89 L |
1457 | q = nr_digits[x_nr_bits - 1].digits1; |
1458 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
1459 | q++; |
1460 | } | |
b2a00c89 | 1461 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 1462 | |
b2a00c89 | 1463 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
1464 | // set invalid flag |
1465 | *pfpsf |= INVALID_EXCEPTION; | |
1466 | // return Integer Indefinite | |
1467 | res = 0x8000000000000000ull; | |
1468 | BID_RETURN (res); | |
b2a00c89 | 1469 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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' | |
b2a00c89 | 1475 | if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1 |
200359e8 L |
1476 | // => set invalid flag |
1477 | *pfpsf |= INVALID_EXCEPTION; | |
1478 | // return Integer Indefinite | |
1479 | res = 0x8000000000000000ull; | |
1480 | BID_RETURN (res); | |
b2a00c89 | 1481 | } else { // if n > 0 and q + exp = 20 |
200359e8 L |
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 | |
b2a00c89 | 1489 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 1500 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
1501 | // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000 |
1502 | // has 21 digits | |
b2a00c89 | 1503 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 1519 | if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) |
200359e8 L |
1520 | // return 0 |
1521 | res = 0x0000000000000000ull; | |
1522 | BID_RETURN (res); | |
b2a00c89 | 1523 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
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 | |
b2a00c89 | 1526 | if (x_sign) { // x <= -1 |
200359e8 L |
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 | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 1542 | // shiftright128[] and maskhigh128[] |
200359e8 | 1543 | // 1 <= x <= 15 |
b2a00c89 | 1544 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
1545 | // C* = C1 * 10^(-x) |
1546 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 1547 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 1548 | Cstar = P128.w[1]; |
b2a00c89 L |
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 | |
200359e8 L |
1551 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
1552 | // correct by Property 1) | |
1553 | // n = C* * 10^(e+x) | |
1554 | ||
b2a00c89 L |
1555 | // shift right C* by Ex-64 = shiftright128[ind] |
1556 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 | 1557 | Cstar = Cstar >> shift; |
b2a00c89 | 1558 | res = Cstar; // the result is positive |
200359e8 L |
1559 | } else if (exp == 0) { |
1560 | // 1 <= q <= 10 | |
1561 | // res = +C (exact) | |
b2a00c89 L |
1562 | res = C1; // the result is positive |
1563 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 1564 | // res = +C * 10^exp (exact) |
b2a00c89 | 1565 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
1566 | } |
1567 | } | |
1568 | BID_RETURN (res); | |
1569 | } | |
1570 | ||
1571 | /***************************************************************************** | |
1572 | * BID64_to_uint64_xint | |
1573 | ****************************************************************************/ | |
1574 | ||
1575 | #if DECIMAL_CALL_BY_REFERENCE | |
1576 | void | |
b2a00c89 | 1577 | bid64_to_uint64_xint (UINT64 * pres, UINT64 * px |
200359e8 L |
1578 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
1579 | _EXC_INFO_PARAM) { | |
1580 | UINT64 x = *px; | |
1581 | #else | |
1582 | UINT64 | |
b2a00c89 | 1583 | bid64_to_uint64_xint (UINT64 x |
200359e8 L |
1584 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
1585 | _EXC_INFO_PARAM) { | |
1586 | #endif | |
1587 | UINT64 res; | |
1588 | UINT64 x_sign; | |
1589 | UINT64 x_exp; | |
b2a00c89 | 1590 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 1597 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 1610 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 1613 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 1614 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 1615 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
1616 | x_exp = 0; |
1617 | C1 = 0; | |
1618 | } | |
1619 | } else { | |
b2a00c89 | 1620 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | |
b2a00c89 | 1634 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 1635 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
1636 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
1637 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
1638 | x_nr_bits = |
1639 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
1640 | } else { // x < 2^32 |
1641 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1642 | x_nr_bits = |
1643 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1644 | } | |
b2a00c89 L |
1645 | } else { // if x < 2^53 |
1646 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1647 | x_nr_bits = |
1648 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1649 | } | |
b2a00c89 | 1650 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 1651 | if (q == 0) { |
b2a00c89 L |
1652 | q = nr_digits[x_nr_bits - 1].digits1; |
1653 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
1654 | q++; |
1655 | } | |
b2a00c89 | 1656 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 1657 | |
b2a00c89 | 1658 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
1659 | // set invalid flag |
1660 | *pfpsf |= INVALID_EXCEPTION; | |
1661 | // return Integer Indefinite | |
1662 | res = 0x8000000000000000ull; | |
1663 | BID_RETURN (res); | |
b2a00c89 | 1664 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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' | |
b2a00c89 | 1670 | if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1 |
200359e8 L |
1671 | // => set invalid flag |
1672 | *pfpsf |= INVALID_EXCEPTION; | |
1673 | // return Integer Indefinite | |
1674 | res = 0x8000000000000000ull; | |
1675 | BID_RETURN (res); | |
b2a00c89 | 1676 | } else { // if n > 0 and q + exp = 20 |
200359e8 L |
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 | |
b2a00c89 | 1684 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 1695 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
1696 | // Note: C * 10^(21-q) has 20 or 21 digits; 0xa0000000000000000 |
1697 | // has 21 digits | |
b2a00c89 | 1698 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 1714 | if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1) |
200359e8 L |
1715 | // set inexact flag |
1716 | *pfpsf |= INEXACT_EXCEPTION; | |
1717 | // return 0 | |
1718 | res = 0x0000000000000000ull; | |
1719 | BID_RETURN (res); | |
b2a00c89 | 1720 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
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 | |
b2a00c89 | 1723 | if (x_sign) { // x <= -1 |
200359e8 L |
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 | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 1739 | // shiftright128[] and maskhigh128[] |
200359e8 | 1740 | // 1 <= x <= 15 |
b2a00c89 | 1741 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
1742 | // C* = C1 * 10^(-x) |
1743 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 1744 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 1745 | Cstar = P128.w[1]; |
b2a00c89 | 1746 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 | 1747 | fstar.w[0] = P128.w[0]; |
b2a00c89 L |
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 | |
200359e8 L |
1750 | // C* = floor(C*) (logical right shift; C has p decimal digits, |
1751 | // correct by Property 1) | |
1752 | // n = C* * 10^(e+x) | |
1753 | ||
b2a00c89 L |
1754 | // shift right C* by Ex-64 = shiftright128[ind] |
1755 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 L |
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 | |
b2a00c89 L |
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] | |
200359e8 L |
1766 | // set the inexact flag |
1767 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
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] | |
200359e8 L |
1773 | // set the inexact flag |
1774 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 | 1775 | } // else the result is exact |
200359e8 L |
1776 | } |
1777 | ||
b2a00c89 | 1778 | res = Cstar; // the result is positive |
200359e8 L |
1779 | } else if (exp == 0) { |
1780 | // 1 <= q <= 10 | |
1781 | // res = +C (exact) | |
b2a00c89 L |
1782 | res = C1; // the result is positive |
1783 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 1784 | // res = +C * 10^exp (exact) |
b2a00c89 | 1785 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
1786 | } |
1787 | } | |
1788 | BID_RETURN (res); | |
1789 | } | |
1790 | ||
1791 | /***************************************************************************** | |
1792 | * BID64_to_uint64_rninta | |
1793 | ****************************************************************************/ | |
1794 | ||
1795 | #if DECIMAL_CALL_BY_REFERENCE | |
1796 | void | |
b2a00c89 | 1797 | bid64_to_uint64_rninta (UINT64 * pres, UINT64 * px |
200359e8 L |
1798 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
1799 | _EXC_INFO_PARAM) { | |
1800 | UINT64 x = *px; | |
1801 | #else | |
1802 | UINT64 | |
b2a00c89 | 1803 | bid64_to_uint64_rninta (UINT64 x |
200359e8 L |
1804 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
1805 | _EXC_INFO_PARAM) { | |
1806 | #endif | |
1807 | UINT64 res; | |
1808 | UINT64 x_sign; | |
1809 | UINT64 x_exp; | |
b2a00c89 | 1810 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 1817 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 1829 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 1832 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 1833 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 1834 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
1835 | x_exp = 0; |
1836 | C1 = 0; | |
1837 | } | |
1838 | } else { | |
b2a00c89 | 1839 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | |
b2a00c89 | 1853 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 1854 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
1855 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
1856 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
1857 | x_nr_bits = |
1858 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
1859 | } else { // x < 2^32 |
1860 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1861 | x_nr_bits = |
1862 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1863 | } | |
b2a00c89 L |
1864 | } else { // if x < 2^53 |
1865 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
1866 | x_nr_bits = |
1867 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
1868 | } | |
b2a00c89 | 1869 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 1870 | if (q == 0) { |
b2a00c89 L |
1871 | q = nr_digits[x_nr_bits - 1].digits1; |
1872 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
1873 | q++; |
1874 | } | |
b2a00c89 | 1875 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 1876 | |
b2a00c89 | 1877 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
1878 | // set invalid flag |
1879 | *pfpsf |= INVALID_EXCEPTION; | |
1880 | // return Integer Indefinite | |
1881 | res = 0x8000000000000000ull; | |
1882 | BID_RETURN (res); | |
b2a00c89 | 1883 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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' | |
b2a00c89 | 1889 | if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1/2 |
200359e8 L |
1890 | // => set invalid flag |
1891 | *pfpsf |= INVALID_EXCEPTION; | |
1892 | // return Integer Indefinite | |
1893 | res = 0x8000000000000000ull; | |
1894 | BID_RETURN (res); | |
b2a00c89 | 1895 | } else { // if n > 0 and q + exp = 20 |
200359e8 L |
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 | |
b2a00c89 | 1903 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 1914 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
1915 | // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb |
1916 | // has 21 digits | |
b2a00c89 | 1917 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 1933 | if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) |
200359e8 L |
1934 | // return 0 |
1935 | res = 0x0000000000000000ull; | |
1936 | BID_RETURN (res); | |
b2a00c89 | 1937 | } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) |
200359e8 L |
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 | |
b2a00c89 L |
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 | |
200359e8 L |
1950 | res = 0x8000000000000000ull; |
1951 | *pfpsf |= INVALID_EXCEPTION; | |
1952 | BID_RETURN (res); | |
1953 | } | |
b2a00c89 | 1954 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
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 | |
b2a00c89 | 1957 | if (x_sign) { // x <= -1 |
200359e8 L |
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 | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 1970 | C1 = C1 + midpoint64[ind - 1]; |
200359e8 L |
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 | |
b2a00c89 | 1974 | // shiftright128[] and maskhigh128[] |
200359e8 | 1975 | // 1 <= x <= 15 |
b2a00c89 | 1976 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
1977 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) |
1978 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 1979 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 1980 | Cstar = P128.w[1]; |
b2a00c89 L |
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 | |
200359e8 L |
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 | ||
b2a00c89 L |
1993 | // shift right C* by Ex-64 = shiftright128[ind] |
1994 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 L |
1995 | Cstar = Cstar >> shift; |
1996 | ||
1997 | // if the result was a midpoint it was rounded away from zero | |
b2a00c89 | 1998 | res = Cstar; // the result is positive |
200359e8 L |
1999 | } else if (exp == 0) { |
2000 | // 1 <= q <= 10 | |
2001 | // res = +C (exact) | |
b2a00c89 L |
2002 | res = C1; // the result is positive |
2003 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 2004 | // res = +C * 10^exp (exact) |
b2a00c89 | 2005 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
2006 | } |
2007 | } | |
2008 | BID_RETURN (res); | |
2009 | } | |
2010 | ||
2011 | /***************************************************************************** | |
2012 | * BID64_to_uint64_xrninta | |
2013 | ****************************************************************************/ | |
2014 | ||
2015 | #if DECIMAL_CALL_BY_REFERENCE | |
2016 | void | |
b2a00c89 | 2017 | bid64_to_uint64_xrninta (UINT64 * pres, UINT64 * px |
200359e8 L |
2018 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
2019 | _EXC_INFO_PARAM) { | |
2020 | UINT64 x = *px; | |
2021 | #else | |
2022 | UINT64 | |
b2a00c89 | 2023 | bid64_to_uint64_xrninta (UINT64 x |
200359e8 L |
2024 | _EXC_FLAGS_PARAM _EXC_MASKS_PARAM |
2025 | _EXC_INFO_PARAM) { | |
2026 | #endif | |
2027 | UINT64 res; | |
2028 | UINT64 x_sign; | |
2029 | UINT64 x_exp; | |
b2a00c89 | 2030 | int exp; // unbiased exponent |
200359e8 L |
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; | |
b2a00c89 | 2038 | UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits |
200359e8 L |
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 | |
b2a00c89 | 2051 | x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative |
200359e8 L |
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) { | |
b2a00c89 | 2054 | x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased |
200359e8 | 2055 | C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2; |
b2a00c89 | 2056 | if (C1 > 9999999999999999ull) { // non-canonical |
200359e8 L |
2057 | x_exp = 0; |
2058 | C1 = 0; | |
2059 | } | |
2060 | } else { | |
b2a00c89 | 2061 | x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased |
200359e8 L |
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 | |
b2a00c89 | 2075 | if (C1 >= 0x0020000000000000ull) { // x >= 2^53 |
200359e8 | 2076 | // split the 64-bit value in two 32-bit halves to avoid rounding errors |
b2a00c89 L |
2077 | if (C1 >= 0x0000000100000000ull) { // x >= 2^32 |
2078 | tmp1.d = (double) (C1 >> 32); // exact conversion | |
200359e8 L |
2079 | x_nr_bits = |
2080 | 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
b2a00c89 L |
2081 | } else { // x < 2^32 |
2082 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
2083 | x_nr_bits = |
2084 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
2085 | } | |
b2a00c89 L |
2086 | } else { // if x < 2^53 |
2087 | tmp1.d = (double) C1; // exact conversion | |
200359e8 L |
2088 | x_nr_bits = |
2089 | 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff); | |
2090 | } | |
b2a00c89 | 2091 | q = nr_digits[x_nr_bits - 1].digits; |
200359e8 | 2092 | if (q == 0) { |
b2a00c89 L |
2093 | q = nr_digits[x_nr_bits - 1].digits1; |
2094 | if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo) | |
200359e8 L |
2095 | q++; |
2096 | } | |
b2a00c89 | 2097 | exp = x_exp - 398; // unbiased exponent |
200359e8 | 2098 | |
b2a00c89 | 2099 | if ((q + exp) > 20) { // x >= 10^20 ~= 2^66.45... (cannot fit in 64 bits) |
200359e8 L |
2100 | // set invalid flag |
2101 | *pfpsf |= INVALID_EXCEPTION; | |
2102 | // return Integer Indefinite | |
2103 | res = 0x8000000000000000ull; | |
2104 | BID_RETURN (res); | |
b2a00c89 | 2105 | } else if ((q + exp) == 20) { // x = c(0)c(1)...c(19).c(20)...c(q-1) |
200359e8 L |
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' | |
b2a00c89 | 2111 | if (x_sign) { // if n < 0 and q + exp = 20 then x is much less than -1/2 |
200359e8 L |
2112 | // => set invalid flag |
2113 | *pfpsf |= INVALID_EXCEPTION; | |
2114 | // return Integer Indefinite | |
2115 | res = 0x8000000000000000ull; | |
2116 | BID_RETURN (res); | |
b2a00c89 | 2117 | } else { // if n > 0 and q + exp = 20 |
200359e8 L |
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 | |
b2a00c89 | 2125 | __mul_128x64_to_128 (C, C1, ten2k128[0]); // 10^20 * C |
200359e8 L |
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' | |
b2a00c89 | 2136 | } else { // if (2 <= q <= 16) => 5 <= 21 - q <= 19 |
200359e8 L |
2137 | // Note: C * 10^(21-q) has 20 or 21 digits; 0x9fffffffffffffffb |
2138 | // has 21 digits | |
b2a00c89 | 2139 | __mul_64x64_to_128MACH (C, C1, ten2k64[21 - q]); |
200359e8 L |
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 | |
b2a00c89 | 2155 | if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1) |
200359e8 L |
2156 | // set inexact flag |
2157 | *pfpsf |= INEXACT_EXCEPTION; | |
2158 | // return 0 | |
2159 | res = 0x0000000000000000ull; | |
2160 | BID_RETURN (res); | |
b2a00c89 | 2161 | } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1) |
200359e8 L |
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 | |
b2a00c89 L |
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 | |
200359e8 L |
2174 | res = 0x8000000000000000ull; |
2175 | *pfpsf |= INVALID_EXCEPTION; | |
2176 | BID_RETURN (res); | |
2177 | } | |
2178 | // set inexact flag | |
2179 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 | 2180 | } else { // if (1 <= q + exp <= 20, 1 <= q <= 16, -15 <= exp <= 19) |
200359e8 L |
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 | |
b2a00c89 | 2183 | if (x_sign) { // x <= -1 |
200359e8 L |
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 | |
b2a00c89 L |
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' | |
200359e8 L |
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 | |
b2a00c89 | 2196 | C1 = C1 + midpoint64[ind - 1]; |
200359e8 L |
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 | |
b2a00c89 | 2200 | // shiftright128[] and maskhigh128[] |
200359e8 | 2201 | // 1 <= x <= 15 |
b2a00c89 | 2202 | // kx = 10^(-x) = ten2mk64[ind - 1] |
200359e8 L |
2203 | // C* = (C1 + 1/2 * 10^x) * 10^(-x) |
2204 | // the approximation of 10^(-x) was rounded up to 54 bits | |
b2a00c89 | 2205 | __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]); |
200359e8 | 2206 | Cstar = P128.w[1]; |
b2a00c89 | 2207 | fstar.w[1] = P128.w[1] & maskhigh128[ind - 1]; |
200359e8 | 2208 | fstar.w[0] = P128.w[0]; |
b2a00c89 L |
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 | |
200359e8 L |
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 | ||
b2a00c89 L |
2221 | // shift right C* by Ex-64 = shiftright128[ind] |
2222 | shift = shiftright128[ind - 1]; // 0 <= shift <= 39 | |
200359e8 L |
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 | |
b2a00c89 | 2229 | if (ind - 1 <= 2) { // fstar.w[1] is 0 |
200359e8 L |
2230 | if (fstar.w[0] > 0x8000000000000000ull) { |
2231 | // f* > 1/2 and the result may be exact | |
b2a00c89 L |
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] | |
200359e8 L |
2236 | // set the inexact flag |
2237 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
2238 | } // else the result is exact |
2239 | } else { // the result is inexact; f2* <= 1/2 | |
200359e8 L |
2240 | // set the inexact flag |
2241 | *pfpsf |= INEXACT_EXCEPTION; | |
2242 | } | |
b2a00c89 L |
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])) { | |
200359e8 L |
2246 | // f2* > 1/2 and the result may be exact |
2247 | // Calculate f2* - 1/2 | |
b2a00c89 L |
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] | |
200359e8 L |
2252 | // set the inexact flag |
2253 | *pfpsf |= INEXACT_EXCEPTION; | |
b2a00c89 L |
2254 | } // else the result is exact |
2255 | } else { // the result is inexact; f2* <= 1/2 | |
200359e8 L |
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 | |
b2a00c89 | 2262 | res = Cstar; // the result is positive |
200359e8 L |
2263 | } else if (exp == 0) { |
2264 | // 1 <= q <= 10 | |
2265 | // res = +C (exact) | |
b2a00c89 L |
2266 | res = C1; // the result is positive |
2267 | } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10 | |
200359e8 | 2268 | // res = +C * 10^exp (exact) |
b2a00c89 | 2269 | res = C1 * ten2k64[exp]; // the result is positive |
200359e8 L |
2270 | } |
2271 | } | |
2272 | BID_RETURN (res); | |
2273 | } |