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

Subversion Repositories openrisc

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

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

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

powered by: WebSVN 2.1.0

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