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

Subversion Repositories openrisc

[/] [openrisc/] [trunk/] [gnu-dev/] [or1k-gcc/] [libgcc/] [config/] [libbid/] [bid_inline_add.h] - Blame information for rev 801

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

Line No. Rev Author Line
1 734 jeremybenn
/* Copyright (C) 2007, 2009  Free Software Foundation, Inc.
2
 
3
This file is part of GCC.
4
 
5
GCC is free software; you can redistribute it and/or modify it under
6
the terms of the GNU General Public License as published by the Free
7
Software Foundation; either version 3, or (at your option) any later
8
version.
9
 
10
GCC is distributed in the hope that it will be useful, but WITHOUT ANY
11
WARRANTY; without even the implied warranty of MERCHANTABILITY or
12
FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
13
for more details.
14
 
15
Under Section 7 of GPL version 3, you are granted additional
16
permissions described in the GCC Runtime Library Exception, version
17
3.1, as published by the Free Software Foundation.
18
 
19
You should have received a copy of the GNU General Public License and
20
a copy of the GCC Runtime Library Exception along with this program;
21
see the files COPYING3 and COPYING.RUNTIME respectively.  If not, see
22
<http://www.gnu.org/licenses/>.  */
23
 
24
/*****************************************************************************
25
 *
26
 *    Helper add functions (for fma)
27
 *
28
 *    __BID_INLINE__ UINT64 get_add64(
29
 *        UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
30
 *        UINT64 sign_y, int exponent_y, UINT64 coefficient_y,
31
 *                                       int rounding_mode)
32
 *
33
 *   __BID_INLINE__ UINT64 get_add128(
34
 *                       UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
35
 *                       UINT64 sign_y, int final_exponent_y, UINT128 CY,
36
 *                       int extra_digits, int rounding_mode)
37
 *
38
 *****************************************************************************
39
 *
40
 *  Algorithm description:
41
 *
42
 *  get_add64:  same as BID64 add, but arguments are unpacked and there
43
 *                                 are no special case checks
44
 *
45
 *  get_add128: add 64-bit coefficient to 128-bit product (which contains
46
 *                                        16+extra_digits decimal digits),
47
 *                         return BID64 result
48
 *              - the exponents are compared and the two coefficients are
49
 *                properly aligned for addition/subtraction
50
 *              - multiple paths are needed
51
 *              - final result exponent is calculated and the lower term is
52
 *                      rounded first if necessary, to avoid manipulating
53
 *                      coefficients longer than 128 bits
54
 *
55
 ****************************************************************************/
56
 
57
#ifndef _INLINE_BID_ADD_H_
58
#define _INLINE_BID_ADD_H_
59
 
60
#include "bid_internal.h"
61
 
62
#define MAX_FORMAT_DIGITS     16
63
#define DECIMAL_EXPONENT_BIAS 398
64
#define MASK_BINARY_EXPONENT  0x7ff0000000000000ull
65
#define BINARY_EXPONENT_BIAS  0x3ff
66
#define UPPER_EXPON_LIMIT     51
67
 
68
///////////////////////////////////////////////////////////////////////
69
//
70
// get_add64() is essentially the same as bid_add(), except that 
71
//             the arguments are unpacked
72
//
73
//////////////////////////////////////////////////////////////////////
74
__BID_INLINE__ UINT64
75
get_add64 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
76
           UINT64 sign_y, int exponent_y, UINT64 coefficient_y,
77
           int rounding_mode, unsigned *fpsc) {
78
  UINT128 CA, CT, CT_new;
79
  UINT64 sign_a, sign_b, coefficient_a, coefficient_b, sign_s, sign_ab,
80
    rem_a;
81
  UINT64 saved_ca, saved_cb, C0_64, C64, remainder_h, T1, carry, tmp,
82
    C64_new;
83
  int_double tempx;
84
  int exponent_a, exponent_b, diff_dec_expon;
85
  int bin_expon_ca, extra_digits, amount, scale_k, scale_ca;
86
  unsigned rmode, status;
87
 
88
  // sort arguments by exponent
89
  if (exponent_x <= exponent_y) {
90
    sign_a = sign_y;
91
    exponent_a = exponent_y;
92
    coefficient_a = coefficient_y;
93
    sign_b = sign_x;
94
    exponent_b = exponent_x;
95
    coefficient_b = coefficient_x;
96
  } else {
97
    sign_a = sign_x;
98
    exponent_a = exponent_x;
99
    coefficient_a = coefficient_x;
100
    sign_b = sign_y;
101
    exponent_b = exponent_y;
102
    coefficient_b = coefficient_y;
103
  }
104
 
105
  // exponent difference
106
  diff_dec_expon = exponent_a - exponent_b;
107
 
108
  /* get binary coefficients of x and y */
109
 
110
  //--- get number of bits in the coefficients of x and y ---
111
 
112
  tempx.d = (double) coefficient_a;
113
  bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
114
 
115
  if (!coefficient_a) {
116
    return get_BID64 (sign_b, exponent_b, coefficient_b, rounding_mode,
117
                      fpsc);
118
  }
119
  if (diff_dec_expon > MAX_FORMAT_DIGITS) {
120
    // normalize a to a 16-digit coefficient
121
 
122
    scale_ca = estimate_decimal_digits[bin_expon_ca];
123
    if (coefficient_a >= power10_table_128[scale_ca].w[0])
124
      scale_ca++;
125
 
126
    scale_k = 16 - scale_ca;
127
 
128
    coefficient_a *= power10_table_128[scale_k].w[0];
129
 
130
    diff_dec_expon -= scale_k;
131
    exponent_a -= scale_k;
132
 
133
    /* get binary coefficients of x and y */
134
 
135
    //--- get number of bits in the coefficients of x and y ---
136
    tempx.d = (double) coefficient_a;
137
    bin_expon_ca = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
138
 
139
    if (diff_dec_expon > MAX_FORMAT_DIGITS) {
140
#ifdef SET_STATUS_FLAGS
141
      if (coefficient_b) {
142
        __set_status_flags (fpsc, INEXACT_EXCEPTION);
143
      }
144
#endif
145
 
146
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
147
#ifndef IEEE_ROUND_NEAREST
148
      if (((rounding_mode) & 3) && coefficient_b)       // not ROUNDING_TO_NEAREST
149
      {
150
        switch (rounding_mode) {
151
        case ROUNDING_DOWN:
152
          if (sign_b) {
153
            coefficient_a -= ((((SINT64) sign_a) >> 63) | 1);
154
            if (coefficient_a < 1000000000000000ull) {
155
              exponent_a--;
156
              coefficient_a = 9999999999999999ull;
157
            } else if (coefficient_a >= 10000000000000000ull) {
158
              exponent_a++;
159
              coefficient_a = 1000000000000000ull;
160
            }
161
          }
162
          break;
163
        case ROUNDING_UP:
164
          if (!sign_b) {
165
            coefficient_a += ((((SINT64) sign_a) >> 63) | 1);
166
            if (coefficient_a < 1000000000000000ull) {
167
              exponent_a--;
168
              coefficient_a = 9999999999999999ull;
169
            } else if (coefficient_a >= 10000000000000000ull) {
170
              exponent_a++;
171
              coefficient_a = 1000000000000000ull;
172
            }
173
          }
174
          break;
175
        default:        // RZ
176
          if (sign_a != sign_b) {
177
            coefficient_a--;
178
            if (coefficient_a < 1000000000000000ull) {
179
              exponent_a--;
180
              coefficient_a = 9999999999999999ull;
181
            }
182
          }
183
          break;
184
        }
185
      } else
186
#endif
187
#endif
188
        // check special case here
189
        if ((coefficient_a == 1000000000000000ull)
190
            && (diff_dec_expon == MAX_FORMAT_DIGITS + 1)
191
            && (sign_a ^ sign_b)
192
            && (coefficient_b > 5000000000000000ull)) {
193
        coefficient_a = 9999999999999999ull;
194
        exponent_a--;
195
      }
196
 
197
      return get_BID64 (sign_a, exponent_a, coefficient_a,
198
                        rounding_mode, fpsc);
199
    }
200
  }
201
  // test whether coefficient_a*10^(exponent_a-exponent_b)  may exceed 2^62
202
  if (bin_expon_ca + estimate_bin_expon[diff_dec_expon] < 60) {
203
    // coefficient_a*10^(exponent_a-exponent_b)<2^63
204
 
205
    // multiply by 10^(exponent_a-exponent_b)
206
    coefficient_a *= power10_table_128[diff_dec_expon].w[0];
207
 
208
    // sign mask
209
    sign_b = ((SINT64) sign_b) >> 63;
210
    // apply sign to coeff. of b
211
    coefficient_b = (coefficient_b + sign_b) ^ sign_b;
212
 
213
    // apply sign to coefficient a
214
    sign_a = ((SINT64) sign_a) >> 63;
215
    coefficient_a = (coefficient_a + sign_a) ^ sign_a;
216
 
217
    coefficient_a += coefficient_b;
218
    // get sign
219
    sign_s = ((SINT64) coefficient_a) >> 63;
220
    coefficient_a = (coefficient_a + sign_s) ^ sign_s;
221
    sign_s &= 0x8000000000000000ull;
222
 
223
    // coefficient_a < 10^16 ?
224
    if (coefficient_a < power10_table_128[MAX_FORMAT_DIGITS].w[0]) {
225
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
226
#ifndef IEEE_ROUND_NEAREST
227
      if (rounding_mode == ROUNDING_DOWN && (!coefficient_a)
228
          && sign_a != sign_b)
229
        sign_s = 0x8000000000000000ull;
230
#endif
231
#endif
232
      return get_BID64 (sign_s, exponent_b, coefficient_a,
233
                        rounding_mode, fpsc);
234
    }
235
    // otherwise rounding is necessary
236
 
237
    // already know coefficient_a<10^19
238
    // coefficient_a < 10^17 ?
239
    if (coefficient_a < power10_table_128[17].w[0])
240
      extra_digits = 1;
241
    else if (coefficient_a < power10_table_128[18].w[0])
242
      extra_digits = 2;
243
    else
244
      extra_digits = 3;
245
 
246
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
247
#ifndef IEEE_ROUND_NEAREST
248
    rmode = rounding_mode;
249
    if (sign_s && (unsigned) (rmode - 1) < 2)
250
      rmode = 3 - rmode;
251
#else
252
    rmode = 0;
253
#endif
254
#else
255
    rmode = 0;
256
#endif
257
    coefficient_a += round_const_table[rmode][extra_digits];
258
 
259
    // get P*(2^M[extra_digits])/10^extra_digits
260
    __mul_64x64_to_128 (CT, coefficient_a,
261
                        reciprocals10_64[extra_digits]);
262
 
263
    // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
264
    amount = short_recip_scale[extra_digits];
265
    C64 = CT.w[1] >> amount;
266
 
267
  } else {
268
    // coefficient_a*10^(exponent_a-exponent_b) is large
269
    sign_s = sign_a;
270
 
271
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
272
#ifndef IEEE_ROUND_NEAREST
273
    rmode = rounding_mode;
274
    if (sign_s && (unsigned) (rmode - 1) < 2)
275
      rmode = 3 - rmode;
276
#else
277
    rmode = 0;
278
#endif
279
#else
280
    rmode = 0;
281
#endif
282
 
283
    // check whether we can take faster path
284
    scale_ca = estimate_decimal_digits[bin_expon_ca];
285
 
286
    sign_ab = sign_a ^ sign_b;
287
    sign_ab = ((SINT64) sign_ab) >> 63;
288
 
289
    // T1 = 10^(16-diff_dec_expon)
290
    T1 = power10_table_128[16 - diff_dec_expon].w[0];
291
 
292
    // get number of digits in coefficient_a
293
    //P_ca = power10_table_128[scale_ca].w[0];
294
    //P_ca_m1 = power10_table_128[scale_ca-1].w[0];
295
    if (coefficient_a >= power10_table_128[scale_ca].w[0]) {
296
      scale_ca++;
297
      //P_ca_m1 = P_ca;
298
      //P_ca = power10_table_128[scale_ca].w[0];
299
    }
300
 
301
    scale_k = 16 - scale_ca;
302
 
303
    // apply sign
304
    //Ts = (T1 + sign_ab) ^ sign_ab;
305
 
306
    // test range of ca
307
    //X = coefficient_a + Ts - P_ca_m1;
308
 
309
    // addition
310
    saved_ca = coefficient_a - T1;
311
    coefficient_a =
312
      (SINT64) saved_ca *(SINT64) power10_table_128[scale_k].w[0];
313
    extra_digits = diff_dec_expon - scale_k;
314
 
315
    // apply sign
316
    saved_cb = (coefficient_b + sign_ab) ^ sign_ab;
317
    // add 10^16 and rounding constant
318
    coefficient_b =
319
      saved_cb + 10000000000000000ull +
320
      round_const_table[rmode][extra_digits];
321
 
322
    // get P*(2^M[extra_digits])/10^extra_digits
323
    __mul_64x64_to_128 (CT, coefficient_b,
324
                        reciprocals10_64[extra_digits]);
325
 
326
    // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
327
    amount = short_recip_scale[extra_digits];
328
    C0_64 = CT.w[1] >> amount;
329
 
330
    // result coefficient 
331
    C64 = C0_64 + coefficient_a;
332
    // filter out difficult (corner) cases
333
    // the following test is equivalent to 
334
    // ( (initial_coefficient_a + Ts) < P_ca && 
335
    //     (initial_coefficient_a + Ts) > P_ca_m1 ), 
336
    // which ensures the number of digits in coefficient_a does not change 
337
    // after adding (the appropriately scaled and rounded) coefficient_b
338
    if ((UINT64) (C64 - 1000000000000000ull - 1) >
339
        9000000000000000ull - 2) {
340
      if (C64 >= 10000000000000000ull) {
341
        // result has more than 16 digits
342
        if (!scale_k) {
343
          // must divide coeff_a by 10
344
          saved_ca = saved_ca + T1;
345
          __mul_64x64_to_128 (CA, saved_ca, 0x3333333333333334ull);
346
          //reciprocals10_64[1]);
347
          coefficient_a = CA.w[1] >> 1;
348
          rem_a =
349
            saved_ca - (coefficient_a << 3) - (coefficient_a << 1);
350
          coefficient_a = coefficient_a - T1;
351
 
352
          saved_cb +=
353
            /*90000000000000000 */ +rem_a *
354
            power10_table_128[diff_dec_expon].w[0];
355
        } else
356
          coefficient_a =
357
            (SINT64) (saved_ca - T1 -
358
                      (T1 << 3)) * (SINT64) power10_table_128[scale_k -
359
                                                              1].w[0];
360
 
361
        extra_digits++;
362
        coefficient_b =
363
          saved_cb + 100000000000000000ull +
364
          round_const_table[rmode][extra_digits];
365
 
366
        // get P*(2^M[extra_digits])/10^extra_digits
367
        __mul_64x64_to_128 (CT, coefficient_b,
368
                            reciprocals10_64[extra_digits]);
369
 
370
        // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
371
        amount = short_recip_scale[extra_digits];
372
        C0_64 = CT.w[1] >> amount;
373
 
374
        // result coefficient 
375
        C64 = C0_64 + coefficient_a;
376
      } else if (C64 <= 1000000000000000ull) {
377
        // less than 16 digits in result
378
        coefficient_a =
379
          (SINT64) saved_ca *(SINT64) power10_table_128[scale_k +
380
                                                        1].w[0];
381
        //extra_digits --;
382
        exponent_b--;
383
        coefficient_b =
384
          (saved_cb << 3) + (saved_cb << 1) + 100000000000000000ull +
385
          round_const_table[rmode][extra_digits];
386
 
387
        // get P*(2^M[extra_digits])/10^extra_digits
388
        __mul_64x64_to_128 (CT_new, coefficient_b,
389
                            reciprocals10_64[extra_digits]);
390
 
391
        // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
392
        amount = short_recip_scale[extra_digits];
393
        C0_64 = CT_new.w[1] >> amount;
394
 
395
        // result coefficient 
396
        C64_new = C0_64 + coefficient_a;
397
        if (C64_new < 10000000000000000ull) {
398
          C64 = C64_new;
399
#ifdef SET_STATUS_FLAGS
400
          CT = CT_new;
401
#endif
402
        } else
403
          exponent_b++;
404
      }
405
 
406
    }
407
 
408
  }
409
 
410
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
411
#ifndef IEEE_ROUND_NEAREST
412
  if (rmode == 0)        //ROUNDING_TO_NEAREST
413
#endif
414
    if (C64 & 1) {
415
      // check whether fractional part of initial_P/10^extra_digits 
416
      // is exactly .5
417
      // this is the same as fractional part of 
418
      //      (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
419
 
420
      // get remainder
421
      remainder_h = CT.w[1] << (64 - amount);
422
 
423
      // test whether fractional part is 0
424
      if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
425
        C64--;
426
      }
427
    }
428
#endif
429
 
430
#ifdef SET_STATUS_FLAGS
431
  status = INEXACT_EXCEPTION;
432
 
433
  // get remainder
434
  remainder_h = CT.w[1] << (64 - amount);
435
 
436
  switch (rmode) {
437
  case ROUNDING_TO_NEAREST:
438
  case ROUNDING_TIES_AWAY:
439
    // test whether fractional part is 0
440
    if ((remainder_h == 0x8000000000000000ull)
441
        && (CT.w[0] < reciprocals10_64[extra_digits]))
442
      status = EXACT_STATUS;
443
    break;
444
  case ROUNDING_DOWN:
445
  case ROUNDING_TO_ZERO:
446
    if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
447
      status = EXACT_STATUS;
448
    break;
449
  default:
450
    // round up
451
    __add_carry_out (tmp, carry, CT.w[0],
452
                     reciprocals10_64[extra_digits]);
453
    if ((remainder_h >> (64 - amount)) + carry >=
454
        (((UINT64) 1) << amount))
455
      status = EXACT_STATUS;
456
    break;
457
  }
458
  __set_status_flags (fpsc, status);
459
 
460
#endif
461
 
462
  return get_BID64 (sign_s, exponent_b + extra_digits, C64,
463
                    rounding_mode, fpsc);
464
}
465
 
466
 
467
///////////////////////////////////////////////////////////////////
468
// round 128-bit coefficient and return result in BID64 format
469
// do not worry about midpoint cases
470
//////////////////////////////////////////////////////////////////
471
static UINT64
472
__bid_simple_round64_sticky (UINT64 sign, int exponent, UINT128 P,
473
                             int extra_digits, int rounding_mode,
474
                             unsigned *fpsc) {
475
  UINT128 Q_high, Q_low, C128;
476
  UINT64 C64;
477
  int amount, rmode;
478
 
479
  rmode = rounding_mode;
480
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
481
#ifndef IEEE_ROUND_NEAREST
482
  if (sign && (unsigned) (rmode - 1) < 2)
483
    rmode = 3 - rmode;
484
#endif
485
#endif
486
  __add_128_64 (P, P, round_const_table[rmode][extra_digits]);
487
 
488
  // get P*(2^M[extra_digits])/10^extra_digits
489
  __mul_128x128_full (Q_high, Q_low, P,
490
                      reciprocals10_128[extra_digits]);
491
 
492
  // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
493
  amount = recip_scale[extra_digits];
494
  __shr_128 (C128, Q_high, amount);
495
 
496
  C64 = __low_64 (C128);
497
 
498
#ifdef SET_STATUS_FLAGS
499
 
500
  __set_status_flags (fpsc, INEXACT_EXCEPTION);
501
 
502
#endif
503
 
504
  return get_BID64 (sign, exponent, C64, rounding_mode, fpsc);
505
}
506
 
507
///////////////////////////////////////////////////////////////////
508
// round 128-bit coefficient and return result in BID64 format
509
///////////////////////////////////////////////////////////////////
510
static UINT64
511
__bid_full_round64 (UINT64 sign, int exponent, UINT128 P,
512
                    int extra_digits, int rounding_mode,
513
                    unsigned *fpsc) {
514
  UINT128 Q_high, Q_low, C128, Stemp, PU;
515
  UINT64 remainder_h, C64, carry, CY;
516
  int amount, amount2, rmode, status = 0;
517
 
518
  if (exponent < 0) {
519
    if (exponent >= -16 && (extra_digits + exponent < 0)) {
520
      extra_digits = -exponent;
521
#ifdef SET_STATUS_FLAGS
522
      if (extra_digits > 0) {
523
        rmode = rounding_mode;
524
        if (sign && (unsigned) (rmode - 1) < 2)
525
          rmode = 3 - rmode;
526
        __add_128_128 (PU, P,
527
                       round_const_table_128[rmode][extra_digits]);
528
        if (__unsigned_compare_gt_128
529
            (power10_table_128[extra_digits + 15], PU))
530
          status = UNDERFLOW_EXCEPTION;
531
      }
532
#endif
533
    }
534
  }
535
 
536
  if (extra_digits > 0) {
537
    exponent += extra_digits;
538
    rmode = rounding_mode;
539
    if (sign && (unsigned) (rmode - 1) < 2)
540
      rmode = 3 - rmode;
541
    __add_128_128 (P, P, round_const_table_128[rmode][extra_digits]);
542
 
543
    // get P*(2^M[extra_digits])/10^extra_digits
544
    __mul_128x128_full (Q_high, Q_low, P,
545
                        reciprocals10_128[extra_digits]);
546
 
547
    // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
548
    amount = recip_scale[extra_digits];
549
    __shr_128_long (C128, Q_high, amount);
550
 
551
    C64 = __low_64 (C128);
552
 
553
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
554
#ifndef IEEE_ROUND_NEAREST
555
    if (rmode == 0)      //ROUNDING_TO_NEAREST
556
#endif
557
      if (C64 & 1) {
558
        // check whether fractional part of initial_P/10^extra_digits 
559
        // is exactly .5
560
 
561
        // get remainder
562
        amount2 = 64 - amount;
563
        remainder_h = 0;
564
        remainder_h--;
565
        remainder_h >>= amount2;
566
        remainder_h = remainder_h & Q_high.w[0];
567
 
568
        if (!remainder_h
569
            && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
570
                || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
571
                    && Q_low.w[0] <
572
                    reciprocals10_128[extra_digits].w[0]))) {
573
          C64--;
574
        }
575
      }
576
#endif
577
 
578
#ifdef SET_STATUS_FLAGS
579
    status |= INEXACT_EXCEPTION;
580
 
581
    // get remainder
582
    remainder_h = Q_high.w[0] << (64 - amount);
583
 
584
    switch (rmode) {
585
    case ROUNDING_TO_NEAREST:
586
    case ROUNDING_TIES_AWAY:
587
      // test whether fractional part is 0
588
      if (remainder_h == 0x8000000000000000ull
589
          && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
590
              || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
591
                  && Q_low.w[0] <
592
                  reciprocals10_128[extra_digits].w[0])))
593
        status = EXACT_STATUS;
594
      break;
595
    case ROUNDING_DOWN:
596
    case ROUNDING_TO_ZERO:
597
      if (!remainder_h
598
          && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
599
              || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
600
                  && Q_low.w[0] <
601
                  reciprocals10_128[extra_digits].w[0])))
602
        status = EXACT_STATUS;
603
      break;
604
    default:
605
      // round up
606
      __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
607
                       reciprocals10_128[extra_digits].w[0]);
608
      __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
609
                          reciprocals10_128[extra_digits].w[1], CY);
610
      if ((remainder_h >> (64 - amount)) + carry >=
611
          (((UINT64) 1) << amount))
612
        status = EXACT_STATUS;
613
    }
614
 
615
    __set_status_flags (fpsc, status);
616
 
617
#endif
618
  } else {
619
    C64 = P.w[0];
620
    if (!C64) {
621
      sign = 0;
622
      if (rounding_mode == ROUNDING_DOWN)
623
        sign = 0x8000000000000000ull;
624
    }
625
  }
626
  return get_BID64 (sign, exponent, C64, rounding_mode, fpsc);
627
}
628
 
629
/////////////////////////////////////////////////////////////////////////////////
630
// round 192-bit coefficient (P, remainder_P) and return result in BID64 format
631
// the lowest 64 bits (remainder_P) are used for midpoint checking only
632
////////////////////////////////////////////////////////////////////////////////
633
static UINT64
634
__bid_full_round64_remainder (UINT64 sign, int exponent, UINT128 P,
635
                              int extra_digits, UINT64 remainder_P,
636
                              int rounding_mode, unsigned *fpsc,
637
                              unsigned uf_status) {
638
  UINT128 Q_high, Q_low, C128, Stemp;
639
  UINT64 remainder_h, C64, carry, CY;
640
  int amount, amount2, rmode, status = uf_status;
641
 
642
  rmode = rounding_mode;
643
  if (sign && (unsigned) (rmode - 1) < 2)
644
    rmode = 3 - rmode;
645
  if (rmode == ROUNDING_UP && remainder_P) {
646
    P.w[0]++;
647
    if (!P.w[0])
648
      P.w[1]++;
649
  }
650
 
651
  if (extra_digits) {
652
    __add_128_64 (P, P, round_const_table[rmode][extra_digits]);
653
 
654
    // get P*(2^M[extra_digits])/10^extra_digits
655
    __mul_128x128_full (Q_high, Q_low, P,
656
                        reciprocals10_128[extra_digits]);
657
 
658
    // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
659
    amount = recip_scale[extra_digits];
660
    __shr_128 (C128, Q_high, amount);
661
 
662
    C64 = __low_64 (C128);
663
 
664
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
665
#ifndef IEEE_ROUND_NEAREST
666
    if (rmode == 0)      //ROUNDING_TO_NEAREST
667
#endif
668
      if (!remainder_P && (C64 & 1)) {
669
        // check whether fractional part of initial_P/10^extra_digits 
670
        // is exactly .5
671
 
672
        // get remainder
673
        amount2 = 64 - amount;
674
        remainder_h = 0;
675
        remainder_h--;
676
        remainder_h >>= amount2;
677
        remainder_h = remainder_h & Q_high.w[0];
678
 
679
        if (!remainder_h
680
            && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
681
                || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
682
                    && Q_low.w[0] <
683
                    reciprocals10_128[extra_digits].w[0]))) {
684
          C64--;
685
        }
686
      }
687
#endif
688
 
689
#ifdef SET_STATUS_FLAGS
690
    status |= INEXACT_EXCEPTION;
691
 
692
    if (!remainder_P) {
693
      // get remainder
694
      remainder_h = Q_high.w[0] << (64 - amount);
695
 
696
      switch (rmode) {
697
      case ROUNDING_TO_NEAREST:
698
      case ROUNDING_TIES_AWAY:
699
        // test whether fractional part is 0
700
        if (remainder_h == 0x8000000000000000ull
701
            && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
702
                || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
703
                    && Q_low.w[0] <
704
                    reciprocals10_128[extra_digits].w[0])))
705
          status = EXACT_STATUS;
706
        break;
707
      case ROUNDING_DOWN:
708
      case ROUNDING_TO_ZERO:
709
        if (!remainder_h
710
            && (Q_low.w[1] < reciprocals10_128[extra_digits].w[1]
711
                || (Q_low.w[1] == reciprocals10_128[extra_digits].w[1]
712
                    && Q_low.w[0] <
713
                    reciprocals10_128[extra_digits].w[0])))
714
          status = EXACT_STATUS;
715
        break;
716
      default:
717
        // round up
718
        __add_carry_out (Stemp.w[0], CY, Q_low.w[0],
719
                         reciprocals10_128[extra_digits].w[0]);
720
        __add_carry_in_out (Stemp.w[1], carry, Q_low.w[1],
721
                            reciprocals10_128[extra_digits].w[1], CY);
722
        if ((remainder_h >> (64 - amount)) + carry >=
723
            (((UINT64) 1) << amount))
724
          status = EXACT_STATUS;
725
      }
726
    }
727
    __set_status_flags (fpsc, status);
728
 
729
#endif
730
  } else {
731
    C64 = P.w[0];
732
#ifdef SET_STATUS_FLAGS
733
    if (remainder_P) {
734
      __set_status_flags (fpsc, uf_status | INEXACT_EXCEPTION);
735
    }
736
#endif
737
  }
738
 
739
  return get_BID64 (sign, exponent + extra_digits, C64, rounding_mode,
740
                    fpsc);
741
}
742
 
743
 
744
///////////////////////////////////////////////////////////////////
745
// get P/10^extra_digits
746
// result fits in 64 bits
747
///////////////////////////////////////////////////////////////////
748
__BID_INLINE__ UINT64
749
__truncate (UINT128 P, int extra_digits)
750
// extra_digits <= 16
751
{
752
  UINT128 Q_high, Q_low, C128;
753
  UINT64 C64;
754
  int amount;
755
 
756
  // get P*(2^M[extra_digits])/10^extra_digits
757
  __mul_128x128_full (Q_high, Q_low, P,
758
                      reciprocals10_128[extra_digits]);
759
 
760
  // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
761
  amount = recip_scale[extra_digits];
762
  __shr_128 (C128, Q_high, amount);
763
 
764
  C64 = __low_64 (C128);
765
 
766
  return C64;
767
}
768
 
769
 
770
///////////////////////////////////////////////////////////////////
771
// return number of decimal digits in 128-bit value X
772
///////////////////////////////////////////////////////////////////
773
__BID_INLINE__ int
774
__get_dec_digits64 (UINT128 X) {
775
  int_double tempx;
776
  int digits_x, bin_expon_cx;
777
 
778
  if (!X.w[1]) {
779
    //--- get number of bits in the coefficients of x and y ---
780
    tempx.d = (double) X.w[0];
781
    bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
782
    // get number of decimal digits in the coeff_x
783
    digits_x = estimate_decimal_digits[bin_expon_cx];
784
    if (X.w[0] >= power10_table_128[digits_x].w[0])
785
      digits_x++;
786
    return digits_x;
787
  }
788
  tempx.d = (double) X.w[1];
789
  bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
790
  // get number of decimal digits in the coeff_x
791
  digits_x = estimate_decimal_digits[bin_expon_cx + 64];
792
  if (__unsigned_compare_ge_128 (X, power10_table_128[digits_x]))
793
    digits_x++;
794
 
795
  return digits_x;
796
}
797
 
798
 
799
////////////////////////////////////////////////////////////////////////////////
800
//
801
// add 64-bit coefficient to 128-bit coefficient, return result in BID64 format
802
//
803
////////////////////////////////////////////////////////////////////////////////
804
__BID_INLINE__ UINT64
805
get_add128 (UINT64 sign_x, int exponent_x, UINT64 coefficient_x,
806
            UINT64 sign_y, int final_exponent_y, UINT128 CY,
807
            int extra_digits, int rounding_mode, unsigned *fpsc) {
808
  UINT128 CY_L, CX, FS, F, CT, ST, T2;
809
  UINT64 CYh, CY0L, T, S, coefficient_y, remainder_y;
810
  SINT64 D = 0;
811
  int_double tempx;
812
  int diff_dec_expon, extra_digits2, exponent_y, status;
813
  int extra_dx, diff_dec2, bin_expon_cx, digits_x, rmode;
814
 
815
  // CY has more than 16 decimal digits
816
 
817
  exponent_y = final_exponent_y - extra_digits;
818
 
819
#ifdef IEEE_ROUND_NEAREST_TIES_AWAY
820
  rounding_mode = 0;
821
#endif
822
#ifdef IEEE_ROUND_NEAREST
823
  rounding_mode = 0;
824
#endif
825
 
826
  if (exponent_x > exponent_y) {
827
    // normalize x
828
    //--- get number of bits in the coefficients of x and y ---
829
    tempx.d = (double) coefficient_x;
830
    bin_expon_cx = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
831
    // get number of decimal digits in the coeff_x
832
    digits_x = estimate_decimal_digits[bin_expon_cx];
833
    if (coefficient_x >= power10_table_128[digits_x].w[0])
834
      digits_x++;
835
 
836
    extra_dx = 16 - digits_x;
837
    coefficient_x *= power10_table_128[extra_dx].w[0];
838
    if ((sign_x ^ sign_y) && (coefficient_x == 1000000000000000ull)) {
839
      extra_dx++;
840
      coefficient_x = 10000000000000000ull;
841
    }
842
    exponent_x -= extra_dx;
843
 
844
    if (exponent_x > exponent_y) {
845
 
846
      // exponent_x > exponent_y
847
      diff_dec_expon = exponent_x - exponent_y;
848
 
849
      if (exponent_x <= final_exponent_y + 1) {
850
        __mul_64x64_to_128 (CX, coefficient_x,
851
                            power10_table_128[diff_dec_expon].w[0]);
852
 
853
        if (sign_x == sign_y) {
854
          __add_128_128 (CT, CY, CX);
855
          if ((exponent_x >
856
               final_exponent_y) /*&& (final_exponent_y>0) */ )
857
            extra_digits++;
858
          if (__unsigned_compare_ge_128
859
              (CT, power10_table_128[16 + extra_digits]))
860
            extra_digits++;
861
        } else {
862
          __sub_128_128 (CT, CY, CX);
863
          if (((SINT64) CT.w[1]) < 0) {
864
            CT.w[0] = 0 - CT.w[0];
865
            CT.w[1] = 0 - CT.w[1];
866
            if (CT.w[0])
867
              CT.w[1]--;
868
            sign_y = sign_x;
869
          } else if (!(CT.w[1] | CT.w[0])) {
870
            sign_y =
871
              (rounding_mode !=
872
               ROUNDING_DOWN) ? 0 : 0x8000000000000000ull;
873
          }
874
          if ((exponent_x + 1 >=
875
               final_exponent_y) /*&& (final_exponent_y>=0) */ ) {
876
            extra_digits = __get_dec_digits64 (CT) - 16;
877
            if (extra_digits <= 0) {
878
              if (!CT.w[0] && rounding_mode == ROUNDING_DOWN)
879
                sign_y = 0x8000000000000000ull;
880
              return get_BID64 (sign_y, exponent_y, CT.w[0],
881
                                rounding_mode, fpsc);
882
            }
883
          } else
884
            if (__unsigned_compare_gt_128
885
                (power10_table_128[15 + extra_digits], CT))
886
            extra_digits--;
887
        }
888
 
889
        return __bid_full_round64 (sign_y, exponent_y, CT, extra_digits,
890
                                   rounding_mode, fpsc);
891
      }
892
      // diff_dec2+extra_digits is the number of digits to eliminate from 
893
      //                           argument CY
894
      diff_dec2 = exponent_x - final_exponent_y;
895
 
896
      if (diff_dec2 >= 17) {
897
#ifndef IEEE_ROUND_NEAREST
898
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
899
        if ((rounding_mode) & 3) {
900
          switch (rounding_mode) {
901
          case ROUNDING_UP:
902
            if (!sign_y) {
903
              D = ((SINT64) (sign_x ^ sign_y)) >> 63;
904
              D = D + D + 1;
905
              coefficient_x += D;
906
            }
907
            break;
908
          case ROUNDING_DOWN:
909
            if (sign_y) {
910
              D = ((SINT64) (sign_x ^ sign_y)) >> 63;
911
              D = D + D + 1;
912
              coefficient_x += D;
913
            }
914
            break;
915
          case ROUNDING_TO_ZERO:
916
            if (sign_y != sign_x) {
917
              D = 0 - 1;
918
              coefficient_x += D;
919
            }
920
            break;
921
          }
922
          if (coefficient_x < 1000000000000000ull) {
923
            coefficient_x -= D;
924
            coefficient_x =
925
              D + (coefficient_x << 1) + (coefficient_x << 3);
926
            exponent_x--;
927
          }
928
        }
929
#endif
930
#endif
931
#ifdef SET_STATUS_FLAGS
932
        if (CY.w[1] | CY.w[0])
933
          __set_status_flags (fpsc, INEXACT_EXCEPTION);
934
#endif
935
        return get_BID64 (sign_x, exponent_x, coefficient_x,
936
                          rounding_mode, fpsc);
937
      }
938
      // here exponent_x <= 16+final_exponent_y
939
 
940
      // truncate CY to 16 dec. digits
941
      CYh = __truncate (CY, extra_digits);
942
 
943
      // get remainder
944
      T = power10_table_128[extra_digits].w[0];
945
      __mul_64x64_to_64 (CY0L, CYh, T);
946
 
947
      remainder_y = CY.w[0] - CY0L;
948
 
949
      // align coeff_x, CYh
950
      __mul_64x64_to_128 (CX, coefficient_x,
951
                          power10_table_128[diff_dec2].w[0]);
952
 
953
      if (sign_x == sign_y) {
954
        __add_128_64 (CT, CX, CYh);
955
        if (__unsigned_compare_ge_128
956
            (CT, power10_table_128[16 + diff_dec2]))
957
          diff_dec2++;
958
      } else {
959
        if (remainder_y)
960
          CYh++;
961
        __sub_128_64 (CT, CX, CYh);
962
        if (__unsigned_compare_gt_128
963
            (power10_table_128[15 + diff_dec2], CT))
964
          diff_dec2--;
965
      }
966
 
967
      return __bid_full_round64_remainder (sign_x, final_exponent_y, CT,
968
                                           diff_dec2, remainder_y,
969
                                           rounding_mode, fpsc, 0);
970
    }
971
  }
972
  // Here (exponent_x <= exponent_y)
973
  {
974
    diff_dec_expon = exponent_y - exponent_x;
975
 
976
    if (diff_dec_expon > MAX_FORMAT_DIGITS) {
977
      rmode = rounding_mode;
978
 
979
      if ((sign_x ^ sign_y)) {
980
        if (!CY.w[0])
981
          CY.w[1]--;
982
        CY.w[0]--;
983
        if (__unsigned_compare_gt_128
984
            (power10_table_128[15 + extra_digits], CY)) {
985
          if (rmode & 3) {
986
            extra_digits--;
987
            final_exponent_y--;
988
          } else {
989
            CY.w[0] = 1000000000000000ull;
990
            CY.w[1] = 0;
991
            extra_digits = 0;
992
          }
993
        }
994
      }
995
      __scale128_10 (CY, CY);
996
      extra_digits++;
997
      CY.w[0] |= 1;
998
 
999
      return __bid_simple_round64_sticky (sign_y, final_exponent_y, CY,
1000
                                          extra_digits, rmode, fpsc);
1001
    }
1002
    // apply sign to coeff_x
1003
    sign_x ^= sign_y;
1004
    sign_x = ((SINT64) sign_x) >> 63;
1005
    CX.w[0] = (coefficient_x + sign_x) ^ sign_x;
1006
    CX.w[1] = sign_x;
1007
 
1008
    // check whether CY (rounded to 16 digits) and CX have 
1009
    //                     any digits in the same position
1010
    diff_dec2 = final_exponent_y - exponent_x;
1011
 
1012
    if (diff_dec2 <= 17) {
1013
      // align CY to 10^ex
1014
      S = power10_table_128[diff_dec_expon].w[0];
1015
      __mul_64x128_short (CY_L, S, CY);
1016
 
1017
      __add_128_128 (ST, CY_L, CX);
1018
      extra_digits2 = __get_dec_digits64 (ST) - 16;
1019
      return __bid_full_round64 (sign_y, exponent_x, ST, extra_digits2,
1020
                                 rounding_mode, fpsc);
1021
    }
1022
    // truncate CY to 16 dec. digits
1023
    CYh = __truncate (CY, extra_digits);
1024
 
1025
    // get remainder
1026
    T = power10_table_128[extra_digits].w[0];
1027
    __mul_64x64_to_64 (CY0L, CYh, T);
1028
 
1029
    coefficient_y = CY.w[0] - CY0L;
1030
    // add rounding constant
1031
    rmode = rounding_mode;
1032
    if (sign_y && (unsigned) (rmode - 1) < 2)
1033
      rmode = 3 - rmode;
1034
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1035
#ifndef IEEE_ROUND_NEAREST
1036
    if (!(rmode & 3))   //ROUNDING_TO_NEAREST
1037
#endif
1038
#endif
1039
    {
1040
      coefficient_y += round_const_table[rmode][extra_digits];
1041
    }
1042
    // align coefficient_y,  coefficient_x
1043
    S = power10_table_128[diff_dec_expon].w[0];
1044
    __mul_64x64_to_128 (F, coefficient_y, S);
1045
 
1046
    // fraction
1047
    __add_128_128 (FS, F, CX);
1048
 
1049
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1050
#ifndef IEEE_ROUND_NEAREST
1051
    if (rmode == 0)      //ROUNDING_TO_NEAREST
1052
#endif
1053
    {
1054
      // rounding code, here RN_EVEN
1055
      // 10^(extra_digits+diff_dec_expon)
1056
      T2 = power10_table_128[diff_dec_expon + extra_digits];
1057
      if (__unsigned_compare_gt_128 (FS, T2)
1058
          || ((CYh & 1) && __test_equal_128 (FS, T2))) {
1059
        CYh++;
1060
        __sub_128_128 (FS, FS, T2);
1061
      }
1062
    }
1063
#endif
1064
#ifndef IEEE_ROUND_NEAREST
1065
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1066
    if (rmode == 4)     //ROUNDING_TO_NEAREST
1067
#endif
1068
    {
1069
      // rounding code, here RN_AWAY
1070
      // 10^(extra_digits+diff_dec_expon)
1071
      T2 = power10_table_128[diff_dec_expon + extra_digits];
1072
      if (__unsigned_compare_ge_128 (FS, T2)) {
1073
        CYh++;
1074
        __sub_128_128 (FS, FS, T2);
1075
      }
1076
    }
1077
#endif
1078
#ifndef IEEE_ROUND_NEAREST
1079
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1080
    switch (rmode) {
1081
    case ROUNDING_DOWN:
1082
    case ROUNDING_TO_ZERO:
1083
      if ((SINT64) FS.w[1] < 0) {
1084
        CYh--;
1085
        if (CYh < 1000000000000000ull) {
1086
          CYh = 9999999999999999ull;
1087
          final_exponent_y--;
1088
        }
1089
      } else {
1090
        T2 = power10_table_128[diff_dec_expon + extra_digits];
1091
        if (__unsigned_compare_ge_128 (FS, T2)) {
1092
          CYh++;
1093
          __sub_128_128 (FS, FS, T2);
1094
        }
1095
      }
1096
      break;
1097
    case ROUNDING_UP:
1098
      if ((SINT64) FS.w[1] < 0)
1099
        break;
1100
      T2 = power10_table_128[diff_dec_expon + extra_digits];
1101
      if (__unsigned_compare_gt_128 (FS, T2)) {
1102
        CYh += 2;
1103
        __sub_128_128 (FS, FS, T2);
1104
      } else if ((FS.w[1] == T2.w[1]) && (FS.w[0] == T2.w[0])) {
1105
        CYh++;
1106
        FS.w[1] = FS.w[0] = 0;
1107
      } else if (FS.w[1] | FS.w[0])
1108
        CYh++;
1109
      break;
1110
    }
1111
#endif
1112
#endif
1113
 
1114
#ifdef SET_STATUS_FLAGS
1115
    status = INEXACT_EXCEPTION;
1116
#ifndef IEEE_ROUND_NEAREST
1117
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1118
    if (!(rmode & 3))
1119
#endif
1120
#endif
1121
    {
1122
      // RN modes
1123
      if ((FS.w[1] ==
1124
           round_const_table_128[0][diff_dec_expon + extra_digits].w[1])
1125
          && (FS.w[0] ==
1126
              round_const_table_128[0][diff_dec_expon +
1127
                                       extra_digits].w[0]))
1128
        status = EXACT_STATUS;
1129
    }
1130
#ifndef IEEE_ROUND_NEAREST
1131
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1132
    else if (!FS.w[1] && !FS.w[0])
1133
      status = EXACT_STATUS;
1134
#endif
1135
#endif
1136
 
1137
    __set_status_flags (fpsc, status);
1138
#endif
1139
 
1140
    return get_BID64 (sign_y, final_exponent_y, CYh, rounding_mode,
1141
                      fpsc);
1142
  }
1143
 
1144
}
1145
 
1146
//////////////////////////////////////////////////////////////////////////
1147
//
1148
//  If coefficient_z is less than 16 digits long, normalize to 16 digits
1149
//
1150
/////////////////////////////////////////////////////////////////////////
1151
static UINT64
1152
BID_normalize (UINT64 sign_z, int exponent_z,
1153
               UINT64 coefficient_z, UINT64 round_dir, int round_flag,
1154
               int rounding_mode, unsigned *fpsc) {
1155
  SINT64 D;
1156
  int_double tempx;
1157
  int digits_z, bin_expon, scale, rmode;
1158
 
1159
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1160
#ifndef IEEE_ROUND_NEAREST
1161
  rmode = rounding_mode;
1162
  if (sign_z && (unsigned) (rmode - 1) < 2)
1163
    rmode = 3 - rmode;
1164
#else
1165
  if (coefficient_z >= power10_table_128[15].w[0])
1166
    return z;
1167
#endif
1168
#endif
1169
 
1170
  //--- get number of bits in the coefficients of x and y ---
1171
  tempx.d = (double) coefficient_z;
1172
  bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
1173
  // get number of decimal digits in the coeff_x
1174
  digits_z = estimate_decimal_digits[bin_expon];
1175
  if (coefficient_z >= power10_table_128[digits_z].w[0])
1176
    digits_z++;
1177
 
1178
  scale = 16 - digits_z;
1179
  exponent_z -= scale;
1180
  if (exponent_z < 0) {
1181
    scale += exponent_z;
1182
    exponent_z = 0;
1183
  }
1184
  coefficient_z *= power10_table_128[scale].w[0];
1185
 
1186
#ifdef SET_STATUS_FLAGS
1187
  if (round_flag) {
1188
    __set_status_flags (fpsc, INEXACT_EXCEPTION);
1189
    if (coefficient_z < 1000000000000000ull)
1190
      __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
1191
    else if ((coefficient_z == 1000000000000000ull) && !exponent_z
1192
             && ((SINT64) (round_dir ^ sign_z) < 0) && round_flag
1193
             && (rmode == ROUNDING_DOWN || rmode == ROUNDING_TO_ZERO))
1194
      __set_status_flags (fpsc, UNDERFLOW_EXCEPTION);
1195
  }
1196
#endif
1197
 
1198
#ifndef IEEE_ROUND_NEAREST_TIES_AWAY
1199
#ifndef IEEE_ROUND_NEAREST
1200
  if (round_flag && (rmode & 3)) {
1201
    D = round_dir ^ sign_z;
1202
 
1203
    if (rmode == ROUNDING_UP) {
1204
      if (D >= 0)
1205
        coefficient_z++;
1206
    } else {
1207
      if (D < 0)
1208
        coefficient_z--;
1209
      if (coefficient_z < 1000000000000000ull && exponent_z) {
1210
        coefficient_z = 9999999999999999ull;
1211
        exponent_z--;
1212
      }
1213
    }
1214
  }
1215
#endif
1216
#endif
1217
 
1218
  return get_BID64 (sign_z, exponent_z, coefficient_z, rounding_mode,
1219
                    fpsc);
1220
}
1221
 
1222
 
1223
//////////////////////////////////////////////////////////////////////////
1224
//
1225
//    0*10^ey + cz*10^ez,   ey<ez  
1226
//
1227
//////////////////////////////////////////////////////////////////////////
1228
 
1229
__BID_INLINE__ UINT64
1230
add_zero64 (int exponent_y, UINT64 sign_z, int exponent_z,
1231
            UINT64 coefficient_z, unsigned *prounding_mode,
1232
            unsigned *fpsc) {
1233
  int_double tempx;
1234
  int bin_expon, scale_k, scale_cz;
1235
  int diff_expon;
1236
 
1237
  diff_expon = exponent_z - exponent_y;
1238
 
1239
  tempx.d = (double) coefficient_z;
1240
  bin_expon = ((tempx.i & MASK_BINARY_EXPONENT) >> 52) - 0x3ff;
1241
  scale_cz = estimate_decimal_digits[bin_expon];
1242
  if (coefficient_z >= power10_table_128[scale_cz].w[0])
1243
    scale_cz++;
1244
 
1245
  scale_k = 16 - scale_cz;
1246
  if (diff_expon < scale_k)
1247
    scale_k = diff_expon;
1248
  coefficient_z *= power10_table_128[scale_k].w[0];
1249
 
1250
  return get_BID64 (sign_z, exponent_z - scale_k, coefficient_z,
1251
                    *prounding_mode, fpsc);
1252
}
1253
#endif

powered by: WebSVN 2.1.0

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