libgcc/config/libbid/bid128_div.c

   1 /* Copyright (C) 2007  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 2, or (at your option) any later
   8 version.
   9
  10 In addition to the permissions in the GNU General Public License, the
  11 Free Software Foundation gives you unlimited permission to link the
  12 compiled version of this file into combinations with other programs,
  13 and to distribute those combinations without any restriction coming
  14 from the use of this file.  (The General Public License restrictions
  15 do apply in other respects; for example, they cover modification of
  16 the file, and distribution when not linked into a combine
  17 executable.)
  18
  19 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
  20 WARRANTY; without even the implied warranty of MERCHANTABILITY or
  21 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  22 for more details.
  23
  24 You should have received a copy of the GNU General Public License
  25 along with GCC; see the file COPYING.  If not, write to the Free
  26 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  27 02110-1301, USA.  */
  28
  29 #define BID_128RES
  30 #include "div_macros.h"
  31
  32 extern UINT32 __bid_convert_table[5][128][2];
  33 extern SINT8 __bid_factors[][2];
  34 extern UINT8 __bid_packed_10000_zeros[];
  35
  36 BID128_FUNCTION_ARG2(__bid128_div, x, y)
  37
  38   UINT256 CA4, CA4r, P256;
  39   UINT128 CX, CY, T128, CQ, CR, CA, TP128, Qh, Ql, res;
  40   UINT64 sign_x, sign_y, T, carry64, D, Q_high, Q_low, QX, X, PD;
  41   int_float fx, fy, f64;
  42   UINT32 QX32, tdigit[3], digit, digit_h, digit_low;
  43   int exponent_x = 0, exponent_y, bin_index, bin_expon, diff_expon, ed2,
  44     digits_q, amount;
  45   int nzeros, i, j, k, d5;
  46   unsigned rmode;
  47
  48
  49   // unpack arguments, check for NaN or Infinity
  50   if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
  51     // test if x is NaN
  52     if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
  53 #ifdef SET_STATUS_FLAGS
  54       if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull || // sNaN
  55           (y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull)
  56         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  57 #endif
  58       res.w[1] = (x.w[1]) & QUIET_MASK64;
  59       res.w[0] = x.w[0];
  60       BID_RETURN (res);
  61     }
  62     // x is Infinity?
  63     if ((x.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
  64       // check if y is Inf.
  65       if (((y.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull))
  66         // return NaN
  67       {
  68 #ifdef SET_STATUS_FLAGS
  69         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  70 #endif
  71         res.w[1] = 0x7c00000000000000ull;
  72         res.w[0] = 0;
  73         BID_RETURN (res);
  74       }
  75       // y is NaN?
  76       if (((y.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull))
  77         // return NaN
  78       {
  79         // return +/-Inf
  80         res.w[1] = ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) |
  81             0x7800000000000000ull;
  82         res.w[0] = 0;
  83         BID_RETURN (res);
  84       }
  85     }
  86     // x is 0
  87     if ((y.w[1] & 0x7800000000000000ull) < 0x7800000000000000ull) {
  88       if ((!y.w[0]) && !(y.w[1] & 0x0001ffffffffffffull)) {
  89 #ifdef SET_STATUS_FLAGS
  90         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  91 #endif
  92         // x=y=0, return NaN
  93         res.w[1] = 0x7c00000000000000ull;
  94         res.w[0] = 0;
  95         BID_RETURN (res);
  96       }
  97       // return 0
  98       res.w[1] = (x.w[1] ^ y.w[1]) & 0x8000000000000000ull;
  99       X = ((y.w[1]) << 1) >> 50;
 100       exponent_x = exponent_x - (int) X + DECIMAL_EXPONENT_BIAS_128;
 101       if (exponent_x > DECIMAL_MAX_EXPON_128)
 102         exponent_x = DECIMAL_MAX_EXPON_128;
 103       else if (exponent_x < 0)
 104         exponent_x = 0;
 105       res.w[1] |= (((UINT64) exponent_x) << 49);
 106       res.w[0] = 0;
 107       BID_RETURN (res);
 108     }
 109   }
 110   if (!unpack_BID128_value (&sign_y, &exponent_y, &CY, y)) {
 111     // y is Inf. or NaN
 112
 113     // test if y is NaN
 114     if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
 115 #ifdef SET_STATUS_FLAGS
 116       if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) // sNaN
 117         __set_status_flags (pfpsf, INVALID_EXCEPTION);
 118 #endif
 119       res.w[1] = y.w[1] & QUIET_MASK64;
 120       res.w[0] = y.w[0];
 121       BID_RETURN (res);
 122     }
 123     // y is Infinity?
 124     if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
 125       // return +/-0
 126       res.w[1] = sign_x ^ sign_y;
 127       res.w[0] = 0;
 128       BID_RETURN (res);
 129     }
 130     // y is 0, return +/-Inf
 131 #ifdef SET_STATUS_FLAGS
 132     __set_status_flags (pfpsf, ZERO_DIVIDE_EXCEPTION);
 133 #endif
 134     res.w[1] =
 135       ((x.w[1] ^ y.w[1]) & 0x8000000000000000ull) | 0x7800000000000000ull;
 136     res.w[0] = 0;
 137     BID_RETURN (res);
 138   }
 139   diff_expon = exponent_x - exponent_y + DECIMAL_EXPONENT_BIAS_128;
 140
 141   if (__unsigned_compare_gt_128 (CY, CX)) {
 142     // CX < CY
 143
 144     // 2^64
 145     f64.i = 0x5f800000;
 146
 147     // fx ~ CX,   fy ~ CY
 148     fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
 149     fy.d = (float) CY.w[1] * f64.d + (float) CY.w[0];
 150     // expon_cy - expon_cx
 151     bin_index = (fy.i - fx.i) >> 23;
 152
 153     if (CX.w[1]) {
 154       T = __bid_power10_index_binexp_128[bin_index].w[0];
 155       __mul_64x128_short (CA, T, CX);
 156     } else {
 157       T128 = __bid_power10_index_binexp_128[bin_index];
 158       __mul_64x128_short (CA, CX.w[0], T128);
 159     }
 160
 161     ed2 = 33;
 162     if (__unsigned_compare_gt_128 (CY, CA))
 163       ed2++;
 164
 165     T128 = __bid_power10_table_128[ed2];
 166     __mul_128x128_to_256 (CA4, CA, T128);
 167
 168     ed2 += __bid_estimate_decimal_digits[bin_index];
 169     CQ.w[0] = CQ.w[1] = 0;
 170     diff_expon = diff_expon - ed2;
 171
 172   } else {
 173     // get CQ = CX/CY
 174     __div_128_by_128 (&CQ, &CR, CX, CY);
 175
 176     if (!CR.w[1] && !CR.w[0]) {
 177       get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
 178                   pfpsf);
 179       BID_RETURN (res);
 180     }
 181     // get number of decimal digits in CQ
 182     // 2^64
 183     f64.i = 0x5f800000;
 184     fx.d = (float) CQ.w[1] * f64.d + (float) CQ.w[0];
 185     // binary expon. of CQ
 186     bin_expon = (fx.i - 0x3f800000) >> 23;
 187
 188     digits_q = __bid_estimate_decimal_digits[bin_expon];
 189     TP128.w[0] = __bid_power10_index_binexp_128[bin_expon].w[0];
 190     TP128.w[1] = __bid_power10_index_binexp_128[bin_expon].w[1];
 191     if (__unsigned_compare_ge_128 (CQ, TP128))
 192       digits_q++;
 193
 194     ed2 = 34 - digits_q;
 195     T128.w[0] = __bid_power10_table_128[ed2].w[0];
 196     T128.w[1] = __bid_power10_table_128[ed2].w[1];
 197     __mul_128x128_to_256 (CA4, CR, T128);
 198     diff_expon = diff_expon - ed2;
 199     __mul_128x128_low (CQ, CQ, T128);
 200
 201   }
 202
 203   __div_256_by_128 (&CQ, &CA4, CY);
 204
 205 #ifdef SET_STATUS_FLAGS
 206   if (CA4.w[0] || CA4.w[1]) {
 207     // set status flags
 208     __set_status_flags (pfpsf, INEXACT_EXCEPTION);
 209   }
 210 #ifndef LEAVE_TRAILING_ZEROS
 211   else
 212 #endif
 213 #else
 214 #ifndef LEAVE_TRAILING_ZEROS
 215   if (!CA4.w[0] && !CA4.w[1])
 216 #endif
 217 #endif
 218 #ifndef LEAVE_TRAILING_ZEROS
 219     // check whether result is exact
 220   {
 221     // check whether CX, CY are short
 222     if (!CX.w[1] && !CY.w[1] && (CX.w[0] <= 1024) && (CY.w[0] <= 1024)) {
 223       i = (int) CY.w[0] - 1;
 224       j = (int) CX.w[0] - 1;
 225       // difference in powers of 2 __bid_factors for Y and X
 226       nzeros = ed2 - __bid_factors[i][0] + __bid_factors[j][0];
 227       // difference in powers of 5 __bid_factors
 228       d5 = ed2 - __bid_factors[i][1] + __bid_factors[j][1];
 229       if (d5 < nzeros)
 230         nzeros = d5;
 231       // get P*(2^M[extra_digits])/10^extra_digits
 232       __mul_128x128_full (Qh, Ql, CQ, __bid_reciprocals10_128[nzeros]);
 233
 234       // now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
 235       amount = __bid_recip_scale[nzeros];
 236       __shr_128_long (CQ, Qh, amount);
 237
 238       diff_expon += nzeros;
 239     } else {
 240       // decompose Q as Qh*10^17 + Ql
 241       //T128 = __bid_reciprocals10_128[17];
 242       T128.w[0] = 0x44909befeb9fad49ull;
 243       T128.w[1] = 0x000b877aa3236a4bull;
 244       __mul_128x128_to_256 (P256, CQ, T128);
 245       //amount = __bid_recip_scale[17];
 246       Q_high = (P256.w[2] >> 44) | (P256.w[3] << (64 - 44));
 247       Q_low = CQ.w[0] - Q_high * 100000000000000000ull;
 248
 249       if (!Q_low) {
 250         diff_expon += 17;
 251
 252         tdigit[0] = Q_high & 0x3ffffff;
 253         tdigit[1] = 0;
 254         QX = Q_high >> 26;
 255         QX32 = QX;
 256         nzeros = 0;
 257
 258         for (j = 0; QX32; j++, QX32 >>= 7) {
 259           k = (QX32 & 127);
 260           tdigit[0] += __bid_convert_table[j][k][0];
 261           tdigit[1] += __bid_convert_table[j][k][1];
 262           if (tdigit[0] >= 100000000) {
 263             tdigit[0] -= 100000000;
 264             tdigit[1]++;
 265           }
 266         }
 267
 268         if (tdigit[1] >= 100000000) {
 269           tdigit[1] -= 100000000;
 270           if (tdigit[1] >= 100000000)
 271             tdigit[1] -= 100000000;
 272         }
 273
 274         digit = tdigit[0];
 275         if (!digit && !tdigit[1])
 276           nzeros += 16;
 277         else {
 278           if (!digit) {
 279             nzeros += 8;
 280             digit = tdigit[1];
 281           }
 282           // decompose digit
 283           PD = (UINT64) digit *0x068DB8BBull;
 284           digit_h = (UINT32) (PD >> 40);
 285           digit_low = digit - digit_h * 10000;
 286
 287           if (!digit_low)
 288             nzeros += 4;
 289           else
 290             digit_h = digit_low;
 291
 292           if (!(digit_h & 1))
 293             nzeros +=
 294               3 & (UINT32) (__bid_packed_10000_zeros[digit_h >> 3] >>
 295                             (digit_h & 7));
 296         }
 297
 298         if (nzeros) {
 299           __mul_64x64_to_128 (CQ, Q_high, __bid_reciprocals10_64[nzeros]);
 300
 301           // now get P/10^extra_digits: shift C64 right by M[extra_digits]-64
 302           amount = __bid_short_recip_scale[nzeros];
 303           CQ.w[0] = CQ.w[1] >> amount;
 304         } else
 305           CQ.w[0] = Q_high;
 306         CQ.w[1] = 0;
 307
 308         diff_expon += nzeros;
 309       } else {
 310         tdigit[0] = Q_low & 0x3ffffff;
 311         tdigit[1] = 0;
 312         QX = Q_low >> 26;
 313         QX32 = QX;
 314         nzeros = 0;
 315
 316         for (j = 0; QX32; j++, QX32 >>= 7) {
 317           k = (QX32 & 127);
 318           tdigit[0] += __bid_convert_table[j][k][0];
 319           tdigit[1] += __bid_convert_table[j][k][1];
 320           if (tdigit[0] >= 100000000) {
 321             tdigit[0] -= 100000000;
 322             tdigit[1]++;
 323           }
 324         }
 325
 326         if (tdigit[1] >= 100000000) {
 327           tdigit[1] -= 100000000;
 328           if (tdigit[1] >= 100000000)
 329             tdigit[1] -= 100000000;
 330         }
 331
 332         digit = tdigit[0];
 333         if (!digit && !tdigit[1])
 334           nzeros += 16;
 335         else {
 336           if (!digit) {
 337             nzeros += 8;
 338             digit = tdigit[1];
 339           }
 340           // decompose digit
 341           PD = (UINT64) digit *0x068DB8BBull;
 342           digit_h = (UINT32) (PD >> 40);
 343           digit_low = digit - digit_h * 10000;
 344
 345           if (!digit_low)
 346             nzeros += 4;
 347           else
 348             digit_h = digit_low;
 349
 350           if (!(digit_h & 1))
 351             nzeros +=
 352               3 & (UINT32) (__bid_packed_10000_zeros[digit_h >> 3] >>
 353                             (digit_h & 7));
 354         }
 355
 356         if (nzeros) {
 357           // get P*(2^M[extra_digits])/10^extra_digits
 358           __mul_128x128_full (Qh, Ql, CQ, __bid_reciprocals10_128[nzeros]);
 359
 360           //now get P/10^extra_digits: shift Q_high right by M[extra_digits]-128
 361           amount = __bid_recip_scale[nzeros];
 362           __shr_128 (CQ, Qh, amount);
 363         }
 364         diff_expon += nzeros;
 365
 366       }
 367     }
 368     get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode,
 369                 pfpsf);
 370     BID_RETURN (res);
 371   }
 372 #endif
 373
 374   if (diff_expon >= 0) {
 375 #ifdef IEEE_ROUND_NEAREST
 376     // rounding
 377     // 2*CA4 - CY
 378     CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
 379     CA4r.w[0] = CA4.w[0] + CA4.w[0];
 380     __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
 381     CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
 382
 383     D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
 384     carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
 385
 386     CQ.w[0] += carry64;
 387     if (CQ.w[0] < carry64)
 388       CQ.w[1]++;
 389 #else
 390 #ifdef IEEE_ROUND_NEAREST_TIES_AWAY
 391     // rounding
 392     // 2*CA4 - CY
 393     CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
 394     CA4r.w[0] = CA4.w[0] + CA4.w[0];
 395     __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
 396     CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
 397
 398     D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
 399     carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
 400
 401     CQ.w[0] += carry64;
 402     if (CQ.w[0] < carry64)
 403       CQ.w[1]++;
 404 #else
 405     rmode = rnd_mode;
 406     if (sign_x ^ sign_y && (unsigned) (rmode - 1) < 2)
 407       rmode = 3 - rmode;
 408     switch (rmode) {
 409     case ROUNDING_TO_NEAREST: // round to nearest code
 410       // rounding
 411       // 2*CA4 - CY
 412       CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
 413       CA4r.w[0] = CA4.w[0] + CA4.w[0];
 414       __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
 415       CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
 416       D = (CA4r.w[1] | CA4r.w[0]) ? 1 : 0;
 417       carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) & ((CQ.w[0]) | D);
 418       CQ.w[0] += carry64;
 419       if (CQ.w[0] < carry64)
 420         CQ.w[1]++;
 421       break;
 422     case ROUNDING_TIES_AWAY:
 423       // rounding
 424       // 2*CA4 - CY
 425       CA4r.w[1] = (CA4.w[1] + CA4.w[1]) | (CA4.w[0] >> 63);
 426       CA4r.w[0] = CA4.w[0] + CA4.w[0];
 427       __sub_borrow_out (CA4r.w[0], carry64, CA4r.w[0], CY.w[0]);
 428       CA4r.w[1] = CA4r.w[1] - CY.w[1] - carry64;
 429       D = (CA4r.w[1] | CA4r.w[0]) ? 0 : 1;
 430       carry64 = (1 + (((SINT64) CA4r.w[1]) >> 63)) | D;
 431       CQ.w[0] += carry64;
 432       if (CQ.w[0] < carry64)
 433         CQ.w[1]++;
 434       break;
 435     case ROUNDING_DOWN:
 436     case ROUNDING_TO_ZERO:
 437       break;
 438     default: // rounding up
 439       CQ.w[0]++;
 440       if (!CQ.w[0])
 441         CQ.w[1]++;
 442       break;
 443     }
 444 #endif
 445 #endif
 446
 447   } else {
 448 #ifdef SET_STATUS_FLAGS
 449     if (CA4.w[0] || CA4.w[1]) {
 450       // set status flags
 451       __set_status_flags (pfpsf, INEXACT_EXCEPTION);
 452     }
 453 #endif
 454
 455     handle_UF_128_rem (&res, sign_x ^ sign_y, diff_expon, CQ,
 456                        CA4.w[1] | CA4.w[0], &rnd_mode, pfpsf);
 457     BID_RETURN (res);
 458
 459   }
 460
 461   get_BID128 (&res, sign_x ^ sign_y, diff_expon, CQ, &rnd_mode, pfpsf);
 462   BID_RETURN (res);
 463
 464 }