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

Subversion Repositories openrisc

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

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

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

powered by: WebSVN 2.1.0

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