libgcc/config/libbid/bid128_quantize.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_internal.h"
  31
  32 BID128_FUNCTION_ARG2 (bid128_quantize, x, y)
  33
  34      UINT256 CT;
  35      UINT128 CX, CY, T, CX2, CR, Stemp, res, REM_H, C2N;
  36      UINT64 sign_x, sign_y, remainder_h, carry, CY64, valid_x;
  37      int_float tempx;
  38      int exponent_x, exponent_y, digits_x, extra_digits, amount;
  39      int expon_diff, total_digits, bin_expon_cx, rmode, status;
  40
  41 valid_x = unpack_BID128_value (&sign_x, &exponent_x, &CX, x);
  42
  43   // unpack arguments, check for NaN or Infinity
  44 if (!unpack_BID128_value (&sign_y, &exponent_y, &CY, y)) {
  45     // y is Inf. or NaN
  46 #ifdef SET_STATUS_FLAGS
  47 if ((x.w[1] & SNAN_MASK64) == SNAN_MASK64)      // y is sNaN
  48   __set_status_flags (pfpsf, INVALID_EXCEPTION);
  49 #endif
  50
  51     // test if y is NaN
  52 if ((y.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
  53 #ifdef SET_STATUS_FLAGS
  54   if ((y.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
  55     // set status flags
  56     __set_status_flags (pfpsf, INVALID_EXCEPTION);
  57   }
  58 #endif
  59   if ((x.w[1] & 0x7c00000000000000ull) != 0x7c00000000000000ull) {
  60     res.w[1] = CY.w[1] & QUIET_MASK64;
  61     res.w[0] = CY.w[0];
  62   } else {
  63     res.w[1] = CX.w[1] & QUIET_MASK64;
  64     res.w[0] = CX.w[0];
  65   }
  66   BID_RETURN (res);
  67 }
  68     // y is Infinity?
  69 if ((y.w[1] & 0x7800000000000000ull) == 0x7800000000000000ull) {
  70   // check if x is not Inf.
  71   if (((x.w[1] & 0x7c00000000000000ull) < 0x7800000000000000ull)) {
  72     // return NaN
  73 #ifdef SET_STATUS_FLAGS
  74     // set status flags
  75     __set_status_flags (pfpsf, INVALID_EXCEPTION);
  76 #endif
  77     res.w[1] = 0x7c00000000000000ull;
  78     res.w[0] = 0;
  79     BID_RETURN (res);
  80   } else
  81     if (((x.w[1] & 0x7c00000000000000ull) <= 0x7800000000000000ull)) {
  82     res.w[1] = CX.w[1] & QUIET_MASK64;
  83     res.w[0] = CX.w[0];
  84     BID_RETURN (res);
  85   }
  86 }
  87
  88 }
  89
  90 if (!valid_x) {
  91   // test if x is NaN or Inf
  92   if ((x.w[1] & 0x7c00000000000000ull) == 0x7800000000000000ull) {
  93 #ifdef SET_STATUS_FLAGS
  94     // set status flags
  95     __set_status_flags (pfpsf, INVALID_EXCEPTION);
  96 #endif
  97     res.w[1] = 0x7c00000000000000ull;
  98     res.w[0] = 0;
  99     BID_RETURN (res);
 100   } else if ((x.w[1] & 0x7c00000000000000ull) == 0x7c00000000000000ull) {
 101     if ((x.w[1] & 0x7e00000000000000ull) == 0x7e00000000000000ull) {
 102 #ifdef SET_STATUS_FLAGS
 103       // set status flags
 104       __set_status_flags (pfpsf, INVALID_EXCEPTION);
 105 #endif
 106     }
 107     res.w[1] = CX.w[1] & QUIET_MASK64;
 108     res.w[0] = CX.w[0];
 109     BID_RETURN (res);
 110   }
 111   if (!CX.w[1] && !CX.w[0]) {
 112     get_BID128_very_fast (&res, sign_x, exponent_y, CX);
 113     BID_RETURN (res);
 114   }
 115 }
 116   // get number of decimal digits in coefficient_x
 117 if (CX.w[1]) {
 118   tempx.d = (float) CX.w[1];
 119   bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f + 64;
 120 } else {
 121   tempx.d = (float) CX.w[0];
 122   bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f;
 123 }
 124
 125 digits_x = estimate_decimal_digits[bin_expon_cx];
 126 if (CX.w[1] > power10_table_128[digits_x].w[1]
 127     || (CX.w[1] == power10_table_128[digits_x].w[1]
 128         && CX.w[0] >= power10_table_128[digits_x].w[0]))
 129   digits_x++;
 130
 131 expon_diff = exponent_x - exponent_y;
 132 total_digits = digits_x + expon_diff;
 133
 134 if ((UINT32) total_digits <= 34) {
 135   if (expon_diff >= 0) {
 136     T = power10_table_128[expon_diff];
 137     __mul_128x128_low (CX2, T, CX);
 138     get_BID128_very_fast (&res, sign_x, exponent_y, CX2);
 139     BID_RETURN (res);
 140   }
 141 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 142 #ifndef IEEE_ROUND_NEAREST
 143   rmode = rnd_mode;
 144   if (sign_x && (unsigned) (rmode - 1) < 2)
 145     rmode = 3 - rmode;
 146 #else
 147   rmode = 0;
 148 #endif
 149 #else
 150   rmode = 0;
 151 #endif
 152   // must round off -expon_diff digits
 153   extra_digits = -expon_diff;
 154   __add_128_128 (CX, CX, round_const_table_128[rmode][extra_digits]);
 155
 156   // get P*(2^M[extra_digits])/10^extra_digits
 157   __mul_128x128_to_256 (CT, CX, reciprocals10_128[extra_digits]);
 158
 159   // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
 160   amount = recip_scale[extra_digits];
 161   CX2.w[0] = CT.w[2];
 162   CX2.w[1] = CT.w[3];
 163   if (amount >= 64) {
 164     CR.w[1] = 0;
 165     CR.w[0] = CX2.w[1] >> (amount - 64);
 166   } else {
 167     __shr_128 (CR, CX2, amount);
 168   }
 169
 170 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 171 #ifndef IEEE_ROUND_NEAREST
 172   if (rnd_mode == 0)
 173 #endif
 174     if (CR.w[0] & 1) {
 175       // check whether fractional part of initial_P/10^extra_digits is
 176       // exactly .5 this is the same as fractional part of
 177       // (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
 178
 179       // get remainder
 180       if (amount >= 64) {
 181         remainder_h = CX2.w[0] | (CX2.w[1] << (128 - amount));
 182       } else
 183         remainder_h = CX2.w[0] << (64 - amount);
 184
 185       // test whether fractional part is 0
 186       if (!remainder_h
 187           && (CT.w[1] < reciprocals10_128[extra_digits].w[1]
 188               || (CT.w[1] == reciprocals10_128[extra_digits].w[1]
 189                   && CT.w[0] < reciprocals10_128[extra_digits].w[0]))) {
 190         CR.w[0]--;
 191       }
 192     }
 193 #endif
 194
 195 #ifdef SET_STATUS_FLAGS
 196   status = INEXACT_EXCEPTION;
 197
 198   // get remainder
 199   if (amount >= 64) {
 200     REM_H.w[1] = (CX2.w[1] << (128 - amount));
 201     REM_H.w[0] = CX2.w[0];
 202   } else {
 203     REM_H.w[1] = CX2.w[0] << (64 - amount);
 204     REM_H.w[0] = 0;
 205   }
 206
 207   switch (rmode) {
 208   case ROUNDING_TO_NEAREST:
 209   case ROUNDING_TIES_AWAY:
 210     // test whether fractional part is 0
 211     if (REM_H.w[1] == 0x8000000000000000ull && !REM_H.w[0]
 212         && (CT.w[1] < reciprocals10_128[extra_digits].w[1]
 213             || (CT.w[1] == reciprocals10_128[extra_digits].w[1]
 214                 && CT.w[0] < reciprocals10_128[extra_digits].w[0])))
 215       status = EXACT_STATUS;
 216     break;
 217   case ROUNDING_DOWN:
 218   case ROUNDING_TO_ZERO:
 219     if (!(REM_H.w[1] | REM_H.w[0])
 220         && (CT.w[1] < reciprocals10_128[extra_digits].w[1]
 221             || (CT.w[1] == reciprocals10_128[extra_digits].w[1]
 222                 && CT.w[0] < reciprocals10_128[extra_digits].w[0])))
 223       status = EXACT_STATUS;
 224     break;
 225   default:
 226     // round up
 227     __add_carry_out (Stemp.w[0], CY64, CT.w[0],
 228                      reciprocals10_128[extra_digits].w[0]);
 229     __add_carry_in_out (Stemp.w[1], carry, CT.w[1],
 230                         reciprocals10_128[extra_digits].w[1], CY64);
 231     if (amount < 64) {
 232       C2N.w[1] = 0;
 233       C2N.w[0] = ((UINT64) 1) << amount;
 234       REM_H.w[0] = REM_H.w[1] >> (64 - amount);
 235       REM_H.w[1] = 0;
 236     } else {
 237       C2N.w[1] = ((UINT64) 1) << (amount - 64);
 238       C2N.w[0] = 0;
 239       REM_H.w[1] >>= (128 - amount);
 240     }
 241     REM_H.w[0] += carry;
 242     if (REM_H.w[0] < carry)
 243       REM_H.w[1]++;
 244     if (__unsigned_compare_ge_128 (REM_H, C2N))
 245       status = EXACT_STATUS;
 246   }
 247
 248   __set_status_flags (pfpsf, status);
 249
 250 #endif
 251   get_BID128_very_fast (&res, sign_x, exponent_y, CR);
 252   BID_RETURN (res);
 253 }
 254 if (total_digits < 0) {
 255   CR.w[1] = CR.w[0] = 0;
 256 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 257 #ifndef IEEE_ROUND_NEAREST
 258   rmode = rnd_mode;
 259   if (sign_x && (unsigned) (rmode - 1) < 2)
 260     rmode = 3 - rmode;
 261   if (rmode == ROUNDING_UP)
 262     CR.w[0] = 1;
 263 #endif
 264 #endif
 265 #ifdef SET_STATUS_FLAGS
 266   __set_status_flags (pfpsf, INEXACT_EXCEPTION);
 267 #endif
 268   get_BID128_very_fast (&res, sign_x, exponent_y, CR);
 269   BID_RETURN (res);
 270 }
 271   // else  more than 34 digits in coefficient
 272 #ifdef SET_STATUS_FLAGS
 273 __set_status_flags (pfpsf, INVALID_EXCEPTION);
 274 #endif
 275 res.w[1] = 0x7c00000000000000ull;
 276 res.w[0] = 0;
 277 BID_RETURN (res);
 278
 279 }