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

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgcc/] [config/] [libbid/] [bid128_round_integral.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
#define BID_128RES
25
 
26
#include "bid_internal.h"
27
 
28
/*****************************************************************************
29
 *  BID128_round_integral_exact
30
 ****************************************************************************/
31
 
32
BID128_FUNCTION_ARG1 (bid128_round_integral_exact, x)
33
 
34
     UINT128 res = { {0xbaddbaddbaddbaddull, 0xbaddbaddbaddbaddull}
35
     };
36
UINT64 x_sign;
37
UINT64 x_exp;
38
int exp;                        // unbiased exponent
39
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
40
UINT64 tmp64;
41
BID_UI64DOUBLE tmp1;
42
unsigned int x_nr_bits;
43
int q, ind, shift;
44
UINT128 C1;
45
UINT256 fstar;
46
UINT256 P256;
47
 
48
  // check for NaN or Infinity
49
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
50
  // x is special
51
  if ((x.w[1] & MASK_NAN) == MASK_NAN) {        // x is NAN
52
    // if x = NaN, then res = Q (x)
53
    // check first for non-canonical NaN payload
54
    if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
55
        (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
56
         (x.w[0] > 0x38c15b09ffffffffull))) {
57
      x.w[1] = x.w[1] & 0xffffc00000000000ull;
58
      x.w[0] = 0x0ull;
59
    }
60
    if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {    // x is SNAN
61
      // set invalid flag
62
      *pfpsf |= INVALID_EXCEPTION;
63
      // return quiet (x)
64
      res.w[1] = x.w[1] & 0xfc003fffffffffffull;        // clear out also G[6]-G[16]
65
      res.w[0] = x.w[0];
66
    } else {    // x is QNaN
67
      // return x
68
      res.w[1] = x.w[1] & 0xfc003fffffffffffull;        // clear out G[6]-G[16]
69
      res.w[0] = x.w[0];
70
    }
71
    BID_RETURN (res)
72
  } else {      // x is not a NaN, so it must be infinity
73
    if ((x.w[1] & MASK_SIGN) == 0x0ull) {       // x is +inf
74
      // return +inf
75
      res.w[1] = 0x7800000000000000ull;
76
      res.w[0] = 0x0000000000000000ull;
77
    } else {    // x is -inf 
78
      // return -inf
79
      res.w[1] = 0xf800000000000000ull;
80
      res.w[0] = 0x0000000000000000ull;
81
    }
82
    BID_RETURN (res);
83
  }
84
}
85
  // unpack x
86
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
87
C1.w[1] = x.w[1] & MASK_COEFF;
88
C1.w[0] = x.w[0];
89
 
90
  // check for non-canonical values (treated as zero)
91
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
92
  // non-canonical
93
  x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
94
  C1.w[1] = 0;   // significand high
95
  C1.w[0] = 0;    // significand low
96
} else {        // G0_G1 != 11
97
  x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
98
  if (C1.w[1] > 0x0001ed09bead87c0ull ||
99
      (C1.w[1] == 0x0001ed09bead87c0ull
100
       && C1.w[0] > 0x378d8e63ffffffffull)) {
101
    // x is non-canonical if coefficient is larger than 10^34 -1
102
    C1.w[1] = 0;
103
    C1.w[0] = 0;
104
  } else {      // canonical
105
    ;
106
  }
107
}
108
 
109
  // test for input equal to zero
110
if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
111
  // x is 0
112
  // return 0 preserving the sign bit and the preferred exponent
113
  // of MAX(Q(x), 0)
114
  if (x_exp <= (0x1820ull << 49)) {
115
    res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
116
  } else {
117
    res.w[1] = x_sign | x_exp;
118
  }
119
  res.w[0] = 0x0000000000000000ull;
120
  BID_RETURN (res);
121
}
122
  // x is not special and is not zero
123
 
124
switch (rnd_mode) {
125
case ROUNDING_TO_NEAREST:
126
case ROUNDING_TIES_AWAY:
127
  // if (exp <= -(p+1)) return 0.0
128
  if (x_exp <= 0x2ffa000000000000ull) { // 0x2ffa000000000000ull == -35
129
    res.w[1] = x_sign | 0x3040000000000000ull;
130
    res.w[0] = 0x0000000000000000ull;
131
    *pfpsf |= INEXACT_EXCEPTION;
132
    BID_RETURN (res);
133
  }
134
  break;
135
case ROUNDING_DOWN:
136
  // if (exp <= -p) return -1.0 or +0.0
137
  if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffa000000000000ull == -34
138
    if (x_sign) {
139
      // if negative, return negative 1, because we know coefficient
140
      // is non-zero (would have been caught above)
141
      res.w[1] = 0xb040000000000000ull;
142
      res.w[0] = 0x0000000000000001ull;
143
    } else {
144
      // if positive, return positive 0, because we know coefficient is
145
      // non-zero (would have been caught above)
146
      res.w[1] = 0x3040000000000000ull;
147
      res.w[0] = 0x0000000000000000ull;
148
    }
149
    *pfpsf |= INEXACT_EXCEPTION;
150
    BID_RETURN (res);
151
  }
152
  break;
153
case ROUNDING_UP:
154
  // if (exp <= -p) return -0.0 or +1.0
155
  if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
156
    if (x_sign) {
157
      // if negative, return negative 0, because we know the coefficient
158
      // is non-zero (would have been caught above)
159
      res.w[1] = 0xb040000000000000ull;
160
      res.w[0] = 0x0000000000000000ull;
161
    } else {
162
      // if positive, return positive 1, because we know coefficient is
163
      // non-zero (would have been caught above)
164
      res.w[1] = 0x3040000000000000ull;
165
      res.w[0] = 0x0000000000000001ull;
166
    }
167
    *pfpsf |= INEXACT_EXCEPTION;
168
    BID_RETURN (res);
169
  }
170
  break;
171
case ROUNDING_TO_ZERO:
172
  // if (exp <= -p) return -0.0 or +0.0
173
  if (x_exp <= 0x2ffc000000000000ull) { // 0x2ffc000000000000ull == -34
174
    res.w[1] = x_sign | 0x3040000000000000ull;
175
    res.w[0] = 0x0000000000000000ull;
176
    *pfpsf |= INEXACT_EXCEPTION;
177
    BID_RETURN (res);
178
  }
179
  break;
180
}
181
 
182
  // q = nr. of decimal digits in x
183
  //  determine first the nr. of bits in x
184
if (C1.w[1] == 0) {
185
  if (C1.w[0] >= 0x0020000000000000ull) {        // x >= 2^53
186
    // split the 64-bit value in two 32-bit halves to avoid rounding errors
187
    if (C1.w[0] >= 0x0000000100000000ull) {      // x >= 2^32
188
      tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
189
      x_nr_bits =
190
        33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
191
    } else {    // x < 2^32
192
      tmp1.d = (double) (C1.w[0]);       // exact conversion
193
      x_nr_bits =
194
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
195
    }
196
  } else {      // if x < 2^53
197
    tmp1.d = (double) C1.w[0];   // exact conversion
198
    x_nr_bits =
199
      1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
200
  }
201
} else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
202
  tmp1.d = (double) C1.w[1];    // exact conversion
203
  x_nr_bits =
204
    65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
205
}
206
 
207
q = nr_digits[x_nr_bits - 1].digits;
208
if (q == 0) {
209
  q = nr_digits[x_nr_bits - 1].digits1;
210
  if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
211
      (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
212
       C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
213
    q++;
214
}
215
exp = (x_exp >> 49) - 6176;
216
if (exp >= 0) {  // -exp <= 0
217
  // the argument is an integer already
218
  res.w[1] = x.w[1];
219
  res.w[0] = x.w[0];
220
  BID_RETURN (res);
221
}
222
  // exp < 0
223
switch (rnd_mode) {
224
case ROUNDING_TO_NEAREST:
225
  if ((q + exp) >= 0) {  // exp < 0 and 1 <= -exp <= q
226
    // need to shift right -exp digits from the coefficient; exp will be 0
227
    ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
228
    // chop off ind digits from the lower part of C1 
229
    // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
230
    tmp64 = C1.w[0];
231
    if (ind <= 19) {
232
      C1.w[0] = C1.w[0] + midpoint64[ind - 1];
233
    } else {
234
      C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
235
      C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
236
    }
237
    if (C1.w[0] < tmp64)
238
      C1.w[1]++;
239
    // calculate C* and f*
240
    // C* is actually floor(C*) in this case
241
    // C* and f* need shifting and masking, as shown by
242
    // shiftright128[] and maskhigh128[]
243
    // 1 <= x <= 34
244
    // kx = 10^(-x) = ten2mk128[ind - 1]
245
    // C* = (C1 + 1/2 * 10^x) * 10^(-x)
246
    // the approximation of 10^(-x) was rounded up to 118 bits
247
    __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
248
    // determine the value of res and fstar
249
 
250
    // determine inexactness of the rounding of C*
251
    // if (0 < f* - 1/2 < 10^(-x)) then
252
    //   the result is exact
253
    // else // if (f* - 1/2 > T*) then
254
    //   the result is inexact
255
    // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[]
256
 
257
    if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
258
      // redundant shift = shiftright128[ind - 1]; // shift = 0
259
      res.w[1] = P256.w[3];
260
      res.w[0] = P256.w[2];
261
      // redundant fstar.w[3] = 0;
262
      // redundant fstar.w[2] = 0;
263
      fstar.w[1] = P256.w[1];
264
      fstar.w[0] = P256.w[0];
265
      // fraction f* < 10^(-x) <=> midpoint
266
      // f* is in the right position to be compared with
267
      // 10^(-x) from ten2mk128[]
268
      // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even)
269
      if ((res.w[0] & 0x0000000000000001ull) &&  // is result odd?
270
          ((fstar.w[1] < (ten2mk128[ind - 1].w[1]))
271
           || ((fstar.w[1] == ten2mk128[ind - 1].w[1])
272
               && (fstar.w[0] < ten2mk128[ind - 1].w[0])))) {
273
        // subract 1 to make even
274
        if (res.w[0]-- == 0) {
275
          res.w[1]--;
276
        }
277
      }
278
      if (fstar.w[1] > 0x8000000000000000ull ||
279
          (fstar.w[1] == 0x8000000000000000ull
280
           && fstar.w[0] > 0x0ull)) {
281
        // f* > 1/2 and the result may be exact
282
        tmp64 = fstar.w[1] - 0x8000000000000000ull;     // f* - 1/2
283
        if (tmp64 > ten2mk128[ind - 1].w[1] ||
284
            (tmp64 == ten2mk128[ind - 1].w[1] &&
285
             fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
286
          // set the inexact flag
287
          *pfpsf |= INEXACT_EXCEPTION;
288
        }       // else the result is exact 
289
      } else {  // the result is inexact; f2* <= 1/2  
290
        // set the inexact flag 
291
        *pfpsf |= INEXACT_EXCEPTION;
292
      }
293
    } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
294
      shift = shiftright128[ind - 1];   // 3 <= shift <= 63
295
      res.w[1] = (P256.w[3] >> shift);
296
      res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
297
      // redundant fstar.w[3] = 0;
298
      fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
299
      fstar.w[1] = P256.w[1];
300
      fstar.w[0] = P256.w[0];
301
      // fraction f* < 10^(-x) <=> midpoint
302
      // f* is in the right position to be compared with
303
      // 10^(-x) from ten2mk128[]
304
      if ((res.w[0] & 0x0000000000000001ull) &&  // is result odd?
305
          fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
306
                              (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
307
                               fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
308
        // subract 1 to make even
309
        if (res.w[0]-- == 0) {
310
          res.w[1]--;
311
        }
312
      }
313
      if (fstar.w[2] > onehalf128[ind - 1] ||
314
          (fstar.w[2] == onehalf128[ind - 1]
315
           && (fstar.w[1] || fstar.w[0]))) {
316
        // f2* > 1/2 and the result may be exact
317
        // Calculate f2* - 1/2
318
        tmp64 = fstar.w[2] - onehalf128[ind - 1];
319
        if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
320
            (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
321
             fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
322
          // set the inexact flag
323
          *pfpsf |= INEXACT_EXCEPTION;
324
        }       // else the result is exact
325
      } else {  // the result is inexact; f2* <= 1/2
326
        // set the inexact flag
327
        *pfpsf |= INEXACT_EXCEPTION;
328
      }
329
    } else {    // 22 <= ind - 1 <= 33
330
      shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
331
      res.w[1] = 0;
332
      res.w[0] = P256.w[3] >> shift;
333
      fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
334
      fstar.w[2] = P256.w[2];
335
      fstar.w[1] = P256.w[1];
336
      fstar.w[0] = P256.w[0];
337
      // fraction f* < 10^(-x) <=> midpoint
338
      // f* is in the right position to be compared with
339
      // 10^(-x) from ten2mk128[]
340
      if ((res.w[0] & 0x0000000000000001ull) &&  // is result odd?
341
          fstar.w[3] == 0 && fstar.w[2] == 0 &&
342
          (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
343
           (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
344
            fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
345
        // subract 1 to make even
346
        if (res.w[0]-- == 0) {
347
          res.w[1]--;
348
        }
349
      }
350
      if (fstar.w[3] > onehalf128[ind - 1] ||
351
          (fstar.w[3] == onehalf128[ind - 1] &&
352
           (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
353
        // f2* > 1/2 and the result may be exact
354
        // Calculate f2* - 1/2
355
        tmp64 = fstar.w[3] - onehalf128[ind - 1];
356
        if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1]
357
            || (fstar.w[1] == ten2mk128[ind - 1].w[1]
358
                && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
359
          // set the inexact flag
360
          *pfpsf |= INEXACT_EXCEPTION;
361
        }       // else the result is exact
362
      } else {  // the result is inexact; f2* <= 1/2
363
        // set the inexact flag
364
        *pfpsf |= INEXACT_EXCEPTION;
365
      }
366
    }
367
    res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
368
    BID_RETURN (res);
369
  } else {      // if ((q + exp) < 0) <=> q < -exp
370
    // the result is +0 or -0
371
    res.w[1] = x_sign | 0x3040000000000000ull;
372
    res.w[0] = 0x0000000000000000ull;
373
    *pfpsf |= INEXACT_EXCEPTION;
374
    BID_RETURN (res);
375
  }
376
  break;
377
case ROUNDING_TIES_AWAY:
378
  if ((q + exp) >= 0) {  // exp < 0 and 1 <= -exp <= q
379
    // need to shift right -exp digits from the coefficient; exp will be 0
380
    ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
381
    // chop off ind digits from the lower part of C1 
382
    // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
383
    tmp64 = C1.w[0];
384
    if (ind <= 19) {
385
      C1.w[0] = C1.w[0] + midpoint64[ind - 1];
386
    } else {
387
      C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
388
      C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
389
    }
390
    if (C1.w[0] < tmp64)
391
      C1.w[1]++;
392
    // calculate C* and f*
393
    // C* is actually floor(C*) in this case
394
    // C* and f* need shifting and masking, as shown by
395
    // shiftright128[] and maskhigh128[]
396
    // 1 <= x <= 34
397
    // kx = 10^(-x) = ten2mk128[ind - 1]
398
    // C* = (C1 + 1/2 * 10^x) * 10^(-x)
399
    // the approximation of 10^(-x) was rounded up to 118 bits
400
    __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
401
    // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
402
    // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
403
    // if (0 < f* < 10^(-x)) then the result is a midpoint
404
    //   if floor(C*) is even then C* = floor(C*) - logical right
405
    //       shift; C* has p decimal digits, correct by Prop. 1)
406
    //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
407
    //       shift; C* has p decimal digits, correct by Pr. 1)
408
    // else
409
    //   C* = floor(C*) (logical right shift; C has p decimal digits,
410
    //       correct by Property 1)
411
    // n = C* * 10^(e+x)
412
 
413
    // determine also the inexactness of the rounding of C*
414
    // if (0 < f* - 1/2 < 10^(-x)) then
415
    //   the result is exact
416
    // else // if (f* - 1/2 > T*) then
417
    //   the result is inexact
418
    // Note: we are going to use ten2mk128[] instead of ten2mk128trunc[]
419
    // shift right C* by Ex-128 = shiftright128[ind]
420
    if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
421
      // redundant shift = shiftright128[ind - 1]; // shift = 0
422
      res.w[1] = P256.w[3];
423
      res.w[0] = P256.w[2];
424
      // redundant fstar.w[3] = 0;
425
      // redundant fstar.w[2] = 0;
426
      fstar.w[1] = P256.w[1];
427
      fstar.w[0] = P256.w[0];
428
      if (fstar.w[1] > 0x8000000000000000ull ||
429
          (fstar.w[1] == 0x8000000000000000ull
430
           && fstar.w[0] > 0x0ull)) {
431
        // f* > 1/2 and the result may be exact
432
        tmp64 = fstar.w[1] - 0x8000000000000000ull;     // f* - 1/2
433
        if ((tmp64 > ten2mk128[ind - 1].w[1] ||
434
             (tmp64 == ten2mk128[ind - 1].w[1] &&
435
              fstar.w[0] >= ten2mk128[ind - 1].w[0]))) {
436
          // set the inexact flag
437
          *pfpsf |= INEXACT_EXCEPTION;
438
        }       // else the result is exact
439
      } else {  // the result is inexact; f2* <= 1/2
440
        // set the inexact flag
441
        *pfpsf |= INEXACT_EXCEPTION;
442
      }
443
    } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
444
      shift = shiftright128[ind - 1];   // 3 <= shift <= 63
445
      res.w[1] = (P256.w[3] >> shift);
446
      res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
447
      // redundant fstar.w[3] = 0;
448
      fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
449
      fstar.w[1] = P256.w[1];
450
      fstar.w[0] = P256.w[0];
451
      if (fstar.w[2] > onehalf128[ind - 1] ||
452
          (fstar.w[2] == onehalf128[ind - 1]
453
           && (fstar.w[1] || fstar.w[0]))) {
454
        // f2* > 1/2 and the result may be exact
455
        // Calculate f2* - 1/2
456
        tmp64 = fstar.w[2] - onehalf128[ind - 1];
457
        if (tmp64 || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
458
            (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
459
             fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
460
          // set the inexact flag
461
          *pfpsf |= INEXACT_EXCEPTION;
462
        }       // else the result is exact
463
      } else {  // the result is inexact; f2* <= 1/2
464
        // set the inexact flag
465
        *pfpsf |= INEXACT_EXCEPTION;
466
      }
467
    } else {    // 22 <= ind - 1 <= 33
468
      shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
469
      res.w[1] = 0;
470
      res.w[0] = P256.w[3] >> shift;
471
      fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
472
      fstar.w[2] = P256.w[2];
473
      fstar.w[1] = P256.w[1];
474
      fstar.w[0] = P256.w[0];
475
      if (fstar.w[3] > onehalf128[ind - 1] ||
476
          (fstar.w[3] == onehalf128[ind - 1] &&
477
           (fstar.w[2] || fstar.w[1] || fstar.w[0]))) {
478
        // f2* > 1/2 and the result may be exact
479
        // Calculate f2* - 1/2
480
        tmp64 = fstar.w[3] - onehalf128[ind - 1];
481
        if (tmp64 || fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1]
482
            || (fstar.w[1] == ten2mk128[ind - 1].w[1]
483
                && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
484
          // set the inexact flag
485
          *pfpsf |= INEXACT_EXCEPTION;
486
        }       // else the result is exact
487
      } else {  // the result is inexact; f2* <= 1/2
488
        // set the inexact flag
489
        *pfpsf |= INEXACT_EXCEPTION;
490
      }
491
    }
492
    // if the result was a midpoint, it was already rounded away from zero
493
    res.w[1] |= x_sign | 0x3040000000000000ull;
494
    BID_RETURN (res);
495
  } else {      // if ((q + exp) < 0) <=> q < -exp
496
    // the result is +0 or -0
497
    res.w[1] = x_sign | 0x3040000000000000ull;
498
    res.w[0] = 0x0000000000000000ull;
499
    *pfpsf |= INEXACT_EXCEPTION;
500
    BID_RETURN (res);
501
  }
502
  break;
503
case ROUNDING_DOWN:
504
  if ((q + exp) > 0) {   // exp < 0 and 1 <= -exp < q
505
    // need to shift right -exp digits from the coefficient; exp will be 0
506
    ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 
507
    // (number of digits to be chopped off)
508
    // chop off ind digits from the lower part of C1 
509
    // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
510
    // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
511
    // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
512
    // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
513
    // tmp64 = C1.w[0];
514
    // if (ind <= 19) {
515
    //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
516
    // } else {
517
    //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
518
    //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
519
    // }
520
    // if (C1.w[0] < tmp64) C1.w[1]++;
521
    // if carry-out from C1.w[0], increment C1.w[1]
522
    // calculate C* and f*
523
    // C* is actually floor(C*) in this case
524
    // C* and f* need shifting and masking, as shown by
525
    // shiftright128[] and maskhigh128[]
526
    // 1 <= x <= 34
527
    // kx = 10^(-x) = ten2mk128[ind - 1]
528
    // C* = (C1 + 1/2 * 10^x) * 10^(-x)
529
    // the approximation of 10^(-x) was rounded up to 118 bits
530
    __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
531
    if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
532
      res.w[1] = P256.w[3];
533
      res.w[0] = P256.w[2];
534
      // redundant fstar.w[3] = 0;
535
      // redundant fstar.w[2] = 0;
536
      // redundant fstar.w[1] = P256.w[1];
537
      // redundant fstar.w[0] = P256.w[0];
538
      // fraction f* > 10^(-x) <=> inexact
539
      // f* is in the right position to be compared with
540
      // 10^(-x) from ten2mk128[]
541
      if ((P256.w[1] > ten2mk128[ind - 1].w[1])
542
          || (P256.w[1] == ten2mk128[ind - 1].w[1]
543
              && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
544
        *pfpsf |= INEXACT_EXCEPTION;
545
        // if positive, the truncated value is already the correct result
546
        if (x_sign) {   // if negative
547
          if (++res.w[0] == 0) {
548
            res.w[1]++;
549
          }
550
        }
551
      }
552
    } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
553
      shift = shiftright128[ind - 1];   // 0 <= shift <= 102
554
      res.w[1] = (P256.w[3] >> shift);
555
      res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
556
      // redundant fstar.w[3] = 0;
557
      fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
558
      fstar.w[1] = P256.w[1];
559
      fstar.w[0] = P256.w[0];
560
      // fraction f* > 10^(-x) <=> inexact
561
      // f* is in the right position to be compared with
562
      // 10^(-x) from ten2mk128[]
563
      if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
564
          (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
565
           fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
566
        *pfpsf |= INEXACT_EXCEPTION;
567
        // if positive, the truncated value is already the correct result
568
        if (x_sign) {   // if negative
569
          if (++res.w[0] == 0) {
570
            res.w[1]++;
571
          }
572
        }
573
      }
574
    } else {    // 22 <= ind - 1 <= 33
575
      shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
576
      res.w[1] = 0;
577
      res.w[0] = P256.w[3] >> shift;
578
      fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
579
      fstar.w[2] = P256.w[2];
580
      fstar.w[1] = P256.w[1];
581
      fstar.w[0] = P256.w[0];
582
      // fraction f* > 10^(-x) <=> inexact
583
      // f* is in the right position to be compared with
584
      // 10^(-x) from ten2mk128[]
585
      if (fstar.w[3] || fstar.w[2]
586
          || fstar.w[1] > ten2mk128[ind - 1].w[1]
587
          || (fstar.w[1] == ten2mk128[ind - 1].w[1]
588
              && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
589
        *pfpsf |= INEXACT_EXCEPTION;
590
        // if positive, the truncated value is already the correct result
591
        if (x_sign) {   // if negative
592
          if (++res.w[0] == 0) {
593
            res.w[1]++;
594
          }
595
        }
596
      }
597
    }
598
    res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
599
    BID_RETURN (res);
600
  } else {      // if exp < 0 and q + exp <= 0
601
    if (x_sign) {       // negative rounds down to -1.0
602
      res.w[1] = 0xb040000000000000ull;
603
      res.w[0] = 0x0000000000000001ull;
604
    } else {    // positive rpunds down to +0.0
605
      res.w[1] = 0x3040000000000000ull;
606
      res.w[0] = 0x0000000000000000ull;
607
    }
608
    *pfpsf |= INEXACT_EXCEPTION;
609
    BID_RETURN (res);
610
  }
611
  break;
612
case ROUNDING_UP:
613
  if ((q + exp) > 0) {   // exp < 0 and 1 <= -exp < q
614
    // need to shift right -exp digits from the coefficient; exp will be 0
615
    ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x' 
616
    // (number of digits to be chopped off)
617
    // chop off ind digits from the lower part of C1 
618
    // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
619
    // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
620
    // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
621
    // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
622
    // tmp64 = C1.w[0];
623
    // if (ind <= 19) {
624
    //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
625
    // } else {
626
    //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
627
    //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
628
    // }
629
    // if (C1.w[0] < tmp64) C1.w[1]++;  
630
    // if carry-out from C1.w[0], increment C1.w[1]
631
    // calculate C* and f*
632
    // C* is actually floor(C*) in this case
633
    // C* and f* need shifting and masking, as shown by
634
    // shiftright128[] and maskhigh128[]
635
    // 1 <= x <= 34
636
    // kx = 10^(-x) = ten2mk128[ind - 1]
637
    // C* = C1 * 10^(-x)
638
    // the approximation of 10^(-x) was rounded up to 118 bits
639
    __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
640
    if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
641
      res.w[1] = P256.w[3];
642
      res.w[0] = P256.w[2];
643
      // redundant fstar.w[3] = 0;
644
      // redundant fstar.w[2] = 0;
645
      // redundant fstar.w[1] = P256.w[1]; 
646
      // redundant fstar.w[0] = P256.w[0];
647
      // fraction f* > 10^(-x) <=> inexact
648
      // f* is in the right position to be compared with 
649
      // 10^(-x) from ten2mk128[]
650
      if ((P256.w[1] > ten2mk128[ind - 1].w[1])
651
          || (P256.w[1] == ten2mk128[ind - 1].w[1]
652
              && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
653
        *pfpsf |= INEXACT_EXCEPTION;
654
        // if negative, the truncated value is already the correct result
655
        if (!x_sign) {  // if positive
656
          if (++res.w[0] == 0) {
657
            res.w[1]++;
658
          }
659
        }
660
      }
661
    } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
662
      shift = shiftright128[ind - 1];   // 3 <= shift <= 63
663
      res.w[1] = (P256.w[3] >> shift);
664
      res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
665
      // redundant fstar.w[3] = 0;
666
      fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
667
      fstar.w[1] = P256.w[1];
668
      fstar.w[0] = P256.w[0];
669
      // fraction f* > 10^(-x) <=> inexact
670
      // f* is in the right position to be compared with 
671
      // 10^(-x) from ten2mk128[]
672
      if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
673
          (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
674
           fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
675
        *pfpsf |= INEXACT_EXCEPTION;
676
        // if negative, the truncated value is already the correct result
677
        if (!x_sign) {  // if positive
678
          if (++res.w[0] == 0) {
679
            res.w[1]++;
680
          }
681
        }
682
      }
683
    } else {    // 22 <= ind - 1 <= 33
684
      shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
685
      res.w[1] = 0;
686
      res.w[0] = P256.w[3] >> shift;
687
      fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
688
      fstar.w[2] = P256.w[2];
689
      fstar.w[1] = P256.w[1];
690
      fstar.w[0] = P256.w[0];
691
      // fraction f* > 10^(-x) <=> inexact
692
      // f* is in the right position to be compared with 
693
      // 10^(-x) from ten2mk128[]
694
      if (fstar.w[3] || fstar.w[2]
695
          || fstar.w[1] > ten2mk128[ind - 1].w[1]
696
          || (fstar.w[1] == ten2mk128[ind - 1].w[1]
697
              && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
698
        *pfpsf |= INEXACT_EXCEPTION;
699
        // if negative, the truncated value is already the correct result
700
        if (!x_sign) {  // if positive
701
          if (++res.w[0] == 0) {
702
            res.w[1]++;
703
          }
704
        }
705
      }
706
    }
707
    res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
708
    BID_RETURN (res);
709
  } else {      // if exp < 0 and q + exp <= 0
710
    if (x_sign) {       // negative rounds up to -0.0
711
      res.w[1] = 0xb040000000000000ull;
712
      res.w[0] = 0x0000000000000000ull;
713
    } else {    // positive rpunds up to +1.0
714
      res.w[1] = 0x3040000000000000ull;
715
      res.w[0] = 0x0000000000000001ull;
716
    }
717
    *pfpsf |= INEXACT_EXCEPTION;
718
    BID_RETURN (res);
719
  }
720
  break;
721
case ROUNDING_TO_ZERO:
722
  if ((q + exp) > 0) {   // exp < 0 and 1 <= -exp < q
723
    // need to shift right -exp digits from the coefficient; exp will be 0
724
    ind = -exp; // 1 <= ind <= 34; ind is a synonym for 'x'
725
    // (number of digits to be chopped off)
726
    // chop off ind digits from the lower part of C1 
727
    // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
728
    // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
729
    // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
730
    // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
731
    //tmp64 = C1.w[0];
732
    // if (ind <= 19) {
733
    //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
734
    // } else {
735
    //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
736
    //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
737
    // }
738
    // if (C1.w[0] < tmp64) C1.w[1]++;  
739
    // if carry-out from C1.w[0], increment C1.w[1]
740
    // calculate C* and f*
741
    // C* is actually floor(C*) in this case
742
    // C* and f* need shifting and masking, as shown by
743
    // shiftright128[] and maskhigh128[]
744
    // 1 <= x <= 34
745
    // kx = 10^(-x) = ten2mk128[ind - 1]
746
    // C* = (C1 + 1/2 * 10^x) * 10^(-x)
747
    // the approximation of 10^(-x) was rounded up to 118 bits
748
    __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
749
    if (ind - 1 <= 2) { // 0 <= ind - 1 <= 2 => shift = 0
750
      res.w[1] = P256.w[3];
751
      res.w[0] = P256.w[2];
752
      // redundant fstar.w[3] = 0;
753
      // redundant fstar.w[2] = 0;
754
      // redundant fstar.w[1] = P256.w[1]; 
755
      // redundant fstar.w[0] = P256.w[0];
756
      // fraction f* > 10^(-x) <=> inexact
757
      // f* is in the right position to be compared with 
758
      // 10^(-x) from ten2mk128[]
759
      if ((P256.w[1] > ten2mk128[ind - 1].w[1])
760
          || (P256.w[1] == ten2mk128[ind - 1].w[1]
761
              && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
762
        *pfpsf |= INEXACT_EXCEPTION;
763
      }
764
    } else if (ind - 1 <= 21) { // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
765
      shift = shiftright128[ind - 1];   // 3 <= shift <= 63
766
      res.w[1] = (P256.w[3] >> shift);
767
      res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
768
      // redundant fstar.w[3] = 0;
769
      fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
770
      fstar.w[1] = P256.w[1];
771
      fstar.w[0] = P256.w[0];
772
      // fraction f* > 10^(-x) <=> inexact
773
      // f* is in the right position to be compared with 
774
      // 10^(-x) from ten2mk128[]
775
      if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
776
          (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
777
           fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
778
        *pfpsf |= INEXACT_EXCEPTION;
779
      }
780
    } else {    // 22 <= ind - 1 <= 33
781
      shift = shiftright128[ind - 1] - 64;      // 2 <= shift <= 38
782
      res.w[1] = 0;
783
      res.w[0] = P256.w[3] >> shift;
784
      fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
785
      fstar.w[2] = P256.w[2];
786
      fstar.w[1] = P256.w[1];
787
      fstar.w[0] = P256.w[0];
788
      // fraction f* > 10^(-x) <=> inexact
789
      // f* is in the right position to be compared with 
790
      // 10^(-x) from ten2mk128[]
791
      if (fstar.w[3] || fstar.w[2]
792
          || fstar.w[1] > ten2mk128[ind - 1].w[1]
793
          || (fstar.w[1] == ten2mk128[ind - 1].w[1]
794
              && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
795
        *pfpsf |= INEXACT_EXCEPTION;
796
      }
797
    }
798
    res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
799
    BID_RETURN (res);
800
  } else {      // if exp < 0 and q + exp <= 0 the result is +0 or -0
801
    res.w[1] = x_sign | 0x3040000000000000ull;
802
    res.w[0] = 0x0000000000000000ull;
803
    *pfpsf |= INEXACT_EXCEPTION;
804
    BID_RETURN (res);
805
  }
806
  break;
807
}
808
 
809
BID_RETURN (res);
810
}
811
 
812
/*****************************************************************************
813
 *  BID128_round_integral_nearest_even
814
 ****************************************************************************/
815
 
816
BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_even, x)
817
 
818
     UINT128 res;
819
     UINT64 x_sign;
820
     UINT64 x_exp;
821
     int exp;                   // unbiased exponent
822
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo (all are UINT64)
823
     UINT64 tmp64;
824
     BID_UI64DOUBLE tmp1;
825
     unsigned int x_nr_bits;
826
     int q, ind, shift;
827
     UINT128 C1;
828
  // UINT128 res is C* at first - represents up to 34 decimal digits ~ 113 bits
829
     UINT256 fstar;
830
     UINT256 P256;
831
 
832
  // check for NaN or Infinity
833
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
834
    // x is special
835
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
836
  // if x = NaN, then res = Q (x)
837
  // check first for non-canonical NaN payload
838
  if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
839
      (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
840
       (x.w[0] > 0x38c15b09ffffffffull))) {
841
    x.w[1] = x.w[1] & 0xffffc00000000000ull;
842
    x.w[0] = 0x0ull;
843
  }
844
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
845
    // set invalid flag
846
    *pfpsf |= INVALID_EXCEPTION;
847
    // return quiet (x)
848
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
849
    res.w[0] = x.w[0];
850
  } else {      // x is QNaN
851
    // return x
852
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
853
    res.w[0] = x.w[0];
854
  }
855
  BID_RETURN (res)
856
} else {        // x is not a NaN, so it must be infinity
857
  if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
858
    // return +inf
859
    res.w[1] = 0x7800000000000000ull;
860
    res.w[0] = 0x0000000000000000ull;
861
  } else {      // x is -inf 
862
    // return -inf
863
    res.w[1] = 0xf800000000000000ull;
864
    res.w[0] = 0x0000000000000000ull;
865
  }
866
  BID_RETURN (res);
867
}
868
}
869
  // unpack x
870
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
871
C1.w[1] = x.w[1] & MASK_COEFF;
872
C1.w[0] = x.w[0];
873
 
874
  // check for non-canonical values (treated as zero)
875
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
876
  // non-canonical
877
  x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
878
  C1.w[1] = 0;   // significand high
879
  C1.w[0] = 0;    // significand low
880
} else {        // G0_G1 != 11
881
  x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
882
  if (C1.w[1] > 0x0001ed09bead87c0ull ||
883
      (C1.w[1] == 0x0001ed09bead87c0ull
884
       && C1.w[0] > 0x378d8e63ffffffffull)) {
885
    // x is non-canonical if coefficient is larger than 10^34 -1
886
    C1.w[1] = 0;
887
    C1.w[0] = 0;
888
  } else {      // canonical
889
    ;
890
  }
891
}
892
 
893
  // test for input equal to zero
894
if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
895
  // x is 0
896
  // return 0 preserving the sign bit and the preferred exponent
897
  // of MAX(Q(x), 0)
898
  if (x_exp <= (0x1820ull << 49)) {
899
    res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
900
  } else {
901
    res.w[1] = x_sign | x_exp;
902
  }
903
  res.w[0] = 0x0000000000000000ull;
904
  BID_RETURN (res);
905
}
906
  // x is not special and is not zero
907
 
908
  // if (exp <= -(p+1)) return 0
909
if (x_exp <= 0x2ffa000000000000ull) {   // 0x2ffa000000000000ull == -35
910
  res.w[1] = x_sign | 0x3040000000000000ull;
911
  res.w[0] = 0x0000000000000000ull;
912
  BID_RETURN (res);
913
}
914
  // q = nr. of decimal digits in x
915
  //  determine first the nr. of bits in x
916
if (C1.w[1] == 0) {
917
  if (C1.w[0] >= 0x0020000000000000ull) {        // x >= 2^53
918
    // split the 64-bit value in two 32-bit halves to avoid rounding errors
919
    if (C1.w[0] >= 0x0000000100000000ull) {      // x >= 2^32
920
      tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
921
      x_nr_bits =
922
        33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
923
    } else {    // x < 2^32
924
      tmp1.d = (double) (C1.w[0]);       // exact conversion
925
      x_nr_bits =
926
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
927
    }
928
  } else {      // if x < 2^53
929
    tmp1.d = (double) C1.w[0];   // exact conversion
930
    x_nr_bits =
931
      1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
932
  }
933
} else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
934
  tmp1.d = (double) C1.w[1];    // exact conversion
935
  x_nr_bits =
936
    65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
937
}
938
 
939
q = nr_digits[x_nr_bits - 1].digits;
940
if (q == 0) {
941
  q = nr_digits[x_nr_bits - 1].digits1;
942
  if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi
943
      || (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
944
          C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
945
    q++;
946
}
947
exp = (x_exp >> 49) - 6176;
948
if (exp >= 0) {  // -exp <= 0
949
  // the argument is an integer already
950
  res.w[1] = x.w[1];
951
  res.w[0] = x.w[0];
952
  BID_RETURN (res);
953
} else if ((q + exp) >= 0) {     // exp < 0 and 1 <= -exp <= q
954
  // need to shift right -exp digits from the coefficient; the exp will be 0
955
  ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x'
956
  // chop off ind digits from the lower part of C1 
957
  // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
958
  tmp64 = C1.w[0];
959
  if (ind <= 19) {
960
    C1.w[0] = C1.w[0] + midpoint64[ind - 1];
961
  } else {
962
    C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
963
    C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
964
  }
965
  if (C1.w[0] < tmp64)
966
    C1.w[1]++;
967
  // calculate C* and f*
968
  // C* is actually floor(C*) in this case
969
  // C* and f* need shifting and masking, as shown by
970
  // shiftright128[] and maskhigh128[]
971
  // 1 <= x <= 34
972
  // kx = 10^(-x) = ten2mk128[ind - 1]
973
  // C* = (C1 + 1/2 * 10^x) * 10^(-x)
974
  // the approximation of 10^(-x) was rounded up to 118 bits
975
  __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
976
  // determine the value of res and fstar
977
  if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
978
    // redundant shift = shiftright128[ind - 1]; // shift = 0
979
    res.w[1] = P256.w[3];
980
    res.w[0] = P256.w[2];
981
    // redundant fstar.w[3] = 0;
982
    // redundant fstar.w[2] = 0;
983
    // redundant fstar.w[1] = P256.w[1];
984
    // redundant fstar.w[0] = P256.w[0];
985
    // fraction f* < 10^(-x) <=> midpoint
986
    // f* is in the right position to be compared with
987
    // 10^(-x) from ten2mk128[]
988
    // if 0 < fstar < 10^(-x), subtract 1 if odd (for rounding to even)
989
    if ((res.w[0] & 0x0000000000000001ull) &&    // is result odd?
990
        ((P256.w[1] < (ten2mk128[ind - 1].w[1]))
991
         || ((P256.w[1] == ten2mk128[ind - 1].w[1])
992
             && (P256.w[0] < ten2mk128[ind - 1].w[0])))) {
993
      // subract 1 to make even
994
      if (res.w[0]-- == 0) {
995
        res.w[1]--;
996
      }
997
    }
998
  } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
999
    shift = shiftright128[ind - 1];     // 3 <= shift <= 63
1000
    res.w[1] = (P256.w[3] >> shift);
1001
    res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1002
    // redundant fstar.w[3] = 0;
1003
    fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1004
    fstar.w[1] = P256.w[1];
1005
    fstar.w[0] = P256.w[0];
1006
    // fraction f* < 10^(-x) <=> midpoint
1007
    // f* is in the right position to be compared with
1008
    // 10^(-x) from ten2mk128[]
1009
    if ((res.w[0] & 0x0000000000000001ull) &&    // is result odd?
1010
        fstar.w[2] == 0 && (fstar.w[1] < ten2mk128[ind - 1].w[1] ||
1011
                            (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1012
                             fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
1013
      // subract 1 to make even
1014
      if (res.w[0]-- == 0) {
1015
        res.w[1]--;
1016
      }
1017
    }
1018
  } else {      // 22 <= ind - 1 <= 33
1019
    shift = shiftright128[ind - 1] - 64;        // 2 <= shift <= 38
1020
    res.w[1] = 0;
1021
    res.w[0] = P256.w[3] >> shift;
1022
    fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1023
    fstar.w[2] = P256.w[2];
1024
    fstar.w[1] = P256.w[1];
1025
    fstar.w[0] = P256.w[0];
1026
    // fraction f* < 10^(-x) <=> midpoint
1027
    // f* is in the right position to be compared with
1028
    // 10^(-x) from ten2mk128[]
1029
    if ((res.w[0] & 0x0000000000000001ull) &&    // is result odd?
1030
        fstar.w[3] == 0 && fstar.w[2] == 0
1031
        && (fstar.w[1] < ten2mk128[ind - 1].w[1]
1032
            || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1033
                && fstar.w[0] < ten2mk128[ind - 1].w[0]))) {
1034
      // subract 1 to make even
1035
      if (res.w[0]-- == 0) {
1036
        res.w[1]--;
1037
      }
1038
    }
1039
  }
1040
  res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1041
  BID_RETURN (res);
1042
} else {        // if ((q + exp) < 0) <=> q < -exp
1043
  // the result is +0 or -0
1044
  res.w[1] = x_sign | 0x3040000000000000ull;
1045
  res.w[0] = 0x0000000000000000ull;
1046
  BID_RETURN (res);
1047
}
1048
}
1049
 
1050
/*****************************************************************************
1051
 *  BID128_round_integral_negative
1052
 ****************************************************************************/
1053
 
1054
BID128_FUNCTION_ARG1_NORND (bid128_round_integral_negative, x)
1055
 
1056
     UINT128 res;
1057
     UINT64 x_sign;
1058
     UINT64 x_exp;
1059
     int exp;                   // unbiased exponent
1060
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo 
1061
  // (all are UINT64)
1062
     BID_UI64DOUBLE tmp1;
1063
     unsigned int x_nr_bits;
1064
     int q, ind, shift;
1065
     UINT128 C1;
1066
  // UINT128 res is C* at first - represents up to 34 decimal digits ~ 
1067
  // 113 bits
1068
     UINT256 fstar;
1069
     UINT256 P256;
1070
 
1071
  // check for NaN or Infinity
1072
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1073
    // x is special
1074
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1075
  // if x = NaN, then res = Q (x)
1076
  // check first for non-canonical NaN payload
1077
  if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1078
      (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1079
       (x.w[0] > 0x38c15b09ffffffffull))) {
1080
    x.w[1] = x.w[1] & 0xffffc00000000000ull;
1081
    x.w[0] = 0x0ull;
1082
  }
1083
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1084
    // set invalid flag
1085
    *pfpsf |= INVALID_EXCEPTION;
1086
    // return quiet (x)
1087
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
1088
    res.w[0] = x.w[0];
1089
  } else {      // x is QNaN
1090
    // return x
1091
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
1092
    res.w[0] = x.w[0];
1093
  }
1094
  BID_RETURN (res)
1095
} else {        // x is not a NaN, so it must be infinity
1096
  if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1097
    // return +inf
1098
    res.w[1] = 0x7800000000000000ull;
1099
    res.w[0] = 0x0000000000000000ull;
1100
  } else {      // x is -inf 
1101
    // return -inf
1102
    res.w[1] = 0xf800000000000000ull;
1103
    res.w[0] = 0x0000000000000000ull;
1104
  }
1105
  BID_RETURN (res);
1106
}
1107
}
1108
  // unpack x
1109
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1110
C1.w[1] = x.w[1] & MASK_COEFF;
1111
C1.w[0] = x.w[0];
1112
 
1113
  // check for non-canonical values (treated as zero)
1114
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
1115
  // non-canonical
1116
  x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
1117
  C1.w[1] = 0;   // significand high
1118
  C1.w[0] = 0;    // significand low
1119
} else {        // G0_G1 != 11
1120
  x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
1121
  if (C1.w[1] > 0x0001ed09bead87c0ull ||
1122
      (C1.w[1] == 0x0001ed09bead87c0ull
1123
       && C1.w[0] > 0x378d8e63ffffffffull)) {
1124
    // x is non-canonical if coefficient is larger than 10^34 -1
1125
    C1.w[1] = 0;
1126
    C1.w[0] = 0;
1127
  } else {      // canonical
1128
    ;
1129
  }
1130
}
1131
 
1132
  // test for input equal to zero
1133
if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1134
  // x is 0
1135
  // return 0 preserving the sign bit and the preferred exponent
1136
  // of MAX(Q(x), 0)
1137
  if (x_exp <= (0x1820ull << 49)) {
1138
    res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1139
  } else {
1140
    res.w[1] = x_sign | x_exp;
1141
  }
1142
  res.w[0] = 0x0000000000000000ull;
1143
  BID_RETURN (res);
1144
}
1145
  // x is not special and is not zero
1146
 
1147
  // if (exp <= -p) return -1.0 or +0.0
1148
if (x_exp <= 0x2ffc000000000000ull) {   // 0x2ffc000000000000ull == -34
1149
  if (x_sign) {
1150
    // if negative, return negative 1, because we know the coefficient
1151
    // is non-zero (would have been caught above)
1152
    res.w[1] = 0xb040000000000000ull;
1153
    res.w[0] = 0x0000000000000001ull;
1154
  } else {
1155
    // if positive, return positive 0, because we know coefficient is
1156
    // non-zero (would have been caught above)
1157
    res.w[1] = 0x3040000000000000ull;
1158
    res.w[0] = 0x0000000000000000ull;
1159
  }
1160
  BID_RETURN (res);
1161
}
1162
  // q = nr. of decimal digits in x
1163
  // determine first the nr. of bits in x
1164
if (C1.w[1] == 0) {
1165
  if (C1.w[0] >= 0x0020000000000000ull) {        // x >= 2^53
1166
    // split the 64-bit value in two 32-bit halves to avoid rounding errors
1167
    if (C1.w[0] >= 0x0000000100000000ull) {      // x >= 2^32
1168
      tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1169
      x_nr_bits =
1170
        33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1171
    } else {    // x < 2^32
1172
      tmp1.d = (double) (C1.w[0]);       // exact conversion
1173
      x_nr_bits =
1174
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1175
    }
1176
  } else {      // if x < 2^53
1177
    tmp1.d = (double) C1.w[0];   // exact conversion
1178
    x_nr_bits =
1179
      1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1180
  }
1181
} else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1182
  tmp1.d = (double) C1.w[1];    // exact conversion
1183
  x_nr_bits =
1184
    65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1185
}
1186
 
1187
q = nr_digits[x_nr_bits - 1].digits;
1188
if (q == 0) {
1189
  q = nr_digits[x_nr_bits - 1].digits1;
1190
  if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1191
      (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1192
       C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1193
    q++;
1194
}
1195
exp = (x_exp >> 49) - 6176;
1196
if (exp >= 0) {  // -exp <= 0
1197
  // the argument is an integer already
1198
  res.w[1] = x.w[1];
1199
  res.w[0] = x.w[0];
1200
  BID_RETURN (res);
1201
} else if ((q + exp) > 0) {      // exp < 0 and 1 <= -exp < q
1202
  // need to shift right -exp digits from the coefficient; the exp will be 0
1203
  ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x' 
1204
  // (number of digits to be chopped off)
1205
  // chop off ind digits from the lower part of C1 
1206
  // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1207
  // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1208
  // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1209
  // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1210
  //tmp64 = C1.w[0];
1211
  // if (ind <= 19) {
1212
  //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1213
  // } else {
1214
  //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1215
  //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1216
  // }
1217
  // if (C1.w[0] < tmp64) C1.w[1]++;
1218
  // if carry-out from C1.w[0], increment C1.w[1]
1219
  // calculate C* and f*
1220
  // C* is actually floor(C*) in this case
1221
  // C* and f* need shifting and masking, as shown by
1222
  // shiftright128[] and maskhigh128[]
1223
  // 1 <= x <= 34
1224
  // kx = 10^(-x) = ten2mk128[ind - 1]
1225
  // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1226
  // the approximation of 10^(-x) was rounded up to 118 bits
1227
  __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1228
  if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
1229
    res.w[1] = P256.w[3];
1230
    res.w[0] = P256.w[2];
1231
    // if positive, the truncated value is already the correct result
1232
    if (x_sign) {       // if negative
1233
      // redundant fstar.w[3] = 0;
1234
      // redundant fstar.w[2] = 0;
1235
      // redundant fstar.w[1] = P256.w[1];
1236
      // redundant fstar.w[0] = P256.w[0];
1237
      // fraction f* > 10^(-x) <=> inexact
1238
      // f* is in the right position to be compared with
1239
      // 10^(-x) from ten2mk128[]
1240
      if ((P256.w[1] > ten2mk128[ind - 1].w[1])
1241
          || (P256.w[1] == ten2mk128[ind - 1].w[1]
1242
              && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
1243
        if (++res.w[0] == 0) {
1244
          res.w[1]++;
1245
        }
1246
      }
1247
    }
1248
  } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1249
    shift = shiftright128[ind - 1];     // 0 <= shift <= 102
1250
    res.w[1] = (P256.w[3] >> shift);
1251
    res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1252
    // if positive, the truncated value is already the correct result
1253
    if (x_sign) {       // if negative
1254
      // redundant fstar.w[3] = 0;
1255
      fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1256
      fstar.w[1] = P256.w[1];
1257
      fstar.w[0] = P256.w[0];
1258
      // fraction f* > 10^(-x) <=> inexact
1259
      // f* is in the right position to be compared with
1260
      // 10^(-x) from ten2mk128[]
1261
      if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
1262
          (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1263
           fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1264
        if (++res.w[0] == 0) {
1265
          res.w[1]++;
1266
        }
1267
      }
1268
    }
1269
  } else {      // 22 <= ind - 1 <= 33
1270
    shift = shiftright128[ind - 1] - 64;        // 2 <= shift <= 38
1271
    res.w[1] = 0;
1272
    res.w[0] = P256.w[3] >> shift;
1273
    // if positive, the truncated value is already the correct result
1274
    if (x_sign) {       // if negative
1275
      fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1276
      fstar.w[2] = P256.w[2];
1277
      fstar.w[1] = P256.w[1];
1278
      fstar.w[0] = P256.w[0];
1279
      // fraction f* > 10^(-x) <=> inexact
1280
      // f* is in the right position to be compared with
1281
      // 10^(-x) from ten2mk128[]
1282
      if (fstar.w[3] || fstar.w[2]
1283
          || fstar.w[1] > ten2mk128[ind - 1].w[1]
1284
          || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1285
              && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1286
        if (++res.w[0] == 0) {
1287
          res.w[1]++;
1288
        }
1289
      }
1290
    }
1291
  }
1292
  res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1293
  BID_RETURN (res);
1294
} else {        // if exp < 0 and q + exp <= 0
1295
  if (x_sign) { // negative rounds down to -1.0
1296
    res.w[1] = 0xb040000000000000ull;
1297
    res.w[0] = 0x0000000000000001ull;
1298
  } else {      // positive rpunds down to +0.0
1299
    res.w[1] = 0x3040000000000000ull;
1300
    res.w[0] = 0x0000000000000000ull;
1301
  }
1302
  BID_RETURN (res);
1303
}
1304
}
1305
 
1306
/*****************************************************************************
1307
 *  BID128_round_integral_positive
1308
 ****************************************************************************/
1309
 
1310
BID128_FUNCTION_ARG1_NORND (bid128_round_integral_positive, x)
1311
 
1312
     UINT128 res;
1313
     UINT64 x_sign;
1314
     UINT64 x_exp;
1315
     int exp;                   // unbiased exponent
1316
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo 
1317
  // (all are UINT64)
1318
     BID_UI64DOUBLE tmp1;
1319
     unsigned int x_nr_bits;
1320
     int q, ind, shift;
1321
     UINT128 C1;
1322
  // UINT128 res is C* at first - represents up to 34 decimal digits ~ 
1323
  // 113 bits
1324
     UINT256 fstar;
1325
     UINT256 P256;
1326
 
1327
  // check for NaN or Infinity
1328
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1329
    // x is special
1330
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1331
  // if x = NaN, then res = Q (x)
1332
  // check first for non-canonical NaN payload
1333
  if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1334
      (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1335
       (x.w[0] > 0x38c15b09ffffffffull))) {
1336
    x.w[1] = x.w[1] & 0xffffc00000000000ull;
1337
    x.w[0] = 0x0ull;
1338
  }
1339
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1340
    // set invalid flag
1341
    *pfpsf |= INVALID_EXCEPTION;
1342
    // return quiet (x)
1343
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
1344
    res.w[0] = x.w[0];
1345
  } else {      // x is QNaN
1346
    // return x
1347
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
1348
    res.w[0] = x.w[0];
1349
  }
1350
  BID_RETURN (res)
1351
} else {        // x is not a NaN, so it must be infinity
1352
  if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1353
    // return +inf
1354
    res.w[1] = 0x7800000000000000ull;
1355
    res.w[0] = 0x0000000000000000ull;
1356
  } else {      // x is -inf 
1357
    // return -inf
1358
    res.w[1] = 0xf800000000000000ull;
1359
    res.w[0] = 0x0000000000000000ull;
1360
  }
1361
  BID_RETURN (res);
1362
}
1363
}
1364
  // unpack x
1365
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1366
C1.w[1] = x.w[1] & MASK_COEFF;
1367
C1.w[0] = x.w[0];
1368
 
1369
  // check for non-canonical values (treated as zero)
1370
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
1371
  // non-canonical
1372
  x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
1373
  C1.w[1] = 0;   // significand high
1374
  C1.w[0] = 0;    // significand low
1375
} else {        // G0_G1 != 11
1376
  x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
1377
  if (C1.w[1] > 0x0001ed09bead87c0ull ||
1378
      (C1.w[1] == 0x0001ed09bead87c0ull
1379
       && C1.w[0] > 0x378d8e63ffffffffull)) {
1380
    // x is non-canonical if coefficient is larger than 10^34 -1
1381
    C1.w[1] = 0;
1382
    C1.w[0] = 0;
1383
  } else {      // canonical
1384
    ;
1385
  }
1386
}
1387
 
1388
  // test for input equal to zero
1389
if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1390
  // x is 0
1391
  // return 0 preserving the sign bit and the preferred exponent 
1392
  // of MAX(Q(x), 0)
1393
  if (x_exp <= (0x1820ull << 49)) {
1394
    res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1395
  } else {
1396
    res.w[1] = x_sign | x_exp;
1397
  }
1398
  res.w[0] = 0x0000000000000000ull;
1399
  BID_RETURN (res);
1400
}
1401
  // x is not special and is not zero
1402
 
1403
  // if (exp <= -p) return -0.0 or +1.0
1404
if (x_exp <= 0x2ffc000000000000ull) {   // 0x2ffc000000000000ull == -34
1405
  if (x_sign) {
1406
    // if negative, return negative 0, because we know the coefficient 
1407
    // is non-zero (would have been caught above)
1408
    res.w[1] = 0xb040000000000000ull;
1409
    res.w[0] = 0x0000000000000000ull;
1410
  } else {
1411
    // if positive, return positive 1, because we know coefficient is 
1412
    // non-zero (would have been caught above)
1413
    res.w[1] = 0x3040000000000000ull;
1414
    res.w[0] = 0x0000000000000001ull;
1415
  }
1416
  BID_RETURN (res);
1417
}
1418
  // q = nr. of decimal digits in x
1419
  // determine first the nr. of bits in x
1420
if (C1.w[1] == 0) {
1421
  if (C1.w[0] >= 0x0020000000000000ull) {        // x >= 2^53
1422
    // split 64-bit value in two 32-bit halves to avoid rounding errors
1423
    if (C1.w[0] >= 0x0000000100000000ull) {      // x >= 2^32
1424
      tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1425
      x_nr_bits =
1426
        33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1427
    } else {    // x < 2^32
1428
      tmp1.d = (double) (C1.w[0]);       // exact conversion
1429
      x_nr_bits =
1430
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1431
    }
1432
  } else {      // if x < 2^53
1433
    tmp1.d = (double) C1.w[0];   // exact conversion
1434
    x_nr_bits =
1435
      1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1436
  }
1437
} else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1438
  tmp1.d = (double) C1.w[1];    // exact conversion
1439
  x_nr_bits =
1440
    65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1441
}
1442
 
1443
q = nr_digits[x_nr_bits - 1].digits;
1444
if (q == 0) {
1445
  q = nr_digits[x_nr_bits - 1].digits1;
1446
  if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1447
      (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1448
       C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1449
    q++;
1450
}
1451
exp = (x_exp >> 49) - 6176;
1452
if (exp >= 0) {  // -exp <= 0
1453
  // the argument is an integer already
1454
  res.w[1] = x.w[1];
1455
  res.w[0] = x.w[0];
1456
  BID_RETURN (res);
1457
} else if ((q + exp) > 0) {      // exp < 0 and 1 <= -exp < q
1458
  // need to shift right -exp digits from the coefficient; exp will be 0
1459
  ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x' 
1460
  // (number of digits to be chopped off)
1461
  // chop off ind digits from the lower part of C1 
1462
  // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1463
  // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1464
  // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1465
  // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1466
  // tmp64 = C1.w[0];
1467
  // if (ind <= 19) {
1468
  //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1469
  // } else {
1470
  //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1471
  //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1472
  // }
1473
  // if (C1.w[0] < tmp64) C1.w[1]++;  
1474
  // if carry-out from C1.w[0], increment C1.w[1]
1475
  // calculate C* and f*
1476
  // C* is actually floor(C*) in this case
1477
  // C* and f* need shifting and masking, as shown by
1478
  // shiftright128[] and maskhigh128[]
1479
  // 1 <= x <= 34
1480
  // kx = 10^(-x) = ten2mk128[ind - 1]
1481
  // C* = C1 * 10^(-x)
1482
  // the approximation of 10^(-x) was rounded up to 118 bits
1483
  __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1484
  if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
1485
    res.w[1] = P256.w[3];
1486
    res.w[0] = P256.w[2];
1487
    // if negative, the truncated value is already the correct result
1488
    if (!x_sign) {      // if positive
1489
      // redundant fstar.w[3] = 0;
1490
      // redundant fstar.w[2] = 0;
1491
      // redundant fstar.w[1] = P256.w[1]; 
1492
      // redundant fstar.w[0] = P256.w[0];
1493
      // fraction f* > 10^(-x) <=> inexact
1494
      // f* is in the right position to be compared with 
1495
      // 10^(-x) from ten2mk128[]
1496
      if ((P256.w[1] > ten2mk128[ind - 1].w[1])
1497
          || (P256.w[1] == ten2mk128[ind - 1].w[1]
1498
              && (P256.w[0] >= ten2mk128[ind - 1].w[0]))) {
1499
        if (++res.w[0] == 0) {
1500
          res.w[1]++;
1501
        }
1502
      }
1503
    }
1504
  } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1505
    shift = shiftright128[ind - 1];     // 3 <= shift <= 63
1506
    res.w[1] = (P256.w[3] >> shift);
1507
    res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1508
    // if negative, the truncated value is already the correct result
1509
    if (!x_sign) {      // if positive
1510
      // redundant fstar.w[3] = 0;
1511
      fstar.w[2] = P256.w[2] & maskhigh128[ind - 1];
1512
      fstar.w[1] = P256.w[1];
1513
      fstar.w[0] = P256.w[0];
1514
      // fraction f* > 10^(-x) <=> inexact
1515
      // f* is in the right position to be compared with 
1516
      // 10^(-x) from ten2mk128[]
1517
      if (fstar.w[2] || fstar.w[1] > ten2mk128[ind - 1].w[1] ||
1518
          (fstar.w[1] == ten2mk128[ind - 1].w[1] &&
1519
           fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1520
        if (++res.w[0] == 0) {
1521
          res.w[1]++;
1522
        }
1523
      }
1524
    }
1525
  } else {      // 22 <= ind - 1 <= 33
1526
    shift = shiftright128[ind - 1] - 64;        // 2 <= shift <= 38
1527
    res.w[1] = 0;
1528
    res.w[0] = P256.w[3] >> shift;
1529
    // if negative, the truncated value is already the correct result
1530
    if (!x_sign) {      // if positive
1531
      fstar.w[3] = P256.w[3] & maskhigh128[ind - 1];
1532
      fstar.w[2] = P256.w[2];
1533
      fstar.w[1] = P256.w[1];
1534
      fstar.w[0] = P256.w[0];
1535
      // fraction f* > 10^(-x) <=> inexact
1536
      // f* is in the right position to be compared with 
1537
      // 10^(-x) from ten2mk128[]
1538
      if (fstar.w[3] || fstar.w[2]
1539
          || fstar.w[1] > ten2mk128[ind - 1].w[1]
1540
          || (fstar.w[1] == ten2mk128[ind - 1].w[1]
1541
              && fstar.w[0] >= ten2mk128[ind - 1].w[0])) {
1542
        if (++res.w[0] == 0) {
1543
          res.w[1]++;
1544
        }
1545
      }
1546
    }
1547
  }
1548
  res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1549
  BID_RETURN (res);
1550
} else {        // if exp < 0 and q + exp <= 0
1551
  if (x_sign) { // negative rounds up to -0.0
1552
    res.w[1] = 0xb040000000000000ull;
1553
    res.w[0] = 0x0000000000000000ull;
1554
  } else {      // positive rpunds up to +1.0
1555
    res.w[1] = 0x3040000000000000ull;
1556
    res.w[0] = 0x0000000000000001ull;
1557
  }
1558
  BID_RETURN (res);
1559
}
1560
}
1561
 
1562
/*****************************************************************************
1563
 *  BID128_round_integral_zero
1564
 ****************************************************************************/
1565
 
1566
BID128_FUNCTION_ARG1_NORND (bid128_round_integral_zero, x)
1567
 
1568
     UINT128 res;
1569
     UINT64 x_sign;
1570
     UINT64 x_exp;
1571
     int exp;                   // unbiased exponent
1572
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo
1573
  // (all are UINT64)
1574
     BID_UI64DOUBLE tmp1;
1575
     unsigned int x_nr_bits;
1576
     int q, ind, shift;
1577
     UINT128 C1;
1578
  // UINT128 res is C* at first - represents up to 34 decimal digits ~
1579
  // 113 bits
1580
     UINT256 P256;
1581
 
1582
  // check for NaN or Infinity
1583
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1584
    // x is special
1585
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1586
  // if x = NaN, then res = Q (x)
1587
  // check first for non-canonical NaN payload
1588
  if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1589
      (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1590
       (x.w[0] > 0x38c15b09ffffffffull))) {
1591
    x.w[1] = x.w[1] & 0xffffc00000000000ull;
1592
    x.w[0] = 0x0ull;
1593
  }
1594
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1595
    // set invalid flag
1596
    *pfpsf |= INVALID_EXCEPTION;
1597
    // return quiet (x)
1598
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
1599
    res.w[0] = x.w[0];
1600
  } else {      // x is QNaN
1601
    // return x
1602
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
1603
    res.w[0] = x.w[0];
1604
  }
1605
  BID_RETURN (res)
1606
} else {        // x is not a NaN, so it must be infinity
1607
  if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1608
    // return +inf
1609
    res.w[1] = 0x7800000000000000ull;
1610
    res.w[0] = 0x0000000000000000ull;
1611
  } else {      // x is -inf 
1612
    // return -inf
1613
    res.w[1] = 0xf800000000000000ull;
1614
    res.w[0] = 0x0000000000000000ull;
1615
  }
1616
  BID_RETURN (res);
1617
}
1618
}
1619
  // unpack x
1620
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1621
C1.w[1] = x.w[1] & MASK_COEFF;
1622
C1.w[0] = x.w[0];
1623
 
1624
  // check for non-canonical values (treated as zero)
1625
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
1626
  // non-canonical
1627
  x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
1628
  C1.w[1] = 0;   // significand high
1629
  C1.w[0] = 0;    // significand low
1630
} else {        // G0_G1 != 11
1631
  x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
1632
  if (C1.w[1] > 0x0001ed09bead87c0ull ||
1633
      (C1.w[1] == 0x0001ed09bead87c0ull
1634
       && C1.w[0] > 0x378d8e63ffffffffull)) {
1635
    // x is non-canonical if coefficient is larger than 10^34 -1
1636
    C1.w[1] = 0;
1637
    C1.w[0] = 0;
1638
  } else {      // canonical
1639
    ;
1640
  }
1641
}
1642
 
1643
  // test for input equal to zero
1644
if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1645
  // x is 0
1646
  // return 0 preserving the sign bit and the preferred exponent
1647
  // of MAX(Q(x), 0)
1648
  if (x_exp <= (0x1820ull << 49)) {
1649
    res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1650
  } else {
1651
    res.w[1] = x_sign | x_exp;
1652
  }
1653
  res.w[0] = 0x0000000000000000ull;
1654
  BID_RETURN (res);
1655
}
1656
  // x is not special and is not zero
1657
 
1658
  // if (exp <= -p) return -0.0 or +0.0
1659
if (x_exp <= 0x2ffc000000000000ull) {   // 0x2ffc000000000000ull == -34
1660
  res.w[1] = x_sign | 0x3040000000000000ull;
1661
  res.w[0] = 0x0000000000000000ull;
1662
  BID_RETURN (res);
1663
}
1664
  // q = nr. of decimal digits in x
1665
  // determine first the nr. of bits in x
1666
if (C1.w[1] == 0) {
1667
  if (C1.w[0] >= 0x0020000000000000ull) {        // x >= 2^53
1668
    // split the 64-bit value in two 32-bit halves to avoid rounding errors
1669
    if (C1.w[0] >= 0x0000000100000000ull) {      // x >= 2^32
1670
      tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1671
      x_nr_bits =
1672
        33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1673
    } else {    // x < 2^32
1674
      tmp1.d = (double) (C1.w[0]);       // exact conversion
1675
      x_nr_bits =
1676
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1677
    }
1678
  } else {      // if x < 2^53
1679
    tmp1.d = (double) C1.w[0];   // exact conversion
1680
    x_nr_bits =
1681
      1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1682
  }
1683
} else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1684
  tmp1.d = (double) C1.w[1];    // exact conversion
1685
  x_nr_bits =
1686
    65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1687
}
1688
 
1689
q = nr_digits[x_nr_bits - 1].digits;
1690
if (q == 0) {
1691
  q = nr_digits[x_nr_bits - 1].digits1;
1692
  if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1693
      (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1694
       C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1695
    q++;
1696
}
1697
exp = (x_exp >> 49) - 6176;
1698
if (exp >= 0) {  // -exp <= 0
1699
  // the argument is an integer already
1700
  res.w[1] = x.w[1];
1701
  res.w[0] = x.w[0];
1702
  BID_RETURN (res);
1703
} else if ((q + exp) > 0) {      // exp < 0 and 1 <= -exp < q
1704
  // need to shift right -exp digits from the coefficient; the exp will be 0
1705
  ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x'
1706
  // (number of digits to be chopped off)
1707
  // chop off ind digits from the lower part of C1 
1708
  // FOR ROUND_TO_NEAREST, WE ADD 1/2 ULP(y) then truncate
1709
  // FOR ROUND_TO_ZERO, WE DON'T NEED TO ADD 1/2 ULP
1710
  // FOR ROUND_TO_POSITIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF POSITIVE
1711
  // FOR ROUND_TO_NEGATIVE_INFINITY, WE TRUNCATE, THEN ADD 1 IF NEGATIVE
1712
  //tmp64 = C1.w[0];
1713
  // if (ind <= 19) {
1714
  //   C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1715
  // } else {
1716
  //   C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1717
  //   C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1718
  // }
1719
  // if (C1.w[0] < tmp64) C1.w[1]++;  
1720
  // if carry-out from C1.w[0], increment C1.w[1]
1721
  // calculate C* and f*
1722
  // C* is actually floor(C*) in this case
1723
  // C* and f* need shifting and masking, as shown by
1724
  // shiftright128[] and maskhigh128[]
1725
  // 1 <= x <= 34
1726
  // kx = 10^(-x) = ten2mk128[ind - 1]
1727
  // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1728
  // the approximation of 10^(-x) was rounded up to 118 bits
1729
  __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1730
  if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
1731
    res.w[1] = P256.w[3];
1732
    res.w[0] = P256.w[2];
1733
  } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1734
    shift = shiftright128[ind - 1];     // 3 <= shift <= 63
1735
    res.w[1] = (P256.w[3] >> shift);
1736
    res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1737
  } else {      // 22 <= ind - 1 <= 33
1738
    shift = shiftright128[ind - 1] - 64;        // 2 <= shift <= 38
1739
    res.w[1] = 0;
1740
    res.w[0] = P256.w[3] >> shift;
1741
  }
1742
  res.w[1] = x_sign | 0x3040000000000000ull | res.w[1];
1743
  BID_RETURN (res);
1744
} else {        // if exp < 0 and q + exp <= 0 the result is +0 or -0
1745
  res.w[1] = x_sign | 0x3040000000000000ull;
1746
  res.w[0] = 0x0000000000000000ull;
1747
  BID_RETURN (res);
1748
}
1749
}
1750
 
1751
/*****************************************************************************
1752
 *  BID128_round_integral_nearest_away
1753
 ****************************************************************************/
1754
 
1755
BID128_FUNCTION_ARG1_NORND (bid128_round_integral_nearest_away, x)
1756
 
1757
     UINT128 res;
1758
     UINT64 x_sign;
1759
     UINT64 x_exp;
1760
     int exp;                   // unbiased exponent
1761
  // Note: C1.w[1], C1.w[0] represent x_signif_hi, x_signif_lo 
1762
  // (all are UINT64)
1763
     UINT64 tmp64;
1764
     BID_UI64DOUBLE tmp1;
1765
     unsigned int x_nr_bits;
1766
     int q, ind, shift;
1767
     UINT128 C1;
1768
  // UINT128 res is C* at first - represents up to 34 decimal digits ~ 
1769
  // 113 bits
1770
  // UINT256 fstar;
1771
     UINT256 P256;
1772
 
1773
  // check for NaN or Infinity
1774
if ((x.w[1] & MASK_SPECIAL) == MASK_SPECIAL) {
1775
    // x is special
1776
if ((x.w[1] & MASK_NAN) == MASK_NAN) {  // x is NAN
1777
  // if x = NaN, then res = Q (x)
1778
  // check first for non-canonical NaN payload
1779
  if (((x.w[1] & 0x00003fffffffffffull) > 0x0000314dc6448d93ull) ||
1780
      (((x.w[1] & 0x00003fffffffffffull) == 0x0000314dc6448d93ull) &&
1781
       (x.w[0] > 0x38c15b09ffffffffull))) {
1782
    x.w[1] = x.w[1] & 0xffffc00000000000ull;
1783
    x.w[0] = 0x0ull;
1784
  }
1785
  if ((x.w[1] & MASK_SNAN) == MASK_SNAN) {      // x is SNAN
1786
    // set invalid flag
1787
    *pfpsf |= INVALID_EXCEPTION;
1788
    // return quiet (x)
1789
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out also G[6]-G[16]
1790
    res.w[0] = x.w[0];
1791
  } else {      // x is QNaN
1792
    // return x
1793
    res.w[1] = x.w[1] & 0xfc003fffffffffffull;  // clear out G[6]-G[16]
1794
    res.w[0] = x.w[0];
1795
  }
1796
  BID_RETURN (res)
1797
} else {        // x is not a NaN, so it must be infinity
1798
  if ((x.w[1] & MASK_SIGN) == 0x0ull) { // x is +inf
1799
    // return +inf
1800
    res.w[1] = 0x7800000000000000ull;
1801
    res.w[0] = 0x0000000000000000ull;
1802
  } else {      // x is -inf 
1803
    // return -inf
1804
    res.w[1] = 0xf800000000000000ull;
1805
    res.w[0] = 0x0000000000000000ull;
1806
  }
1807
  BID_RETURN (res);
1808
}
1809
}
1810
  // unpack x
1811
x_sign = x.w[1] & MASK_SIGN;    // 0 for positive, MASK_SIGN for negative
1812
C1.w[1] = x.w[1] & MASK_COEFF;
1813
C1.w[0] = x.w[0];
1814
 
1815
  // check for non-canonical values (treated as zero)
1816
if ((x.w[1] & 0x6000000000000000ull) == 0x6000000000000000ull) {        // G0_G1=11
1817
  // non-canonical
1818
  x_exp = (x.w[1] << 2) & MASK_EXP;     // biased and shifted left 49 bits
1819
  C1.w[1] = 0;   // significand high
1820
  C1.w[0] = 0;    // significand low
1821
} else {        // G0_G1 != 11
1822
  x_exp = x.w[1] & MASK_EXP;    // biased and shifted left 49 bits
1823
  if (C1.w[1] > 0x0001ed09bead87c0ull ||
1824
      (C1.w[1] == 0x0001ed09bead87c0ull
1825
       && C1.w[0] > 0x378d8e63ffffffffull)) {
1826
    // x is non-canonical if coefficient is larger than 10^34 -1
1827
    C1.w[1] = 0;
1828
    C1.w[0] = 0;
1829
  } else {      // canonical
1830
    ;
1831
  }
1832
}
1833
 
1834
  // test for input equal to zero
1835
if ((C1.w[1] == 0x0ull) && (C1.w[0] == 0x0ull)) {
1836
  // x is 0
1837
  // return 0 preserving the sign bit and the preferred exponent
1838
  // of MAX(Q(x), 0)
1839
  if (x_exp <= (0x1820ull << 49)) {
1840
    res.w[1] = (x.w[1] & 0x8000000000000000ull) | 0x3040000000000000ull;
1841
  } else {
1842
    res.w[1] = x_sign | x_exp;
1843
  }
1844
  res.w[0] = 0x0000000000000000ull;
1845
  BID_RETURN (res);
1846
}
1847
  // x is not special and is not zero
1848
 
1849
  // if (exp <= -(p+1)) return 0.0
1850
if (x_exp <= 0x2ffa000000000000ull) {   // 0x2ffa000000000000ull == -35
1851
  res.w[1] = x_sign | 0x3040000000000000ull;
1852
  res.w[0] = 0x0000000000000000ull;
1853
  BID_RETURN (res);
1854
}
1855
  // q = nr. of decimal digits in x
1856
  //  determine first the nr. of bits in x
1857
if (C1.w[1] == 0) {
1858
  if (C1.w[0] >= 0x0020000000000000ull) {        // x >= 2^53
1859
    // split the 64-bit value in two 32-bit halves to avoid rounding errors
1860
    if (C1.w[0] >= 0x0000000100000000ull) {      // x >= 2^32
1861
      tmp1.d = (double) (C1.w[0] >> 32); // exact conversion
1862
      x_nr_bits =
1863
        33 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1864
    } else {    // x < 2^32
1865
      tmp1.d = (double) (C1.w[0]);       // exact conversion
1866
      x_nr_bits =
1867
        1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1868
    }
1869
  } else {      // if x < 2^53
1870
    tmp1.d = (double) C1.w[0];   // exact conversion
1871
    x_nr_bits =
1872
      1 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1873
  }
1874
} else {        // C1.w[1] != 0 => nr. bits = 64 + nr_bits (C1.w[1])
1875
  tmp1.d = (double) C1.w[1];    // exact conversion
1876
  x_nr_bits =
1877
    65 + ((((unsigned int) (tmp1.ui64 >> 52)) & 0x7ff) - 0x3ff);
1878
}
1879
 
1880
q = nr_digits[x_nr_bits - 1].digits;
1881
if (q == 0) {
1882
  q = nr_digits[x_nr_bits - 1].digits1;
1883
  if (C1.w[1] > nr_digits[x_nr_bits - 1].threshold_hi ||
1884
      (C1.w[1] == nr_digits[x_nr_bits - 1].threshold_hi &&
1885
       C1.w[0] >= nr_digits[x_nr_bits - 1].threshold_lo))
1886
    q++;
1887
}
1888
exp = (x_exp >> 49) - 6176;
1889
if (exp >= 0) {  // -exp <= 0
1890
  // the argument is an integer already
1891
  res.w[1] = x.w[1];
1892
  res.w[0] = x.w[0];
1893
  BID_RETURN (res);
1894
} else if ((q + exp) >= 0) {     // exp < 0 and 1 <= -exp <= q
1895
  // need to shift right -exp digits from the coefficient; the exp will be 0
1896
  ind = -exp;   // 1 <= ind <= 34; ind is a synonym for 'x'
1897
  // chop off ind digits from the lower part of C1 
1898
  // C1 = C1 + 1/2 * 10^x where the result C1 fits in 127 bits
1899
  tmp64 = C1.w[0];
1900
  if (ind <= 19) {
1901
    C1.w[0] = C1.w[0] + midpoint64[ind - 1];
1902
  } else {
1903
    C1.w[0] = C1.w[0] + midpoint128[ind - 20].w[0];
1904
    C1.w[1] = C1.w[1] + midpoint128[ind - 20].w[1];
1905
  }
1906
  if (C1.w[0] < tmp64)
1907
    C1.w[1]++;
1908
  // calculate C* and f*
1909
  // C* is actually floor(C*) in this case
1910
  // C* and f* need shifting and masking, as shown by
1911
  // shiftright128[] and maskhigh128[]
1912
  // 1 <= x <= 34
1913
  // kx = 10^(-x) = ten2mk128[ind - 1]
1914
  // C* = (C1 + 1/2 * 10^x) * 10^(-x)
1915
  // the approximation of 10^(-x) was rounded up to 118 bits
1916
  __mul_128x128_to_256 (P256, C1, ten2mk128[ind - 1]);
1917
  // the top Ex bits of 10^(-x) are T* = ten2mk128trunc[ind], e.g.
1918
  // if x=1, T*=ten2mk128trunc[0]=0x19999999999999999999999999999999
1919
  // if (0 < f* < 10^(-x)) then the result is a midpoint
1920
  //   if floor(C*) is even then C* = floor(C*) - logical right
1921
  //       shift; C* has p decimal digits, correct by Prop. 1)
1922
  //   else if floor(C*) is odd C* = floor(C*)-1 (logical right
1923
  //       shift; C* has p decimal digits, correct by Pr. 1)
1924
  // else
1925
  //   C* = floor(C*) (logical right shift; C has p decimal digits,
1926
  //       correct by Property 1)
1927
  // n = C* * 10^(e+x)
1928
 
1929
  // shift right C* by Ex-128 = shiftright128[ind]
1930
  if (ind - 1 <= 2) {   // 0 <= ind - 1 <= 2 => shift = 0
1931
    res.w[1] = P256.w[3];
1932
    res.w[0] = P256.w[2];
1933
  } else if (ind - 1 <= 21) {   // 3 <= ind - 1 <= 21 => 3 <= shift <= 63
1934
    shift = shiftright128[ind - 1];     // 3 <= shift <= 63
1935
    res.w[0] = (P256.w[3] << (64 - shift)) | (P256.w[2] >> shift);
1936
    res.w[1] = (P256.w[3] >> shift);
1937
  } else {      // 22 <= ind - 1 <= 33
1938
    shift = shiftright128[ind - 1];     // 2 <= shift <= 38
1939
    res.w[1] = 0;
1940
    res.w[0] = (P256.w[3] >> (shift - 64));      // 2 <= shift - 64 <= 38
1941
  }
1942
  // if the result was a midpoint, it was already rounded away from zero
1943
  res.w[1] |= x_sign | 0x3040000000000000ull;
1944
  BID_RETURN (res);
1945
} else {        // if ((q + exp) < 0) <=> q < -exp
1946
  // the result is +0 or -0
1947
  res.w[1] = x_sign | 0x3040000000000000ull;
1948
  res.w[0] = 0x0000000000000000ull;
1949
  BID_RETURN (res);
1950
}
1951
}

powered by: WebSVN 2.1.0

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