|
0
|
1 /* Copyright (C) 2007, 2009 Free Software Foundation, Inc.
|
|
|
2
|
|
|
3 This file is part of GCC.
|
|
|
4
|
|
|
5 GCC is free software; you can redistribute it and/or modify it under
|
|
|
6 the terms of the GNU General Public License as published by the Free
|
|
|
7 Software Foundation; either version 3, or (at your option) any later
|
|
|
8 version.
|
|
|
9
|
|
|
10 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
|
11 WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
|
12 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
|
13 for more details.
|
|
|
14
|
|
|
15 Under Section 7 of GPL version 3, you are granted additional
|
|
|
16 permissions described in the GCC Runtime Library Exception, version
|
|
|
17 3.1, as published by the Free Software Foundation.
|
|
|
18
|
|
|
19 You should have received a copy of the GNU General Public License and
|
|
|
20 a copy of the GCC Runtime Library Exception along with this program;
|
|
|
21 see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
|
|
|
22 <http://www.gnu.org/licenses/>. */
|
|
|
23
|
|
|
24 #include "bid_internal.h"
|
|
|
25
|
|
|
26 /*****************************************************************************
|
|
|
27 * BID64_to_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 }
|
