libgcc/config/libbid/bid64_rem.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 /*****************************************************************************
  30  *    BID64 remainder
  31  *****************************************************************************
  32  *
  33  *  Algorithm description:
  34  *
  35  *  if(exponent_x < exponent_y)
  36  *    scale coefficient_y so exponents are aligned
  37  *    perform coefficient divide (64-bit integer divide), unless
  38  *            coefficient_y is longer than 64 bits (clearly larger
  39  *                                               than coefficient_x)
  40  *  else  // exponent_x > exponent_y
  41  *     use a loop to scale coefficient_x to 18_digits, divide by
  42  *         coefficient_y (64-bit integer divide), calculate remainder
  43  *         as new_coefficient_x and repeat until final remainder is obtained
  44  *         (when new_exponent_x < exponent_y)
  45  *
  46  ****************************************************************************/
  47
  48 #include "bid_internal.h"
  49
  50 #define MAX_FORMAT_DIGITS     16
  51 #define DECIMAL_EXPONENT_BIAS 398
  52 #define MASK_BINARY_EXPONENT  0x7ff0000000000000ull
  53 #define BINARY_EXPONENT_BIAS  0x3ff
  54 #define UPPER_EXPON_LIMIT     51
  55
  56 #if DECIMAL_CALL_BY_REFERENCE
  57
  58 void
  59 bid64_rem (UINT64 * pres, UINT64 * px,
  60            UINT64 *
  61            py _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  62   UINT64 x, y;
  63 #else
  64
  65 UINT64
  66 bid64_rem (UINT64 x,
  67            UINT64 y _EXC_FLAGS_PARAM _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  68 #endif
  69   UINT128 CY;
  70   UINT64 sign_x, sign_y, coefficient_x, coefficient_y, res;
  71   UINT64 Q, R, R2, T, valid_y, valid_x;
  72   int_float tempx;
  73   int exponent_x, exponent_y, bin_expon, e_scale;
  74   int digits_x, diff_expon;
  75
  76 #if DECIMAL_CALL_BY_REFERENCE
  77   x = *px;
  78   y = *py;
  79 #endif
  80
  81   valid_y = unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y);
  82   valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
  83
  84   // unpack arguments, check for NaN or Infinity
  85   if (!valid_x) {
  86     // x is Inf. or NaN or 0
  87 #ifdef SET_STATUS_FLAGS
  88     if ((y & SNAN_MASK64) == SNAN_MASK64)       // y is sNaN
  89       __set_status_flags (pfpsf, INVALID_EXCEPTION);
  90 #endif
  91
  92     // test if x is NaN
  93     if ((x & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
  94 #ifdef SET_STATUS_FLAGS
  95       if (((x & SNAN_MASK64) == SNAN_MASK64))
  96         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  97 #endif
  98       res = coefficient_x & QUIET_MASK64;;
  99       BID_RETURN (res);
 100     }
 101     // x is Infinity?
 102     if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
 103       if (((y & NAN_MASK64) != NAN_MASK64)) {
 104 #ifdef SET_STATUS_FLAGS
 105         __set_status_flags (pfpsf, INVALID_EXCEPTION);
 106 #endif
 107         // return NaN
 108         res = 0x7c00000000000000ull;
 109         BID_RETURN (res);
 110       }
 111     }
 112     // x is 0
 113     // return x if y != 0
 114     if (((y & 0x7800000000000000ull) < 0x7800000000000000ull) &&
 115         coefficient_y) {
 116       if ((y & 0x6000000000000000ull) == 0x6000000000000000ull)
 117         exponent_y = (y >> 51) & 0x3ff;
 118       else
 119         exponent_y = (y >> 53) & 0x3ff;
 120
 121       if (exponent_y < exponent_x)
 122         exponent_x = exponent_y;
 123
 124       x = exponent_x;
 125       x <<= 53;
 126
 127       res = x | sign_x;
 128       BID_RETURN (res);
 129     }
 130
 131   }
 132   if (!valid_y) {
 133     // y is Inf. or NaN
 134
 135     // test if y is NaN
 136     if ((y & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
 137 #ifdef SET_STATUS_FLAGS
 138       if (((y & SNAN_MASK64) == SNAN_MASK64))
 139         __set_status_flags (pfpsf, INVALID_EXCEPTION);
 140 #endif
 141       res = coefficient_y & QUIET_MASK64;;
 142       BID_RETURN (res);
 143     }
 144     // y is Infinity?
 145     if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) {
 146       res = very_fast_get_BID64 (sign_x, exponent_x, coefficient_x);
 147       BID_RETURN (res);
 148     }
 149     // y is 0, return NaN
 150     {
 151 #ifdef SET_STATUS_FLAGS
 152       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 153 #endif
 154       res = 0x7c00000000000000ull;
 155       BID_RETURN (res);
 156     }
 157   }
 158
 159
 160   diff_expon = exponent_x - exponent_y;
 161   if (diff_expon <= 0) {
 162     diff_expon = -diff_expon;
 163
 164     if (diff_expon > 16) {
 165       // |x|<|y| in this case
 166       res = x;
 167       BID_RETURN (res);
 168     }
 169     // set exponent of y to exponent_x, scale coefficient_y
 170     T = power10_table_128[diff_expon].w[0];
 171     __mul_64x64_to_128 (CY, coefficient_y, T);
 172
 173     if (CY.w[1] || CY.w[0] > (coefficient_x << 1)) {
 174       res = x;
 175       BID_RETURN (res);
 176     }
 177
 178     Q = coefficient_x / CY.w[0];
 179     R = coefficient_x - Q * CY.w[0];
 180
 181     R2 = R + R;
 182     if (R2 > CY.w[0] || (R2 == CY.w[0] && (Q & 1))) {
 183       R = CY.w[0] - R;
 184       sign_x ^= 0x8000000000000000ull;
 185     }
 186
 187     res = very_fast_get_BID64 (sign_x, exponent_x, R);
 188     BID_RETURN (res);
 189   }
 190
 191
 192   while (diff_expon > 0) {
 193     // get number of digits in coeff_x
 194     tempx.d = (float) coefficient_x;
 195     bin_expon = ((tempx.i >> 23) & 0xff) - 0x7f;
 196     digits_x = estimate_decimal_digits[bin_expon];
 197     // will not use this test, dividend will have 18 or 19 digits
 198     //if(coefficient_x >= power10_table_128[digits_x].w[0])
 199     //      digits_x++;
 200
 201     e_scale = 18 - digits_x;
 202     if (diff_expon >= e_scale) {
 203       diff_expon -= e_scale;
 204     } else {
 205       e_scale = diff_expon;
 206       diff_expon = 0;
 207     }
 208
 209     // scale dividend to 18 or 19 digits
 210     coefficient_x *= power10_table_128[e_scale].w[0];
 211
 212     // quotient
 213     Q = coefficient_x / coefficient_y;
 214     // remainder
 215     coefficient_x -= Q * coefficient_y;
 216
 217     // check for remainder == 0
 218     if (!coefficient_x) {
 219       res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, 0);
 220       BID_RETURN (res);
 221     }
 222   }
 223
 224   R2 = coefficient_x + coefficient_x;
 225   if (R2 > coefficient_y || (R2 == coefficient_y && (Q & 1))) {
 226     coefficient_x = coefficient_y - coefficient_x;
 227     sign_x ^= 0x8000000000000000ull;
 228   }
 229
 230   res = very_fast_get_BID64 (sign_x, exponent_y, coefficient_x);
 231   BID_RETURN (res);
 232
 233 }