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