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

Subversion Repositories openrisc_me

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

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

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

powered by: WebSVN 2.1.0

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