libgcc/config/libbid/bid128_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 #define BID_128RES
  30 #include "bid_div_macros.h"
  31
  32
  33 BID128_FUNCTION_ARG2_NORND_CUSTOMRESTYPE (UINT128, bid128_rem, x, y)
  34
  35      UINT256 P256;
  36      UINT128 CX, CY, CX2, CQ, CR, T, CXS, P128, res;
  37      UINT64 sign_x, sign_y, valid_y;
  38      SINT64 D;
  39      int_float f64, fx;
  40      int exponent_x, exponent_y, diff_expon, bin_expon_cx, scale,
  41        scale0;
  42
  43   // unpack arguments, check for NaN or Infinity
  44
  45 valid_y = unpack_BID128_value (&sign_y, &exponent_y, &CY, y);
  46
  47 if (!unpack_BID128_value (&sign_x, &exponent_x, &CX, x)) {
  48 #ifdef SET_STATUS_FLAGS
  49 if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64)      // y is sNaN
  50   __set_status_flags (pfpsf, INVALID_EXCEPTION);
  51 #endif
  52     // test if x is NaN
  53 if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
  54 #ifdef SET_STATUS_FLAGS
  55   if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64)    // y is sNaN
  56     __set_status_flags (pfpsf, INVALID_EXCEPTION);
  57 #endif
  58   res.w[1] = CX.w[1] & QUIET_MASK64;
  59   res.w[0] = CX.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) != 0x7c00000000000000ull))
  66     // return NaN
  67   {
  68 #ifdef SET_STATUS_FLAGS
  69     // set status flags
  70     __set_status_flags (pfpsf, INVALID_EXCEPTION);
  71 #endif
  72     res.w[1] = 0x7c00000000000000ull;
  73     res.w[0] = 0;
  74     BID_RETURN (res);
  75   }
  76
  77 }
  78     // x is 0
  79 if ((!CY.w[1]) && (!CY.w[0])) {
  80 #ifdef SET_STATUS_FLAGS
  81   // set status flags
  82   __set_status_flags (pfpsf, INVALID_EXCEPTION);
  83 #endif
  84   // x=y=0, return NaN
  85   res.w[1] = 0x7c00000000000000ull;
  86   res.w[0] = 0;
  87   BID_RETURN (res);
  88 }
  89 if (valid_y || ((y.w[1] & NAN_MASK64) == INFINITY_MASK64)) {
  90   // return 0
  91   if ((exponent_x > exponent_y)
  92       && ((y.w[1] & NAN_MASK64) != INFINITY_MASK64))
  93     exponent_x = exponent_y;
  94
  95   res.w[1] = sign_x | (((UINT64) exponent_x) << 49);
  96   res.w[0] = 0;
  97   BID_RETURN (res);
  98 }
  99 }
 100 if (!valid_y) {
 101   // y is Inf. or NaN
 102
 103   // test if y is NaN
 104   if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
 105 #ifdef SET_STATUS_FLAGS
 106     if ((y.w[1] & SNAN_MASK64) == SNAN_MASK64)  // y is sNaN
 107       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 108 #endif
 109     res.w[1] = CY.w[1] & QUIET_MASK64;
 110     res.w[0] = CY.w[0];
 111     BID_RETURN (res);
 112   }
 113   // y is Infinity?
 114   if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
 115     // return x
 116     res.w[1] = x.w[1];
 117     res.w[0] = x.w[0];
 118     BID_RETURN (res);
 119   }
 120   // y is 0
 121 #ifdef SET_STATUS_FLAGS
 122   // set status flags
 123   __set_status_flags (pfpsf, INVALID_EXCEPTION);
 124 #endif
 125   res.w[1] = 0x7c00000000000000ull;
 126   res.w[0] = 0;
 127   BID_RETURN (res);
 128 }
 129
 130 diff_expon = exponent_x - exponent_y;
 131
 132 if (diff_expon <= 0) {
 133   diff_expon = -diff_expon;
 134
 135   if (diff_expon > 34) {
 136     // |x|<|y| in this case
 137     res = x;
 138     BID_RETURN (res);
 139   }
 140   // set exponent of y to exponent_x, scale coefficient_y
 141   T = power10_table_128[diff_expon];
 142   __mul_128x128_to_256 (P256, CY, T);
 143
 144   if (P256.w[2] || P256.w[3]) {
 145     // |x|<|y| in this case
 146     res = x;
 147     BID_RETURN (res);
 148   }
 149
 150   CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63);
 151   CX2.w[0] = CX.w[0] << 1;
 152   if (__unsigned_compare_ge_128 (P256, CX2)) {
 153     // |x|<|y| in this case
 154     res = x;
 155     BID_RETURN (res);
 156   }
 157
 158   P128.w[0] = P256.w[0];
 159   P128.w[1] = P256.w[1];
 160   __div_128_by_128 (&CQ, &CR, CX, P128);
 161
 162   CX2.w[1] = (CR.w[1] << 1) | (CR.w[0] >> 63);
 163   CX2.w[0] = CR.w[0] << 1;
 164   if ((__unsigned_compare_gt_128 (CX2, P256))
 165       || (CX2.w[1] == P256.w[1] && CX2.w[0] == P256.w[0]
 166           && (CQ.w[0] & 1))) {
 167     __sub_128_128 (CR, P256, CR);
 168     sign_x ^= 0x8000000000000000ull;
 169   }
 170
 171   get_BID128_very_fast (&res, sign_x, exponent_x, CR);
 172   BID_RETURN (res);
 173 }
 174   // 2^64
 175 f64.i = 0x5f800000;
 176
 177 scale0 = 38;
 178 if (!CY.w[1])
 179   scale0 = 34;
 180
 181 while (diff_expon > 0) {
 182   // get number of digits in CX and scale=38-digits
 183   // fx ~ CX
 184   fx.d = (float) CX.w[1] * f64.d + (float) CX.w[0];
 185   bin_expon_cx = ((fx.i >> 23) & 0xff) - 0x7f;
 186   scale = scale0 - estimate_decimal_digits[bin_expon_cx];
 187   // scale = 38-estimate_decimal_digits[bin_expon_cx];
 188   D = CX.w[1] - power10_index_binexp_128[bin_expon_cx].w[1];
 189   if (D > 0
 190       || (!D && CX.w[0] >= power10_index_binexp_128[bin_expon_cx].w[0]))
 191     scale--;
 192
 193   if (diff_expon >= scale)
 194     diff_expon -= scale;
 195   else {
 196     scale = diff_expon;
 197     diff_expon = 0;
 198   }
 199
 200   T = power10_table_128[scale];
 201   __mul_128x128_low (CXS, CX, T);
 202
 203   __div_128_by_128 (&CQ, &CX, CXS, CY);
 204
 205   // check for remainder == 0
 206   if (!CX.w[1] && !CX.w[0]) {
 207     get_BID128_very_fast (&res, sign_x, exponent_y, CX);
 208     BID_RETURN (res);
 209   }
 210 }
 211
 212 CX2.w[1] = (CX.w[1] << 1) | (CX.w[0] >> 63);
 213 CX2.w[0] = CX.w[0] << 1;
 214 if ((__unsigned_compare_gt_128 (CX2, CY))
 215     || (CX2.w[1] == CY.w[1] && CX2.w[0] == CY.w[0] && (CQ.w[0] & 1))) {
 216   __sub_128_128 (CX, CY, CX);
 217   sign_x ^= 0x8000000000000000ull;
 218 }
 219
 220 get_BID128_very_fast (&res, sign_x, exponent_y, CX);
 221 BID_RETURN (res);
 222 }