OpenCores
URL https://opencores.org/ocsvn/openrisc/openrisc/trunk

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgcc/] [config/] [libbid/] [bid128_to_uint32.c] - Blame information for rev 741

Go to most recent revision | Details | Compare with Previous | View Log

Line No. Rev Author Line
1 734 jeremybenn 
/* 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 
 *  BID128_to_uint32_rnint
28 
 ****************************************************************************/
29 
 
30 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
31 
                                          bid128_to_uint32_rnint, x)
32 
 
33 
     unsigned int res;
34 
     UINT64 x_sign;
35 
     UINT64 x_exp;
36 
     int exp;                   // unbiased exponent
37 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
38 
     UINT64 tmp64;
39 
     BID_UI64DOUBLE tmp1;
40 
     unsigned int x_nr_bits;
41 
     int q, ind, shift;
42 
     UINT128 C1, C;
43 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
44 
     UINT256 fstar;
45 
     UINT256 P256;
46 
 
47 
  // unpack x
48 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
49 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
50 
C1.w[1] = x.w[1] & MASK_COEFF;
51 
C1.w[0] = x.w[0];
52 
 
53 
  // check for NaN or Infinity
54 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
55 
    // x is special
56 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
57 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
58 
    // set invalid flag
59 
    *pfpsf |= INVALID_EXCEPTION;
60 
    // return Integer Indefinite
61 
    res = 0x80000000;
62 
  } else {      // x is QNaN
63 
    // set invalid flag
64 
    *pfpsf |= INVALID_EXCEPTION;
65 
    // return Integer Indefinite
66 
    res = 0x80000000;
67 
  }
68 
  BID_RETURN (res);
69 
} else {        // x is not a NaN, so it must be infinity
70 
  if (!x_sign) {        // x is +inf
71 
    // set invalid flag
72 
    *pfpsf |= INVALID_EXCEPTION;
73 
    // return Integer Indefinite
74 
    res = 0x80000000;
75 
  } else {      // x is -inf
76 
    // set invalid flag
77 
    *pfpsf |= INVALID_EXCEPTION;
78 
    // return Integer Indefinite
79 
    res = 0x80000000;
80 
  }
81 
  BID_RETURN (res);
82 
}
83 
}
84 
  // check for non-canonical values (after the check for special values)
85 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
86 
    || (C1.w[1] == 0x0001ed09bead87c0ull
87 
        && (C1.w[0] > 0x378d8e63ffffffffull))
88 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
89 
  res = 0x00000000;
90 
  BID_RETURN (res);
91 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
92 
  // x is 0
93 
  res = 0x00000000;
94 
  BID_RETURN (res);
95 
} else {        // x is not special and is not zero
96 
 
97 
  // q = nr. of decimal digits in x
98 
  //  determine first the nr. of bits in x
99 
  if (C1.w[1] == 0) {
100 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
101 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
102 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
103 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
104 
        x_nr_bits =
105 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
106 
      } else {  // x < 2^32
107 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
108 
        x_nr_bits =
109 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
110 
      }
111 
    } else {    // if x < 2^53
112 
      tmp1.d = (double) C1.w[0]; // exact conversion
113 
      x_nr_bits =
114 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
115 
    }
116 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
117 
    tmp1.d = (double) C1.w[1];  // exact conversion
118 
    x_nr_bits =
119 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
120 
  }
121 
  q = nr_digits[x_nr_bits - 1].digits;
122 
  if (q == 0) {
123 
    q = nr_digits[x_nr_bits - 1].digits1;
124 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
125 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
126 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
127 
      q++;
128 
  }
129 
  exp = (x_exp >> 49) - 6176;
130 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
131 
    // set invalid flag
132 
    *pfpsf |= INVALID_EXCEPTION;
133 
    // return Integer Indefinite
134 
    res = 0x80000000;
135 
    BID_RETURN (res);
136 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
137 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
138 
    // so x rounded to an integer may or may not fit in an unsigned 32-bit int
139 
    // the cases that do not fit are identified here; the ones that fit
140 
    // fall through and will be handled with other cases further,
141 
    // under '1 <= q + exp <= 10'
142 
    if (x_sign) {       // if n < 0 and q + exp = 10
143 
      // if n < -1/2 then n cannot be converted to uint32 with RN
144 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 1/2
145 
      // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x05, 1<=q<=34
146 
      if (q <= 11) {
147 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
148 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
149 
        if (tmp64 > 0x05ull) {
150 
          // set invalid flag
151 
          *pfpsf |= INVALID_EXCEPTION;
152 
          // return Integer Indefinite
153 
          res = 0x80000000;
154 
          BID_RETURN (res);
155 
        }
156 
        // else cases that can be rounded to a 32-bit int fall through
157 
        // to '1 <= q + exp <= 10'
158 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
159 
        // 0.c(0)c(1)...c(q-1) * 10^11 > 0x05 <=>
160 
        // C > 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
161 
        // (scale 1/2 up)
162 
        tmp64 = 0x05ull;
163 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
164 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
165 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
166 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
167 
        }
168 
        if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
169 
          // set invalid flag
170 
          *pfpsf |= INVALID_EXCEPTION;
171 
          // return Integer Indefinite
172 
          res = 0x80000000;
173 
          BID_RETURN (res);
174 
        }
175 
        // else cases that can be rounded to a 32-bit int fall through
176 
        // to '1 <= q + exp <= 10'
177 
      }
178 
    } else {    // if n > 0 and q + exp = 10
179 
      // if n >= 2^32 - 1/2 then n is too large
180 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
181 
      // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
182 
      if (q <= 11) {
183 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
184 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
185 
        if (tmp64 >= 0x9fffffffbull) {
186 
          // set invalid flag
187 
          *pfpsf |= INVALID_EXCEPTION;
188 
          // return Integer Indefinite
189 
          res = 0x80000000;
190 
          BID_RETURN (res);
191 
        }
192 
        // else cases that can be rounded to a 32-bit int fall through
193 
        // to '1 <= q + exp <= 10'
194 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
195 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
196 
        // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
197 
        // (scale 2^32-1/2 up)
198 
        tmp64 = 0x9fffffffbull;
199 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
200 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
201 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
202 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
203 
        }
204 
        if (C1.w[1] > C.w[1]
205 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
206 
          // set invalid flag
207 
          *pfpsf |= INVALID_EXCEPTION;
208 
          // return Integer Indefinite
209 
          res = 0x80000000;
210 
          BID_RETURN (res);
211 
        }
212 
        // else cases that can be rounded to a 32-bit int fall through
213 
        // to '1 <= q + exp <= 10'
214 
      }
215 
    }
216 
  }
217 
  // n is not too large to be converted to int32: -1/2 <= n < 2^32 - 1/2
218 
  // Note: some of the cases tested for above fall through to this point
219 
  if ((q + exp) < 0) {   // n = +/-0.0...c(0)c(1)...c(q-1)
220 
    // return 0
221 
    res = 0x00000000;
222 
    BID_RETURN (res);
223 
  } else if ((q + exp) == 0) {   // n = +/-0.c(0)c(1)...c(q-1)
224 
    // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
225 
    //   res = 0
226 
    // else if x > 0
227 
    //   res = +1
228 
    // else // if x < 0
229 
    //   invalid exc
230 
    ind = q - 1;
231 
    if (ind <= 18) {    // 0 <= ind <= 18
232 
      if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
233 
        res = 0x00000000;       // return 0
234 
      } else if (!x_sign) {     // n > 0
235 
        res = 0x00000001;       // return +1
236 
      } else {
237 
        res = 0x80000000;
238 
        *pfpsf |= INVALID_EXCEPTION;
239 
      }
240 
    } else {    // 19 <= ind <= 33
241 
      if ((C1.w[1] < midpoint128[ind - 19].w[1])
242 
          || ((C1.w[1] == midpoint128[ind - 19].w[1])
243 
              && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
244 
        res = 0x00000000;       // return 0
245 
      } else if (!x_sign) {     // n > 0
246 
        res = 0x00000001;       // return +1
247 
      } else {
248 
        res = 0x80000000;
249 
        *pfpsf |= INVALID_EXCEPTION;
250 
      }
251 
    }
252 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
253 
    if (x_sign) {       // x <= -1
254 
      // set invalid flag
255 
      *pfpsf |= INVALID_EXCEPTION;
256 
      // return Integer Indefinite
257 
      res = 0x80000000;
258 
      BID_RETURN (res);
259 
    }
260 
    // 1 <= x < 2^32-1/2 so x can be rounded
261 
    // to nearest to a 32-bit unsigned integer
262 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
263 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
264 
      // chop off ind digits from the lower part of C1
265 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
266 
      tmp64 = C1.w[0];
267 
      if (ind <= 19) {
268 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
269 
      } else {
270 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
271 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
272 
      }
273 
      if (C1.w[0] < tmp64)
274 
        C1.w[1]++;
275 
      // calculate C* and f*
276 
      // C* is actually floor(C*) in this case
277 
      // C* and f* need shifting and masking, as shown by
278 
      // shiftright128[] and maskhigh128[]
279 
      // 1 <= x <= 33
280 
      // kx = 10^(-x) = ten2mk128[ind - 1]
281 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
282 
      // the approximation of 10^(-x) was rounded up to 118 bits
283 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
284 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
285 
        Cstar.w[1] = P256.w[3];
286 
        Cstar.w[0] = P256.w[2];
287 
        fstar.w[3] = 0;
288 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
289 
        fstar.w[1] = P256.w[1];
290 
        fstar.w[0] = P256.w[0];
291 
      } else {  // 22 <= ind - 1 <= 33
292 
        Cstar.w[1] = 0;
293 
        Cstar.w[0] = P256.w[3];
294 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
295 
        fstar.w[2] = P256.w[2];
296 
        fstar.w[1] = P256.w[1];
297 
        fstar.w[0] = P256.w[0];
298 
      }
299 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
300 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
301 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
302 
      //   if floor(C*) is even then C* = floor(C*) - logical right
303 
      //       shift; C* has p decimal digits, correct by Prop. 1)
304 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
305 
      //       shift; C* has p decimal digits, correct by Pr. 1)
306 
      // else
307 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
308 
      //       correct by Property 1)
309 
      // n = C* * 10^(e+x)
310 
 
311 
      // shift right C* by Ex-128 = shiftright128[ind]
312 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
313 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
314 
        Cstar.w[0] =
315 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
316 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
317 
      } else {  // 22 <= ind - 1 <= 33
318 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
319 
      }
320 
      // if the result was a midpoint it was rounded away from zero, so
321 
      // it will need a correction
322 
      // check for midpoints
323 
      if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
324 
          && (fstar.w[1] || fstar.w[0])
325 
          && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
326 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
327 
                  && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
328 
        // the result is a midpoint; round to nearest
329 
        if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
330 
          // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
331 
          Cstar.w[0]--;  // Cstar.w[0] is now even
332 
        }       // else MP in [ODD, EVEN]
333 
      }
334 
      res = Cstar.w[0];  // the result is positive
335 
    } else if (exp == 0) {
336 
      // 1 <= q <= 10
337 
      // res = C (exact)
338 
      res = C1.w[0];
339 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
340 
      // res = C * 10^exp (exact)
341 
      res = C1.w[0] * ten2k64[exp];
342 
    }
343 
  }
344 
}
345 
 
346 
BID_RETURN (res);
347 
}
348 
 
349 
/*****************************************************************************
350 
 *  BID128_to_uint32_xrnint
351 
 ****************************************************************************/
352 
 
353 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
354 
                                          bid128_to_uint32_xrnint, x)
355 
 
356 
     unsigned int res;
357 
     UINT64 x_sign;
358 
     UINT64 x_exp;
359 
     int exp;                   // unbiased exponent
360 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
361 
     UINT64 tmp64, tmp64A;
362 
     BID_UI64DOUBLE tmp1;
363 
     unsigned int x_nr_bits;
364 
     int q, ind, shift;
365 
     UINT128 C1, C;
366 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
367 
     UINT256 fstar;
368 
     UINT256 P256;
369 
     unsigned int tmp_inexact = 0;
370 
 
371 
  // unpack x
372 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
373 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
374 
C1.w[1] = x.w[1] & MASK_COEFF;
375 
C1.w[0] = x.w[0];
376 
 
377 
  // check for NaN or Infinity
378 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
379 
    // x is special
380 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
381 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
382 
    // set invalid flag
383 
    *pfpsf |= INVALID_EXCEPTION;
384 
    // return Integer Indefinite
385 
    res = 0x80000000;
386 
  } else {      // x is QNaN
387 
    // set invalid flag
388 
    *pfpsf |= INVALID_EXCEPTION;
389 
    // return Integer Indefinite
390 
    res = 0x80000000;
391 
  }
392 
  BID_RETURN (res);
393 
} else {        // x is not a NaN, so it must be infinity
394 
  if (!x_sign) {        // x is +inf
395 
    // set invalid flag
396 
    *pfpsf |= INVALID_EXCEPTION;
397 
    // return Integer Indefinite
398 
    res = 0x80000000;
399 
  } else {      // x is -inf
400 
    // set invalid flag
401 
    *pfpsf |= INVALID_EXCEPTION;
402 
    // return Integer Indefinite
403 
    res = 0x80000000;
404 
  }
405 
  BID_RETURN (res);
406 
}
407 
}
408 
  // check for non-canonical values (after the check for special values)
409 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
410 
    || (C1.w[1] == 0x0001ed09bead87c0ull
411 
        && (C1.w[0] > 0x378d8e63ffffffffull))
412 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
413 
  res = 0x00000000;
414 
  BID_RETURN (res);
415 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
416 
  // x is 0
417 
  res = 0x00000000;
418 
  BID_RETURN (res);
419 
} else {        // x is not special and is not zero
420 
 
421 
  // q = nr. of decimal digits in x
422 
  //  determine first the nr. of bits in x
423 
  if (C1.w[1] == 0) {
424 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
425 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
426 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
427 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
428 
        x_nr_bits =
429 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
430 
      } else {  // x < 2^32
431 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
432 
        x_nr_bits =
433 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
434 
      }
435 
    } else {    // if x < 2^53
436 
      tmp1.d = (double) C1.w[0]; // exact conversion
437 
      x_nr_bits =
438 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
439 
    }
440 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
441 
    tmp1.d = (double) C1.w[1];  // exact conversion
442 
    x_nr_bits =
443 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
444 
  }
445 
  q = nr_digits[x_nr_bits - 1].digits;
446 
  if (q == 0) {
447 
    q = nr_digits[x_nr_bits - 1].digits1;
448 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
449 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
450 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
451 
      q++;
452 
  }
453 
  exp = (x_exp >> 49) - 6176;
454 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
455 
    // set invalid flag
456 
    *pfpsf |= INVALID_EXCEPTION;
457 
    // return Integer Indefinite
458 
    res = 0x80000000;
459 
    BID_RETURN (res);
460 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
461 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
462 
    // so x rounded to an integer may or may not fit in an unsigned 32-bit int
463 
    // the cases that do not fit are identified here; the ones that fit
464 
    // fall through and will be handled with other cases further,
465 
    // under '1 <= q + exp <= 10'
466 
    if (x_sign) {       // if n < 0 and q + exp = 10
467 
      // if n < -1/2 then n cannot be converted to uint32 with RN
468 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 1/2
469 
      // <=> 0.c(0)c(1)...c(q-1) * 10^11 > 0x05, 1<=q<=34
470 
      if (q <= 11) {
471 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
472 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
473 
        if (tmp64 > 0x05ull) {
474 
          // set invalid flag
475 
          *pfpsf |= INVALID_EXCEPTION;
476 
          // return Integer Indefinite
477 
          res = 0x80000000;
478 
          BID_RETURN (res);
479 
        }
480 
        // else cases that can be rounded to a 32-bit int fall through
481 
        // to '1 <= q + exp <= 10'
482 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
483 
        // 0.c(0)c(1)...c(q-1) * 10^11 > 0x05 <=>
484 
        // C > 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
485 
        // (scale 1/2 up)
486 
        tmp64 = 0x05ull;
487 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
488 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
489 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
490 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
491 
        }
492 
        if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
493 
          // set invalid flag
494 
          *pfpsf |= INVALID_EXCEPTION;
495 
          // return Integer Indefinite
496 
          res = 0x80000000;
497 
          BID_RETURN (res);
498 
        }
499 
        // else cases that can be rounded to a 32-bit int fall through
500 
        // to '1 <= q + exp <= 10'
501 
      }
502 
    } else {    // if n > 0 and q + exp = 10
503 
      // if n >= 2^32 - 1/2 then n is too large
504 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
505 
      // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
506 
      if (q <= 11) {
507 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
508 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
509 
        if (tmp64 >= 0x9fffffffbull) {
510 
          // set invalid flag
511 
          *pfpsf |= INVALID_EXCEPTION;
512 
          // return Integer Indefinite
513 
          res = 0x80000000;
514 
          BID_RETURN (res);
515 
        }
516 
        // else cases that can be rounded to a 32-bit int fall through
517 
        // to '1 <= q + exp <= 10'
518 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
519 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
520 
        // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
521 
        // (scale 2^32-1/2 up)
522 
        tmp64 = 0x9fffffffbull;
523 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
524 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
525 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
526 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
527 
        }
528 
        if (C1.w[1] > C.w[1]
529 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
530 
          // set invalid flag
531 
          *pfpsf |= INVALID_EXCEPTION;
532 
          // return Integer Indefinite
533 
          res = 0x80000000;
534 
          BID_RETURN (res);
535 
        }
536 
        // else cases that can be rounded to a 32-bit int fall through
537 
        // to '1 <= q + exp <= 10'
538 
      }
539 
    }
540 
  }
541 
  // n is not too large to be converted to int32: -1/2 <= n < 2^32 - 1/2
542 
  // Note: some of the cases tested for above fall through to this point
543 
  if ((q + exp) < 0) {   // n = +/-0.0...c(0)c(1)...c(q-1)
544 
    // set inexact flag
545 
    *pfpsf |= INEXACT_EXCEPTION;
546 
    // return 0
547 
    res = 0x00000000;
548 
    BID_RETURN (res);
549 
  } else if ((q + exp) == 0) {   // n = +/-0.c(0)c(1)...c(q-1)
550 
    // if 0.c(0)c(1)...c(q-1) <= 0.5 <=> c(0)c(1)...c(q-1) <= 5 * 10^(q-1)
551 
    //   res = 0
552 
    // else if x > 0
553 
    //   res = +1
554 
    // else // if x < 0
555 
    //   invalid exc
556 
    ind = q - 1;
557 
    if (ind <= 18) {    // 0 <= ind <= 18
558 
      if ((C1.w[1] == 0) && (C1.w[0] <= midpoint64[ind])) {
559 
        res = 0x00000000;       // return 0
560 
      } else if (!x_sign) {     // n > 0
561 
        res = 0x00000001;       // return +1
562 
      } else {
563 
        res = 0x80000000;
564 
        *pfpsf |= INVALID_EXCEPTION;
565 
        BID_RETURN (res);
566 
      }
567 
    } else {    // 19 <= ind <= 33
568 
      if ((C1.w[1] < midpoint128[ind - 19].w[1])
569 
          || ((C1.w[1] == midpoint128[ind - 19].w[1])
570 
              && (C1.w[0] <= midpoint128[ind - 19].w[0]))) {
571 
        res = 0x00000000;       // return 0
572 
      } else if (!x_sign) {     // n > 0
573 
        res = 0x00000001;       // return +1
574 
      } else {
575 
        res = 0x80000000;
576 
        *pfpsf |= INVALID_EXCEPTION;
577 
        BID_RETURN (res);
578 
      }
579 
    }
580 
    // set inexact flag
581 
    *pfpsf |= INEXACT_EXCEPTION;
582 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
583 
    if (x_sign) {       // x <= -1
584 
      // set invalid flag
585 
      *pfpsf |= INVALID_EXCEPTION;
586 
      // return Integer Indefinite
587 
      res = 0x80000000;
588 
      BID_RETURN (res);
589 
    }
590 
    // 1 <= x < 2^32-1/2 so x can be rounded
591 
    // to nearest to a 32-bit unsigned integer
592 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
593 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
594 
      // chop off ind digits from the lower part of C1
595 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
596 
      tmp64 = C1.w[0];
597 
      if (ind <= 19) {
598 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
599 
      } else {
600 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
601 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
602 
      }
603 
      if (C1.w[0] < tmp64)
604 
        C1.w[1]++;
605 
      // calculate C* and f*
606 
      // C* is actually floor(C*) in this case
607 
      // C* and f* need shifting and masking, as shown by
608 
      // shiftright128[] and maskhigh128[]
609 
      // 1 <= x <= 33
610 
      // kx = 10^(-x) = ten2mk128[ind - 1]
611 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
612 
      // the approximation of 10^(-x) was rounded up to 118 bits
613 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
614 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
615 
        Cstar.w[1] = P256.w[3];
616 
        Cstar.w[0] = P256.w[2];
617 
        fstar.w[3] = 0;
618 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
619 
        fstar.w[1] = P256.w[1];
620 
        fstar.w[0] = P256.w[0];
621 
      } else {  // 22 <= ind - 1 <= 33
622 
        Cstar.w[1] = 0;
623 
        Cstar.w[0] = P256.w[3];
624 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
625 
        fstar.w[2] = P256.w[2];
626 
        fstar.w[1] = P256.w[1];
627 
        fstar.w[0] = P256.w[0];
628 
      }
629 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
630 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
631 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
632 
      //   if floor(C*) is even then C* = floor(C*) - logical right
633 
      //       shift; C* has p decimal digits, correct by Prop. 1)
634 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
635 
      //       shift; C* has p decimal digits, correct by Pr. 1)
636 
      // else
637 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
638 
      //       correct by Property 1)
639 
      // n = C* * 10^(e+x)
640 
 
641 
      // shift right C* by Ex-128 = shiftright128[ind]
642 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
643 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
644 
        Cstar.w[0] =
645 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
646 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
647 
      } else {  // 22 <= ind - 1 <= 33
648 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
649 
      }
650 
      // determine inexactness of the rounding of C*
651 
      // if (0 < f* - 1/2 < 10^(-x)) then
652 
      //   the result is exact
653 
      // else // if (f* - 1/2 > T*) then
654 
      //   the result is inexact
655 
      if (ind - 1 <= 2) {
656 
        if (fstar.w[1] > 0x8000000000000000ull ||
657 
            (fstar.w[1] == 0x8000000000000000ull
658 
             && fstar.w[0] > 0x0ull)) {
659 
          // f* > 1/2 and the result may be exact
660 
          tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
661 
          if (tmp64 > ten2mk128trunc[ind - 1].w[1]
662 
              || (tmp64 == ten2mk128trunc[ind - 1].w[1]
663 
                  && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
664 
            // set the inexact flag
665 
            // *pfpsf |= INEXACT_EXCEPTION;
666 
            tmp_inexact = 1;
667 
          }     // else the result is exact
668 
        } else {        // the result is inexact; f2* <= 1/2
669 
          // set the inexact flag
670 
          // *pfpsf |= INEXACT_EXCEPTION;
671 
          tmp_inexact = 1;
672 
        }
673 
      } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
674 
        if (fstar.w[3] > 0x0 ||
675 
            (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
676 
            (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
677 
             (fstar.w[1] || fstar.w[0]))) {
678 
          // f2* > 1/2 and the result may be exact
679 
          // Calculate f2* - 1/2
680 
          tmp64 = fstar.w[2] - onehalf128[ind - 1];
681 
          tmp64A = fstar.w[3];
682 
          if (tmp64 > fstar.w[2])
683 
            tmp64A--;
684 
          if (tmp64A || tmp64
685 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
686 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
687 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
688 
            // set the inexact flag
689 
            // *pfpsf |= INEXACT_EXCEPTION;
690 
            tmp_inexact = 1;
691 
          }     // else the result is exact
692 
        } else {        // the result is inexact; f2* <= 1/2
693 
          // set the inexact flag
694 
          // *pfpsf |= INEXACT_EXCEPTION;
695 
          tmp_inexact = 1;
696 
        }
697 
      } else {  // if 22 <= ind <= 33
698 
        if (fstar.w[3] > onehalf128[ind - 1] ||
699 
            (fstar.w[3] == onehalf128[ind - 1] &&
700 
             (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
701 
          // f2* > 1/2 and the result may be exact
702 
          // Calculate f2* - 1/2
703 
          tmp64 = fstar.w[3] - onehalf128[ind - 1];
704 
          if (tmp64 || fstar.w[2]
705 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
706 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
707 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
708 
            // set the inexact flag
709 
            // *pfpsf |= INEXACT_EXCEPTION;
710 
            tmp_inexact = 1;
711 
          }     // else the result is exact
712 
        } else {        // the result is inexact; f2* <= 1/2
713 
          // set the inexact flag
714 
          // *pfpsf |= INEXACT_EXCEPTION;
715 
          tmp_inexact = 1;
716 
        }
717 
      }
718 
 
719 
      // if the result was a midpoint it was rounded away from zero, so
720 
      // it will need a correction
721 
      // check for midpoints
722 
      if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
723 
          && (fstar.w[1] || fstar.w[0])
724 
          && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
725 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
726 
                  && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
727 
        // the result is a midpoint; round to nearest
728 
        if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
729 
          // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
730 
          Cstar.w[0]--;  // Cstar.w[0] is now even
731 
        }       // else MP in [ODD, EVEN]
732 
      }
733 
      res = Cstar.w[0];  // the result is positive
734 
      if (tmp_inexact)
735 
        *pfpsf |= INEXACT_EXCEPTION;
736 
    } else if (exp == 0) {
737 
      // 1 <= q <= 10
738 
      // res = C (exact)
739 
      res = C1.w[0];
740 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
741 
      // res = C * 10^exp (exact)
742 
      res = C1.w[0] * ten2k64[exp];
743 
    }
744 
  }
745 
}
746 
 
747 
BID_RETURN (res);
748 
}
749 
 
750 
/*****************************************************************************
751 
 *  BID128_to_uint32_floor
752 
 ****************************************************************************/
753 
 
754 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
755 
                                          bid128_to_uint32_floor, x)
756 
 
757 
     unsigned int res;
758 
     UINT64 x_sign;
759 
     UINT64 x_exp;
760 
     int exp;                   // unbiased exponent
761 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
762 
     UINT64 tmp64, tmp64A;
763 
     BID_UI64DOUBLE tmp1;
764 
     unsigned int x_nr_bits;
765 
     int q, ind, shift;
766 
     UINT128 C1, C;
767 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
768 
     UINT256 fstar;
769 
     UINT256 P256;
770 
     int is_inexact_gt_midpoint = 0;
771 
     int is_midpoint_lt_even = 0;
772 
 
773 
  // unpack x
774 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
775 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
776 
C1.w[1] = x.w[1] & MASK_COEFF;
777 
C1.w[0] = x.w[0];
778 
 
779 
  // check for NaN or Infinity
780 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
781 
    // x is special
782 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
783 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
784 
    // set invalid flag
785 
    *pfpsf |= INVALID_EXCEPTION;
786 
    // return Integer Indefinite
787 
    res = 0x80000000;
788 
  } else {      // x is QNaN
789 
    // set invalid flag
790 
    *pfpsf |= INVALID_EXCEPTION;
791 
    // return Integer Indefinite
792 
    res = 0x80000000;
793 
  }
794 
  BID_RETURN (res);
795 
} else {        // x is not a NaN, so it must be infinity
796 
  if (!x_sign) {        // x is +inf
797 
    // set invalid flag
798 
    *pfpsf |= INVALID_EXCEPTION;
799 
    // return Integer Indefinite
800 
    res = 0x80000000;
801 
  } else {      // x is -inf
802 
    // set invalid flag
803 
    *pfpsf |= INVALID_EXCEPTION;
804 
    // return Integer Indefinite
805 
    res = 0x80000000;
806 
  }
807 
  BID_RETURN (res);
808 
}
809 
}
810 
  // check for non-canonical values (after the check for special values)
811 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
812 
    || (C1.w[1] == 0x0001ed09bead87c0ull
813 
        && (C1.w[0] > 0x378d8e63ffffffffull))
814 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
815 
  res = 0x00000000;
816 
  BID_RETURN (res);
817 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
818 
  // x is 0
819 
  res = 0x00000000;
820 
  BID_RETURN (res);
821 
} else {        // x is not special and is not zero
822 
  // x < 0 is invalid
823 
  if (x_sign) {
824 
    // set invalid flag
825 
    *pfpsf |= INVALID_EXCEPTION;
826 
    // return Integer Indefinite
827 
    res = 0x80000000;
828 
    BID_RETURN (res);
829 
  }
830 
  // x > 0 from this point on
831 
  // q = nr. of decimal digits in x
832 
  //  determine first the nr. of bits in x
833 
  if (C1.w[1] == 0) {
834 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
835 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
836 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
837 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
838 
        x_nr_bits =
839 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
840 
      } else {  // x < 2^32
841 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
842 
        x_nr_bits =
843 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
844 
      }
845 
    } else {    // if x < 2^53
846 
      tmp1.d = (double) C1.w[0]; // exact conversion
847 
      x_nr_bits =
848 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
849 
    }
850 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
851 
    tmp1.d = (double) C1.w[1];  // exact conversion
852 
    x_nr_bits =
853 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
854 
  }
855 
  q = nr_digits[x_nr_bits - 1].digits;
856 
  if (q == 0) {
857 
    q = nr_digits[x_nr_bits - 1].digits1;
858 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
859 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
860 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
861 
      q++;
862 
  }
863 
  exp = (x_exp >> 49) - 6176;
864 
 
865 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
866 
    // set invalid flag
867 
    *pfpsf |= INVALID_EXCEPTION;
868 
    // return Integer Indefinite
869 
    res = 0x80000000;
870 
    BID_RETURN (res);
871 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
872 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
873 
    // so x rounded to an integer may or may not fit in a signed 32-bit int
874 
    // the cases that do not fit are identified here; the ones that fit
875 
    // fall through and will be handled with other cases further,
876 
    // under '1 <= q + exp <= 10'
877 
    // n > 0 and q + exp = 10
878 
    // if n >= 2^32 then n is too large
879 
    // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
880 
    // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
881 
    if (q <= 11) {
882 
      tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
883 
      // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
884 
      if (tmp64 >= 0xa00000000ull) {
885 
        // set invalid flag
886 
        *pfpsf |= INVALID_EXCEPTION;
887 
        // return Integer Indefinite
888 
        res = 0x80000000;
889 
        BID_RETURN (res);
890 
      }
891 
      // else cases that can be rounded to a 32-bit int fall through
892 
      // to '1 <= q + exp <= 10'
893 
    } else {    // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
894 
      // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
895 
      // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
896 
      // (scale 2^32 up)
897 
      tmp64 = 0xa00000000ull;
898 
      if (q - 11 <= 19) {       // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
899 
        __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
900 
      } else {  // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
901 
        __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
902 
      }
903 
      if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
904 
        // set invalid flag
905 
        *pfpsf |= INVALID_EXCEPTION;
906 
        // return Integer Indefinite
907 
        res = 0x80000000;
908 
        BID_RETURN (res);
909 
      }
910 
      // else cases that can be rounded to a 32-bit int fall through
911 
      // to '1 <= q + exp <= 10'
912 
    }
913 
  }
914 
  // n is not too large to be converted to int32: 0 <= n < 2^31
915 
  // Note: some of the cases tested for above fall through to this point
916 
  if ((q + exp) <= 0) {
917 
    // n = +0.0...c(0)c(1)...c(q-1) or n = +0.c(0)c(1)...c(q-1)
918 
    // return 0
919 
    res = 0x00000000;
920 
    BID_RETURN (res);
921 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
922 
    // 1 <= x < 2^32 so x can be rounded down to a 32-bit unsigned integer
923 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
924 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
925 
      // chop off ind digits from the lower part of C1
926 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
927 
      tmp64 = C1.w[0];
928 
      if (ind <= 19) {
929 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
930 
      } else {
931 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
932 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
933 
      }
934 
      if (C1.w[0] < tmp64)
935 
        C1.w[1]++;
936 
      // calculate C* and f*
937 
      // C* is actually floor(C*) in this case
938 
      // C* and f* need shifting and masking, as shown by
939 
      // shiftright128[] and maskhigh128[]
940 
      // 1 <= x <= 33
941 
      // kx = 10^(-x) = ten2mk128[ind - 1]
942 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
943 
      // the approximation of 10^(-x) was rounded up to 118 bits
944 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
945 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
946 
        Cstar.w[1] = P256.w[3];
947 
        Cstar.w[0] = P256.w[2];
948 
        fstar.w[3] = 0;
949 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
950 
        fstar.w[1] = P256.w[1];
951 
        fstar.w[0] = P256.w[0];
952 
      } else {  // 22 <= ind - 1 <= 33
953 
        Cstar.w[1] = 0;
954 
        Cstar.w[0] = P256.w[3];
955 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
956 
        fstar.w[2] = P256.w[2];
957 
        fstar.w[1] = P256.w[1];
958 
        fstar.w[0] = P256.w[0];
959 
      }
960 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
961 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
962 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
963 
      //   if floor(C*) is even then C* = floor(C*) - logical right
964 
      //       shift; C* has p decimal digits, correct by Prop. 1)
965 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
966 
      //       shift; C* has p decimal digits, correct by Pr. 1)
967 
      // else
968 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
969 
      //       correct by Property 1)
970 
      // n = C* * 10^(e+x)
971 
 
972 
      // shift right C* by Ex-128 = shiftright128[ind]
973 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
974 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
975 
        Cstar.w[0] =
976 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
977 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
978 
      } else {  // 22 <= ind - 1 <= 33
979 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
980 
      }
981 
      // determine inexactness of the rounding of C*
982 
      // if (0 < f* - 1/2 < 10^(-x)) then
983 
      //   the result is exact
984 
      // else // if (f* - 1/2 > T*) then
985 
      //   the result is inexact
986 
      if (ind - 1 <= 2) {
987 
        if (fstar.w[1] > 0x8000000000000000ull ||
988 
            (fstar.w[1] == 0x8000000000000000ull
989 
             && fstar.w[0] > 0x0ull)) {
990 
          // f* > 1/2 and the result may be exact
991 
          tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
992 
          if (tmp64 > ten2mk128trunc[ind - 1].w[1]
993 
              || (tmp64 == ten2mk128trunc[ind - 1].w[1]
994 
                  && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
995 
          }     // else the result is exact
996 
        } else {        // the result is inexact; f2* <= 1/2
997 
          is_inexact_gt_midpoint = 1;
998 
        }
999 
      } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1000 
        if (fstar.w[3] > 0x0 ||
1001 
            (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
1002 
            (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
1003 
             (fstar.w[1] || fstar.w[0]))) {
1004 
          // f2* > 1/2 and the result may be exact
1005 
          // Calculate f2* - 1/2
1006 
          tmp64 = fstar.w[2] - onehalf128[ind - 1];
1007 
          tmp64A = fstar.w[3];
1008 
          if (tmp64 > fstar.w[2])
1009 
            tmp64A--;
1010 
          if (tmp64A || tmp64
1011 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1012 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1013 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1014 
          }     // else the result is exact
1015 
        } else {        // the result is inexact; f2* <= 1/2
1016 
          is_inexact_gt_midpoint = 1;
1017 
        }
1018 
      } else {  // if 22 <= ind <= 33
1019 
        if (fstar.w[3] > onehalf128[ind - 1] ||
1020 
            (fstar.w[3] == onehalf128[ind - 1] &&
1021 
             (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
1022 
          // f2* > 1/2 and the result may be exact
1023 
          // Calculate f2* - 1/2
1024 
          tmp64 = fstar.w[3] - onehalf128[ind - 1];
1025 
          if (tmp64 || fstar.w[2]
1026 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1027 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1028 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1029 
          }     // else the result is exact
1030 
        } else {        // the result is inexact; f2* <= 1/2
1031 
          is_inexact_gt_midpoint = 1;
1032 
        }
1033 
      }
1034 
 
1035 
      // if the result was a midpoint it was rounded away from zero, so
1036 
      // it will need a correction
1037 
      // check for midpoints
1038 
      if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
1039 
          && (fstar.w[1] || fstar.w[0])
1040 
          && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
1041 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1042 
                  && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
1043 
        // the result is a midpoint; round to nearest
1044 
        if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1045 
          // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1046 
          Cstar.w[0]--;  // Cstar.w[0] is now even
1047 
          is_inexact_gt_midpoint = 0;
1048 
        } else {        // else MP in [ODD, EVEN]
1049 
          is_midpoint_lt_even = 1;
1050 
          is_inexact_gt_midpoint = 0;
1051 
        }
1052 
      }
1053 
      // general correction for RM
1054 
      if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
1055 
        Cstar.w[0] = Cstar.w[0] - 1;
1056 
      } else {
1057 
        ;       // the result is already correct
1058 
      }
1059 
      res = Cstar.w[0];
1060 
    } else if (exp == 0) {
1061 
      // 1 <= q <= 10
1062 
      // res = +C (exact)
1063 
      res = C1.w[0];
1064 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1065 
      // res = +C * 10^exp (exact)
1066 
      res = C1.w[0] * ten2k64[exp];
1067 
    }
1068 
  }
1069 
}
1070 
 
1071 
BID_RETURN (res);
1072 
}
1073 
 
1074 
/*****************************************************************************
1075 
 *  BID128_to_uint32_xfloor
1076 
 ****************************************************************************/
1077 
 
1078 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1079 
                                          bid128_to_uint32_xfloor, x)
1080 
 
1081 
     unsigned int res;
1082 
     UINT64 x_sign;
1083 
     UINT64 x_exp;
1084 
     int exp;                   // unbiased exponent
1085 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1086 
     UINT64 tmp64, tmp64A;
1087 
     BID_UI64DOUBLE tmp1;
1088 
     unsigned int x_nr_bits;
1089 
     int q, ind, shift;
1090 
     UINT128 C1, C;
1091 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1092 
     UINT256 fstar;
1093 
     UINT256 P256;
1094 
     int is_inexact_gt_midpoint = 0;
1095 
     int is_midpoint_lt_even = 0;
1096 
 
1097 
  // unpack x
1098 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1099 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1100 
C1.w[1] = x.w[1] & MASK_COEFF;
1101 
C1.w[0] = x.w[0];
1102 
 
1103 
  // check for NaN or Infinity
1104 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1105 
    // x is special
1106 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1107 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1108 
    // set invalid flag
1109 
    *pfpsf |= INVALID_EXCEPTION;
1110 
    // return Integer Indefinite
1111 
    res = 0x80000000;
1112 
  } else {      // x is QNaN
1113 
    // set invalid flag
1114 
    *pfpsf |= INVALID_EXCEPTION;
1115 
    // return Integer Indefinite
1116 
    res = 0x80000000;
1117 
  }
1118 
  BID_RETURN (res);
1119 
} else {        // x is not a NaN, so it must be infinity
1120 
  if (!x_sign) {        // x is +inf
1121 
    // set invalid flag
1122 
    *pfpsf |= INVALID_EXCEPTION;
1123 
    // return Integer Indefinite
1124 
    res = 0x80000000;
1125 
  } else {      // x is -inf
1126 
    // set invalid flag
1127 
    *pfpsf |= INVALID_EXCEPTION;
1128 
    // return Integer Indefinite
1129 
    res = 0x80000000;
1130 
  }
1131 
  BID_RETURN (res);
1132 
}
1133 
}
1134 
  // check for non-canonical values (after the check for special values)
1135 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
1136 
    || (C1.w[1] == 0x0001ed09bead87c0ull
1137 
        && (C1.w[0] > 0x378d8e63ffffffffull))
1138 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1139 
  res = 0x00000000;
1140 
  BID_RETURN (res);
1141 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1142 
  // x is 0
1143 
  res = 0x00000000;
1144 
  BID_RETURN (res);
1145 
} else {        // x is not special and is not zero
1146 
  // x < 0 is invalid
1147 
  if (x_sign) {
1148 
    // set invalid flag
1149 
    *pfpsf |= INVALID_EXCEPTION;
1150 
    // return Integer Indefinite
1151 
    res = 0x80000000;
1152 
    BID_RETURN (res);
1153 
  }
1154 
  // x > 0 from this point on
1155 
  // q = nr. of decimal digits in x
1156 
  //  determine first the nr. of bits in x
1157 
  if (C1.w[1] == 0) {
1158 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
1159 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
1160 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
1161 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
1162 
        x_nr_bits =
1163 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1164 
      } else {  // x < 2^32
1165 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
1166 
        x_nr_bits =
1167 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1168 
      }
1169 
    } else {    // if x < 2^53
1170 
      tmp1.d = (double) C1.w[0]; // exact conversion
1171 
      x_nr_bits =
1172 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1173 
    }
1174 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1175 
    tmp1.d = (double) C1.w[1];  // exact conversion
1176 
    x_nr_bits =
1177 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1178 
  }
1179 
  q = nr_digits[x_nr_bits - 1].digits;
1180 
  if (q == 0) {
1181 
    q = nr_digits[x_nr_bits - 1].digits1;
1182 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1183 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1184 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1185 
      q++;
1186 
  }
1187 
  exp = (x_exp >> 49) - 6176;
1188 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1189 
    // set invalid flag
1190 
    *pfpsf |= INVALID_EXCEPTION;
1191 
    // return Integer Indefinite
1192 
    res = 0x80000000;
1193 
    BID_RETURN (res);
1194 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1195 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1196 
    // so x rounded to an integer may or may not fit in a signed 32-bit int
1197 
    // the cases that do not fit are identified here; the ones that fit
1198 
    // fall through and will be handled with other cases further,
1199 
    // under '1 <= q + exp <= 10'
1200 
    // n > 0 and q + exp = 10
1201 
    // if n >= 2^32 then n is too large
1202 
    // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
1203 
    // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
1204 
    if (q <= 11) {
1205 
      tmp64 = C1.w[0] * ten2k64[11 - q]; // C scaled up to 11-digit int
1206 
      // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1207 
      if (tmp64 >= 0xa00000000ull) {
1208 
        // set invalid flag
1209 
        *pfpsf |= INVALID_EXCEPTION;
1210 
        // return Integer Indefinite
1211 
        res = 0x80000000;
1212 
        BID_RETURN (res);
1213 
      }
1214 
      // else cases that can be rounded to a 32-bit int fall through
1215 
      // to '1 <= q + exp <= 10'
1216 
    } else {    // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1217 
      // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
1218 
      // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
1219 
      // (scale 2^32 up)
1220 
      tmp64 = 0xa00000000ull;
1221 
      if (q - 11 <= 19) {       // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1222 
        __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1223 
      } else {  // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1224 
        __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1225 
      }
1226 
      if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1227 
        // set invalid flag
1228 
        *pfpsf |= INVALID_EXCEPTION;
1229 
        // return Integer Indefinite
1230 
        res = 0x80000000;
1231 
        BID_RETURN (res);
1232 
      }
1233 
      // else cases that can be rounded to a 32-bit int fall through
1234 
      // to '1 <= q + exp <= 10'
1235 
    }
1236 
  }
1237 
  // n is not too large to be converted to int32: 0 <= n < 2^31
1238 
  // Note: some of the cases tested for above fall through to this point
1239 
  if ((q + exp) <= 0) {
1240 
    // n = +0.0...c(0)c(1)...c(q-1) or n = +0.c(0)c(1)...c(q-1)
1241 
    // set inexact flag
1242 
    *pfpsf |= INEXACT_EXCEPTION;
1243 
    // return 0
1244 
    res = 0x00000000;
1245 
    BID_RETURN (res);
1246 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1247 
    // 1 <= x < 2^32 so x can be rounded down to a 32-bit unsigned integer
1248 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1249 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1250 
      // chop off ind digits from the lower part of C1
1251 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1252 
      tmp64 = C1.w[0];
1253 
      if (ind <= 19) {
1254 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1255 
      } else {
1256 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1257 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1258 
      }
1259 
      if (C1.w[0] < tmp64)
1260 
        C1.w[1]++;
1261 
      // calculate C* and f*
1262 
      // C* is actually floor(C*) in this case
1263 
      // C* and f* need shifting and masking, as shown by
1264 
      // shiftright128[] and maskhigh128[]
1265 
      // 1 <= x <= 33
1266 
      // kx = 10^(-x) = ten2mk128[ind - 1]
1267 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1268 
      // the approximation of 10^(-x) was rounded up to 118 bits
1269 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1270 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1271 
        Cstar.w[1] = P256.w[3];
1272 
        Cstar.w[0] = P256.w[2];
1273 
        fstar.w[3] = 0;
1274 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1275 
        fstar.w[1] = P256.w[1];
1276 
        fstar.w[0] = P256.w[0];
1277 
      } else {  // 22 <= ind - 1 <= 33
1278 
        Cstar.w[1] = 0;
1279 
        Cstar.w[0] = P256.w[3];
1280 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1281 
        fstar.w[2] = P256.w[2];
1282 
        fstar.w[1] = P256.w[1];
1283 
        fstar.w[0] = P256.w[0];
1284 
      }
1285 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1286 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1287 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
1288 
      //   if floor(C*) is even then C* = floor(C*) - logical right
1289 
      //       shift; C* has p decimal digits, correct by Prop. 1)
1290 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
1291 
      //       shift; C* has p decimal digits, correct by Pr. 1)
1292 
      // else
1293 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
1294 
      //       correct by Property 1)
1295 
      // n = C* * 10^(e+x)
1296 
 
1297 
      // shift right C* by Ex-128 = shiftright128[ind]
1298 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1299 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1300 
        Cstar.w[0] =
1301 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1302 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1303 
      } else {  // 22 <= ind - 1 <= 33
1304 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
1305 
      }
1306 
      // determine inexactness of the rounding of C*
1307 
      // if (0 < f* - 1/2 < 10^(-x)) then
1308 
      //   the result is exact
1309 
      // else // if (f* - 1/2 > T*) then
1310 
      //   the result is inexact
1311 
      if (ind - 1 <= 2) {
1312 
        if (fstar.w[1] > 0x8000000000000000ull ||
1313 
            (fstar.w[1] == 0x8000000000000000ull
1314 
             && fstar.w[0] > 0x0ull)) {
1315 
          // f* > 1/2 and the result may be exact
1316 
          tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
1317 
          if (tmp64 > ten2mk128trunc[ind - 1].w[1]
1318 
              || (tmp64 == ten2mk128trunc[ind - 1].w[1]
1319 
                  && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
1320 
            // set the inexact flag
1321 
            *pfpsf |= INEXACT_EXCEPTION;
1322 
          }     // else the result is exact
1323 
        } else {        // the result is inexact; f2* <= 1/2
1324 
          // set the inexact flag
1325 
          *pfpsf |= INEXACT_EXCEPTION;
1326 
          is_inexact_gt_midpoint = 1;
1327 
        }
1328 
      } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1329 
        if (fstar.w[3] > 0x0 ||
1330 
            (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
1331 
            (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
1332 
             (fstar.w[1] || fstar.w[0]))) {
1333 
          // f2* > 1/2 and the result may be exact
1334 
          // Calculate f2* - 1/2
1335 
          tmp64 = fstar.w[2] - onehalf128[ind - 1];
1336 
          tmp64A = fstar.w[3];
1337 
          if (tmp64 > fstar.w[2])
1338 
            tmp64A--;
1339 
          if (tmp64A || tmp64
1340 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1341 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1342 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1343 
            // set the inexact flag
1344 
            *pfpsf |= INEXACT_EXCEPTION;
1345 
          }     // else the result is exact
1346 
        } else {        // the result is inexact; f2* <= 1/2
1347 
          // set the inexact flag
1348 
          *pfpsf |= INEXACT_EXCEPTION;
1349 
          is_inexact_gt_midpoint = 1;
1350 
        }
1351 
      } else {  // if 22 <= ind <= 33
1352 
        if (fstar.w[3] > onehalf128[ind - 1] ||
1353 
            (fstar.w[3] == onehalf128[ind - 1] &&
1354 
             (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
1355 
          // f2* > 1/2 and the result may be exact
1356 
          // Calculate f2* - 1/2
1357 
          tmp64 = fstar.w[3] - onehalf128[ind - 1];
1358 
          if (tmp64 || fstar.w[2]
1359 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1360 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1361 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1362 
            // set the inexact flag
1363 
            *pfpsf |= INEXACT_EXCEPTION;
1364 
          }     // else the result is exact
1365 
        } else {        // the result is inexact; f2* <= 1/2
1366 
          // set the inexact flag
1367 
          *pfpsf |= INEXACT_EXCEPTION;
1368 
          is_inexact_gt_midpoint = 1;
1369 
        }
1370 
      }
1371 
 
1372 
      // if the result was a midpoint it was rounded away from zero, so
1373 
      // it will need a correction
1374 
      // check for midpoints
1375 
      if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
1376 
          && (fstar.w[1] || fstar.w[0])
1377 
          && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
1378 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1379 
                  && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
1380 
        // the result is a midpoint; round to nearest
1381 
        if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1382 
          // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1383 
          Cstar.w[0]--;  // Cstar.w[0] is now even
1384 
          is_inexact_gt_midpoint = 0;
1385 
        } else {        // else MP in [ODD, EVEN]
1386 
          is_midpoint_lt_even = 1;
1387 
          is_inexact_gt_midpoint = 0;
1388 
        }
1389 
      }
1390 
      // general correction for RM
1391 
      if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
1392 
        Cstar.w[0] = Cstar.w[0] - 1;
1393 
      } else {
1394 
        ;       // the result is already correct
1395 
      }
1396 
      res = Cstar.w[0];
1397 
    } else if (exp == 0) {
1398 
      // 1 <= q <= 10
1399 
      // res = +C (exact)
1400 
      res = C1.w[0];
1401 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1402 
      // res = +C * 10^exp (exact)
1403 
      res = C1.w[0] * ten2k64[exp];
1404 
    }
1405 
  }
1406 
}
1407 
 
1408 
BID_RETURN (res);
1409 
}
1410 
 
1411 
/*****************************************************************************
1412 
 *  BID128_to_uint32_ceil
1413 
 ****************************************************************************/
1414 
 
1415 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1416 
                                          bid128_to_uint32_ceil, x)
1417 
 
1418 
     unsigned int res;
1419 
     UINT64 x_sign;
1420 
     UINT64 x_exp;
1421 
     int exp;                   // unbiased exponent
1422 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1423 
     UINT64 tmp64, tmp64A;
1424 
     BID_UI64DOUBLE tmp1;
1425 
     unsigned int x_nr_bits;
1426 
     int q, ind, shift;
1427 
     UINT128 C1, C;
1428 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1429 
     UINT256 fstar;
1430 
     UINT256 P256;
1431 
     int is_inexact_lt_midpoint = 0;
1432 
     int is_midpoint_gt_even = 0;
1433 
 
1434 
  // unpack x
1435 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1436 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1437 
C1.w[1] = x.w[1] & MASK_COEFF;
1438 
C1.w[0] = x.w[0];
1439 
 
1440 
  // check for NaN or Infinity
1441 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1442 
    // x is special
1443 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1444 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1445 
    // set invalid flag
1446 
    *pfpsf |= INVALID_EXCEPTION;
1447 
    // return Integer Indefinite
1448 
    res = 0x80000000;
1449 
  } else {      // x is QNaN
1450 
    // set invalid flag
1451 
    *pfpsf |= INVALID_EXCEPTION;
1452 
    // return Integer Indefinite
1453 
    res = 0x80000000;
1454 
  }
1455 
  BID_RETURN (res);
1456 
} else {        // x is not a NaN, so it must be infinity
1457 
  if (!x_sign) {        // x is +inf
1458 
    // set invalid flag
1459 
    *pfpsf |= INVALID_EXCEPTION;
1460 
    // return Integer Indefinite
1461 
    res = 0x80000000;
1462 
  } else {      // x is -inf
1463 
    // set invalid flag
1464 
    *pfpsf |= INVALID_EXCEPTION;
1465 
    // return Integer Indefinite
1466 
    res = 0x80000000;
1467 
  }
1468 
  BID_RETURN (res);
1469 
}
1470 
}
1471 
  // check for non-canonical values (after the check for special values)
1472 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
1473 
    || (C1.w[1] == 0x0001ed09bead87c0ull
1474 
        && (C1.w[0] > 0x378d8e63ffffffffull))
1475 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1476 
  res = 0x00000000;
1477 
  BID_RETURN (res);
1478 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1479 
  // x is 0
1480 
  res = 0x00000000;
1481 
  BID_RETURN (res);
1482 
} else {        // x is not special and is not zero
1483 
 
1484 
  // q = nr. of decimal digits in x
1485 
  //  determine first the nr. of bits in x
1486 
  if (C1.w[1] == 0) {
1487 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
1488 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
1489 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
1490 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
1491 
        x_nr_bits =
1492 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1493 
      } else {  // x < 2^32
1494 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
1495 
        x_nr_bits =
1496 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1497 
      }
1498 
    } else {    // if x < 2^53
1499 
      tmp1.d = (double) C1.w[0]; // exact conversion
1500 
      x_nr_bits =
1501 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1502 
    }
1503 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1504 
    tmp1.d = (double) C1.w[1];  // exact conversion
1505 
    x_nr_bits =
1506 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1507 
  }
1508 
  q = nr_digits[x_nr_bits - 1].digits;
1509 
  if (q == 0) {
1510 
    q = nr_digits[x_nr_bits - 1].digits1;
1511 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1512 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1513 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1514 
      q++;
1515 
  }
1516 
  exp = (x_exp >> 49) - 6176;
1517 
 
1518 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1519 
    // set invalid flag
1520 
    *pfpsf |= INVALID_EXCEPTION;
1521 
    // return Integer Indefinite
1522 
    res = 0x80000000;
1523 
    BID_RETURN (res);
1524 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1525 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1526 
    // so x rounded to an integer may or may not fit in a signed 32-bit int
1527 
    // the cases that do not fit are identified here; the ones that fit
1528 
    // fall through and will be handled with other cases further,
1529 
    // under '1 <= q + exp <= 10'
1530 
    if (x_sign) {       // if n < 0 and q + exp = 10
1531 
      // if n <= -1 then n is too large
1532 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
1533 
      // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
1534 
      if (q <= 11) {
1535 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
1536 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1537 
        if (tmp64 >= 0x0aull) {
1538 
          // set invalid flag
1539 
          *pfpsf |= INVALID_EXCEPTION;
1540 
          // return Integer Indefinite
1541 
          res = 0x80000000;
1542 
          BID_RETURN (res);
1543 
        }
1544 
        // else cases that can be rounded to a 32-bit int fall through
1545 
        // to '1 <= q + exp <= 10'
1546 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1547 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
1548 
        // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
1549 
        // (scale 1 up)
1550 
        tmp64 = 0x0aull;
1551 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1552 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1553 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1554 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1555 
        }
1556 
        if (C1.w[1] > C.w[1]
1557 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
1558 
          // set invalid flag
1559 
          *pfpsf |= INVALID_EXCEPTION;
1560 
          // return Integer Indefinite
1561 
          res = 0x80000000;
1562 
          BID_RETURN (res);
1563 
        }
1564 
        // else cases that can be rounded to a 32-bit int fall through
1565 
        // to '1 <= q + exp <= 10'
1566 
      }
1567 
    } else {    // if n > 0 and q + exp = 10
1568 
      // if n > 2^32 - 1 then n is too large
1569 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
1570 
      // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=34
1571 
      if (q <= 11) {
1572 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
1573 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1574 
        if (tmp64 > 0x9fffffff6ull) {
1575 
          // set invalid flag
1576 
          *pfpsf |= INVALID_EXCEPTION;
1577 
          // return Integer Indefinite
1578 
          res = 0x80000000;
1579 
          BID_RETURN (res);
1580 
        }
1581 
        // else cases that can be rounded to a 32-bit int fall through
1582 
        // to '1 <= q + exp <= 10'
1583 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1584 
        // 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6 <=>
1585 
        // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
1586 
        // (scale 2^32 up)
1587 
        tmp64 = 0x9fffffff6ull;
1588 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1589 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1590 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1591 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1592 
        }
1593 
        if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1594 
          // set invalid flag
1595 
          *pfpsf |= INVALID_EXCEPTION;
1596 
          // return Integer Indefinite
1597 
          res = 0x80000000;
1598 
          BID_RETURN (res);
1599 
        }
1600 
        // else cases that can be rounded to a 32-bit int fall through
1601 
        // to '1 <= q + exp <= 10'
1602 
      }
1603 
    }
1604 
  }
1605 
  // n is not too large to be converted to int32: -2^32-1 < n <= 2^32-1
1606 
  // Note: some of the cases tested for above fall through to this point
1607 
  if ((q + exp) <= 0) {
1608 
    // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
1609 
    // return 0
1610 
    if (x_sign)
1611 
      res = 0x00000000;
1612 
    else
1613 
      res = 0x00000001;
1614 
    BID_RETURN (res);
1615 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1616 
    // -2^32-1 < x <= -1 or 1 <= x <= 2^32-1 so x can be rounded
1617 
    // toward positive infinity to a 32-bit signed integer
1618 
    if (x_sign) {       // x <= -1 is invalid
1619 
      // set invalid flag
1620 
      *pfpsf |= INVALID_EXCEPTION;
1621 
      // return Integer Indefinite
1622 
      res = 0x80000000;
1623 
      BID_RETURN (res);
1624 
    }
1625 
    // x > 0 from this point on
1626 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1627 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1628 
      // chop off ind digits from the lower part of C1
1629 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1630 
      tmp64 = C1.w[0];
1631 
      if (ind <= 19) {
1632 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1633 
      } else {
1634 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1635 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1636 
      }
1637 
      if (C1.w[0] < tmp64)
1638 
        C1.w[1]++;
1639 
      // calculate C* and f*
1640 
      // C* is actually floor(C*) in this case
1641 
      // C* and f* need shifting and masking, as shown by
1642 
      // shiftright128[] and maskhigh128[]
1643 
      // 1 <= x <= 33
1644 
      // kx = 10^(-x) = ten2mk128[ind - 1]
1645 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1646 
      // the approximation of 10^(-x) was rounded up to 118 bits
1647 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1648 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1649 
        Cstar.w[1] = P256.w[3];
1650 
        Cstar.w[0] = P256.w[2];
1651 
        fstar.w[3] = 0;
1652 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1653 
        fstar.w[1] = P256.w[1];
1654 
        fstar.w[0] = P256.w[0];
1655 
      } else {  // 22 <= ind - 1 <= 33
1656 
        Cstar.w[1] = 0;
1657 
        Cstar.w[0] = P256.w[3];
1658 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1659 
        fstar.w[2] = P256.w[2];
1660 
        fstar.w[1] = P256.w[1];
1661 
        fstar.w[0] = P256.w[0];
1662 
      }
1663 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1664 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1665 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
1666 
      //   if floor(C*) is even then C* = floor(C*) - logical right
1667 
      //       shift; C* has p decimal digits, correct by Prop. 1)
1668 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
1669 
      //       shift; C* has p decimal digits, correct by Pr. 1)
1670 
      // else
1671 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
1672 
      //       correct by Property 1)
1673 
      // n = C* * 10^(e+x)
1674 
 
1675 
      // shift right C* by Ex-128 = shiftright128[ind]
1676 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
1677 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
1678 
        Cstar.w[0] =
1679 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
1680 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
1681 
      } else {  // 22 <= ind - 1 <= 33
1682 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
1683 
      }
1684 
      // determine inexactness of the rounding of C*
1685 
      // if (0 < f* - 1/2 < 10^(-x)) then
1686 
      //   the result is exact
1687 
      // else // if (f* - 1/2 > T*) then
1688 
      //   the result is inexact
1689 
      if (ind - 1 <= 2) {
1690 
        if (fstar.w[1] > 0x8000000000000000ull ||
1691 
            (fstar.w[1] == 0x8000000000000000ull
1692 
             && fstar.w[0] > 0x0ull)) {
1693 
          // f* > 1/2 and the result may be exact
1694 
          tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
1695 
          if (tmp64 > ten2mk128trunc[ind - 1].w[1]
1696 
              || (tmp64 == ten2mk128trunc[ind - 1].w[1]
1697 
                  && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
1698 
            is_inexact_lt_midpoint = 1;
1699 
          }     // else the result is exact
1700 
        } else {        // the result is inexact; f2* <= 1/2
1701 
          ;
1702 
        }
1703 
      } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
1704 
        if (fstar.w[3] > 0x0 ||
1705 
            (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
1706 
            (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
1707 
             (fstar.w[1] || fstar.w[0]))) {
1708 
          // f2* > 1/2 and the result may be exact
1709 
          // Calculate f2* - 1/2
1710 
          tmp64 = fstar.w[2] - onehalf128[ind - 1];
1711 
          tmp64A = fstar.w[3];
1712 
          if (tmp64 > fstar.w[2])
1713 
            tmp64A--;
1714 
          if (tmp64A || tmp64
1715 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1716 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1717 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1718 
            is_inexact_lt_midpoint = 1;
1719 
          }     // else the result is exact
1720 
        } else {        // the result is inexact; f2* <= 1/2
1721 
        }
1722 
      } else {  // if 22 <= ind <= 33
1723 
        if (fstar.w[3] > onehalf128[ind - 1] ||
1724 
            (fstar.w[3] == onehalf128[ind - 1] &&
1725 
             (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
1726 
          // f2* > 1/2 and the result may be exact
1727 
          // Calculate f2* - 1/2
1728 
          tmp64 = fstar.w[3] - onehalf128[ind - 1];
1729 
          if (tmp64 || fstar.w[2]
1730 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
1731 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1732 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
1733 
            is_inexact_lt_midpoint = 1;
1734 
          }     // else the result is exact
1735 
        } else {        // the result is inexact; f2* <= 1/2
1736 
          ;
1737 
        }
1738 
      }
1739 
 
1740 
      // if the result was a midpoint it was rounded away from zero, so
1741 
      // it will need a correction
1742 
      // check for midpoints
1743 
      if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
1744 
          && (fstar.w[1] || fstar.w[0])
1745 
          && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
1746 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
1747 
                  && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
1748 
        // the result is a midpoint; round to nearest
1749 
        if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
1750 
          // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
1751 
          Cstar.w[0]--;  // Cstar.w[0] is now even
1752 
          is_midpoint_gt_even = 1;
1753 
          is_inexact_lt_midpoint = 0;
1754 
        } else {        // else MP in [ODD, EVEN]
1755 
          is_inexact_lt_midpoint = 0;
1756 
        }
1757 
      }
1758 
      // general correction for RM
1759 
      if (is_midpoint_gt_even || is_inexact_lt_midpoint) {
1760 
        Cstar.w[0] = Cstar.w[0] + 1;
1761 
      } else {
1762 
        ;       // the result is already correct
1763 
      }
1764 
      res = Cstar.w[0];
1765 
    } else if (exp == 0) {
1766 
      // 1 <= q <= 10
1767 
      // res = +C (exact)
1768 
      res = C1.w[0];
1769 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
1770 
      // res = +C * 10^exp (exact)
1771 
      res = C1.w[0] * ten2k64[exp];
1772 
    }
1773 
  }
1774 
}
1775 
 
1776 
BID_RETURN (res);
1777 
}
1778 
 
1779 
/*****************************************************************************
1780 
 *  BID128_to_uint32_xceil
1781 
 ****************************************************************************/
1782 
 
1783 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
1784 
                                          bid128_to_uint32_xceil, x)
1785 
 
1786 
     unsigned int res;
1787 
     UINT64 x_sign;
1788 
     UINT64 x_exp;
1789 
     int exp;                   // unbiased exponent
1790 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
1791 
     UINT64 tmp64, tmp64A;
1792 
     BID_UI64DOUBLE tmp1;
1793 
     unsigned int x_nr_bits;
1794 
     int q, ind, shift;
1795 
     UINT128 C1, C;
1796 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
1797 
     UINT256 fstar;
1798 
     UINT256 P256;
1799 
     int is_inexact_lt_midpoint = 0;
1800 
     int is_midpoint_gt_even = 0;
1801 
 
1802 
  // unpack x
1803 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1804 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
1805 
C1.w[1] = x.w[1] & MASK_COEFF;
1806 
C1.w[0] = x.w[0];
1807 
 
1808 
  // check for NaN or Infinity
1809 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1810 
    // x is special
1811 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1812 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1813 
    // set invalid flag
1814 
    *pfpsf |= INVALID_EXCEPTION;
1815 
    // return Integer Indefinite
1816 
    res = 0x80000000;
1817 
  } else {      // x is QNaN
1818 
    // set invalid flag
1819 
    *pfpsf |= INVALID_EXCEPTION;
1820 
    // return Integer Indefinite
1821 
    res = 0x80000000;
1822 
  }
1823 
  BID_RETURN (res);
1824 
} else {        // x is not a NaN, so it must be infinity
1825 
  if (!x_sign) {        // x is +inf
1826 
    // set invalid flag
1827 
    *pfpsf |= INVALID_EXCEPTION;
1828 
    // return Integer Indefinite
1829 
    res = 0x80000000;
1830 
  } else {      // x is -inf
1831 
    // set invalid flag
1832 
    *pfpsf |= INVALID_EXCEPTION;
1833 
    // return Integer Indefinite
1834 
    res = 0x80000000;
1835 
  }
1836 
  BID_RETURN (res);
1837 
}
1838 
}
1839 
  // check for non-canonical values (after the check for special values)
1840 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
1841 
    || (C1.w[1] == 0x0001ed09bead87c0ull
1842 
        && (C1.w[0] > 0x378d8e63ffffffffull))
1843 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
1844 
  res = 0x00000000;
1845 
  BID_RETURN (res);
1846 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1847 
  // x is 0
1848 
  res = 0x00000000;
1849 
  BID_RETURN (res);
1850 
} else {        // x is not special and is not zero
1851 
 
1852 
  // q = nr. of decimal digits in x
1853 
  //  determine first the nr. of bits in x
1854 
  if (C1.w[1] == 0) {
1855 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
1856 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
1857 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
1858 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
1859 
        x_nr_bits =
1860 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1861 
      } else {  // x < 2^32
1862 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
1863 
        x_nr_bits =
1864 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1865 
      }
1866 
    } else {    // if x < 2^53
1867 
      tmp1.d = (double) C1.w[0]; // exact conversion
1868 
      x_nr_bits =
1869 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1870 
    }
1871 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1872 
    tmp1.d = (double) C1.w[1];  // exact conversion
1873 
    x_nr_bits =
1874 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1875 
  }
1876 
  q = nr_digits[x_nr_bits - 1].digits;
1877 
  if (q == 0) {
1878 
    q = nr_digits[x_nr_bits - 1].digits1;
1879 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
1880 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
1881 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1882 
      q++;
1883 
  }
1884 
  exp = (x_exp >> 49) - 6176;
1885 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
1886 
    // set invalid flag
1887 
    *pfpsf |= INVALID_EXCEPTION;
1888 
    // return Integer Indefinite
1889 
    res = 0x80000000;
1890 
    BID_RETURN (res);
1891 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
1892 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
1893 
    // so x rounded to an integer may or may not fit in a signed 32-bit int
1894 
    // the cases that do not fit are identified here; the ones that fit
1895 
    // fall through and will be handled with other cases further,
1896 
    // under '1 <= q + exp <= 10'
1897 
    if (x_sign) {       // if n < 0 and q + exp = 10
1898 
      // if n <= -1 then n is too large
1899 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
1900 
      // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x50000000a, 1<=q<=34
1901 
      if (q <= 11) {
1902 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
1903 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1904 
        if (tmp64 >= 0x0aull) {
1905 
          // set invalid flag
1906 
          *pfpsf |= INVALID_EXCEPTION;
1907 
          // return Integer Indefinite
1908 
          res = 0x80000000;
1909 
          BID_RETURN (res);
1910 
        }
1911 
        // else cases that can be rounded to a 32-bit int fall through
1912 
        // to '1 <= q + exp <= 10'
1913 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1914 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
1915 
        // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
1916 
        // (scale 1 up)
1917 
        tmp64 = 0x0aull;
1918 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1919 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1920 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1921 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1922 
        }
1923 
        if (C1.w[1] > C.w[1]
1924 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
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 
    } else {    // if n > 0 and q + exp = 10
1935 
      // if n > 2^32 - 1 then n is too large
1936 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) > 2^32 - 1
1937 
      // too large if 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6, 1<=q<=34
1938 
      if (q <= 11) {
1939 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
1940 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
1941 
        if (tmp64 > 0x9fffffff6ull) {
1942 
          // set invalid flag
1943 
          *pfpsf |= INVALID_EXCEPTION;
1944 
          // return Integer Indefinite
1945 
          res = 0x80000000;
1946 
          BID_RETURN (res);
1947 
        }
1948 
        // else cases that can be rounded to a 32-bit int fall through
1949 
        // to '1 <= q + exp <= 10'
1950 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
1951 
        // 0.c(0)c(1)...c(q-1) * 10^11 > 0x9fffffff6 <=>
1952 
        // C > 0x9fffffff6 * 10^(q-11) where 1 <= q - 11 <= 23
1953 
        // (scale 2^32 up)
1954 
        tmp64 = 0x9fffffff6ull;
1955 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
1956 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
1957 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
1958 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
1959 
        }
1960 
        if (C1.w[1] > C.w[1] || (C1.w[1] == C.w[1] && C1.w[0] > C.w[0])) {
1961 
          // set invalid flag
1962 
          *pfpsf |= INVALID_EXCEPTION;
1963 
          // return Integer Indefinite
1964 
          res = 0x80000000;
1965 
          BID_RETURN (res);
1966 
        }
1967 
        // else cases that can be rounded to a 32-bit int fall through
1968 
        // to '1 <= q + exp <= 10'
1969 
      }
1970 
    }
1971 
  }
1972 
  // n is not too large to be converted to int32: -2^32-1 < n <= 2^32-1
1973 
  // Note: some of the cases tested for above fall through to this point
1974 
  if ((q + exp) <= 0) {
1975 
    // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
1976 
    // set inexact flag
1977 
    *pfpsf |= INEXACT_EXCEPTION;
1978 
    // return 0
1979 
    if (x_sign)
1980 
      res = 0x00000000;
1981 
    else
1982 
      res = 0x00000001;
1983 
    BID_RETURN (res);
1984 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
1985 
    // -2^32-1 < x <= -1 or 1 <= x <= 2^32-1 so x can be rounded
1986 
    // toward positive infinity to a 32-bit signed integer
1987 
    if (x_sign) {       // x <= -1 is invalid
1988 
      // set invalid flag
1989 
      *pfpsf |= INVALID_EXCEPTION;
1990 
      // return Integer Indefinite
1991 
      res = 0x80000000;
1992 
      BID_RETURN (res);
1993 
    }
1994 
    // x > 0 from this point on
1995 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
1996 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
1997 
      // chop off ind digits from the lower part of C1
1998 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
1999 
      tmp64 = C1.w[0];
2000 
      if (ind <= 19) {
2001 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2002 
      } else {
2003 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2004 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2005 
      }
2006 
      if (C1.w[0] < tmp64)
2007 
        C1.w[1]++;
2008 
      // calculate C* and f*
2009 
      // C* is actually floor(C*) in this case
2010 
      // C* and f* need shifting and masking, as shown by
2011 
      // shiftright128[] and maskhigh128[]
2012 
      // 1 <= x <= 33
2013 
      // kx = 10^(-x) = ten2mk128[ind - 1]
2014 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2015 
      // the approximation of 10^(-x) was rounded up to 118 bits
2016 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2017 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2018 
        Cstar.w[1] = P256.w[3];
2019 
        Cstar.w[0] = P256.w[2];
2020 
        fstar.w[3] = 0;
2021 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
2022 
        fstar.w[1] = P256.w[1];
2023 
        fstar.w[0] = P256.w[0];
2024 
      } else {  // 22 <= ind - 1 <= 33
2025 
        Cstar.w[1] = 0;
2026 
        Cstar.w[0] = P256.w[3];
2027 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
2028 
        fstar.w[2] = P256.w[2];
2029 
        fstar.w[1] = P256.w[1];
2030 
        fstar.w[0] = P256.w[0];
2031 
      }
2032 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2033 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2034 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
2035 
      //   if floor(C*) is even then C* = floor(C*) - logical right
2036 
      //       shift; C* has p decimal digits, correct by Prop. 1)
2037 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
2038 
      //       shift; C* has p decimal digits, correct by Pr. 1)
2039 
      // else
2040 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
2041 
      //       correct by Property 1)
2042 
      // n = C* * 10^(e+x)
2043 
 
2044 
      // shift right C* by Ex-128 = shiftright128[ind]
2045 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2046 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2047 
        Cstar.w[0] =
2048 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2049 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2050 
      } else {  // 22 <= ind - 1 <= 33
2051 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
2052 
      }
2053 
      // determine inexactness of the rounding of C*
2054 
      // if (0 < f* - 1/2 < 10^(-x)) then
2055 
      //   the result is exact
2056 
      // else // if (f* - 1/2 > T*) then
2057 
      //   the result is inexact
2058 
      if (ind - 1 <= 2) {
2059 
        if (fstar.w[1] > 0x8000000000000000ull ||
2060 
            (fstar.w[1] == 0x8000000000000000ull
2061 
             && fstar.w[0] > 0x0ull)) {
2062 
          // f* > 1/2 and the result may be exact
2063 
          tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
2064 
          if (tmp64 > ten2mk128trunc[ind - 1].w[1]
2065 
              || (tmp64 == ten2mk128trunc[ind - 1].w[1]
2066 
                  && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
2067 
            // set the inexact flag
2068 
            *pfpsf |= INEXACT_EXCEPTION;
2069 
            is_inexact_lt_midpoint = 1;
2070 
          }     // else the result is exact
2071 
        } else {        // the result is inexact; f2* <= 1/2
2072 
          // set the inexact flag
2073 
          *pfpsf |= INEXACT_EXCEPTION;
2074 
        }
2075 
      } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
2076 
        if (fstar.w[3] > 0x0 ||
2077 
            (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
2078 
            (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
2079 
             (fstar.w[1] || fstar.w[0]))) {
2080 
          // f2* > 1/2 and the result may be exact
2081 
          // Calculate f2* - 1/2
2082 
          tmp64 = fstar.w[2] - onehalf128[ind - 1];
2083 
          tmp64A = fstar.w[3];
2084 
          if (tmp64 > fstar.w[2])
2085 
            tmp64A--;
2086 
          if (tmp64A || tmp64
2087 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2088 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2089 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2090 
            // set the inexact flag
2091 
            *pfpsf |= INEXACT_EXCEPTION;
2092 
            is_inexact_lt_midpoint = 1;
2093 
          }     // else the result is exact
2094 
        } else {        // the result is inexact; f2* <= 1/2
2095 
          // set the inexact flag
2096 
          *pfpsf |= INEXACT_EXCEPTION;
2097 
        }
2098 
      } else {  // if 22 <= ind <= 33
2099 
        if (fstar.w[3] > onehalf128[ind - 1] ||
2100 
            (fstar.w[3] == onehalf128[ind - 1] &&
2101 
             (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
2102 
          // f2* > 1/2 and the result may be exact
2103 
          // Calculate f2* - 1/2
2104 
          tmp64 = fstar.w[3] - onehalf128[ind - 1];
2105 
          if (tmp64 || fstar.w[2]
2106 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2107 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2108 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2109 
            // set the inexact flag
2110 
            *pfpsf |= INEXACT_EXCEPTION;
2111 
            is_inexact_lt_midpoint = 1;
2112 
          }     // else the result is exact
2113 
        } else {        // the result is inexact; f2* <= 1/2
2114 
          // set the inexact flag
2115 
          *pfpsf |= INEXACT_EXCEPTION;
2116 
        }
2117 
      }
2118 
 
2119 
      // if the result was a midpoint it was rounded away from zero, so
2120 
      // it will need a correction
2121 
      // check for midpoints
2122 
      if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
2123 
          && (fstar.w[1] || fstar.w[0])
2124 
          && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
2125 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2126 
                  && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
2127 
        // the result is a midpoint; round to nearest
2128 
        if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2129 
          // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2130 
          Cstar.w[0]--;  // Cstar.w[0] is now even
2131 
          is_midpoint_gt_even = 1;
2132 
          is_inexact_lt_midpoint = 0;
2133 
        } else {        // else MP in [ODD, EVEN]
2134 
          is_inexact_lt_midpoint = 0;
2135 
        }
2136 
      }
2137 
      // general correction for RM
2138 
      if (is_midpoint_gt_even || is_inexact_lt_midpoint) {
2139 
        Cstar.w[0] = Cstar.w[0] + 1;
2140 
      } else {
2141 
        ;       // the result is already correct
2142 
      }
2143 
      res = Cstar.w[0];
2144 
    } else if (exp == 0) {
2145 
      // 1 <= q <= 10
2146 
      // res = +C (exact)
2147 
      res = C1.w[0];
2148 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2149 
      // res = +C * 10^exp (exact)
2150 
      res = C1.w[0] * ten2k64[exp];
2151 
    }
2152 
  }
2153 
}
2154 
 
2155 
BID_RETURN (res);
2156 
}
2157 
 
2158 
/*****************************************************************************
2159 
 *  BID128_to_uint32_int
2160 
 ****************************************************************************/
2161 
 
2162 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2163 
                                          bid128_to_uint32_int, x)
2164 
 
2165 
     int res;
2166 
     UINT64 x_sign;
2167 
     UINT64 x_exp;
2168 
     int exp;                   // unbiased exponent
2169 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2170 
     UINT64 tmp64, tmp64A;
2171 
     BID_UI64DOUBLE tmp1;
2172 
     unsigned int x_nr_bits;
2173 
     int q, ind, shift;
2174 
     UINT128 C1, C;
2175 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2176 
     UINT256 fstar;
2177 
     UINT256 P256;
2178 
     int is_inexact_gt_midpoint = 0;
2179 
     int is_midpoint_lt_even = 0;
2180 
 
2181 
  // unpack x
2182 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2183 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2184 
C1.w[1] = x.w[1] & MASK_COEFF;
2185 
C1.w[0] = x.w[0];
2186 
 
2187 
  // check for NaN or Infinity
2188 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2189 
    // x is special
2190 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2191 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2192 
    // set invalid flag
2193 
    *pfpsf |= INVALID_EXCEPTION;
2194 
    // return Integer Indefinite
2195 
    res = 0x80000000;
2196 
  } else {      // x is QNaN
2197 
    // set invalid flag
2198 
    *pfpsf |= INVALID_EXCEPTION;
2199 
    // return Integer Indefinite
2200 
    res = 0x80000000;
2201 
  }
2202 
  BID_RETURN (res);
2203 
} else {        // x is not a NaN, so it must be infinity
2204 
  if (!x_sign) {        // x is +inf
2205 
    // set invalid flag
2206 
    *pfpsf |= INVALID_EXCEPTION;
2207 
    // return Integer Indefinite
2208 
    res = 0x80000000;
2209 
  } else {      // x is -inf
2210 
    // set invalid flag
2211 
    *pfpsf |= INVALID_EXCEPTION;
2212 
    // return Integer Indefinite
2213 
    res = 0x80000000;
2214 
  }
2215 
  BID_RETURN (res);
2216 
}
2217 
}
2218 
  // check for non-canonical values (after the check for special values)
2219 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
2220 
    || (C1.w[1] == 0x0001ed09bead87c0ull
2221 
        && (C1.w[0] > 0x378d8e63ffffffffull))
2222 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2223 
  res = 0x00000000;
2224 
  BID_RETURN (res);
2225 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2226 
  // x is 0
2227 
  res = 0x00000000;
2228 
  BID_RETURN (res);
2229 
} else {        // x is not special and is not zero
2230 
 
2231 
  // q = nr. of decimal digits in x
2232 
  //  determine first the nr. of bits in x
2233 
  if (C1.w[1] == 0) {
2234 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
2235 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
2236 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
2237 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
2238 
        x_nr_bits =
2239 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2240 
      } else {  // x < 2^32
2241 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
2242 
        x_nr_bits =
2243 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2244 
      }
2245 
    } else {    // if x < 2^53
2246 
      tmp1.d = (double) C1.w[0]; // exact conversion
2247 
      x_nr_bits =
2248 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2249 
    }
2250 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2251 
    tmp1.d = (double) C1.w[1];  // exact conversion
2252 
    x_nr_bits =
2253 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2254 
  }
2255 
  q = nr_digits[x_nr_bits - 1].digits;
2256 
  if (q == 0) {
2257 
    q = nr_digits[x_nr_bits - 1].digits1;
2258 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2259 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2260 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2261 
      q++;
2262 
  }
2263 
  exp = (x_exp >> 49) - 6176;
2264 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2265 
    // set invalid flag
2266 
    *pfpsf |= INVALID_EXCEPTION;
2267 
    // return Integer Indefinite
2268 
    res = 0x80000000;
2269 
    BID_RETURN (res);
2270 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2271 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2272 
    // so x rounded to an integer may or may not fit in a signed 32-bit int
2273 
    // the cases that do not fit are identified here; the ones that fit
2274 
    // fall through and will be handled with other cases further,
2275 
    // under '1 <= q + exp <= 10'
2276 
    if (x_sign) {       // if n < 0 and q + exp = 10
2277 
      // if n <= -1 then n is too large
2278 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
2279 
      // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a, 1<=q<=34
2280 
      if (q <= 11) {
2281 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
2282 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2283 
        if (tmp64 >= 0x0aull) {
2284 
          // set invalid flag
2285 
          *pfpsf |= INVALID_EXCEPTION;
2286 
          // return Integer Indefinite
2287 
          res = 0x80000000;
2288 
          BID_RETURN (res);
2289 
        }
2290 
        // else cases that can be rounded to a 32-bit uint fall through
2291 
        // to '1 <= q + exp <= 10'
2292 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2293 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
2294 
        // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
2295 
        // (scale 1 up)
2296 
        tmp64 = 0x0aull;
2297 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2298 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
2299 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2300 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
2301 
        }
2302 
        if (C1.w[1] > C.w[1]
2303 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2304 
          // set invalid flag
2305 
          *pfpsf |= INVALID_EXCEPTION;
2306 
          // return Integer Indefinite
2307 
          res = 0x80000000;
2308 
          BID_RETURN (res);
2309 
        }
2310 
        // else cases that can be rounded to a 32-bit int fall through
2311 
        // to '1 <= q + exp <= 10'
2312 
      }
2313 
    } else {    // if n > 0 and q + exp = 10
2314 
      // if n >= 2^32 then n is too large
2315 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
2316 
      // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
2317 
      if (q <= 11) {
2318 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
2319 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2320 
        if (tmp64 >= 0xa00000000ull) {
2321 
          // set invalid flag
2322 
          *pfpsf |= INVALID_EXCEPTION;
2323 
          // return Integer Indefinite
2324 
          res = 0x80000000;
2325 
          BID_RETURN (res);
2326 
        }
2327 
        // else cases that can be rounded to a 32-bit uint fall through
2328 
        // to '1 <= q + exp <= 10'
2329 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2330 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
2331 
        // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
2332 
        // (scale 2^32 up)
2333 
        tmp64 = 0xa00000000ull;
2334 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2335 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
2336 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2337 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
2338 
        }
2339 
        if (C1.w[1] > C.w[1]
2340 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2341 
          // set invalid flag
2342 
          *pfpsf |= INVALID_EXCEPTION;
2343 
          // return Integer Indefinite
2344 
          res = 0x80000000;
2345 
          BID_RETURN (res);
2346 
        }
2347 
        // else cases that can be rounded to a 32-bit int fall through
2348 
        // to '1 <= q + exp <= 10'
2349 
      }
2350 
    }
2351 
  }
2352 
  // n is not too large to be converted to uint32: -2^32 < n < 2^32
2353 
  // Note: some of the cases tested for above fall through to this point
2354 
  if ((q + exp) <= 0) {
2355 
    // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2356 
    // return 0
2357 
    res = 0x00000000;
2358 
    BID_RETURN (res);
2359 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2360 
    // x = d(0)...d(k).d(k+1)..., k >= 0, d(0) != 0
2361 
    if (x_sign) {       // x <= -1
2362 
      // set invalid flag
2363 
      *pfpsf |= INVALID_EXCEPTION;
2364 
      // return Integer Indefinite
2365 
      res = 0x80000000;
2366 
      BID_RETURN (res);
2367 
    }
2368 
    // x > 0 from this point on
2369 
    // 1 <= x < 2^32 so x can be rounded to zero to a 32-bit unsigned integer
2370 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2371 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
2372 
      // chop off ind digits from the lower part of C1
2373 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2374 
      tmp64 = C1.w[0];
2375 
      if (ind <= 19) {
2376 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2377 
      } else {
2378 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2379 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2380 
      }
2381 
      if (C1.w[0] < tmp64)
2382 
        C1.w[1]++;
2383 
      // calculate C* and f*
2384 
      // C* is actually floor(C*) in this case
2385 
      // C* and f* need shifting and masking, as shown by
2386 
      // shiftright128[] and maskhigh128[]
2387 
      // 1 <= x <= 33
2388 
      // kx = 10^(-x) = ten2mk128[ind - 1]
2389 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2390 
      // the approximation of 10^(-x) was rounded up to 118 bits
2391 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2392 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2393 
        Cstar.w[1] = P256.w[3];
2394 
        Cstar.w[0] = P256.w[2];
2395 
        fstar.w[3] = 0;
2396 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
2397 
        fstar.w[1] = P256.w[1];
2398 
        fstar.w[0] = P256.w[0];
2399 
      } else {  // 22 <= ind - 1 <= 33
2400 
        Cstar.w[1] = 0;
2401 
        Cstar.w[0] = P256.w[3];
2402 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
2403 
        fstar.w[2] = P256.w[2];
2404 
        fstar.w[1] = P256.w[1];
2405 
        fstar.w[0] = P256.w[0];
2406 
      }
2407 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2408 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2409 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
2410 
      //   if floor(C*) is even then C* = floor(C*) - logical right
2411 
      //       shift; C* has p decimal digits, correct by Prop. 1)
2412 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
2413 
      //       shift; C* has p decimal digits, correct by Pr. 1)
2414 
      // else
2415 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
2416 
      //       correct by Property 1)
2417 
      // n = C* * 10^(e+x)
2418 
 
2419 
      // shift right C* by Ex-128 = shiftright128[ind]
2420 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2421 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2422 
        Cstar.w[0] =
2423 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2424 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2425 
      } else {  // 22 <= ind - 1 <= 33
2426 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
2427 
      }
2428 
      // determine inexactness of the rounding of C*
2429 
      // if (0 < f* - 1/2 < 10^(-x)) then
2430 
      //   the result is exact
2431 
      // else // if (f* - 1/2 > T*) then
2432 
      //   the result is inexact
2433 
      if (ind - 1 <= 2) {
2434 
        if (fstar.w[1] > 0x8000000000000000ull ||
2435 
            (fstar.w[1] == 0x8000000000000000ull
2436 
             && fstar.w[0] > 0x0ull)) {
2437 
          // f* > 1/2 and the result may be exact
2438 
          tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
2439 
          if (tmp64 > ten2mk128trunc[ind - 1].w[1]
2440 
              || (tmp64 == ten2mk128trunc[ind - 1].w[1]
2441 
                  && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
2442 
          }     // else the result is exact
2443 
        } else {        // the result is inexact; f2* <= 1/2
2444 
          is_inexact_gt_midpoint = 1;
2445 
        }
2446 
      } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
2447 
        if (fstar.w[3] > 0x0 ||
2448 
            (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
2449 
            (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
2450 
             (fstar.w[1] || fstar.w[0]))) {
2451 
          // f2* > 1/2 and the result may be exact
2452 
          // Calculate f2* - 1/2
2453 
          tmp64 = fstar.w[2] - onehalf128[ind - 1];
2454 
          tmp64A = fstar.w[3];
2455 
          if (tmp64 > fstar.w[2])
2456 
            tmp64A--;
2457 
          if (tmp64A || tmp64
2458 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2459 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2460 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2461 
          }     // else the result is exact
2462 
        } else {        // the result is inexact; f2* <= 1/2
2463 
          is_inexact_gt_midpoint = 1;
2464 
        }
2465 
      } else {  // if 22 <= ind <= 33
2466 
        if (fstar.w[3] > onehalf128[ind - 1] ||
2467 
            (fstar.w[3] == onehalf128[ind - 1] &&
2468 
             (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
2469 
          // f2* > 1/2 and the result may be exact
2470 
          // Calculate f2* - 1/2
2471 
          tmp64 = fstar.w[3] - onehalf128[ind - 1];
2472 
          if (tmp64 || fstar.w[2]
2473 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2474 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2475 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2476 
          }     // else the result is exact
2477 
        } else {        // the result is inexact; f2* <= 1/2
2478 
          is_inexact_gt_midpoint = 1;
2479 
        }
2480 
      }
2481 
 
2482 
      // if the result was a midpoint it was rounded away from zero, so
2483 
      // it will need a correction
2484 
      // check for midpoints
2485 
      if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
2486 
          && (fstar.w[1] || fstar.w[0])
2487 
          && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
2488 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2489 
                  && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
2490 
        // the result is a midpoint; round to nearest
2491 
        if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2492 
          // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2493 
          Cstar.w[0]--;  // Cstar.w[0] is now even
2494 
          is_inexact_gt_midpoint = 0;
2495 
        } else {        // else MP in [ODD, EVEN]
2496 
          is_midpoint_lt_even = 1;
2497 
          is_inexact_gt_midpoint = 0;
2498 
        }
2499 
      }
2500 
      // general correction for RZ
2501 
      if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
2502 
        Cstar.w[0] = Cstar.w[0] - 1;
2503 
      } else {
2504 
        ;       // exact, the result is already correct
2505 
      }
2506 
      res = Cstar.w[0];
2507 
    } else if (exp == 0) {
2508 
      // 1 <= q <= 10
2509 
      // res = +C (exact)
2510 
      res = C1.w[0];
2511 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2512 
      // res = +C * 10^exp (exact)
2513 
      res = C1.w[0] * ten2k64[exp];
2514 
    }
2515 
  }
2516 
}
2517 
 
2518 
BID_RETURN (res);
2519 
}
2520 
 
2521 
/*****************************************************************************
2522 
 *  BID128_to_uint32_xint
2523 
 ****************************************************************************/
2524 
 
2525 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2526 
                                          bid128_to_uint32_xint, x)
2527 
 
2528 
     int res;
2529 
     UINT64 x_sign;
2530 
     UINT64 x_exp;
2531 
     int exp;                   // unbiased exponent
2532 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2533 
     UINT64 tmp64, tmp64A;
2534 
     BID_UI64DOUBLE tmp1;
2535 
     unsigned int x_nr_bits;
2536 
     int q, ind, shift;
2537 
     UINT128 C1, C;
2538 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2539 
     UINT256 fstar;
2540 
     UINT256 P256;
2541 
     int is_inexact_gt_midpoint = 0;
2542 
     int is_midpoint_lt_even = 0;
2543 
 
2544 
  // unpack x
2545 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2546 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2547 
C1.w[1] = x.w[1] & MASK_COEFF;
2548 
C1.w[0] = x.w[0];
2549 
 
2550 
  // check for NaN or Infinity
2551 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2552 
    // x is special
2553 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2554 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2555 
    // set invalid flag
2556 
    *pfpsf |= INVALID_EXCEPTION;
2557 
    // return Integer Indefinite
2558 
    res = 0x80000000;
2559 
  } else {      // x is QNaN
2560 
    // set invalid flag
2561 
    *pfpsf |= INVALID_EXCEPTION;
2562 
    // return Integer Indefinite
2563 
    res = 0x80000000;
2564 
  }
2565 
  BID_RETURN (res);
2566 
} else {        // x is not a NaN, so it must be infinity
2567 
  if (!x_sign) {        // x is +inf
2568 
    // set invalid flag
2569 
    *pfpsf |= INVALID_EXCEPTION;
2570 
    // return Integer Indefinite
2571 
    res = 0x80000000;
2572 
  } else {      // x is -inf
2573 
    // set invalid flag
2574 
    *pfpsf |= INVALID_EXCEPTION;
2575 
    // return Integer Indefinite
2576 
    res = 0x80000000;
2577 
  }
2578 
  BID_RETURN (res);
2579 
}
2580 
}
2581 
  // check for non-canonical values (after the check for special values)
2582 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
2583 
    || (C1.w[1] == 0x0001ed09bead87c0ull
2584 
        && (C1.w[0] > 0x378d8e63ffffffffull))
2585 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2586 
  res = 0x00000000;
2587 
  BID_RETURN (res);
2588 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2589 
  // x is 0
2590 
  res = 0x00000000;
2591 
  BID_RETURN (res);
2592 
} else {        // x is not special and is not zero
2593 
 
2594 
  // q = nr. of decimal digits in x
2595 
  //  determine first the nr. of bits in x
2596 
  if (C1.w[1] == 0) {
2597 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
2598 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
2599 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
2600 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
2601 
        x_nr_bits =
2602 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2603 
      } else {  // x < 2^32
2604 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
2605 
        x_nr_bits =
2606 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2607 
      }
2608 
    } else {    // if x < 2^53
2609 
      tmp1.d = (double) C1.w[0]; // exact conversion
2610 
      x_nr_bits =
2611 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2612 
    }
2613 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2614 
    tmp1.d = (double) C1.w[1];  // exact conversion
2615 
    x_nr_bits =
2616 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2617 
  }
2618 
  q = nr_digits[x_nr_bits - 1].digits;
2619 
  if (q == 0) {
2620 
    q = nr_digits[x_nr_bits - 1].digits1;
2621 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2622 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2623 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2624 
      q++;
2625 
  }
2626 
  exp = (x_exp >> 49) - 6176;
2627 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
2628 
    // set invalid flag
2629 
    *pfpsf |= INVALID_EXCEPTION;
2630 
    // return Integer Indefinite
2631 
    res = 0x80000000;
2632 
    BID_RETURN (res);
2633 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
2634 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
2635 
    // so x rounded to an integer may or may not fit in a signed 32-bit int
2636 
    // the cases that do not fit are identified here; the ones that fit
2637 
    // fall through and will be handled with other cases further,
2638 
    // under '1 <= q + exp <= 10'
2639 
    if (x_sign) {       // if n < 0 and q + exp = 10
2640 
      // if n <= -1 then n is too large
2641 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1
2642 
      // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a, 1<=q<=34
2643 
      if (q <= 11) {
2644 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
2645 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2646 
        if (tmp64 >= 0x0aull) {
2647 
          // set invalid flag
2648 
          *pfpsf |= INVALID_EXCEPTION;
2649 
          // return Integer Indefinite
2650 
          res = 0x80000000;
2651 
          BID_RETURN (res);
2652 
        }
2653 
        // else cases that can be rounded to a 32-bit uint fall through
2654 
        // to '1 <= q + exp <= 10'
2655 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2656 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x0a <=>
2657 
        // C >= 0x0a * 10^(q-11) where 1 <= q - 11 <= 23
2658 
        // (scale 1 up)
2659 
        tmp64 = 0x0aull;
2660 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2661 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
2662 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2663 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
2664 
        }
2665 
        if (C1.w[1] > C.w[1]
2666 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2667 
          // set invalid flag
2668 
          *pfpsf |= INVALID_EXCEPTION;
2669 
          // return Integer Indefinite
2670 
          res = 0x80000000;
2671 
          BID_RETURN (res);
2672 
        }
2673 
        // else cases that can be rounded to a 32-bit int fall through
2674 
        // to '1 <= q + exp <= 10'
2675 
      }
2676 
    } else {    // if n > 0 and q + exp = 10
2677 
      // if n >= 2^32 then n is too large
2678 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32
2679 
      // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000, 1<=q<=34
2680 
      if (q <= 11) {
2681 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
2682 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
2683 
        if (tmp64 >= 0xa00000000ull) {
2684 
          // set invalid flag
2685 
          *pfpsf |= INVALID_EXCEPTION;
2686 
          // return Integer Indefinite
2687 
          res = 0x80000000;
2688 
          BID_RETURN (res);
2689 
        }
2690 
        // else cases that can be rounded to a 32-bit uint fall through
2691 
        // to '1 <= q + exp <= 10'
2692 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
2693 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0xa00000000 <=>
2694 
        // C >= 0xa00000000 * 10^(q-11) where 1 <= q - 11 <= 23
2695 
        // (scale 2^32 up)
2696 
        tmp64 = 0xa00000000ull;
2697 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
2698 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
2699 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
2700 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
2701 
        }
2702 
        if (C1.w[1] > C.w[1]
2703 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
2704 
          // set invalid flag
2705 
          *pfpsf |= INVALID_EXCEPTION;
2706 
          // return Integer Indefinite
2707 
          res = 0x80000000;
2708 
          BID_RETURN (res);
2709 
        }
2710 
        // else cases that can be rounded to a 32-bit int fall through
2711 
        // to '1 <= q + exp <= 10'
2712 
      }
2713 
    }
2714 
  }
2715 
  // n is not too large to be converted to uint32: -2^32 < n < 2^32
2716 
  // Note: some of the cases tested for above fall through to this point
2717 
  if ((q + exp) <= 0) {
2718 
    // n = +/-0.0...c(0)c(1)...c(q-1) or n = +/-0.c(0)c(1)...c(q-1)
2719 
    // set inexact flag
2720 
    *pfpsf |= INEXACT_EXCEPTION;
2721 
    // return 0
2722 
    res = 0x00000000;
2723 
    BID_RETURN (res);
2724 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
2725 
    // x = d(0)...d(k).d(k+1)..., k >= 0, d(0) != 0
2726 
    if (x_sign) {       // x <= -1
2727 
      // set invalid flag
2728 
      *pfpsf |= INVALID_EXCEPTION;
2729 
      // return Integer Indefinite
2730 
      res = 0x80000000;
2731 
      BID_RETURN (res);
2732 
    }
2733 
    // x > 0 from this point on
2734 
    // 1 <= x < 2^32 so x can be rounded to zero to a 32-bit unsigned integer
2735 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
2736 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
2737 
      // chop off ind digits from the lower part of C1
2738 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
2739 
      tmp64 = C1.w[0];
2740 
      if (ind <= 19) {
2741 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
2742 
      } else {
2743 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
2744 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
2745 
      }
2746 
      if (C1.w[0] < tmp64)
2747 
        C1.w[1]++;
2748 
      // calculate C* and f*
2749 
      // C* is actually floor(C*) in this case
2750 
      // C* and f* need shifting and masking, as shown by
2751 
      // shiftright128[] and maskhigh128[]
2752 
      // 1 <= x <= 33
2753 
      // kx = 10^(-x) = ten2mk128[ind - 1]
2754 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
2755 
      // the approximation of 10^(-x) was rounded up to 118 bits
2756 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
2757 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2758 
        Cstar.w[1] = P256.w[3];
2759 
        Cstar.w[0] = P256.w[2];
2760 
        fstar.w[3] = 0;
2761 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
2762 
        fstar.w[1] = P256.w[1];
2763 
        fstar.w[0] = P256.w[0];
2764 
      } else {  // 22 <= ind - 1 <= 33
2765 
        Cstar.w[1] = 0;
2766 
        Cstar.w[0] = P256.w[3];
2767 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
2768 
        fstar.w[2] = P256.w[2];
2769 
        fstar.w[1] = P256.w[1];
2770 
        fstar.w[0] = P256.w[0];
2771 
      }
2772 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
2773 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
2774 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
2775 
      //   if floor(C*) is even then C* = floor(C*) - logical right
2776 
      //       shift; C* has p decimal digits, correct by Prop. 1)
2777 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
2778 
      //       shift; C* has p decimal digits, correct by Pr. 1)
2779 
      // else
2780 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
2781 
      //       correct by Property 1)
2782 
      // n = C* * 10^(e+x)
2783 
 
2784 
      // shift right C* by Ex-128 = shiftright128[ind]
2785 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
2786 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
2787 
        Cstar.w[0] =
2788 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
2789 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
2790 
      } else {  // 22 <= ind - 1 <= 33
2791 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
2792 
      }
2793 
      // determine inexactness of the rounding of C*
2794 
      // if (0 < f* - 1/2 < 10^(-x)) then
2795 
      //   the result is exact
2796 
      // else // if (f* - 1/2 > T*) then
2797 
      //   the result is inexact
2798 
      if (ind - 1 <= 2) {
2799 
        if (fstar.w[1] > 0x8000000000000000ull ||
2800 
            (fstar.w[1] == 0x8000000000000000ull
2801 
             && fstar.w[0] > 0x0ull)) {
2802 
          // f* > 1/2 and the result may be exact
2803 
          tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
2804 
          if (tmp64 > ten2mk128trunc[ind - 1].w[1]
2805 
              || (tmp64 == ten2mk128trunc[ind - 1].w[1]
2806 
                  && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
2807 
            // set the inexact flag
2808 
            *pfpsf |= INEXACT_EXCEPTION;
2809 
          }     // else the result is exact
2810 
        } else {        // the result is inexact; f2* <= 1/2
2811 
          // set the inexact flag
2812 
          *pfpsf |= INEXACT_EXCEPTION;
2813 
          is_inexact_gt_midpoint = 1;
2814 
        }
2815 
      } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
2816 
        if (fstar.w[3] > 0x0 ||
2817 
            (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
2818 
            (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
2819 
             (fstar.w[1] || fstar.w[0]))) {
2820 
          // f2* > 1/2 and the result may be exact
2821 
          // Calculate f2* - 1/2
2822 
          tmp64 = fstar.w[2] - onehalf128[ind - 1];
2823 
          tmp64A = fstar.w[3];
2824 
          if (tmp64 > fstar.w[2])
2825 
            tmp64A--;
2826 
          if (tmp64A || tmp64
2827 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2828 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2829 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2830 
            // set the inexact flag
2831 
            *pfpsf |= INEXACT_EXCEPTION;
2832 
          }     // else the result is exact
2833 
        } else {        // the result is inexact; f2* <= 1/2
2834 
          // set the inexact flag
2835 
          *pfpsf |= INEXACT_EXCEPTION;
2836 
          is_inexact_gt_midpoint = 1;
2837 
        }
2838 
      } else {  // if 22 <= ind <= 33
2839 
        if (fstar.w[3] > onehalf128[ind - 1] ||
2840 
            (fstar.w[3] == onehalf128[ind - 1] &&
2841 
             (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
2842 
          // f2* > 1/2 and the result may be exact
2843 
          // Calculate f2* - 1/2
2844 
          tmp64 = fstar.w[3] - onehalf128[ind - 1];
2845 
          if (tmp64 || fstar.w[2]
2846 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
2847 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2848 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
2849 
            // set the inexact flag
2850 
            *pfpsf |= INEXACT_EXCEPTION;
2851 
          }     // else the result is exact
2852 
        } else {        // the result is inexact; f2* <= 1/2
2853 
          // set the inexact flag
2854 
          *pfpsf |= INEXACT_EXCEPTION;
2855 
          is_inexact_gt_midpoint = 1;
2856 
        }
2857 
      }
2858 
 
2859 
      // if the result was a midpoint it was rounded away from zero, so
2860 
      // it will need a correction
2861 
      // check for midpoints
2862 
      if ((fstar.w[3] == 0) && (fstar.w[2] == 0)
2863 
          && (fstar.w[1] || fstar.w[0])
2864 
          && (fstar.w[1] < ten2mk128trunc[ind - 1].w[1]
2865 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
2866 
                  && fstar.w[0] <= ten2mk128trunc[ind - 1].w[0]))) {
2867 
        // the result is a midpoint; round to nearest
2868 
        if (Cstar.w[0] & 0x01) { // Cstar.w[0] is odd; MP in [EVEN, ODD]
2869 
          // if floor(C*) is odd C = floor(C*) - 1; the result >= 1
2870 
          Cstar.w[0]--;  // Cstar.w[0] is now even
2871 
          is_inexact_gt_midpoint = 0;
2872 
        } else {        // else MP in [ODD, EVEN]
2873 
          is_midpoint_lt_even = 1;
2874 
          is_inexact_gt_midpoint = 0;
2875 
        }
2876 
      }
2877 
      // general correction for RZ
2878 
      if (is_midpoint_lt_even || is_inexact_gt_midpoint) {
2879 
        Cstar.w[0] = Cstar.w[0] - 1;
2880 
      } else {
2881 
        ;       // exact, the result is already correct
2882 
      }
2883 
      res = Cstar.w[0];
2884 
    } else if (exp == 0) {
2885 
      // 1 <= q <= 10
2886 
      // res = +C (exact)
2887 
      res = C1.w[0];
2888 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
2889 
      // res = +C * 10^exp (exact)
2890 
      res = C1.w[0] * ten2k64[exp];
2891 
    }
2892 
  }
2893 
}
2894 
 
2895 
BID_RETURN (res);
2896 
}
2897 
 
2898 
/*****************************************************************************
2899 
 *  BID128_to_uint32_rninta
2900 
 ****************************************************************************/
2901 
 
2902 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
2903 
                                          bid128_to_uint32_rninta, x)
2904 
 
2905 
     unsigned int res;
2906 
     UINT64 x_sign;
2907 
     UINT64 x_exp;
2908 
     int exp;                   // unbiased exponent
2909 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
2910 
     UINT64 tmp64;
2911 
     BID_UI64DOUBLE tmp1;
2912 
     unsigned int x_nr_bits;
2913 
     int q, ind, shift;
2914 
     UINT128 C1, C;
2915 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
2916 
     UINT256 P256;
2917 
 
2918 
  // unpack x
2919 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
2920 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
2921 
C1.w[1] = x.w[1] & MASK_COEFF;
2922 
C1.w[0] = x.w[0];
2923 
 
2924 
  // check for NaN or Infinity
2925 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
2926 
    // x is special
2927 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
2928 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
2929 
    // set invalid flag
2930 
    *pfpsf |= INVALID_EXCEPTION;
2931 
    // return Integer Indefinite
2932 
    res = 0x80000000;
2933 
  } else {      // x is QNaN
2934 
    // set invalid flag
2935 
    *pfpsf |= INVALID_EXCEPTION;
2936 
    // return Integer Indefinite
2937 
    res = 0x80000000;
2938 
  }
2939 
  BID_RETURN (res);
2940 
} else {        // x is not a NaN, so it must be infinity
2941 
  if (!x_sign) {        // x is +inf
2942 
    // set invalid flag
2943 
    *pfpsf |= INVALID_EXCEPTION;
2944 
    // return Integer Indefinite
2945 
    res = 0x80000000;
2946 
  } else {      // x is -inf
2947 
    // set invalid flag
2948 
    *pfpsf |= INVALID_EXCEPTION;
2949 
    // return Integer Indefinite
2950 
    res = 0x80000000;
2951 
  }
2952 
  BID_RETURN (res);
2953 
}
2954 
}
2955 
  // check for non-canonical values (after the check for special values)
2956 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
2957 
    || (C1.w[1] == 0x0001ed09bead87c0ull
2958 
        && (C1.w[0] > 0x378d8e63ffffffffull))
2959 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
2960 
  res = 0x00000000;
2961 
  BID_RETURN (res);
2962 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
2963 
  // x is 0
2964 
  res = 0x00000000;
2965 
  BID_RETURN (res);
2966 
} else {        // x is not special and is not zero
2967 
 
2968 
  // q = nr. of decimal digits in x
2969 
  //  determine first the nr. of bits in x
2970 
  if (C1.w[1] == 0) {
2971 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
2972 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
2973 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
2974 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
2975 
        x_nr_bits =
2976 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2977 
      } else {  // x < 2^32
2978 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
2979 
        x_nr_bits =
2980 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2981 
      }
2982 
    } else {    // if x < 2^53
2983 
      tmp1.d = (double) C1.w[0]; // exact conversion
2984 
      x_nr_bits =
2985 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2986 
    }
2987 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
2988 
    tmp1.d = (double) C1.w[1];  // exact conversion
2989 
    x_nr_bits =
2990 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
2991 
  }
2992 
  q = nr_digits[x_nr_bits - 1].digits;
2993 
  if (q == 0) {
2994 
    q = nr_digits[x_nr_bits - 1].digits1;
2995 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
2996 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
2997 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
2998 
      q++;
2999 
  }
3000 
  exp = (x_exp >> 49) - 6176;
3001 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3002 
    // set invalid flag
3003 
    *pfpsf |= INVALID_EXCEPTION;
3004 
    // return Integer Indefinite
3005 
    res = 0x80000000;
3006 
    BID_RETURN (res);
3007 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
3008 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
3009 
    // so x rounded to an integer may or may not fit in a signed 32-bit int
3010 
    // the cases that do not fit are identified here; the ones that fit
3011 
    // fall through and will be handled with other cases further,
3012 
    // under '1 <= q + exp <= 10'
3013 
    if (x_sign) {       // if n < 0 and q + exp = 10
3014 
      // if n <= -1/2 then n is too large
3015 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1/2
3016 
      // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05, 1<=q<=34
3017 
      if (q <= 11) {
3018 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
3019 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3020 
        if (tmp64 >= 0x05ull) {
3021 
          // set invalid flag
3022 
          *pfpsf |= INVALID_EXCEPTION;
3023 
          // return Integer Indefinite
3024 
          res = 0x80000000;
3025 
          BID_RETURN (res);
3026 
        }
3027 
        // else cases that can be rounded to a 32-bit int fall through
3028 
        // to '1 <= q + exp <= 10'
3029 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3030 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05 <=>
3031 
        // C >= 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
3032 
        // (scale 1/2 up)
3033 
        tmp64 = 0x05ull;
3034 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3035 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
3036 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3037 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
3038 
        }
3039 
        if (C1.w[1] > C.w[1]
3040 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3041 
          // set invalid flag
3042 
          *pfpsf |= INVALID_EXCEPTION;
3043 
          // return Integer Indefinite
3044 
          res = 0x80000000;
3045 
          BID_RETURN (res);
3046 
        }
3047 
        // else cases that can be rounded to a 32-bit int fall through
3048 
        // to '1 <= q + exp <= 10'
3049 
      }
3050 
    } else {    // if n > 0 and q + exp = 10
3051 
      // if n >= 2^32 - 1/2 then n is too large
3052 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
3053 
      // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
3054 
      if (q <= 11) {
3055 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
3056 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3057 
        if (tmp64 >= 0x9fffffffbull) {
3058 
          // set invalid flag
3059 
          *pfpsf |= INVALID_EXCEPTION;
3060 
          // return Integer Indefinite
3061 
          res = 0x80000000;
3062 
          BID_RETURN (res);
3063 
        }
3064 
        // else cases that can be rounded to a 32-bit int fall through
3065 
        // to '1 <= q + exp <= 10'
3066 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3067 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
3068 
        // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
3069 
        // (scale 2^32-1/2 up)
3070 
        tmp64 = 0x9fffffffbull;
3071 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3072 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
3073 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3074 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
3075 
        }
3076 
        if (C1.w[1] > C.w[1]
3077 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3078 
          // set invalid flag
3079 
          *pfpsf |= INVALID_EXCEPTION;
3080 
          // return Integer Indefinite
3081 
          res = 0x80000000;
3082 
          BID_RETURN (res);
3083 
        }
3084 
        // else cases that can be rounded to a 32-bit int fall through
3085 
        // to '1 <= q + exp <= 10'
3086 
      }
3087 
    }
3088 
  }
3089 
  // n is not too large to be converted to int32: -1/2 < n < 2^32 - 1/2
3090 
  // Note: some of the cases tested for above fall through to this point
3091 
  if ((q + exp) < 0) {   // n = +/-0.0...c(0)c(1)...c(q-1)
3092 
    // return 0
3093 
    res = 0x00000000;
3094 
    BID_RETURN (res);
3095 
  } else if ((q + exp) == 0) {   // n = +/-0.c(0)c(1)...c(q-1)
3096 
    // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3097 
    //   res = 0
3098 
    // else if x > 0
3099 
    //   res = +1
3100 
    // else // if x < 0
3101 
    //   invalid exc
3102 
    ind = q - 1;
3103 
    if (ind <= 18) {    // 0 <= ind <= 18
3104 
      if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
3105 
        res = 0x00000000;       // return 0
3106 
      } else if (!x_sign) {     // n > 0
3107 
        res = 0x00000001;       // return +1
3108 
      } else {
3109 
        res = 0x80000000;
3110 
        *pfpsf |= INVALID_EXCEPTION;
3111 
      }
3112 
    } else {    // 19 <= ind <= 33
3113 
      if ((C1.w[1] < midpoint128[ind - 19].w[1])
3114 
          || ((C1.w[1] == midpoint128[ind - 19].w[1])
3115 
              && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
3116 
        res = 0x00000000;       // return 0
3117 
      } else if (!x_sign) {     // n > 0
3118 
        res = 0x00000001;       // return +1
3119 
      } else {
3120 
        res = 0x80000000;
3121 
        *pfpsf |= INVALID_EXCEPTION;
3122 
      }
3123 
    }
3124 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
3125 
    if (x_sign) {       // x <= -1
3126 
      // set invalid flag
3127 
      *pfpsf |= INVALID_EXCEPTION;
3128 
      // return Integer Indefinite
3129 
      res = 0x80000000;
3130 
      BID_RETURN (res);
3131 
    }
3132 
    // 1 <= x < 2^31-1/2 so x can be rounded
3133 
    // to nearest-away to a 32-bit signed integer
3134 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
3135 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
3136 
      // chop off ind digits from the lower part of C1
3137 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3138 
      tmp64 = C1.w[0];
3139 
      if (ind <= 19) {
3140 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
3141 
      } else {
3142 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
3143 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
3144 
      }
3145 
      if (C1.w[0] < tmp64)
3146 
        C1.w[1]++;
3147 
      // calculate C* and f*
3148 
      // C* is actually floor(C*) in this case
3149 
      // C* and f* need shifting and masking, as shown by
3150 
      // shiftright128[] and maskhigh128[]
3151 
      // 1 <= x <= 33
3152 
      // kx = 10^(-x) = ten2mk128[ind - 1]
3153 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3154 
      // the approximation of 10^(-x) was rounded up to 118 bits
3155 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
3156 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
3157 
        Cstar.w[1] = P256.w[3];
3158 
        Cstar.w[0] = P256.w[2];
3159 
      } else {  // 22 <= ind - 1 <= 33
3160 
        Cstar.w[1] = 0;
3161 
        Cstar.w[0] = P256.w[3];
3162 
      }
3163 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3164 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3165 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
3166 
      //   if floor(C*) is even then C* = floor(C*) - logical right
3167 
      //       shift; C* has p decimal digits, correct by Prop. 1)
3168 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
3169 
      //       shift; C* has p decimal digits, correct by Pr. 1)
3170 
      // else
3171 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
3172 
      //       correct by Property 1)
3173 
      // n = C* * 10^(e+x)
3174 
 
3175 
      // shift right C* by Ex-128 = shiftright128[ind]
3176 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
3177 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
3178 
        Cstar.w[0] =
3179 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
3180 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3181 
      } else {  // 22 <= ind - 1 <= 33
3182 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
3183 
      }
3184 
      // if the result was a midpoint, it was already rounded away from zero
3185 
      res = Cstar.w[0];  // always positive
3186 
      // no need to check for midpoints - already rounded away from zero!
3187 
    } else if (exp == 0) {
3188 
      // 1 <= q <= 10
3189 
      // res = +C (exact)
3190 
      res = C1.w[0];
3191 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3192 
      // res = +C * 10^exp (exact)
3193 
      res = C1.w[0] * ten2k64[exp];
3194 
    }
3195 
  }
3196 
}
3197 
 
3198 
BID_RETURN (res);
3199 
}
3200 
 
3201 
/*****************************************************************************
3202 
 *  BID128_to_uint32_xrninta
3203 
 ****************************************************************************/
3204 
 
3205 
BID128_FUNCTION_ARG1_NORND_CUSTOMRESTYPE (unsigned int,
3206 
                                          bid128_to_uint32_xrninta, x)
3207 
 
3208 
     unsigned int res;
3209 
     UINT64 x_sign;
3210 
     UINT64 x_exp;
3211 
     int exp;                   // unbiased exponent
3212 
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
3213 
     UINT64 tmp64, tmp64A;
3214 
     BID_UI64DOUBLE tmp1;
3215 
     unsigned int x_nr_bits;
3216 
     int q, ind, shift;
3217 
     UINT128 C1, C;
3218 
     UINT128 Cstar;             // C* represents up to 34 decimal digits ~ 113 bits
3219 
     UINT256 fstar;
3220 
     UINT256 P256;
3221 
     unsigned int tmp_inexact = 0;
3222 
 
3223 
  // unpack x
3224 
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
3225 
x_exp = x.w[1] & MASK_EXP;      // biased and shifted left 49 bit positions
3226 
C1.w[1] = x.w[1] & MASK_COEFF;
3227 
C1.w[0] = x.w[0];
3228 
 
3229 
  // check for NaN or Infinity
3230 
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
3231 
    // x is special
3232 
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
3233 
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
3234 
    // set invalid flag
3235 
    *pfpsf |= INVALID_EXCEPTION;
3236 
    // return Integer Indefinite
3237 
    res = 0x80000000;
3238 
  } else {      // x is QNaN
3239 
    // set invalid flag
3240 
    *pfpsf |= INVALID_EXCEPTION;
3241 
    // return Integer Indefinite
3242 
    res = 0x80000000;
3243 
  }
3244 
  BID_RETURN (res);
3245 
} else {        // x is not a NaN, so it must be infinity
3246 
  if (!x_sign) {        // x is +inf
3247 
    // set invalid flag
3248 
    *pfpsf |= INVALID_EXCEPTION;
3249 
    // return Integer Indefinite
3250 
    res = 0x80000000;
3251 
  } else {      // x is -inf
3252 
    // set invalid flag
3253 
    *pfpsf |= INVALID_EXCEPTION;
3254 
    // return Integer Indefinite
3255 
    res = 0x80000000;
3256 
  }
3257 
  BID_RETURN (res);
3258 
}
3259 
}
3260 
  // check for non-canonical values (after the check for special values)
3261 
if ((C1.w[1] > 0x0001ed09bead87c0ull)
3262 
    || (C1.w[1] == 0x0001ed09bead87c0ull
3263 
        && (C1.w[0] > 0x378d8e63ffffffffull))
3264 
    || ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull)) {
3265 
  res = 0x00000000;
3266 
  BID_RETURN (res);
3267 
} else if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
3268 
  // x is 0
3269 
  res = 0x00000000;
3270 
  BID_RETURN (res);
3271 
} else {        // x is not special and is not zero
3272 
 
3273 
  // q = nr. of decimal digits in x
3274 
  //  determine first the nr. of bits in x
3275 
  if (C1.w[1] == 0) {
3276 
    if (C1.w[0] >= 0x0020000000000000ull) {      // x >= 2^53
3277 
      // split the 64-bit value in two 32-bit halves to avoid rounding errors
3278 
      if (C1.w[0] >= 0x0000000100000000ull) {    // x >= 2^32
3279 
        tmp1.d = (double) (C1.w[0] >> 32);       // exact conversion
3280 
        x_nr_bits =
3281 
          33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3282 
      } else {  // x < 2^32
3283 
        tmp1.d = (double) (C1.w[0]);     // exact conversion
3284 
        x_nr_bits =
3285 
          1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3286 
      }
3287 
    } else {    // if x < 2^53
3288 
      tmp1.d = (double) C1.w[0]; // exact conversion
3289 
      x_nr_bits =
3290 
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3291 
    }
3292 
  } else {      // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
3293 
    tmp1.d = (double) C1.w[1];  // exact conversion
3294 
    x_nr_bits =
3295 
      65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
3296 
  }
3297 
  q = nr_digits[x_nr_bits - 1].digits;
3298 
  if (q == 0) {
3299 
    q = nr_digits[x_nr_bits - 1].digits1;
3300 
    if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
3301 
        || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi
3302 
            && C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
3303 
      q++;
3304 
  }
3305 
  exp = (x_exp >> 49) - 6176;
3306 
  if ((q + exp) > 10) { // x >= 10^10 ~= 2^33.2... (cannot fit in 32 bits)
3307 
    // set invalid flag
3308 
    *pfpsf |= INVALID_EXCEPTION;
3309 
    // return Integer Indefinite
3310 
    res = 0x80000000;
3311 
    BID_RETURN (res);
3312 
  } else if ((q + exp) == 10) { // x = c(0)c(1)...c(9).c(10)...c(q-1)
3313 
    // in this case 2^29.89... ~= 10^9 <= x < 10^10 ~= 2^33.2...
3314 
    // so x rounded to an integer may or may not fit in a signed 32-bit int
3315 
    // the cases that do not fit are identified here; the ones that fit
3316 
    // fall through and will be handled with other cases further,
3317 
    // under '1 <= q + exp <= 10'
3318 
    if (x_sign) {       // if n < 0 and q + exp = 10
3319 
      // if n <= -1/2 then n is too large
3320 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 1/2
3321 
      // <=> 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05, 1<=q<=34
3322 
      if (q <= 11) {
3323 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
3324 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3325 
        if (tmp64 >= 0x05ull) {
3326 
          // set invalid flag
3327 
          *pfpsf |= INVALID_EXCEPTION;
3328 
          // return Integer Indefinite
3329 
          res = 0x80000000;
3330 
          BID_RETURN (res);
3331 
        }
3332 
        // else cases that can be rounded to a 32-bit int fall through
3333 
        // to '1 <= q + exp <= 10'
3334 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3335 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x05 <=>
3336 
        // C >= 0x05 * 10^(q-11) where 1 <= q - 11 <= 23
3337 
        // (scale 1/2 up)
3338 
        tmp64 = 0x05ull;
3339 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3340 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
3341 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3342 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
3343 
        }
3344 
        if (C1.w[1] > C.w[1]
3345 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3346 
          // set invalid flag
3347 
          *pfpsf |= INVALID_EXCEPTION;
3348 
          // return Integer Indefinite
3349 
          res = 0x80000000;
3350 
          BID_RETURN (res);
3351 
        }
3352 
        // else cases that can be rounded to a 32-bit int fall through
3353 
        // to '1 <= q + exp <= 10'
3354 
      }
3355 
    } else {    // if n > 0 and q + exp = 10
3356 
      // if n >= 2^32 - 1/2 then n is too large
3357 
      // too large if c(0)c(1)...c(9).c(10)...c(q-1) >= 2^32-1/2
3358 
      // too large if 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb, 1<=q<=34
3359 
      if (q <= 11) {
3360 
        tmp64 = C1.w[0] * ten2k64[11 - q];       // C scaled up to 11-digit int
3361 
        // c(0)c(1)...c(9)c(10) or c(0)c(1)...c(q-1)0...0 (11 digits)
3362 
        if (tmp64 >= 0x9fffffffbull) {
3363 
          // set invalid flag
3364 
          *pfpsf |= INVALID_EXCEPTION;
3365 
          // return Integer Indefinite
3366 
          res = 0x80000000;
3367 
          BID_RETURN (res);
3368 
        }
3369 
        // else cases that can be rounded to a 32-bit int fall through
3370 
        // to '1 <= q + exp <= 10'
3371 
      } else {  // if (q > 11), i.e. 12 <= q <= 34 and so -24 <= exp <= -2
3372 
        // 0.c(0)c(1)...c(q-1) * 10^11 >= 0x9fffffffb <=>
3373 
        // C >= 0x9fffffffb * 10^(q-11) where 1 <= q - 11 <= 23
3374 
        // (scale 2^32-1/2 up)
3375 
        tmp64 = 0x9fffffffbull;
3376 
        if (q - 11 <= 19) {     // 1 <= q - 11 <= 19; 10^(q-11) requires 64 bits
3377 
          __mul_64x64_to_128MACH (C, tmp64, ten2k64[q - 11]);
3378 
        } else {        // 20 <= q - 11 <= 23, and 10^(q-11) requires 128 bits
3379 
          __mul_128x64_to_128 (C, tmp64, ten2k128[q - 31]);
3380 
        }
3381 
        if (C1.w[1] > C.w[1]
3382 
            || (C1.w[1] == C.w[1] && C1.w[0] >= C.w[0])) {
3383 
          // set invalid flag
3384 
          *pfpsf |= INVALID_EXCEPTION;
3385 
          // return Integer Indefinite
3386 
          res = 0x80000000;
3387 
          BID_RETURN (res);
3388 
        }
3389 
        // else cases that can be rounded to a 32-bit int fall through
3390 
        // to '1 <= q + exp <= 10'
3391 
      }
3392 
    }
3393 
  }
3394 
  // n is not too large to be converted to int32: -1/2 < n < 2^32 - 1/2
3395 
  // Note: some of the cases tested for above fall through to this point
3396 
  if ((q + exp) < 0) {   // n = +/-0.0...c(0)c(1)...c(q-1)
3397 
    // set inexact flag
3398 
    *pfpsf |= INEXACT_EXCEPTION;
3399 
    // return 0
3400 
    res = 0x00000000;
3401 
    BID_RETURN (res);
3402 
  } else if ((q + exp) == 0) {   // n = +/-0.c(0)c(1)...c(q-1)
3403 
    // if 0.c(0)c(1)...c(q-1) < 0.5 <=> c(0)c(1)...c(q-1) < 5 * 10^(q-1)
3404 
    //   res = 0
3405 
    // else if x > 0
3406 
    //   res = +1
3407 
    // else // if x < 0
3408 
    //   invalid exc
3409 
    ind = q - 1;
3410 
    if (ind <= 18) {    // 0 <= ind <= 18
3411 
      if ((C1.w[1] == 0) && (C1.w[0] < midpoint64[ind])) {
3412 
        res = 0x00000000;       // return 0
3413 
      } else if (!x_sign) {     // n > 0
3414 
        res = 0x00000001;       // return +1
3415 
      } else {
3416 
        res = 0x80000000;
3417 
        *pfpsf |= INVALID_EXCEPTION;
3418 
        BID_RETURN (res);
3419 
      }
3420 
    } else {    // 19 <= ind <= 33
3421 
      if ((C1.w[1] < midpoint128[ind - 19].w[1])
3422 
          || ((C1.w[1] == midpoint128[ind - 19].w[1])
3423 
              && (C1.w[0] < midpoint128[ind - 19].w[0]))) {
3424 
        res = 0x00000000;       // return 0
3425 
      } else if (!x_sign) {     // n > 0
3426 
        res = 0x00000001;       // return +1
3427 
      } else {
3428 
        res = 0x80000000;
3429 
        *pfpsf |= INVALID_EXCEPTION;
3430 
        BID_RETURN (res);
3431 
      }
3432 
    }
3433 
    // set inexact flag
3434 
    *pfpsf |= INEXACT_EXCEPTION;
3435 
  } else {      // if (1 <= q + exp <= 10, 1 <= q <= 34, -33 <= exp <= 9)
3436 
    if (x_sign) {       // x <= -1
3437 
      // set invalid flag
3438 
      *pfpsf |= INVALID_EXCEPTION;
3439 
      // return Integer Indefinite
3440 
      res = 0x80000000;
3441 
      BID_RETURN (res);
3442 
    }
3443 
    // 1 <= x < 2^31-1/2 so x can be rounded
3444 
    // to nearest-away to a 32-bit signed integer
3445 
    if (exp < 0) {       // 2 <= q <= 34, -33 <= exp <= -1, 1 <= q + exp <= 10
3446 
      ind = -exp;       // 1 <= ind <= 33; ind is a synonym for 'x'
3447 
      // chop off ind digits from the lower part of C1
3448 
      // C1 = C1 + 1/2 * 10^ind where the result C1 fits in 127 bits
3449 
      tmp64 = C1.w[0];
3450 
      if (ind <= 19) {
3451 
        C1.w[0] = C1.w[0] + midpoint64[ind - 1];
3452 
      } else {
3453 
        C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
3454 
        C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
3455 
      }
3456 
      if (C1.w[0] < tmp64)
3457 
        C1.w[1]++;
3458 
      // calculate C* and f*
3459 
      // C* is actually floor(C*) in this case
3460 
      // C* and f* need shifting and masking, as shown by
3461 
      // shiftright128[] and maskhigh128[]
3462 
      // 1 <= x <= 33
3463 
      // kx = 10^(-x) = ten2mk128[ind - 1]
3464 
      // C* = (C1 + 1/2 * 10^x) * 10^(-x)
3465 
      // the approximation of 10^(-x) was rounded up to 118 bits
3466 
      __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
3467 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
3468 
        Cstar.w[1] = P256.w[3];
3469 
        Cstar.w[0] = P256.w[2];
3470 
        fstar.w[3] = 0;
3471 
        fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
3472 
        fstar.w[1] = P256.w[1];
3473 
        fstar.w[0] = P256.w[0];
3474 
      } else {  // 22 <= ind - 1 <= 33
3475 
        Cstar.w[1] = 0;
3476 
        Cstar.w[0] = P256.w[3];
3477 
        fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
3478 
        fstar.w[2] = P256.w[2];
3479 
        fstar.w[1] = P256.w[1];
3480 
        fstar.w[0] = P256.w[0];
3481 
      }
3482 
      // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
3483 
      // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
3484 
      // if (0 < f* < 10^(-x)) then the result is a midpoint
3485 
      //   if floor(C*) is even then C* = floor(C*) - logical right
3486 
      //       shift; C* has p decimal digits, correct by Prop. 1)
3487 
      //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
3488 
      //       shift; C* has p decimal digits, correct by Pr. 1)
3489 
      // else
3490 
      //   C* = floor(C*) (logical right shift; C has p decimal digits,
3491 
      //       correct by Property 1)
3492 
      // n = C* * 10^(e+x)
3493 
 
3494 
      // shift right C* by Ex-128 = shiftright128[ind]
3495 
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
3496 
      if (ind - 1 <= 21) {      // 0 <= ind - 1 <= 21
3497 
        Cstar.w[0] =
3498 
          (Cstar.w[0] >> shift) | (Cstar.w[1] << (64 - shift));
3499 
        // redundant, it will be 0! Cstar.w[1] = (Cstar.w[1] >> shift);
3500 
      } else {  // 22 <= ind - 1 <= 33
3501 
        Cstar.w[0] = (Cstar.w[0] >> (shift - 64));        // 2 <= shift - 64 <= 38
3502 
      }
3503 
      // if the result was a midpoint, it was already rounded away from zero
3504 
      // determine inexactness of the rounding of C*
3505 
      // if (0 < f* - 1/2 < 10^(-x)) then
3506 
      //   the result is exact
3507 
      // else // if (f* - 1/2 > T*) then
3508 
      //   the result is inexact
3509 
      if (ind - 1 <= 2) {
3510 
        if (fstar.w[1] > 0x8000000000000000ull ||
3511 
            (fstar.w[1] == 0x8000000000000000ull
3512 
             && fstar.w[0] > 0x0ull)) {
3513 
          // f* > 1/2 and the result may be exact
3514 
          tmp64 = fstar.w[1] - 0x8000000000000000ull;   // f* - 1/2
3515 
          if (tmp64 > ten2mk128trunc[ind - 1].w[1]
3516 
              || (tmp64 == ten2mk128trunc[ind - 1].w[1]
3517 
                  && fstar.w[0] >= ten2mk128trunc[ind - 1].w[0])) {
3518 
            // set the inexact flag
3519 
            // *pfpsf |= INEXACT_EXCEPTION;
3520 
            tmp_inexact = 1;
3521 
          }     // else the result is exact
3522 
        } else {        // the result is inexact; f2* <= 1/2
3523 
          // set the inexact flag
3524 
          // *pfpsf |= INEXACT_EXCEPTION;
3525 
          tmp_inexact = 1;
3526 
        }
3527 
      } else if (ind - 1 <= 21) {       // if 3 <= ind <= 21
3528 
        if (fstar.w[3] > 0x0 ||
3529 
            (fstar.w[3] == 0x0 && fstar.w[2] > onehalf128[ind - 1]) ||
3530 
            (fstar.w[3] == 0x0 && fstar.w[2] == onehalf128[ind - 1] &&
3531 
             (fstar.w[1] || fstar.w[0]))) {
3532 
          // f2* > 1/2 and the result may be exact
3533 
          // Calculate f2* - 1/2
3534 
          tmp64 = fstar.w[2] - onehalf128[ind - 1];
3535 
          tmp64A = fstar.w[3];
3536 
          if (tmp64 > fstar.w[2])
3537 
            tmp64A--;
3538 
          if (tmp64A || tmp64
3539 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
3540 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
3541 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
3542 
            // set the inexact flag
3543 
            // *pfpsf |= INEXACT_EXCEPTION;
3544 
            tmp_inexact = 1;
3545 
          }     // else the result is exact
3546 
        } else {        // the result is inexact; f2* <= 1/2
3547 
          // set the inexact flag
3548 
          // *pfpsf |= INEXACT_EXCEPTION;
3549 
          tmp_inexact = 1;
3550 
        }
3551 
      } else {  // if 22 <= ind <= 33
3552 
        if (fstar.w[3] > onehalf128[ind - 1] ||
3553 
            (fstar.w[3] == onehalf128[ind - 1] &&
3554 
             (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
3555 
          // f2* > 1/2 and the result may be exact
3556 
          // Calculate f2* - 1/2
3557 
          tmp64 = fstar.w[3] - onehalf128[ind - 1];
3558 
          if (tmp64 || fstar.w[2]
3559 
              || fstar.w[1] > ten2mk128trunc[ind - 1].w[1]
3560 
              || (fstar.w[1] == ten2mk128trunc[ind - 1].w[1]
3561 
                  && fstar.w[0] > ten2mk128trunc[ind - 1].w[0])) {
3562 
            // set the inexact flag
3563 
            // *pfpsf |= INEXACT_EXCEPTION;
3564 
            tmp_inexact = 1;
3565 
          }     // else the result is exact
3566 
        } else {        // the result is inexact; f2* <= 1/2
3567 
          // set the inexact flag
3568 
          // *pfpsf |= INEXACT_EXCEPTION;
3569 
          tmp_inexact = 1;
3570 
        }
3571 
      }
3572 
      // no need to check for midpoints - already rounded away from zero!
3573 
      res = Cstar.w[0];  // the result is positive
3574 
      if (tmp_inexact)
3575 
        *pfpsf |= INEXACT_EXCEPTION;
3576 
    } else if (exp == 0) {
3577 
      // 1 <= q <= 10
3578 
      // res = +C (exact)
3579 
      res = C1.w[0];
3580 
    } else {    // if (exp > 0) => 1 <= exp <= 9, 1 <= q < 9, 2 <= q + exp <= 10
3581 
      // res = +C * 10^exp (exact)
3582 
      res = C1.w[0] * ten2k64[exp];
3583 
    }
3584 
  }
3585 
}
3586 
 
3587 
BID_RETURN (res);
3588 
}

powered by: WebSVN 2.1.0

© copyright 1999-2024 OpenCores.org, equivalent to Oliscience, all rights reserved. OpenCores®, registered trademark.