libgcc/config/libbid/bid64_rem.c

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