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

Subversion Repositories openrisc

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

Details | Compare with Previous | View Log

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

powered by: WebSVN 2.1.0

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