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

Subversion Repositories openrisc

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

Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

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