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