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

Subversion Repositories openrisc

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

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

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

powered by: WebSVN 2.1.0

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