]> git.ipfire.org Git - thirdparty/gcc.git/blame - libgcc/config/libbid/bid64_to_uint32.c
Update copyright years.
[thirdparty/gcc.git] / libgcc / config / libbid / bid64_to_uint32.c
CommitLineData
8d9254fc 1/* Copyright (C) 2007-2020 Free Software Foundation, Inc.
200359e8
L
2
3This file is part of GCC.
4
5GCC is free software; you can redistribute it and/or modify it under
6the terms of the GNU General Public License as published by the Free
748086b7 7Software Foundation; either version 3, or (at your option) any later
200359e8
L
8version.
9
200359e8
L
10GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11WARRANTY; without even the implied warranty of MERCHANTABILITY or
12FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
13for more details.
14
748086b7
JJ
15Under Section 7 of GPL version 3, you are granted additional
16permissions described in the GCC Runtime Library Exception, version
173.1, as published by the Free Software Foundation.
18
19You should have received a copy of the GNU General Public License and
20a copy of the GCC Runtime Library Exception along with this program;
21see the files COPYING3 and COPYING.RUNTIME respectively. If not, see
22<http://www.gnu.org/licenses/>. */
200359e8
L
23
24#include "bid_internal.h"
25
26/*****************************************************************************
27 * BID64_to_uint32_rnint
28 ****************************************************************************/
29
30#if DECIMAL_CALL_BY_REFERENCE
31void
b2a00c89 32bid64_to_uint32_rnint (unsigned int *pres, UINT64 * px
200359e8
L
33 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
34 _EXC_INFO_PARAM) {
35 UINT64 x = *px;
36#else
37unsigned int
b2a00c89 38bid64_to_uint32_rnint (UINT64 x
200359e8
L
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;
b2a00c89 45 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 52 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
53 UINT128 fstar;
54 UINT128 P128;
55
56 // check for NaN or Infinity
57 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
58 // set invalid flag
59 *pfpsf |= INVALID_EXCEPTION;
60 // return Integer Indefinite
61 res = 0x80000000;
62 BID_RETURN (res);
63 }
64 // unpack x
b2a00c89 65 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
66 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
67 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
b2a00c89 68 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 69 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 70 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
71 x_exp = 0;
72 C1 = 0;
73 }
74 } else {
b2a00c89 75 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
76 C1 = x & MASK_BINARY_SIG1;
77 }
78
79 // check for zeros (possibly from non-canonical values)
80 if (C1 == 0x0ull) {
81 // x is 0
82 res = 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
b2a00c89 89 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 90 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
91 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
92 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
93 x_nr_bits =
94 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
95 } else { // x < 2^32
96 tmp1.d = (double) C1; // exact conversion
200359e8
L
97 x_nr_bits =
98 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
99 }
b2a00c89
L
100 } else { // if x < 2^53
101 tmp1.d = (double) C1; // exact conversion
200359e8
L
102 x_nr_bits =
103 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
104 }
b2a00c89 105 q = nr_digits[x_nr_bits - 1].digits;
200359e8 106 if (q == 0) {
b2a00c89
L
107 q = nr_digits[x_nr_bits - 1].digits1;
108 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
109 q++;
110 }
b2a00c89 111 exp = x_exp - 398; // unbiased exponent
200359e8 112
b2a00c89 113 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
114 // set invalid flag
115 *pfpsf |= INVALID_EXCEPTION;
116 // return Integer Indefinite
117 res = 0x80000000;
118 BID_RETURN (res);
b2a00c89 119 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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'
b2a00c89 125 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1/2
200359e8
L
126 // => set invalid flag
127 *pfpsf |= INVALID_EXCEPTION;
128 // return Integer Indefinite
129 res = 0x80000000;
130 BID_RETURN (res);
b2a00c89 131 } else { // if n > 0 and q + exp = 10
200359e8
L
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
b2a00c89 138 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 149 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 154 tmp64 = 0x9fffffffbull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 169 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
200359e8
L
170 // return 0
171 res = 0x00000000;
172 BID_RETURN (res);
b2a00c89 173 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
200359e8
L
174 // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
175 // res = 0
176 // else if x > 0
177 // res = +1
178 // else // if x < 0
179 // invalid exc
180 ind = q - 1;
b2a00c89
L
181 if (C1 <= midpoint64[ind]) {
182 res = 0x00000000; // return 0
183 } else if (x_sign) { // n < 0
200359e8
L
184 // set invalid flag
185 *pfpsf |= INVALID_EXCEPTION;
186 // return Integer Indefinite
187 res = 0x80000000;
188 BID_RETURN (res);
b2a00c89
L
189 } else { // n > 0
190 res = 0x00000001; // return +1
200359e8 191 }
b2a00c89 192 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
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
b2a00c89 195 if (x_sign) { // x <= -1
200359e8
L
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
b2a00c89
L
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'
200359e8
L
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
b2a00c89 208 C1 = C1 + midpoint64[ind - 1];
200359e8
L
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
b2a00c89 212 // shiftright128[] and maskhigh128[]
200359e8 213 // 1 <= x <= 15
b2a00c89 214 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
215 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
216 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 217 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 218 Cstar = P128.w[1];
b2a00c89 219 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
200359e8 220 fstar.w[0] = P128.w[0];
b2a00c89
L
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
200359e8
L
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
b2a00c89
L
233 // shift right C* by Ex-64 = shiftright128[ind]
234 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
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] &&
b2a00c89
L
241 (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) {
242 // ten2mk128trunc[ind -1].w[1] is identical to
243 // ten2mk128[ind -1].w[1]
200359e8 244 // the result is a midpoint; round to nearest
b2a00c89 245 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
200359e8 246 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
b2a00c89
L
247 Cstar--; // Cstar is now even
248 } // else MP in [ODD, EVEN]
200359e8 249 }
b2a00c89 250 res = Cstar; // the result is positive
200359e8
L
251 } else if (exp == 0) {
252 // 1 <= q <= 10
253 // res = +C (exact)
b2a00c89
L
254 res = C1; // the result is positive
255 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 256 // res = +C * 10^exp (exact)
b2a00c89 257 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
258 }
259 }
260 BID_RETURN (res);
261}
262
263/*****************************************************************************
264 * BID64_to_uint32_xrnint
265 ****************************************************************************/
266
267#if DECIMAL_CALL_BY_REFERENCE
268void
b2a00c89 269bid64_to_uint32_xrnint (unsigned int *pres, UINT64 * px
200359e8
L
270 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
271 _EXC_INFO_PARAM) {
272 UINT64 x = *px;
273#else
274unsigned int
b2a00c89 275bid64_to_uint32_xrnint (UINT64 x
200359e8
L
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;
b2a00c89 282 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 289 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
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
b2a00c89 302 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
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) {
b2a00c89 305 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 306 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 307 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
308 x_exp = 0;
309 C1 = 0;
310 }
311 } else {
b2a00c89 312 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
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
b2a00c89 326 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 327 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
328 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
329 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
330 x_nr_bits =
331 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
332 } else { // x < 2^32
333 tmp1.d = (double) C1; // exact conversion
200359e8
L
334 x_nr_bits =
335 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
336 }
b2a00c89
L
337 } else { // if x < 2^53
338 tmp1.d = (double) C1; // exact conversion
200359e8
L
339 x_nr_bits =
340 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
341 }
b2a00c89 342 q = nr_digits[x_nr_bits - 1].digits;
200359e8 343 if (q == 0) {
b2a00c89
L
344 q = nr_digits[x_nr_bits - 1].digits1;
345 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
346 q++;
347 }
b2a00c89 348 exp = x_exp - 398; // unbiased exponent
200359e8 349
b2a00c89 350 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
351 // set invalid flag
352 *pfpsf |= INVALID_EXCEPTION;
353 // return Integer Indefinite
354 res = 0x80000000;
355 BID_RETURN (res);
b2a00c89 356 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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'
b2a00c89 362 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1/2
200359e8
L
363 // => set invalid flag
364 *pfpsf |= INVALID_EXCEPTION;
365 // return Integer Indefinite
366 res = 0x80000000;
367 BID_RETURN (res);
b2a00c89 368 } else { // if n > 0 and q + exp = 10
200359e8
L
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
b2a00c89 375 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 386 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 391 tmp64 = 0x9fffffffbull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 406 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
200359e8
L
407 // set inexact flag
408 *pfpsf |= INEXACT_EXCEPTION;
409 // return 0
410 res = 0x00000000;
411 BID_RETURN (res);
b2a00c89 412 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
200359e8
L
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;
b2a00c89
L
420 if (C1 <= midpoint64[ind]) {
421 res = 0x00000000; // return 0
422 } else if (x_sign) { // n < 0
200359e8
L
423 // set invalid flag
424 *pfpsf |= INVALID_EXCEPTION;
425 // return Integer Indefinite
426 res = 0x80000000;
427 BID_RETURN (res);
b2a00c89
L
428 } else { // n > 0
429 res = 0x00000001; // return +1
200359e8
L
430 }
431 // set inexact flag
432 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89 433 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
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
b2a00c89 436 if (x_sign) { // x <= -1
200359e8
L
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
b2a00c89
L
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'
200359e8
L
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
b2a00c89 449 C1 = C1 + midpoint64[ind - 1];
200359e8
L
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
b2a00c89 453 // shiftright128[] and maskhigh128[]
200359e8 454 // 1 <= x <= 15
b2a00c89 455 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
456 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
457 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 458 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 459 Cstar = P128.w[1];
b2a00c89 460 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
200359e8 461 fstar.w[0] = P128.w[0];
b2a00c89
L
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
200359e8
L
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
b2a00c89
L
474 // shift right C* by Ex-64 = shiftright128[ind]
475 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
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
b2a00c89 482 if (ind - 1 <= 2) { // fstar.w[1] is 0
200359e8
L
483 if (fstar.w[0] > 0x8000000000000000ull) {
484 // f* > 1/2 and the result may be exact
b2a00c89
L
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]
200359e8
L
489 // set the inexact flag
490 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L
491 } // else the result is exact
492 } else { // the result is inexact; f2* <= 1/2
200359e8
L
493 // set the inexact flag
494 *pfpsf |= INEXACT_EXCEPTION;
495 }
b2a00c89
L
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])) {
200359e8
L
499 // f2* > 1/2 and the result may be exact
500 // Calculate f2* - 1/2
b2a00c89
L
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]
200359e8
L
505 // set the inexact flag
506 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L
507 } // else the result is exact
508 } else { // the result is inexact; f2* <= 1/2
200359e8
L
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] &&
b2a00c89
L
518 (fstar.w[0] <= ten2mk128trunc[ind - 1].w[1])) {
519 // ten2mk128trunc[ind -1].w[1] is identical to
520 // ten2mk128[ind -1].w[1]
200359e8 521 // the result is a midpoint; round to nearest
b2a00c89 522 if (Cstar & 0x01) { // Cstar is odd; MP in [EVEN, ODD]
200359e8 523 // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
b2a00c89
L
524 Cstar--; // Cstar is now even
525 } // else MP in [ODD, EVEN]
200359e8 526 }
b2a00c89 527 res = Cstar; // the result is positive
200359e8
L
528 } else if (exp == 0) {
529 // 1 <= q <= 10
530 // res = +C (exact)
b2a00c89
L
531 res = C1; // the result is positive
532 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 533 // res = +C * 10^exp (exact)
b2a00c89 534 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
535 }
536 }
537 BID_RETURN (res);
538}
539
540/*****************************************************************************
541 * BID64_to_uint32_floor
542 ****************************************************************************/
543
544#if DECIMAL_CALL_BY_REFERENCE
545void
b2a00c89 546bid64_to_uint32_floor (unsigned int *pres, UINT64 * px
200359e8
L
547 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
548 _EXC_INFO_PARAM) {
549 UINT64 x = *px;
550#else
551unsigned int
b2a00c89 552bid64_to_uint32_floor (UINT64 x
200359e8
L
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;
b2a00c89 559 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 566 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
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
b2a00c89 578 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
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) {
b2a00c89 581 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 582 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 583 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
584 x_exp = 0;
585 C1 = 0;
586 }
587 } else {
b2a00c89 588 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
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
b2a00c89 600 if (x_sign) { // if n < 0 the conversion is invalid
200359e8
L
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
b2a00c89 609 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 610 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
611 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
612 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
613 x_nr_bits =
614 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
615 } else { // x < 2^32
616 tmp1.d = (double) C1; // exact conversion
200359e8
L
617 x_nr_bits =
618 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
619 }
b2a00c89
L
620 } else { // if x < 2^53
621 tmp1.d = (double) C1; // exact conversion
200359e8
L
622 x_nr_bits =
623 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
624 }
b2a00c89 625 q = nr_digits[x_nr_bits - 1].digits;
200359e8 626 if (q == 0) {
b2a00c89
L
627 q = nr_digits[x_nr_bits - 1].digits1;
628 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
629 q++;
630 }
b2a00c89 631 exp = x_exp - 398; // unbiased exponent
200359e8 632
b2a00c89 633 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
634 // set invalid flag
635 *pfpsf |= INVALID_EXCEPTION;
636 // return Integer Indefinite
637 res = 0x80000000;
638 BID_RETURN (res);
b2a00c89 639 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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
b2a00c89 652 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 663 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 668 tmp64 = 0xa00000000ull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 682 if ((q + exp) <= 0) { // n = +0.[0...0]c(0)c(1)...c(q-1)
200359e8
L
683 // return 0
684 res = 0x00000000;
685 BID_RETURN (res);
b2a00c89 686 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
687 // 1 <= x < 2^32 so x can be rounded
688 // to nearest to a 32-bit unsigned integer
b2a00c89
L
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'
200359e8
L
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
b2a00c89 696 // shiftright128[] and maskhigh128[]
200359e8 697 // 1 <= x <= 15
b2a00c89 698 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
699 // C* = C1 * 10^(-x)
700 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 701 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 702 Cstar = P128.w[1];
b2a00c89
L
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
200359e8
L
705 // C* = floor(C*) (logical right shift; C has p decimal digits,
706 // correct by Property 1)
707 // n = C* * 10^(e+x)
708
b2a00c89
L
709 // shift right C* by Ex-64 = shiftright128[ind]
710 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
711 Cstar = Cstar >> shift;
712
b2a00c89 713 res = Cstar; // the result is positive
200359e8
L
714 } else if (exp == 0) {
715 // 1 <= q <= 10
716 // res = +C (exact)
b2a00c89
L
717 res = C1; // the result is positive
718 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 719 // res = +C * 10^exp (exact)
b2a00c89 720 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
721 }
722 }
723 BID_RETURN (res);
724}
725
726/*****************************************************************************
727 * BID64_to_uint32_xfloor
728 ****************************************************************************/
729
730#if DECIMAL_CALL_BY_REFERENCE
731void
b2a00c89 732bid64_to_uint32_xfloor (unsigned int *pres, UINT64 * px
200359e8
L
733 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
734 _EXC_INFO_PARAM) {
735 UINT64 x = *px;
736#else
737unsigned int
b2a00c89 738bid64_to_uint32_xfloor (UINT64 x
200359e8
L
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;
b2a00c89 745 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 752 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
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
b2a00c89 765 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
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) {
b2a00c89 768 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 769 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 770 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
771 x_exp = 0;
772 C1 = 0;
773 }
774 } else {
b2a00c89 775 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
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
b2a00c89 787 if (x_sign) { // if n < 0 the conversion is invalid
200359e8
L
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
b2a00c89 796 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 797 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
798 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
799 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
800 x_nr_bits =
801 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
802 } else { // x < 2^32
803 tmp1.d = (double) C1; // exact conversion
200359e8
L
804 x_nr_bits =
805 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
806 }
b2a00c89
L
807 } else { // if x < 2^53
808 tmp1.d = (double) C1; // exact conversion
200359e8
L
809 x_nr_bits =
810 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
811 }
b2a00c89 812 q = nr_digits[x_nr_bits - 1].digits;
200359e8 813 if (q == 0) {
b2a00c89
L
814 q = nr_digits[x_nr_bits - 1].digits1;
815 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
816 q++;
817 }
b2a00c89 818 exp = x_exp - 398; // unbiased exponent
200359e8 819
b2a00c89 820 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
821 // set invalid flag
822 *pfpsf |= INVALID_EXCEPTION;
823 // return Integer Indefinite
824 res = 0x80000000;
825 BID_RETURN (res);
b2a00c89 826 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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
b2a00c89 839 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 850 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 855 tmp64 = 0xa00000000ull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 869 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
200359e8
L
870 // set inexact flag
871 *pfpsf |= INEXACT_EXCEPTION;
872 // return 0
873 res = 0x00000000;
874 BID_RETURN (res);
b2a00c89 875 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
876 // 1 <= x < 2^32 so x can be rounded
877 // to nearest to a 32-bit unsigned integer
b2a00c89
L
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'
200359e8
L
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
b2a00c89 885 // shiftright128[] and maskhigh128[]
200359e8 886 // 1 <= x <= 15
b2a00c89 887 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
888 // C* = C1 * 10^(-x)
889 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 890 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 891 Cstar = P128.w[1];
b2a00c89 892 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
200359e8 893 fstar.w[0] = P128.w[0];
b2a00c89
L
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
200359e8
L
896 // C* = floor(C*) (logical right shift; C has p decimal digits,
897 // correct by Property 1)
898 // n = C* * 10^(e+x)
899
b2a00c89
L
900 // shift right C* by Ex-64 = shiftright128[ind]
901 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
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) {
b2a00c89
L
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]
200359e8
L
912 // set the inexact flag
913 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L
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]
200359e8
L
919 // set the inexact flag
920 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89 921 } // else the result is exact
200359e8
L
922 }
923
b2a00c89 924 res = Cstar; // the result is positive
200359e8
L
925 } else if (exp == 0) {
926 // 1 <= q <= 10
927 // res = +C (exact)
b2a00c89
L
928 res = C1; // the result is positive
929 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 930 // res = +C * 10^exp (exact)
b2a00c89 931 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
932 }
933 }
934 BID_RETURN (res);
935}
936
937/*****************************************************************************
938 * BID64_to_uint32_ceil
939 ****************************************************************************/
940
941#if DECIMAL_CALL_BY_REFERENCE
942void
b2a00c89 943bid64_to_uint32_ceil (unsigned int *pres, UINT64 * px
200359e8
L
944 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
945 _EXC_INFO_PARAM) {
946 UINT64 x = *px;
947#else
948unsigned int
b2a00c89 949bid64_to_uint32_ceil (UINT64 x
200359e8
L
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;
b2a00c89 956 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 963 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
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
b2a00c89 976 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
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) {
b2a00c89 979 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 980 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 981 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
982 x_exp = 0;
983 C1 = 0;
984 }
985 } else {
b2a00c89 986 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
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
b2a00c89 1000 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 1001 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
1002 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1003 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
1004 x_nr_bits =
1005 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
1006 } else { // x < 2^32
1007 tmp1.d = (double) C1; // exact conversion
200359e8
L
1008 x_nr_bits =
1009 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1010 }
b2a00c89
L
1011 } else { // if x < 2^53
1012 tmp1.d = (double) C1; // exact conversion
200359e8
L
1013 x_nr_bits =
1014 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1015 }
b2a00c89 1016 q = nr_digits[x_nr_bits - 1].digits;
200359e8 1017 if (q == 0) {
b2a00c89
L
1018 q = nr_digits[x_nr_bits - 1].digits1;
1019 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
1020 q++;
1021 }
b2a00c89 1022 exp = x_exp - 398; // unbiased exponent
200359e8 1023
b2a00c89 1024 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
1025 // set invalid flag
1026 *pfpsf |= INVALID_EXCEPTION;
1027 // return Integer Indefinite
1028 res = 0x80000000;
1029 BID_RETURN (res);
b2a00c89 1030 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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'
b2a00c89 1036 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1
200359e8
L
1037 // => set invalid flag
1038 *pfpsf |= INVALID_EXCEPTION;
1039 // return Integer Indefinite
1040 res = 0x80000000;
1041 BID_RETURN (res);
b2a00c89 1042 } else { // if n > 0 and q + exp = 10
200359e8
L
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
b2a00c89 1049 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 1060 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 1065 tmp64 = 0x9fffffff6ull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 1080 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
200359e8
L
1081 // return 0 or 1
1082 if (x_sign)
1083 res = 0x00000000;
1084 else
1085 res = 0x00000001;
1086 BID_RETURN (res);
b2a00c89 1087 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
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
b2a00c89 1090 if (x_sign) { // x <= -1
200359e8
L
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
b2a00c89
L
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'
200359e8
L
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
b2a00c89 1106 // shiftright128[] and maskhigh128[]
200359e8 1107 // 1 <= x <= 15
b2a00c89 1108 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
1109 // C* = C1 * 10^(-x)
1110 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 1111 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 1112 Cstar = P128.w[1];
b2a00c89 1113 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
200359e8 1114 fstar.w[0] = P128.w[0];
b2a00c89
L
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
200359e8
L
1117 // C* = floor(C*) (logical right shift; C has p decimal digits,
1118 // correct by Property 1)
1119 // n = C* * 10^(e+x)
1120
b2a00c89
L
1121 // shift right C* by Ex-64 = shiftright128[ind]
1122 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
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
b2a00c89
L
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]
200359e8 1133 Cstar++;
b2a00c89
L
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]
200359e8 1139 Cstar++;
b2a00c89 1140 } // else the result is exact
200359e8
L
1141 }
1142
b2a00c89 1143 res = Cstar; // the result is positive
200359e8
L
1144 } else if (exp == 0) {
1145 // 1 <= q <= 10
1146 // res = +C (exact)
b2a00c89
L
1147 res = C1; // the result is positive
1148 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 1149 // res = +C * 10^exp (exact)
b2a00c89 1150 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
1151 }
1152 }
1153 BID_RETURN (res);
1154}
1155
1156/*****************************************************************************
1157 * BID64_to_uint32_xceil
1158 ****************************************************************************/
1159
1160#if DECIMAL_CALL_BY_REFERENCE
1161void
b2a00c89 1162bid64_to_uint32_xceil (unsigned int *pres, UINT64 * px
200359e8
L
1163 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1164 _EXC_INFO_PARAM) {
1165 UINT64 x = *px;
1166#else
1167unsigned int
b2a00c89 1168bid64_to_uint32_xceil (UINT64 x
200359e8
L
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;
b2a00c89 1175 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 1182 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
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
b2a00c89 1195 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
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) {
b2a00c89 1198 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 1199 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 1200 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
1201 x_exp = 0;
1202 C1 = 0;
1203 }
1204 } else {
b2a00c89 1205 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
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
b2a00c89 1219 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 1220 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
1221 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1222 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
1223 x_nr_bits =
1224 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
1225 } else { // x < 2^32
1226 tmp1.d = (double) C1; // exact conversion
200359e8
L
1227 x_nr_bits =
1228 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1229 }
b2a00c89
L
1230 } else { // if x < 2^53
1231 tmp1.d = (double) C1; // exact conversion
200359e8
L
1232 x_nr_bits =
1233 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1234 }
b2a00c89 1235 q = nr_digits[x_nr_bits - 1].digits;
200359e8 1236 if (q == 0) {
b2a00c89
L
1237 q = nr_digits[x_nr_bits - 1].digits1;
1238 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
1239 q++;
1240 }
b2a00c89 1241 exp = x_exp - 398; // unbiased exponent
200359e8 1242
b2a00c89 1243 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
1244 // set invalid flag
1245 *pfpsf |= INVALID_EXCEPTION;
1246 // return Integer Indefinite
1247 res = 0x80000000;
1248 BID_RETURN (res);
b2a00c89 1249 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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'
b2a00c89 1255 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1
200359e8
L
1256 // => set invalid flag
1257 *pfpsf |= INVALID_EXCEPTION;
1258 // return Integer Indefinite
1259 res = 0x80000000;
1260 BID_RETURN (res);
b2a00c89 1261 } else { // if n > 0 and q + exp = 10
200359e8
L
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
b2a00c89 1268 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 1279 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 1284 tmp64 = 0x9fffffff6ull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 1299 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
200359e8
L
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);
b2a00c89 1308 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
1309 // x <= -1 or 1 <= x < 2^32 so if positive, x can be
1310 // rounded to nearest to a 32-bit unsigned integer
b2a00c89 1311 if (x_sign) { // x <= -1
200359e8
L
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
b2a00c89
L
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'
200359e8
L
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
b2a00c89 1327 // shiftright128[] and maskhigh128[]
200359e8 1328 // 1 <= x <= 15
b2a00c89 1329 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
1330 // C* = C1 * 10^(-x)
1331 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 1332 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 1333 Cstar = P128.w[1];
b2a00c89 1334 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
200359e8 1335 fstar.w[0] = P128.w[0];
b2a00c89
L
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
200359e8
L
1338 // C* = floor(C*) (logical right shift; C has p decimal digits,
1339 // correct by Property 1)
1340 // n = C* * 10^(e+x)
1341
b2a00c89
L
1342 // shift right C* by Ex-64 = shiftright128[ind]
1343 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
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
b2a00c89
L
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]
200359e8
L
1354 Cstar++;
1355 // set the inexact flag
1356 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L
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]
200359e8
L
1362 Cstar++;
1363 // set the inexact flag
1364 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89 1365 } // else the result is exact
200359e8
L
1366 }
1367
b2a00c89 1368 res = Cstar; // the result is positive
200359e8
L
1369 } else if (exp == 0) {
1370 // 1 <= q <= 10
1371 // res = +C (exact)
b2a00c89
L
1372 res = C1; // the result is positive
1373 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 1374 // res = +C * 10^exp (exact)
b2a00c89 1375 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
1376 }
1377 }
1378 BID_RETURN (res);
1379}
1380
1381/*****************************************************************************
1382 * BID64_to_uint32_int
1383 ****************************************************************************/
1384
1385#if DECIMAL_CALL_BY_REFERENCE
1386void
b2a00c89 1387bid64_to_uint32_int (unsigned int *pres, UINT64 * px
200359e8
L
1388 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM)
1389{
1390 UINT64 x = *px;
1391#else
1392unsigned int
b2a00c89 1393bid64_to_uint32_int (UINT64 x
200359e8
L
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;
b2a00c89 1400 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 1407 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
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
b2a00c89 1419 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
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) {
b2a00c89 1422 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 1423 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 1424 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
1425 x_exp = 0;
1426 C1 = 0;
1427 }
1428 } else {
b2a00c89 1429 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
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
b2a00c89 1443 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 1444 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
1445 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1446 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
1447 x_nr_bits =
1448 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
1449 } else { // x < 2^32
1450 tmp1.d = (double) C1; // exact conversion
200359e8
L
1451 x_nr_bits =
1452 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1453 }
b2a00c89
L
1454 } else { // if x < 2^53
1455 tmp1.d = (double) C1; // exact conversion
200359e8
L
1456 x_nr_bits =
1457 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1458 }
b2a00c89 1459 q = nr_digits[x_nr_bits - 1].digits;
200359e8 1460 if (q == 0) {
b2a00c89
L
1461 q = nr_digits[x_nr_bits - 1].digits1;
1462 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
1463 q++;
1464 }
b2a00c89 1465 exp = x_exp - 398; // unbiased exponent
200359e8 1466
b2a00c89 1467 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
1468 // set invalid flag
1469 *pfpsf |= INVALID_EXCEPTION;
1470 // return Integer Indefinite
1471 res = 0x80000000;
1472 BID_RETURN (res);
b2a00c89 1473 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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'
b2a00c89 1479 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1
200359e8
L
1480 // => set invalid flag
1481 *pfpsf |= INVALID_EXCEPTION;
1482 // return Integer Indefinite
1483 res = 0x80000000;
1484 BID_RETURN (res);
b2a00c89 1485 } else { // if n > 0 and q + exp = 10
200359e8
L
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
b2a00c89 1492 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 1503 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 1508 tmp64 = 0xa00000000ull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 1523 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
200359e8
L
1524 // return 0
1525 res = 0x00000000;
1526 BID_RETURN (res);
b2a00c89 1527 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
1528 // x <= -1 or 1 <= x < 2^32 so if positive, x can be
1529 // rounded to nearest to a 32-bit unsigned integer
b2a00c89 1530 if (x_sign) { // x <= -1
200359e8
L
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
b2a00c89
L
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'
200359e8
L
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
b2a00c89 1546 // shiftright128[] and maskhigh128[]
200359e8 1547 // 1 <= x <= 15
b2a00c89 1548 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
1549 // C* = C1 * 10^(-x)
1550 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 1551 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 1552 Cstar = P128.w[1];
b2a00c89
L
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
200359e8
L
1555 // C* = floor(C*) (logical right shift; C has p decimal digits,
1556 // correct by Property 1)
1557 // n = C* * 10^(e+x)
1558
b2a00c89
L
1559 // shift right C* by Ex-64 = shiftright128[ind]
1560 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
1561 Cstar = Cstar >> shift;
1562
b2a00c89 1563 res = Cstar; // the result is positive
200359e8
L
1564 } else if (exp == 0) {
1565 // 1 <= q <= 10
1566 // res = +C (exact)
b2a00c89
L
1567 res = C1; // the result is positive
1568 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 1569 // res = +C * 10^exp (exact)
b2a00c89 1570 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
1571 }
1572 }
1573 BID_RETURN (res);
1574}
1575
1576/*****************************************************************************
1577 * BID64_to_uint32_xint
1578 ****************************************************************************/
1579
1580#if DECIMAL_CALL_BY_REFERENCE
1581void
b2a00c89 1582bid64_to_uint32_xint (unsigned int *pres, UINT64 * px
200359e8
L
1583 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1584 _EXC_INFO_PARAM) {
1585 UINT64 x = *px;
1586#else
1587unsigned int
b2a00c89 1588bid64_to_uint32_xint (UINT64 x
200359e8
L
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;
b2a00c89 1595 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 1602 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
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
b2a00c89 1615 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
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) {
b2a00c89 1618 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 1619 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 1620 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
1621 x_exp = 0;
1622 C1 = 0;
1623 }
1624 } else {
b2a00c89 1625 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
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
b2a00c89 1639 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 1640 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
1641 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1642 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
1643 x_nr_bits =
1644 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
1645 } else { // x < 2^32
1646 tmp1.d = (double) C1; // exact conversion
200359e8
L
1647 x_nr_bits =
1648 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1649 }
b2a00c89
L
1650 } else { // if x < 2^53
1651 tmp1.d = (double) C1; // exact conversion
200359e8
L
1652 x_nr_bits =
1653 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1654 }
b2a00c89 1655 q = nr_digits[x_nr_bits - 1].digits;
200359e8 1656 if (q == 0) {
b2a00c89
L
1657 q = nr_digits[x_nr_bits - 1].digits1;
1658 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
1659 q++;
1660 }
b2a00c89 1661 exp = x_exp - 398; // unbiased exponent
200359e8 1662
b2a00c89 1663 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
1664 // set invalid flag
1665 *pfpsf |= INVALID_EXCEPTION;
1666 // return Integer Indefinite
1667 res = 0x80000000;
1668 BID_RETURN (res);
b2a00c89 1669 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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'
b2a00c89 1675 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1
200359e8
L
1676 // => set invalid flag
1677 *pfpsf |= INVALID_EXCEPTION;
1678 // return Integer Indefinite
1679 res = 0x80000000;
1680 BID_RETURN (res);
b2a00c89 1681 } else { // if n > 0 and q + exp = 10
200359e8
L
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
b2a00c89 1688 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 1699 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 1704 tmp64 = 0xa00000000ull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 1719 if ((q + exp) <= 0) { // n = +/-0.[0...0]c(0)c(1)...c(q-1)
200359e8
L
1720 // set inexact flag
1721 *pfpsf |= INEXACT_EXCEPTION;
1722 // return 0
1723 res = 0x00000000;
1724 BID_RETURN (res);
b2a00c89 1725 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
1726 // x <= -1 or 1 <= x < 2^32 so if positive, x can be
1727 // rounded to nearest to a 32-bit unsigned integer
b2a00c89 1728 if (x_sign) { // x <= -1
200359e8
L
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
b2a00c89
L
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'
200359e8
L
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
b2a00c89 1744 // shiftright128[] and maskhigh128[]
200359e8 1745 // 1 <= x <= 15
b2a00c89 1746 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
1747 // C* = C1 * 10^(-x)
1748 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 1749 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 1750 Cstar = P128.w[1];
b2a00c89 1751 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
200359e8 1752 fstar.w[0] = P128.w[0];
b2a00c89
L
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
200359e8
L
1755 // C* = floor(C*) (logical right shift; C has p decimal digits,
1756 // correct by Property 1)
1757 // n = C* * 10^(e+x)
1758
b2a00c89
L
1759 // shift right C* by Ex-64 = shiftright128[ind]
1760 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
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
b2a00c89
L
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]
200359e8
L
1771 // set the inexact flag
1772 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L
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]
200359e8
L
1778 // set the inexact flag
1779 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89 1780 } // else the result is exact
200359e8
L
1781 }
1782
b2a00c89 1783 res = Cstar; // the result is positive
200359e8
L
1784 } else if (exp == 0) {
1785 // 1 <= q <= 10
1786 // res = +C (exact)
b2a00c89
L
1787 res = C1; // the result is positive
1788 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 1789 // res = +C * 10^exp (exact)
b2a00c89 1790 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
1791 }
1792 }
1793 BID_RETURN (res);
1794}
1795
1796/*****************************************************************************
1797 * BID64_to_uint32_rninta
1798 ****************************************************************************/
1799
1800#if DECIMAL_CALL_BY_REFERENCE
1801void
b2a00c89 1802bid64_to_uint32_rninta (unsigned int *pres, UINT64 * px
200359e8
L
1803 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
1804 _EXC_INFO_PARAM) {
1805 UINT64 x = *px;
1806#else
1807unsigned int
b2a00c89 1808bid64_to_uint32_rninta (UINT64 x
200359e8
L
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;
b2a00c89 1815 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 1822 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
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
b2a00c89 1834 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
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) {
b2a00c89 1837 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 1838 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 1839 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
1840 x_exp = 0;
1841 C1 = 0;
1842 }
1843 } else {
b2a00c89 1844 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
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
b2a00c89 1858 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 1859 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
1860 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
1861 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
1862 x_nr_bits =
1863 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
1864 } else { // x < 2^32
1865 tmp1.d = (double) C1; // exact conversion
200359e8
L
1866 x_nr_bits =
1867 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1868 }
b2a00c89
L
1869 } else { // if x < 2^53
1870 tmp1.d = (double) C1; // exact conversion
200359e8
L
1871 x_nr_bits =
1872 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1873 }
b2a00c89 1874 q = nr_digits[x_nr_bits - 1].digits;
200359e8 1875 if (q == 0) {
b2a00c89
L
1876 q = nr_digits[x_nr_bits - 1].digits1;
1877 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
1878 q++;
1879 }
b2a00c89 1880 exp = x_exp - 398; // unbiased exponent
200359e8 1881
b2a00c89 1882 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
1883 // set invalid flag
1884 *pfpsf |= INVALID_EXCEPTION;
1885 // return Integer Indefinite
1886 res = 0x80000000;
1887 BID_RETURN (res);
b2a00c89 1888 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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'
b2a00c89 1894 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1/2
200359e8
L
1895 // => set invalid flag
1896 *pfpsf |= INVALID_EXCEPTION;
1897 // return Integer Indefinite
1898 res = 0x80000000;
1899 BID_RETURN (res);
b2a00c89 1900 } else { // if n > 0 and q + exp = 10
200359e8
L
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
b2a00c89 1907 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 1918 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 1923 tmp64 = 0x9fffffffbull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 1938 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
200359e8
L
1939 // return 0
1940 res = 0x00000000;
1941 BID_RETURN (res);
b2a00c89 1942 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
200359e8
L
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;
b2a00c89
L
1950 if (C1 < midpoint64[ind]) {
1951 res = 0x00000000; // return 0
1952 } else if (x_sign) { // n < 0
200359e8
L
1953 // set invalid flag
1954 *pfpsf |= INVALID_EXCEPTION;
1955 // return Integer Indefinite
1956 res = 0x80000000;
1957 BID_RETURN (res);
b2a00c89
L
1958 } else { // n > 0
1959 res = 0x00000001; // return +1
200359e8 1960 }
b2a00c89 1961 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
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
b2a00c89 1964 if (x_sign) { // x <= -1
200359e8
L
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
b2a00c89
L
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'
200359e8
L
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
b2a00c89 1977 C1 = C1 + midpoint64[ind - 1];
200359e8
L
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
b2a00c89 1981 // shiftright128[] and maskhigh128[]
200359e8 1982 // 1 <= x <= 15
b2a00c89 1983 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
1984 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1985 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 1986 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 1987 Cstar = P128.w[1];
b2a00c89
L
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
200359e8
L
1990 // C* = floor(C*) (logical right shift; C has p decimal digits,
1991 // correct by Property 1)
1992 // n = C* * 10^(e+x)
1993
b2a00c89
L
1994 // shift right C* by Ex-64 = shiftright128[ind]
1995 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
1996 Cstar = Cstar >> shift;
1997
1998 // if the result was a midpoint it was rounded away from zero
b2a00c89 1999 res = Cstar; // the result is positive
200359e8
L
2000 } else if (exp == 0) {
2001 // 1 <= q <= 10
2002 // res = +C (exact)
b2a00c89
L
2003 res = C1; // the result is positive
2004 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 2005 // res = +C * 10^exp (exact)
b2a00c89 2006 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
2007 }
2008 }
2009 BID_RETURN (res);
2010}
2011
2012/*****************************************************************************
2013 * BID64_to_uint32_xrninta
2014 ****************************************************************************/
2015
2016#if DECIMAL_CALL_BY_REFERENCE
2017void
b2a00c89 2018bid64_to_uint32_xrninta (unsigned int *pres, UINT64 * px
200359e8
L
2019 _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
2020 _EXC_INFO_PARAM) {
2021 UINT64 x = *px;
2022#else
2023unsigned int
b2a00c89 2024bid64_to_uint32_xrninta (UINT64 x
200359e8
L
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;
b2a00c89 2031 int exp; // unbiased exponent
200359e8
L
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;
b2a00c89 2038 UINT64 Cstar; // C* represents up to 16 decimal digits ~ 54 bits
200359e8
L
2039 UINT128 fstar;
2040 UINT128 P128;
2041
2042 // check for NaN or Infinity
2043 if ((x & MASK_NAN) == MASK_NAN || (x & MASK_INF) == MASK_INF) {
2044 // set invalid flag
2045 *pfpsf |= INVALID_EXCEPTION;
2046 // return Integer Indefinite
2047 res = 0x80000000;
2048 BID_RETURN (res);
2049 }
2050 // unpack x
b2a00c89 2051 x_sign = x & MASK_SIGN; // 0 for positive, MASK_SIGN for negative
200359e8
L
2052 // if steering bits are 11 (condition will be 0), then exponent is G[0:w+1] =>
2053 if ((x & MASK_STEERING_BITS) == MASK_STEERING_BITS) {
b2a00c89 2054 x_exp = (x & MASK_BINARY_EXPONENT2) >> 51; // biased
200359e8 2055 C1 = (x & MASK_BINARY_SIG2) | MASK_BINARY_OR2;
b2a00c89 2056 if (C1 > 9999999999999999ull) { // non-canonical
200359e8
L
2057 x_exp = 0;
2058 C1 = 0;
2059 }
2060 } else {
b2a00c89 2061 x_exp = (x & MASK_BINARY_EXPONENT1) >> 53; // biased
200359e8
L
2062 C1 = x & MASK_BINARY_SIG1;
2063 }
2064
2065 // check for zeros (possibly from non-canonical values)
2066 if (C1 == 0x0ull) {
2067 // x is 0
2068 res = 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
b2a00c89 2075 if (C1 >= 0x0020000000000000ull) { // x >= 2^53
200359e8 2076 // split the 64-bit value in two 32-bit halves to avoid rounding errors
b2a00c89
L
2077 if (C1 >= 0x0000000100000000ull) { // x >= 2^32
2078 tmp1.d = (double) (C1 >> 32); // exact conversion
200359e8
L
2079 x_nr_bits =
2080 33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
b2a00c89
L
2081 } else { // x < 2^32
2082 tmp1.d = (double) C1; // exact conversion
200359e8
L
2083 x_nr_bits =
2084 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2085 }
b2a00c89
L
2086 } else { // if x < 2^53
2087 tmp1.d = (double) C1; // exact conversion
200359e8
L
2088 x_nr_bits =
2089 1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2090 }
b2a00c89 2091 q = nr_digits[x_nr_bits - 1].digits;
200359e8 2092 if (q == 0) {
b2a00c89
L
2093 q = nr_digits[x_nr_bits - 1].digits1;
2094 if (C1 >= nr_digits[x_nr_bits - 1].threshold_lo)
200359e8
L
2095 q++;
2096 }
b2a00c89 2097 exp = x_exp - 398; // unbiased exponent
200359e8 2098
b2a00c89 2099 if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
200359e8
L
2100 // set invalid flag
2101 *pfpsf |= INVALID_EXCEPTION;
2102 // return Integer Indefinite
2103 res = 0x80000000;
2104 BID_RETURN (res);
b2a00c89 2105 } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
200359e8
L
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'
b2a00c89 2111 if (x_sign) { // if n < 0 and q + exp = 10 then x is much less than -1/2
200359e8
L
2112 // => set invalid flag
2113 *pfpsf |= INVALID_EXCEPTION;
2114 // return Integer Indefinite
2115 res = 0x80000000;
2116 BID_RETURN (res);
b2a00c89 2117 } else { // if n > 0 and q + exp = 10
200359e8
L
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
b2a00c89 2124 tmp64 = C1 * ten2k64[11 - q]; // C scaled up to 11-digit int
200359e8
L
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'
b2a00c89 2135 } else { // if (q > 11), i.e. 12 <= q <= 16 and so -15 <= exp <= -2
200359e8
L
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
b2a00c89 2140 tmp64 = 0x9fffffffbull * ten2k64[q - 11];
200359e8
L
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
b2a00c89 2155 if ((q + exp) < 0) { // n = +/-0.0...c(0)c(1)...c(q-1)
200359e8
L
2156 // set inexact flag
2157 *pfpsf |= INEXACT_EXCEPTION;
2158 // return 0
2159 res = 0x00000000;
2160 BID_RETURN (res);
b2a00c89 2161 } else if ((q + exp) == 0) { // n = +/-0.c(0)c(1)...c(q-1)
200359e8
L
2162 // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
2163 // res = 0
2164 // else if x > 0
2165 // res = +1
2166 // else // if x < 0
2167 // invalid exc
2168 ind = q - 1;
b2a00c89
L
2169 if (C1 < midpoint64[ind]) {
2170 res = 0x00000000; // return 0
2171 } else if (x_sign) { // n < 0
200359e8
L
2172 // set invalid flag
2173 *pfpsf |= INVALID_EXCEPTION;
2174 // return Integer Indefinite
2175 res = 0x80000000;
2176 BID_RETURN (res);
b2a00c89
L
2177 } else { // n > 0
2178 res = 0x00000001; // return +1
200359e8
L
2179 }
2180 // set inexact flag
2181 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89 2182 } else { // if (1 <= q + exp <= 10, 1 <= q <= 16, -15 <= exp <= 9)
200359e8
L
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
b2a00c89 2185 if (x_sign) { // x <= -1
200359e8
L
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
b2a00c89
L
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'
200359e8
L
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
b2a00c89 2198 C1 = C1 + midpoint64[ind - 1];
200359e8
L
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
b2a00c89 2202 // shiftright128[] and maskhigh128[]
200359e8 2203 // 1 <= x <= 15
b2a00c89 2204 // kx = 10^(-x) = ten2mk64[ind - 1]
200359e8
L
2205 // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2206 // the approximation of 10^(-x) was rounded up to 54 bits
b2a00c89 2207 __mul_64x64_to_128MACH (P128, C1, ten2mk64[ind - 1]);
200359e8 2208 Cstar = P128.w[1];
b2a00c89 2209 fstar.w[1] = P128.w[1] & maskhigh128[ind - 1];
200359e8 2210 fstar.w[0] = P128.w[0];
b2a00c89
L
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
200359e8
L
2213 // C* = floor(C*) (logical right shift; C has p decimal digits,
2214 // correct by Property 1)
2215 // n = C* * 10^(e+x)
2216
b2a00c89
L
2217 // shift right C* by Ex-64 = shiftright128[ind]
2218 shift = shiftright128[ind - 1]; // 0 <= shift <= 39
200359e8
L
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
b2a00c89 2226 if (ind - 1 <= 2) { // fstar.w[1] is 0
200359e8
L
2227 if (fstar.w[0] > 0x8000000000000000ull) {
2228 // f* > 1/2 and the result may be exact
b2a00c89
L
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]
200359e8
L
2233 // set the inexact flag
2234 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L
2235 } // else the result is exact
2236 } else { // the result is inexact; f2* <= 1/2
200359e8
L
2237 // set the inexact flag
2238 *pfpsf |= INEXACT_EXCEPTION;
2239 }
b2a00c89
L
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])) {
200359e8
L
2243 // f2* > 1/2 and the result may be exact
2244 // Calculate f2* - 1/2
b2a00c89
L
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]
200359e8
L
2249 // set the inexact flag
2250 *pfpsf |= INEXACT_EXCEPTION;
b2a00c89
L
2251 } // else the result is exact
2252 } else { // the result is inexact; f2* <= 1/2
200359e8
L
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
b2a00c89 2259 res = Cstar; // the result is positive
200359e8
L
2260 } else if (exp == 0) {
2261 // 1 <= q <= 10
2262 // res = +C (exact)
b2a00c89
L
2263 res = C1; // the result is positive
2264 } else { // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
200359e8 2265 // res = +C * 10^exp (exact)
b2a00c89 2266 res = C1 * ten2k64[exp]; // the result is positive
200359e8
L
2267 }
2268 }
2269 BID_RETURN (res);
2270}