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

Subversion Repositories openrisc_me

[/] [openrisc/] [trunk/] [gnu-src/] [gcc-4.5.1/] [libgcc/] [config/] [libbid/] [bid128_to_uint64.c] - Blame information for rev 280

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

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

powered by: WebSVN 2.1.0

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