libgcc/config/libbid/bid64_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 #include "bid_internal.h"
  30
  31 #define MAX_FORMAT_DIGITS     16
  32 #define DECIMAL_EXPONENT_BIAS 398
  33 #define MAX_DECIMAL_EXPONENT  767
  34
  35 #if DECIMAL_CALL_BY_REFERENCE
  36
  37 void
  38 bid64_quantize (UINT64 * pres, UINT64 * px,
  39                 UINT64 *
  40                 py _RND_MODE_PARAM _EXC_FLAGS_PARAM _EXC_MASKS_PARAM
  41                 _EXC_INFO_PARAM) {
  42   UINT64 x, y;
  43 #else
  44
  45 UINT64
  46 bid64_quantize (UINT64 x,
  47                 UINT64 y _RND_MODE_PARAM _EXC_FLAGS_PARAM
  48                 _EXC_MASKS_PARAM _EXC_INFO_PARAM) {
  49 #endif
  50   UINT128 CT;
  51   UINT64 sign_x, sign_y, coefficient_x, coefficient_y, remainder_h, C64,
  52     valid_x;
  53   UINT64 tmp, carry, res;
  54   int_float tempx;
  55   int exponent_x, exponent_y, digits_x, extra_digits, amount, amount2;
  56   int expon_diff, total_digits, bin_expon_cx;
  57   unsigned rmode, status;
  58
  59 #if DECIMAL_CALL_BY_REFERENCE
  60 #if !DECIMAL_GLOBAL_ROUNDING
  61   _IDEC_round rnd_mode = *prnd_mode;
  62 #endif
  63   x = *px;
  64   y = *py;
  65 #endif
  66
  67   valid_x = unpack_BID64 (&sign_x, &exponent_x, &coefficient_x, x);
  68   // unpack arguments, check for NaN or Infinity
  69   if (!unpack_BID64 (&sign_y, &exponent_y, &coefficient_y, y)) {
  70     // Inf. or NaN or 0
  71 #ifdef SET_STATUS_FLAGS
  72     if ((x & SNAN_MASK64) == SNAN_MASK64)       // y is sNaN
  73       __set_status_flags (pfpsf, INVALID_EXCEPTION);
  74 #endif
  75
  76     // x=Inf, y=Inf?
  77     if (((coefficient_x << 1) == 0xf000000000000000ull)
  78         && ((coefficient_y << 1) == 0xf000000000000000ull)) {
  79       res = coefficient_x;
  80       BID_RETURN (res);
  81     }
  82     // Inf or NaN?
  83     if ((y & 0x7800000000000000ull) == 0x7800000000000000ull) {
  84 #ifdef SET_STATUS_FLAGS
  85       if (((y & 0x7e00000000000000ull) == 0x7e00000000000000ull)        // sNaN
  86           || (((y & 0x7c00000000000000ull) == 0x7800000000000000ull) && //Inf
  87               ((x & 0x7c00000000000000ull) < 0x7800000000000000ull)))
  88         __set_status_flags (pfpsf, INVALID_EXCEPTION);
  89 #endif
  90       if ((y & NAN_MASK64) != NAN_MASK64)
  91         coefficient_y = 0;
  92       if ((x & NAN_MASK64) != NAN_MASK64) {
  93         res = 0x7c00000000000000ull | (coefficient_y & QUIET_MASK64);
  94         if (((y & NAN_MASK64) != NAN_MASK64) && ((x & NAN_MASK64) == 0x7800000000000000ull))
  95                 res = x;
  96         BID_RETURN (res);
  97       }
  98     }
  99   }
 100   // unpack arguments, check for NaN or Infinity
 101   if (!valid_x) {
 102     // x is Inf. or NaN or 0
 103
 104     // Inf or NaN?
 105     if ((x & 0x7800000000000000ull) == 0x7800000000000000ull) {
 106 #ifdef SET_STATUS_FLAGS
 107       if (((x & 0x7e00000000000000ull) == 0x7e00000000000000ull)        // sNaN
 108           || ((x & 0x7c00000000000000ull) == 0x7800000000000000ull))    //Inf
 109         __set_status_flags (pfpsf, INVALID_EXCEPTION);
 110 #endif
 111       if ((x & NAN_MASK64) != NAN_MASK64)
 112         coefficient_x = 0;
 113       res = 0x7c00000000000000ull | (coefficient_x & QUIET_MASK64);
 114       BID_RETURN (res);
 115     }
 116
 117     res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, 0);
 118     BID_RETURN (res);
 119   }
 120   // get number of decimal digits in coefficient_x
 121   tempx.d = (float) coefficient_x;
 122   bin_expon_cx = ((tempx.i >> 23) & 0xff) - 0x7f;
 123   digits_x = estimate_decimal_digits[bin_expon_cx];
 124   if (coefficient_x >= power10_table_128[digits_x].w[0])
 125     digits_x++;
 126
 127   expon_diff = exponent_x - exponent_y;
 128   total_digits = digits_x + expon_diff;
 129
 130   // check range of scaled coefficient
 131   if ((UINT32) (total_digits + 1) <= 17) {
 132     if (expon_diff >= 0) {
 133       coefficient_x *= power10_table_128[expon_diff].w[0];
 134       res = very_fast_get_BID64 (sign_x, exponent_y, coefficient_x);
 135       BID_RETURN (res);
 136     }
 137     // must round off -expon_diff digits
 138     extra_digits = -expon_diff;
 139 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 140 #ifndef IEEE_ROUND_NEAREST
 141     rmode = rnd_mode;
 142     if (sign_x && (unsigned) (rmode - 1) < 2)
 143       rmode = 3 - rmode;
 144 #else
 145     rmode = 0;
 146 #endif
 147 #else
 148     rmode = 0;
 149 #endif
 150     coefficient_x += round_const_table[rmode][extra_digits];
 151
 152     // get P*(2^M[extra_digits])/10^extra_digits
 153     __mul_64x64_to_128 (CT, coefficient_x,
 154                         reciprocals10_64[extra_digits]);
 155
 156     // now get P/10^extra_digits: shift C64 right by M[extra_digits]-128
 157     amount = short_recip_scale[extra_digits];
 158     C64 = CT.w[1] >> amount;
 159 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 160 #ifndef IEEE_ROUND_NEAREST
 161     if (rnd_mode == 0)
 162 #endif
 163       if (C64 & 1) {
 164         // check whether fractional part of initial_P/10^extra_digits
 165         // is exactly .5
 166         // this is the same as fractional part of
 167         //   (initial_P + 0.5*10^extra_digits)/10^extra_digits is exactly zero
 168
 169         // get remainder
 170         amount2 = 64 - amount;
 171         remainder_h = 0;
 172         remainder_h--;
 173         remainder_h >>= amount2;
 174         remainder_h = remainder_h & CT.w[1];
 175
 176         // test whether fractional part is 0
 177         if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits])) {
 178           C64--;
 179         }
 180       }
 181 #endif
 182
 183 #ifdef SET_STATUS_FLAGS
 184     status = INEXACT_EXCEPTION;
 185     // get remainder
 186     remainder_h = CT.w[1] << (64 - amount);
 187     switch (rmode) {
 188     case ROUNDING_TO_NEAREST:
 189     case ROUNDING_TIES_AWAY:
 190       // test whether fractional part is 0
 191       if ((remainder_h == 0x8000000000000000ull)
 192           && (CT.w[0] < reciprocals10_64[extra_digits]))
 193         status = EXACT_STATUS;
 194       break;
 195     case ROUNDING_DOWN:
 196     case ROUNDING_TO_ZERO:
 197       if (!remainder_h && (CT.w[0] < reciprocals10_64[extra_digits]))
 198         status = EXACT_STATUS;
 199       //if(!C64 && rmode==ROUNDING_DOWN) sign_s=sign_y;
 200       break;
 201     default:
 202       // round up
 203       __add_carry_out (tmp, carry, CT.w[0],
 204                        reciprocals10_64[extra_digits]);
 205       if ((remainder_h >> (64 - amount)) + carry >=
 206           (((UINT64) 1) << amount))
 207         status = EXACT_STATUS;
 208       break;
 209     }
 210     __set_status_flags (pfpsf, status);
 211 #endif
 212
 213     res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, C64);
 214     BID_RETURN (res);
 215   }
 216
 217   if (total_digits < 0) {
 218 #ifdef SET_STATUS_FLAGS
 219     __set_status_flags (pfpsf, INEXACT_EXCEPTION);
 220 #endif
 221     C64 = 0;
 222 #ifndef IEEE_ROUND_NEAREST_TIES_AWAY
 223 #ifndef IEEE_ROUND_NEAREST
 224     rmode = rnd_mode;
 225     if (sign_x && (unsigned) (rmode - 1) < 2)
 226       rmode = 3 - rmode;
 227     if (rmode == ROUNDING_UP)
 228       C64 = 1;
 229 #endif
 230 #endif
 231     res = very_fast_get_BID64_small_mantissa (sign_x, exponent_y, C64);
 232     BID_RETURN (res);
 233   }
 234   // else  more than 16 digits in coefficient
 235 #ifdef SET_STATUS_FLAGS
 236   __set_status_flags (pfpsf, INVALID_EXCEPTION);
 237 #endif
 238   res = 0x7c00000000000000ull;
 239   BID_RETURN (res);
 240
 241 }