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

Subversion Repositories openrisc

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

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

powered by: WebSVN 2.1.0

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