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